A computational Eulerian–Lagrangian model (ORSA2D_WT) is used for modelling the movement of floating rigid bodies on the water surface. The two-dimensional transport is computed with a dynamic approach, modifying existing formulations for the transport of bodies within fluid flows for the case of floating bodies, by adopting suitable added mass, drag and side coefficients. An original formulation for planar rotation is proposed, which includes the effect of the hydrodynamic torque and a resistance term, named added inertia, based on the difference between the angular velocity of the flow and that of the body. The value of the added inertia coefficient is calibrated against experiments made on purpose, involving the transport of a cylinder in a flume with two side obstacles. The calibrated code is applied to a slightly larger set of experiments for its preliminary evaluation. The outcome of the simulations shows that the streamwise and transversal displacements are well modelled, while some inaccuracies arise when considering the cylinder orientation. The effects of the initial conditions on the cylinders' trajectory and rotation are discussed, showing their influence on the evolution of the rotation angles.

Most of the scientific literature dealing with floating body transport considers it as a complementary aspect of river dynamics (Abbe & Montgomery 1996; Bocchiola 2011; Gurnell 2012; Picco et al. 2017), monitoring wood motion during floods (e.g. Ravazzolo et al. 2014) and considering its effect on river bed forms or on the aquatic ecosystem.

Less attention is paid to the relation between large wood (LW) and flood risk. However, it has been proven that transported LW may significantly intensify the drawback of a flood (Comiti et al. 2008), especially if natural or artificial channel narrowing is present, as in urban areas. LW can also gather in reservoirs and at spillways, provoking backwater rise (e.g. Schalko et al. 2018) and representing an extra issue for reservoirs managers.

In general, the problem of protection against floating transport is solved with practical measures (Uchiogi et al. 1996; Bradley et al. 2005) like those employed for debris flow. However, while for debris flows numerical models are quite popular (Silva et al. 2016; Canelas et al. 2017), the presence of LW is generally not included in the numerical simulations of flood risk or dam breach (e.g. Costabile & Macchione 2015; Macchione et al. 2016). In recent years, debris management strategies in case of floods (including LW) have been developed, also with the aid of computer intelligence methods (e.g. Fotovatikhah et al. 2018), to reduce the additional hazard related to the presence of floating bodies during floods. Some attempts to achieve a physically-based design of safety structures have been reported (Denk et al. 2008; Comiti et al. 2012; Schmocker & Weitbrecht 2013), and a few ground-breaking applications of innovative software which considers LW effects can be found (e.g. Ruiz-Villanueva et al. 2013).

The prediction of LW-related flood hazard would benefit from the use of numerical models that include the transport of floating debris in the hydraulic simulations. Some numerical models do already exist, although it is still debated which should be the most reliable and effective way to couple the dynamics of the two phases, namely the discrete floating elements and the continuous water flow.

One possible approach is to consider a relatively large amount of debris entrained in the flood, and to try to predict its final position following the streamlines along the entire river basin (Mazzorana et al. 2011). The planar displacement of large volumes of wood is hence computed, disregarding the physical response of the single floating object on the water.

A different method resides in the adoption of a Lagrangian–Lagrangian approach, usually by applying the Smoothed Particles Hydrodynamics (SPH) technique to the equations of fluid motion. By considering water as a group of particles and floating objects as single large particles or rigid shells, SPH allows the one- or two-way coupling of the water and wood dynamics, computing their reciprocal influence with high precision, both in two- and three-dimensions (2D and 3D, Solenthaler et al. 2011; Prakash et al. 2014; Amicarelli et al. 2015).

Another option is the use of hybrid 2D methods, which couple two different techniques for the solution of the two phases. These methods usually estimate the flow velocity by solving the Shallow-Water Equations (SWE) and evaluate the motion of the rigid body through an Eulerian–Lagrangian approach, applying either a kinematic model (Ruiz-Villanueva et al. 2014) or a dynamic one (Alonso 2004; Stockstill et al. 2009; Persi et al. 2018a; Petaccia et al. 2018) to compute the motion of the discontinuous phase (i.e. LW). In kinematic models the LW velocity is assigned directly from the speed of the fluid phase, while dynamic ones require the computation of the forces and torques exerted by the flow on the rigid body and are expected to be more consistent with the physics of the phenomenon.

Overall, Eulerian–Lagrangian methods are halfway between faster cell-by-cell models based on a simplified wood transport along streamlines (Mazzorana et al. 2011) and the more accurate, but computationally expensive, description of the detailed interaction between flow and logs, such as the one obtained through the SPH technique. The critical aspect of the Eulerian–Lagrangian approach is the evaluation of the correct forces acting on the body and their computation, which can be integrated on the body length or approximated. There is no shared opinion on the most appropriate dynamic description for floating wood transport, and existing approaches do not always include the same forces. In general, only the drag force is included, disregarding the transversal component of the hydrodynamic force (i.e. side force). Stockstill et al. (2009) also neglect the added mass force. Furthermore, these forces are often computed under the hypothesis of complete submergence, while floating LW are in general semi-submerged objects, thus requiring a modification of the coefficients. Regarding the rotation formulation, different approaches can be found in the literature: Stockstill et al. (2009) implement the Euler equation while Alonso (2004) uses the angular momentum equation, but none of them provide a systematic validation of their method.

The aim of this paper is to present the model ORSA2D_WT, a dynamic Eulerian–Lagrangian model developed by the authors, which includes all the relevant forces for LW translation and adopts an original formulation for the computation of rotation, proving its validity through the comparison with laboratory experiments.

Since the methodology proposed in existing dynamic LW transport models is not univocal, the forces ruling the motion are obtained from a more general description, which is then adapted to the specific case of floating bodies transport. The translation equation, already drafted in Persi et al. (2018a), is based on the Basset–Boussinesque–Oseen (BBO) equation (e.g. Corrsin & Lumley 1956), which gives a general physical interpretation of the unsteady forces on spherical particles settling in a fluid at rest. The equation, which has been proven to be valid to solve problems related to the motion of small particles at low Reynolds numbers (Magnaudet & Eames 2000; Tagawa et al. 2013), is extended to higher Reynolds number flows and to the case of non-spherical particles (Yin et al. 2003; Mandø & Rosendahl 2010). The shape-dependent forces (drag, lift, side force and added mass force) can be adjusted thanks to the choice of proper coefficients.

In literature, dynamic approaches for the computation of rotation depend on the body shape. The torque on spheres is generally computed by applying the conservation of linear and angular momentum (Bagchi & Balachandar 2002), while for non-spherical bodies, it is computed as the sum of three components (offset torque, resistance torque and cross term, which are meaningful for 3D modelling, as from Mandø & Rosendahl (2010)). Other authors dealing with LW transport adopted different formulations but do not exhaustively prove their validity. In Persi et al. (2018a) a formulation similar to the one proposed by Mandø & Rosendahl (2010) is adapted for 2D rotation, but the lack of comparison with experimental data did not allow its validation. As will be shown in the present paper, that formulation does not match well with the laboratory experiments. For this reason, this paper proposes an alternative formulation for rotation in 2D modelling, which includes the centre of mass torque and a resistance term, named added inertia, which is proportional to the relative angular acceleration through a coefficient, defined as the added inertia coefficient.

The resulting Eulerian–Lagrangian model is calibrated against experiments carried out on purpose, to determine the proper value for the added inertia coefficient. Since the initial conditions appear to be a key parameter for simulation of floating body transport, both for the acquisition inaccuracy and for their unpredictability in real events, a sensitivity analysis is carried out, verifying their effect on the simulation.

Numerical model of wood transport in shallow water flows

The Eulerian–Lagrangian model, named ORSA2D_WT, is obtained by including a Lagrangian routine in ORSA2D, a 2D finite volume Eulerian solver of the SWE (Petaccia et al. 2010; Petaccia et al. 2016), which implements a Roe's Riemann solver and is first-order accurate in time and space (Roe 1981).

Although nowadays 3D models are also employed for flood risk simulations, the choice of a 2D hydraulic model is consistent with the characteristics of the physical problem here considered, which is the planar transportation of LW debris. In addition, 2D models are still more affordable from the computational point of view. However, in the case of 3D flows and obstruction formation with vertical accumulation, both the hydraulic and the DEM method should shift to 3D, in order to provide a more realistic description of this phenomenon.

The SWE are first solved, to obtain the water depth and the linear and angular velocities. Then, a localization algorithm, which is a combination of the nearest neighbour and the random walk techniques (Petaccia et al. 2018), allows the identification of the cell (or cells) of the mesh in which the body is placed. The algorithm requires few iterations and does not increase the computational time.

The assignment of the correct flow variables to each body, independently from the relative body-cell dimensions, is performed with a two-step linear interpolation from cell-centred values as detailed in Petaccia et al. (2018).

Once the body location is identified and the hydraulic variables are known, the flotation condition is verified with a force balance (Braudrick & Grant 2001) and the Lagrangian transport of LW is modelled with a Discrete Element Method (DEM), which includes the solution of planar translation and planar rotation equations, as well as those for the body trajectory and orientation:
(1)
(2)
(3)
Here, the suffixes b and f refer to the body and to the fluid respectively, m is the mass, A is the reference area (LD, cylinder length by diameter), V=(u, v) is the velocity vector, ω is the angular velocity, CD, CS and CAM are the drag, side and added mass coefficients, ρ is the density, I is the cylinder moment of inertia, r is the distance between the body centre of mass and the application point of each force F acting on the body, CAI is the added inertia coefficient, X=(x, y) is the position vector and ϑ is the body orientation.

The translation equation (Equation (1)) is adapted from the Maxey–Riley equation (Maxey & Riley 1983), which is an extension of the BBO equation to the case of non-uniform flows and is valid for rigid small spheres in non-uniform creeping flows (i.e. with a particle Reynolds number ). It includes the drag force, the side force (which corresponds to the lift force on a horizontal plane), the added mass force and the pressure gradient force, in which the total derivative of the flow velocity appears (Corrsin & Lumley 1956).

Equation (2) is an original formulation for the computation of body rotation on the water surface, which differs from the formulation which is generally employed and can be found in Mandø & Rosendahl (2010). Their formulation includes an offset torque, due to the application of the hydrodynamic forces in the centre of pressure, which does not coincide with the body centre of mass, and a resistance torque, which originates from the integration of resistance terms on the body main dimension. It was adapted to the case of floating bodies and presented in Persi et al. (2018a). In order to derive Equation (2), a different approach is followed. The first term on the right-hand side takes into account the distribution of the forces (drag, side, added mass, pressure gradient) on the main body length, and refers to the centre of mass of the body, not to the centre of pressure. This choice is justified by the fact that considered bodies are not totally submerged in water, hence the pressure does not present the same distribution as the elongated objects studied by Mandø & Rosendahl (2010). Here, the focus is rather on the force variation along the body, which is a consequence of the presence of velocity gradients (Figure 1).

Figure 1

(a) Velocity gradients on large bodies; (b) subdivision of the diameter (for a sphere) and of the axis (for a cylinder) in four subsections. The body centre of mass (CM) is also shown.

Figure 1

(a) Velocity gradients on large bodies; (b) subdivision of the diameter (for a sphere) and of the axis (for a cylinder) in four subsections. The body centre of mass (CM) is also shown.

Close modal
To properly represent this phenomenon, which is particularly important in non-uniform flows, forces are evaluated not only in the centre of mass but in four different sections along the body (Figure 1). For the solution of Equation (1) the forces are applied at the centre of each of the four segments, i.e. points 1–4 in Figure 1, and, for rotation, r becomes the distance between each of the four points and the centre of mass (CM). This allows the computation of the torque with respect to the centre of mass. Note that the added mass term, included in the forces in the first term of Equation (2), is computed by considering the relative linear acceleration:
(4)
where Vol is the body volume.

As suggested by Mandø & Rosendahl (2010), a second term, which acts as a resistance to rotation, is needed. This resistance term (second term on the right-hand side of Equation (2)) is named added inertia torque and presents a formulation analogue to the one of the added mass force in Equation (1): the rate of change of the body angular velocity is connected to the difference between the body and the flow angular accelerations, with a proportionality coefficient CAI which is named added inertia coefficient. The local flow angular acceleration is expressed as the total derivative of the flow vorticity, in order to take into account both the temporal and spatial variations.

The total derivatives of the flow linear and angular velocities require the computation of the partial derivatives in x and y, which are computed with the trapezoidal interpolation suggested by Hirsch (2007) for the case of a triangular mesh:
(5)
(6)
(7)
Here, N is the number of nodes of the cell of interest, Ψf represents the considered variable (either uf, vf, ωf), Ar is the cell area, x and y are the streamwise and transversal coordinates.

The computation of gradients takes advantage of the interpolation from cell central values to nodal values (Petaccia et al. 2018) and is performed only for those cells in which the computational points are found.

The SWE solver and the DEM are one way coupled, i.e. the effect of LW on the water flow is currently not included. The final accelerations (linear and angular) are computed and the body position and orientation are then adjusted. Finally, the time step, which is the same for both parts of the code, is evaluated according to the Courant–Friedrichs–Lewy (CFL) condition (Morales-Hernández et al. 2013).

Within one time step, during which the flow field is kept unchanged, the DEM is solved with an explicit fourth order Runge–Kutta (RK) time-integration scheme (see e.g. Chow 1979). The position and orientation vector ξ= (x, y, ϑ) and the linear and angular velocity vector Ψb = (ub, vb, ωb) are computed according to the following equations:
(8)
(9)
where N is the time step index, h is the time step and Δ and Δ are the increments at each kth RK-step for the computation of the weighted average. Within the main time loop, the increments of the position vector ξ are computed according to Equation (3) by considering linear and angular velocities from the previous step, while the increments of the velocity vector Ψ are computed according to Equations (1) and (2) by applying a convenient displacement based on the solution at the previous step.
The sequence of the RK strategy is summarized in the following equations, where F() represents the right side of Equations (1) and (2):
(10a)
(10b)
(10c)
(10d)

Note that t is the instant considered for the computation, which is kept constant within the four steps of the RK cycle, as well as the flow field. The choice of the RK scheme does not play a major role in the approximation of the partial differential equation, since the time step computed with the CFL condition is already small (of the order of 3 × 10–4 seconds) due to the reduced cells dimension (sides of 1 cm approximately).

For the correct computation of the hydrodynamic forces, an accurate estimation of the corresponding coefficient is required. The drag and side coefficients in Equations (1) and (2) are computed according to the body shape. For spheres, the drag coefficient is set as a function of the particle Reynolds number, while the side coefficient is set as a function of literature value C’ (e.g. Truscott & Techet 2009) times the sign of the product among the relative velocity and the relative vorticity (Table 1). For cylinders, both the drag and side coefficients vary with body orientation with the flow, as shown in the literature (Gippel et al. 1996; Hoang et al. 2015). The values implemented in ORSA2D_WT (Table 1) are the interpolant curves obtained from the results of a specific experimental campaign performed on semi-submerged cylinders (Persi et al. 2018b).

Table 1

Drag and side coefficients for a sphere and a cylinder

SphereCylinder
CD    
  
  
  
CS   
SphereCylinder
CD    
  
  
  
CS   

Suffix i stands for each body subsection; ϑ is the relative angle among the cylinder axis and the relative velocity, in sexagesimal degrees.

The added mass coefficient is assumed constant, since no data exist on its variation for semi-submerged cylindrical bodies. Literature values refer to totally submerged bodies and need to be adapted to the case of floating objects. For submerged bodies, the value 1 is proposed for cylinders (Dean & Dalrymple 1991) and 0.5 for spheres. Note that in the expression of the added mass force in Equation (1), the added mass coefficient should be equal to 2 for cylinders and 1 for spheres, since the force is halved with respect to the standard equation. Overall, the value of the added mass coefficient for a cylinder means that the added mass corresponds to a mass of fluid of the same volume as the sample.

For the case of semi-submerged bodies, the added mass coefficient is computed by dividing the submerged added mass volume (dashed line in Figure 2(b)) by the body volume (solid line in Figure 2(a)). The value obtained, 1.41, refers to a wooden cylinder with density 774 kg m–3. The same procedure can be repeated for the case of semi-submerged spheres, resulting in the value of about 0.69, for sphere density equal to 694 kg m–3.

Figure 2

Diagram for the estimation of the added mass coefficient for a cylinder. (a) Volume of the cylinder (solid line) and added volume (dashed lines); (b) volumes are adjusted for sinking and the submerged height for the added volume can be estimated.

Figure 2

Diagram for the estimation of the added mass coefficient for a cylinder. (a) Volume of the cylinder (solid line) and added volume (dashed lines); (b) volumes are adjusted for sinking and the submerged height for the added volume can be estimated.

Close modal

Experimental set-up for model validation

Elongated bodies and rotating spheres require the joint calibration of the translation and rotation formulations, since the two phenomena are strictly connected. In the literature, analytical, numerical and laboratory experiments can be found to calibrate the translation formulation, as shown in Persi et al. (2018a). In some cases, the rotation of spheres is considered (Truscott & Techet 2009) while, to the authors' knowledge, no detailed analysis exists with reference to the rotation of elongated bodies. In order to fill the gap of joint calibration, 16 experiments with cylinders floating on the water surface were performed in a laboratory flume at the University of Zaragoza.

The prismatic straight flume is 3.25 m long, 0.24 m wide and has transparent poly methyl methacrylate (PMMA) walls with a height of 0.16 m. It is equipped with two rectangular side obstacles (0.08 × 0.07 m, Figure 3(a)) the first placed on the right-hand side at 2.00 m and the second on the left-hand side, at 2.50 m from the inlet, and has a horizontal bottom. The obstacles are made of plastic and present a smooth surface comparable to that of the acrylic. Steady-state tests were performed, with the pumping system providing a constant discharge of 15.3 m3 h–1, measured by an electromagnetic flow meter (COPA-XE DE43F by ABB, which has an accuracy of 0.5% of the maximum rate, 60 m3 h–1).

Figure 3

(a) Sketch of the flume with two alternate side obstacles; (b) vertical section of the device for cylindrical objects release. The insertion of the cylinder, its motion, as well as the motion of the piston, are shown by arrows; (c) orthorectified image of the portion of the flume framed by the camera (dashed line in (a)).

Figure 3

(a) Sketch of the flume with two alternate side obstacles; (b) vertical section of the device for cylindrical objects release. The insertion of the cylinder, its motion, as well as the motion of the piston, are shown by arrows; (c) orthorectified image of the portion of the flume framed by the camera (dashed line in (a)).

Close modal

Cylindrical wooden samples (L = 0.073 m, d = 0.01 m, ρ = 774 kg m–3) are released about 1.25 m downstream the inlet section, perpendicular to the flow direction thanks to an ad hoc built device (Figure 3(b)). The cylinder is inserted from the top in the vertical box, and then it is slowly pushed forward in a horizontal box thanks to a pneumatic cylinder, until it reaches a slot and falls in the water. The device is placed at about 0.02 m above the water surface. In this way, cylinders are dropped into the water, flow initially totally submerged and then re-emerge on the water surface in an un-controlled manner, not influenced by the operator manual dexterity. However, the device does not allow replication of the exact initial condition for each test, since the push of the piston and the water entering do not occur exactly in the same conditions. Note that cylinders are released in sequence, in order to avoid any interference and to focus on the description of the motion of the singular element, in agreement with calibration purposes.

The experiments are recorded from top view (Figure 3(c)), hanging a Nikon camera (Nikon D810, with a Nikon 24–70 mm f/2.8G lens, set at a focal distance of 24 mm) which provides videos with a resolution of 30 fps. Due to the distortion introduced by the lens, the orthorectification of the images is performed with a MATLAB code (The MathWorks, Inc., Natick, Massachusetts, USA). The log tracking is performed frame by frame with the software TRACKER, which allows a semi-automatic analysis (Brown & Christian 2011). The time evolution of the body position and orientation are thus obtained.

Hydraulic simulation

The hydraulic simulation, with an unstructured triangular mesh of about 8,500 elements, is performed with a constant inlet discharge of 4.25 L s–1 as upstream boundary condition and critical flow as downstream boundary condition. The resistance to the flow is modelled by the Manning coefficient n = 0.01 s m–1/3. The side obstacles are represented as side walls, which has been shown by Petaccia et al. (2018) to be the method that provides the best performances.

Measured and simulated water levels are compared in Figure 4 for the right and left side of the flume. The comparison of depth-averaged velocities is shown in Figure 5. Velocity measurements were performed with a digital flowmeter (MiniAri20, with the probe Mini 95.0004 by PCE Instruments) at two points along each vertical (each point, in Figure 5(a)), at approximately one-third and two-thirds of the water depth. The final velocity is computed as the average of the two measurements. According to the measurements, the Froude number, computed with the water level and velocity upstream the obstacles, is about 0.35, and the particle Reynolds number is 3,000.

Figure 4

Simulated and measured water level for the configuration with two rectangular side obstacles. (a) Right side, with the first obstacle at x = 2.00m; (b) left side, with second obstacle at x = 2.50m.

Figure 4

Simulated and measured water level for the configuration with two rectangular side obstacles. (a) Right side, with the first obstacle at x = 2.00m; (b) left side, with second obstacle at x = 2.50m.

Close modal
Figure 5

Comparison of the measured and simulated depth-averaged flow velocities, for the case with two side obstacles. Beware different axes limits. (a) Planar sketch of the flume with points of measure; (b) streamwise velocity component u; (c) transversal velocity component v.

Figure 5

Comparison of the measured and simulated depth-averaged flow velocities, for the case with two side obstacles. Beware different axes limits. (a) Planar sketch of the flume with points of measure; (b) streamwise velocity component u; (c) transversal velocity component v.

Close modal

Measured and simulated water levels are well comparable, with a determination coefficient higher than 0.99 for both sides. The largest differences are observed at x = 2.25 m, with an overestimation of the water level of about 0.006 m for both sides, corresponding to an error of 8% over the maximum measured level.

Regarding the velocity field, in the first part of the flume the transversal component is zero and the streamwise component is nearly constant in the entire section, as shown by the measured and simulated values. Then, both components increase, and vary across the section. In general, the simulation of the longitudinal velocity component is more accurate than that of the transversal one, especially at the downstream sections.

The determination coefficients, obtained by comparing the simulated and measured values at the points shown in Figure 5 are: 0.793, 0.939 and 0.666 for the streamwise component; 0.875 (correlation performed on only six values over 10, since in the other points the measure was not possible), 0.846 and 0.172 for the transversal one. The values are provided for the right section (y = 0.06 m), axis (y = 0.12 m) and left section (y = 0.12 m) respectively.

The correlation between measured and simulated velocity appears very low for the left side, especially for the transversal component, probably due to some 3D effects which are not well reproduced by the 2D model. In fact, the presence of rectangular obstacles causes an abrupt deviation of the flow and a recirculation downstream from the obstacles, introducing vertical variations of the velocity which cannot be simulated by a two-dimensional model. In addition, local turbulence and dissipation tend to reduce the velocity: such a reduction is, in percentage, more significant for the transversal velocity than for the streamwise component, which presents higher values. The differences among the numerical and experimental transversal velocity are thus increased.

Overall, the hydraulic simulation shows that the flow is in general well reproduced, although some inaccuracies are observed for the velocity field on the left side. Such mismatching should be considered when analysing the results of the coupled solution between the flow and the floating bodies model.

Wood transport simulation

The simulation of wooden cylinders transport is performed with ORSA2D_WT, with added mass, drag and side coefficients set at the previously mentioned values. The initial conditions for cylindrical samples (position, orientation, linear and angular velocity) are estimated from the experimental ones at 0.5 s from the first observation of the cylinder, in order to avoid the bias due to the effect of the cylinder inflow and the following transition from totally submerged conditions to flotation. The flow conditions are those reported in the previous paragraph.

As a first step, the added inertia coefficient, CAI, is calibrated. Three values are tested, showing the comparison of trajectories and angles in Figure 6 for one cylinder released in the flume with two side obstacles. By increasing the added inertia coefficient, the trajectory tends to the left side of the flume and the angle is more similar to the experimental data. The largest differences for trajectory and orientation are observed in the downstream part of the flume where, for the lower values of CAI, a strong increase of the angle is observed. It is worth highlighting that in this part turbulence is most relevant, and 3D effects are encountered as observed from velocity measurements. Figure 6 shows also the results, in term of displacement, rotation and angular velocity, obtained with the formulation presented in Persi et al. (2018a), based on the expression by Mandø & Rosendahl (2010). Although in this case the trajectory is slightly nearer to the experimental data, the orientation is totally missed, as well as the increment of angular velocity (Figure 6(c)).

Figure 6

(a) Comparison of the measured and simulated trajectory for one cylinder; (b) simulated and measured angle versus time for the same cylinder; (c) simulated and experimental angular velocity.

Figure 6

(a) Comparison of the measured and simulated trajectory for one cylinder; (b) simulated and measured angle versus time for the same cylinder; (c) simulated and experimental angular velocity.

Close modal

In order to evaluate the effect of the variation of the added inertia coefficient on x, y and ϑ, confidence intervals of width ±5% of the range of each variable (which is 3.25 m in x, 0.24 m in y and 360° for ϑ) are considered. The percentage of observations which are included in the confidence interval is used as a measure of the accuracy of the simulation. Such an approach is preferred to the computation of the correlation between experiments and simulation, since a more deterministic comparison would not be in harmony with the uncertainty in the experimental data and the inherent randomness of the phenomenon. Table 2 reports the percentages of occurrence within the interval for eight simulations, showing that the added inertia coefficient does not greatly influence the computation of cylinder trajectories, while the effects on rotation are well visible. It appears that the higher the value of the added inertia, the higher the accuracy of computed orientation. For the maximum value tested, CAI=4.0, almost 80% of the simulated orientations are inside the confidence interval, although a slight reduction in the accuracy in y direction is observed. However, such a value produces a strong delay in the adaptation of the cylinder behaviour to the flow conditions. As an example, consider the variation of the body angular velocity in Figure 6(c): smooth variations of angular velocity are well simulated independently from the added inertia coefficient, while the abrupt variation which starts around t = 3 s in the experiments is delayed and its peak is reduced, in particular for the highest CAI. The intermediate value of 1.8 is thus preferred, since it guarantees a value of the cylinder time response closer to the experimental observation.

Table 2

Comparison of the percentage of data in the confidence interval for different values of the added inertia coefficient and for the previous formulation

CAIx (%)y (%)ϑ (%)
0.0 92 70 65 
1.8 92 70 77 
4.0 92 69 79 
Persi et al. (2018a)  93 71 58 
CAIx (%)y (%)ϑ (%)
0.0 92 70 65 
1.8 92 70 77 
4.0 92 69 79 
Persi et al. (2018a)  93 71 58 

The percentages obtained with the torque plus added inertia formulation are compared with those obtained with the formulation in Persi et al. (2018a): the latter presents slightly higher accuracy for linear displacements, but much lower data are included in the confidence interval for cylinders orientation. Such results are in agreement with the trends shown in Figure 6 and confirm that the new formulation is more suitable for the simulation of floating cylinders transport.

ORSA2D_WT is then applied to the entire set of experiments (16 tests) performed in the two-side obstacles flume at the University of Zaragoza and the results in terms of linear and angular displacements and velocities are shown in Figure 7, together with the corresponding experimental results. The initial conditions (position and rotation angle) of the experiments are extracted both visually or automatically with the software TRACKER, while the initial linear and angular velocities are obtained by averaging the rate of displacement in the range 0.5–1.5 s.

Figure 7

(a) Experimental displacement (top) and velocities (bottom), obtained from videos analysis; (b) simulated displacement (top) and velocities (bottom), obtained with ORSA2D_WT: (c) comparison of time averages in x, y and ϑ.

Figure 7

(a) Experimental displacement (top) and velocities (bottom), obtained from videos analysis; (b) simulated displacement (top) and velocities (bottom), obtained with ORSA2D_WT: (c) comparison of time averages in x, y and ϑ.

Close modal

The experimental results appear more variable and the velocities present abrupt oscillations due to the acquisition process and to flow turbulence. However, the trends are clear and are well reproduced by ORSA2D_WT. Streamwise displacements (x) are fairly grouped with respect to the experimental ones being slightly underestimated: at the maximum considered time (4 s) the maximum simulated travelled distance is 10 cm lower than the experimental one. This is shown also by the streamwise velocity, which strongly diminishes around t = 3 s, slowing down the logs translation. As regards the transversal displacement (y), wider oscillations from the left to the right bank are observed in the numerical results, with transversal simulated velocities showing higher minimum and maximum values than the experimental ones. The cylinder orientation presents the largest differences, with experimental data being more variable than the numerical ones. The logs are initially perpendicular to the flow, both in the experiments and in the simulation. Then, experimental angles tend to align with the flow, dividing into two groups around 5 rad (seven cylinders, average angle 293 ± 17°) and below 2 rad (nine cylinders, average angle 74 ± 35°), while numerical orientations are gradually distributed in the range 1.2–5.7 rad, not showing a particular tendency to align with the current. Overall, the correlation among time averages (Figure 7(c)) computed from the experiments and from the simulations is 0.999 in x, 0.916 in y and 0.830 in ϑ. The single log behaviour is however not perfectly replicated: in some case the trajectories are similar and the angle is missed, while in other cases the opposite condition occurs. An example of these two conditions, together with a fairly well reproduced test, is shown in Figure 8. However, the comparison of numerical and simulated results shows that the original formulation proposed to compute cylinder rotation provides good results in the upstream part of the flume, where the flow velocity variation is smooth, while largest differences are observed downstream of the second obstacle.

Figure 8

Contour map of the flow field with experimental (black line) and numerical (blue line) log axis. The lines show the position and orientation of the cylinders at increasing times for tests (a) number 3; (b) number 9: (c) number 16. Please refer to the online version of this paper to see this figure in colour: http://dx.doi.org/10.2166/hydro.2018.085.

Figure 8

Contour map of the flow field with experimental (black line) and numerical (blue line) log axis. The lines show the position and orientation of the cylinders at increasing times for tests (a) number 3; (b) number 9: (c) number 16. Please refer to the online version of this paper to see this figure in colour: http://dx.doi.org/10.2166/hydro.2018.085.

Close modal

Numerical results are compared with the experimental ones in Figure 9, where the confidence intervals for each variable are also shown. The data in the confidence intervals are 98% in x, 58% in y and 44% in ϑ, respectively. Lower percentages are obtained with respect to Table 2, mainly due to the reduced accuracy in the assignment of the initial conditions, which have a large influence on the numerical results.

Figure 9

Comparison of measured and experimental data. (a) Streamwise displacement; (b) transversal displacement; (c) angle.

Figure 9

Comparison of measured and experimental data. (a) Streamwise displacement; (b) transversal displacement; (c) angle.

Close modal

Effect of the initial conditions

The assignment of the initial conditions plays a major role in the outcome of the simulation. The analysis of the orthorectified videos may introduce random errors during the extraction of data, due to the difficulties in evaluating the exact water level (which is deduced by locating the water surface profile along the flume walls) for orthorectification. A maximum error of 0.005 m in the log positioning has been estimated, especially for the transversal positioning. Higher errors in angular displacement and velocity may be introduced, due to the uncertainties in maintaining the correct alignment of the two ends of the cylinders and to the video resolution.

In order to analyse how the initial conditions affect the simulation, sensitivity analysis is performed, by varying the initial values of 10% of their initial range. The range, computed as the difference between the maximum and minimum values observed for the 16 considered experiments, is shown in Table 3.

Table 3

Initial range for each considered variable and its variation adopted for the sensitivity analysis

Initial conditionsUnitInitial rangeVariation
y 0.041 ±0.004 
u m s–1 0.072 ±0.007 
v m s–1 0.040 ±0.004 
ϑ rad 1.404 ±0.140 
ω rad s–1 0.907 ±0.091 
Initial conditionsUnitInitial rangeVariation
y 0.041 ±0.004 
u m s–1 0.072 ±0.007 
v m s–1 0.040 ±0.004 
ϑ rad 1.404 ±0.140 
ω rad s–1 0.907 ±0.091 

In Table 4 the correlation coefficients among the simulation performed with the original initial conditions and those with the modified ones are reported. The orientation is the most affected by the variation of the initial conditions, with correlation coefficients that diminish up to 0.367, while the x position is the most stable variable. The variation of the transversal position does not significantly affect the simulation outcome while the largest differences are observed by decreasing the angular velocity, with a strong variation of both the transversal coordinate and angular correlation coefficient.

Table 4

Correlation coefficients computed for different initial conditions for the simulation of the transport of one cylinder

y+y–u+u–v+v–ϑ+ϑ–ω+ω
Rx2 1.000 1.000 1.000 1.000 1.000 1.000 0.999 0.999 1.000 0.997 
Ry2 0.992 0.993 0.992 0.966 0.968 0.974 0.966 0.947 0.987 0.738 
Rϑ2 0.974 0.926 0.930 0.849 0.876 0.937 0.874 0.819 0.904 0.367 
y+y–u+u–v+v–ϑ+ϑ–ω+ω
Rx2 1.000 1.000 1.000 1.000 1.000 1.000 0.999 0.999 1.000 0.997 
Ry2 0.992 0.993 0.992 0.966 0.968 0.974 0.966 0.947 0.987 0.738 
Rϑ2 0.974 0.926 0.930 0.849 0.876 0.937 0.874 0.819 0.904 0.367 

Figure 10 shows the effect of the initial angular velocity variation on the trajectory and orientation of a cylinder. The trajectory is not much affected, and similar trends can be observed for the simulation with standard values and for the two simulations with varied angular velocity. The orientation does not change much, and is very near to the experimental values, up to t = 3.5 s. Then, the effect of the angular velocity is clear: if it is increased, the angle tends to increase with respect to the simulation with standard values. The opposite is observed for a reduction of the initial angular velocity. The final difference between the two angles with modified initial condition is about 1.48 rad, which is considerably higher than the initial distribution of angles observed during the experimental campaign.

Figure 10

(a) Comparison of the experimental (solid square) and simulated trajectories. (b) Comparison of the experimental and numerical angles as a function of time. Solid lines show the results obtained with different values of the angular velocity for log 10.

Figure 10

(a) Comparison of the experimental (solid square) and simulated trajectories. (b) Comparison of the experimental and numerical angles as a function of time. Solid lines show the results obtained with different values of the angular velocity for log 10.

Close modal

The proposed dynamic DEM is one-way coupled with the Eulerian solution of the SWE for the computation of the 2D displacement and of the planar rotation of wooden cylinders floating on the water surface. The code, named ORSA2D_WT, includes the computation of drag force, side force, added mass force and pressure gradient force to estimate wood translation, while the rotation of cylindrical bodies is calculated from the hydrodynamic force balance on the horizontal plane (i.e. the water surface) plus an added inertia term which represents a resistance to rotation due to differences in flow and body angular acceleration. The added inertia coefficient is calibrated against ad hoc laboratory experiments, showing that a constant value can be implemented to obtain acceptable agreement among experimental and numerical data. The original formulation is then applied to 16 tests performed in a prismatic flume with two rectangular side obstacles. The numerical and experimental trajectories are comparable, as well as the time evolution of the average displacements, while the final orientations of numerical logs do not fit the expected range.

Moving the attention to the simulation of specific experiments, some errors in the trajectories, and in particular in the final angle of the cylinders, are observed. Such mismatch is mainly attributable to the dependency on the initial conditions, whose determination from the experiments is not always free from uncertainties. These errors appear to depend mostly on the initial values of the orientation and of the angular velocity, which strongly influence the results with a major effect on cylinders orientation. Free surface roughness and turbulence introduce additional randomness, which increases the dispersion of the results, especially with reference to the orientation, which cannot be caught by a deterministic model as ORSA2D_WT.

Note that the effect of the flow velocity may also contribute to increase the differences between the experiments and the simulation, especially in the downstream part of the flume where the hydraulic simulation is slightly inexact.

For a wider applicability, the calibration should be performed on different flow conditions, to verify the value of the added inertia coefficient and to assess the general accuracy of the resistance term.

The considerable differences observed among the simulations with standard values and those with varied initial conditions, especially for the cylinder orientation, make the initial values an extremely important parameter for future application of the code. They not only require attention in the estimation of possible range but, most significantly, demand a scenario-based approach in the numerical modelling, in order to investigate the entire set of conditions.

Overall, ORSA2D_WT gives satisfactory results for the simulation of floating cylinders transport, resulting in a promising tool for the numerical modelling of wood transport in real rivers. A preliminary application of ORSA2D_WT to a real-scale experiment on an alpine river (Persi et al. 2018c) has given promising results confirming that the one-way coupling between the DEM and the SWE solver is the correct approach when focusing on wood transport of single elements. On the other hand, the simulation of obstructions, which involves a large amount of wooden debris, necessitates the modelling of body interactions and raft formation and is thus out of the purposes of this contribution, and would possibly require the application of a two-way coupled model or, eventually, a complete 3D analysis.

Abbe
T. B.
&
Montgomery
D. R.
1996
Large woody debris jams, channel hydraulics and habitat formation in large rivers
.
Regul. Rivers
12
(
2–3
),
201
221
.
Alonso
C. V.
2004
Transport mechanics of stream-borne logs
. In:
Riparian Vegetation and Fluvial Geomorphology
(
Bennet
S. J.
&
Simon
A.
, eds).
Wiley
,
Hoboken, New Jersey
, pp.
59
69
.
Amicarelli
A.
,
Albano
R.
,
Mirauda
D.
,
Agate
G.
,
Sole
A.
&
Guandalini
R.
2015
A smoothed particle hydrodynamics model for 3D solid body transport in free surface flows
.
Comput. Fluids
116
,
205
228
.
Bradley
J. B.
,
Richards
D. L.
&
Bahner
C. C.
2005
Debris Control Structures: Evaluation and Countermeasures
.
US Department of Transportation, Federal Highway Administration, Washington, DC
.
Brown
D.
&
Christian
W.
2011
Simulating what you see: combining computer modeling with video analysis
. In
Proceedings of MPTL16
,
15–17 September
,
Ljubljana
.
Canelas
R. B.
,
Domínguez
J. M.
,
Crespo
A. J. C.
,
Gómez-Ggesteira
M.
&
Ferreira
R. M. L.
2017
Resolved simulation of a granular-fluid flow with a coupled SH-DCDEM model
.
J. Hydr. Eng.
143
(
9
),
06017012
.
Chow
C. Y.
1979
An Introduction to Computational Fluid Mechanics
.
Seminole Publishing Company
,
Boulder
,
Colorado
.
Comiti
F.
,
Mao
L.
,
Preciso
E.
,
Picco
L.
,
Marchi
L.
&
Borg
M.
2008
Large wood and flash floods: evidence from the 2007 event in the Davča basin (Slovenia)
.
WIT Trans. Eng. Sci.
60
,
173
182
.
Comiti
F.
,
Agostino
V. D.
,
Moser
M.
,
Lenzi
M. A.
,
Bettella
F.
,
Agnese
A. D.
,
Rigon
E.
,
Gius
S.
&
Mazzorana
B.
2012
Preventing wood-related hazards in mountain basins: from wood load estimation to designing retention structures
. In:
Proceedings of the 12th Congress INTERPRAEVENT
,
Grenoble, France
, pp.
23
26
.
Corrsin
S.
&
Lumley
J.
1956
On the equation of motion for a particle in turbulent fluid
.
Appl. Sci. Res.
6
,
114
116
.
Costabile
P.
&
Macchione
F.
2015
Enhancing river model set-up for 2-D dynamic flood modelling
.
Environ. Model. Softw.
67
,
89
107
.
Dean
R. G.
&
Dalrymple
R. A.
1991
Water Wave Mechanics for Engineers and Scientists
.
Advanced Series on Ocean Engineering. World Scientific Pub.
,
Singapore
.
Denk
M.
,
Rimböck
A.
&
Wendeler
C.
2008
Schwemmholz-,geschiebe- und murgangrückhalt mit flexiblen ringentzsperren (Driftwood retention, bedload retention, and debris flow retention with flexible ring net barriers)
.
Wasser Energie Luft
100
(
4
),
317
320
(in German)
.
Fotovatikhah
F.
,
Herrera
M.
,
Shamshirband
S.
,
Chau
K. W.
,
Faizollahzadeh Ardabili
S.
&
Md. Piran
J.
2018
Survey of computational intelligence as basis to big flood management: challenges, research directions and future work
.
Eng. Appl. Comput. Fluid Mech.
12
(
1
),
411
437
.
Gippel
C. J.
,
O'Neill
I. C.
,
Finlayson
B. L.
&
Schnatz
I. N. G. O.
1996
Hydraulic guidelines for the re-introduction and management of large woody debris in lowland rivers
.
Regul. Rivers Res. Manage.
12
(
23
),
223
236
.
Hirsch
C.
2007
Numerical Computation of Internal and External Flows
.
Butterworth-Heinmann
,
Oxford
.
Hoang
M. C.
,
Laneville
A.
&
Légeron
F.
2015
Experimental study on aerodynamic coefficients of yawed cylinders
.
J. Fluids Struct.
54
,
597
611
.
Macchione
F.
,
Costabile
P.
,
Costanzo
C.
,
Lorenzo
G. D.
&
Razdar
B.
2016
Dam breach modelling: influence on downstream water levels and a proposal of a physically based module for flood propagation software
.
J. Hydroinform.
18
(
4
),
615
633
.
Magnaudet
J.
&
Eames
I.
2000
The motion of high-reynolds-number bubbles in inhomogeneous flows
.
Annu. Rev. Fluid Mech.
32
,
659
708
.
Mandø
M.
&
Rosendahl
L.
2010
On the motion of non-spherical particles at high Reynolds number
.
Powder Technol.
202
(
1–3
),
1
13
.
Maxey
M. R.
&
Riley
J. J.
1983
Equation of motion for a small rigid sphere in a nonuniform flow
.
Phys. Fluids
26
(
4
),
883
889
.
Mazzorana
B.
,
Hübl
J.
,
Zischg
A.
&
Largiader
A.
2011
Modelling woody material transport and deposition in alpine rivers
.
Nat. Hazards
56
(
2
),
425
449
.
Morales-Hernández
M.
,
García-Navarro
P.
,
Burguete
J.
&
Brufau
P.
2013
A conservative strategy to couple 1D and 2D models for shallow water flow simulation
.
Comput. Fluids
81
,
26
44
.
Persi
E.
,
Petaccia
G.
&
Sibilla
S.
2018a
Large wood transport modelling by a coupled Eulerian–Lagrangian approach
.
Nat. Hazards
91
(
1
),
59
74
.
Persi
E.
,
Petaccia
G.
,
Fenocchi
A.
,
Manenti
S.
,
Ghilardi
P.
&
Sibilla
S.
2018b
(Submitted)
Hydrodynamic coefficients of yawed cylinders in open-channel flow
.
Flow Meas. Instrum.
Persi
E.
,
Petaccia
G.
,
Sibilla
S.
,
Lucía
A.
,
Andreoli
A.
&
Comiti
F.
2018c
Modelling the displacement of large wood in the Rienz River
. In:
Proc. of the 5th IAHR Europe Congress
,
12–14 June
,
Trento, Italy
, pp.
637
638
.
Petaccia
G.
,
Soares-Frazão
S.
,
Savi
F.
,
Natale
L.
&
Zech
Y.
2010
Simplified versus detailed two-dimensional approaches to transient flow modeling in urban areas
.
J. Hydr. Eng. ASCE
136
(
4
),
262
266
.
Petaccia
G.
,
Leporati
F.
&
Torti
E.
2016
OpenMP and CUDA simulations of Sella Zerbino Dam break on unstructured grids
.
Comput. Geosci.
20
(
5
),
1123
1132
.
Petaccia
G.
,
Persi
E.
,
Sibilla
S.
,
Brufau
P.
&
García-Navarro
P.
2018
Enhanced one-way coupled SWE-DE model for floating body transport
.
Ital. J. Eng. Geol. Environ.
1
.
doi: 10.4408/IJEGE.2018-01-S-15
.
Picco
L.
,
Bertoldi
W.
&
Comiti
F.
2017
Dynamics and ecology of wood in world rivers
.
Geomorphology
279
,
10
11
.
Prakash
M.
,
Rothauge
K.
&
Cleary
P. W.
2014
Modelling the impact of dam failure scenarios on flood inundation using SPH
.
Appl. Math. Model.
38
(
23
),
5515
5534
.
Ruiz-Villanueva
V.
,
Bodeque
J. M.
,
Díez-Herrero
A.
,
Eguibar
M. A.
&
Pardo-Igúzquiza
E.
2013
Reconstruction of a flash flood with large wood transport and its influence on hazard patterns in an ungauged mountain basin
.
Hydrol. Process.
27
,
3424
3437
.
Ruiz-Villanueva
V.
,
Bladé
E.
,
Sánchez-Juny
M.
,
Marti-Cardona
B.
,
Díez-Herrero
A.
&
Bodoque
J. M.
2014
Two-dimensional numerical modeling of wood transport
.
J. Hydroinform.
16
(
5
),
1077
1096
.
Schalko
I.
,
Schmocker
L.
,
Weitbrech
V.
&
Boes
R. M.
2018
Backwater rise due to large wood accumulations
.
J. Hydr. Eng.
144
(
9
),
04018056
.
Schmocker
L.
&
Weitbrecht
V.
2013
Driftwood: risk analysis and engineering measures
.
J. Hydr. Eng. ASCE
139
(
7
),
683
695
.
Silva
M.
,
Costa
S.
,
Canelas
R. B.
,
Pinheiro
A. N.
&
Cardoso
A. H.
2016
Experimental and numerical study of Slit-Check dams
.
Int. J. Sustain. Dev. Plan.
11
(
2
),
107
118
.
Solenthaler
B.
,
Bucher
P.
,
Chentanez
N.
,
Müller
M.
,
Gross
M.
2011
SPH based shallow water simulation
. In:
Proceedings of the 8th Workshop on Virtual Reality Interactions and Physical Simulations
(
Bender
J.
,
Erleben
K.
&
Galin
E.
, eds).
Eurographics Association, Lyon
,
France
, pp.
39
46
.
Stockstill
R. L.
,
Daly
S. F.
&
Hopkins
M. A.
2009
Modeling floating objects at river structures
.
J. Hydr. Eng. ASCE
135
(
5
),
403
414
.
Tagawa
Y.
,
van der Molen
J.
,
van Wijngaarden
L.
&
Sun
C.
2013
Wall forces on a sphere in a rotating liquid-filled cylinder
.
Phys. Fluids
25
(
6
),
063302
.
Truscott
T. T.
&
Techet
A. H.
2009
Water entry of spinning spheres
.
J. Fluid Mech.
625
,
135
165
.
Uchiogi
T.
,
Shima
J.
,
Tajima
H.
&
Ishikawa
Y.
1996
Design methods for wood-debris entrapment
. In:
Proceedings of INTERPRAEVENT Conference
.
Garmisch-Partenkirchen
,
Germany
, pp.
279
288
.
Yin
C.
,
Rosendahl
L.
,
Knudsen
K. S.
&
Sørensen
H.
2003
Modelling the motion of cylindrical particles in a nonuniform flow
.
Chem. Eng. Sci.
58
(
15
),
3489
3498
.