Calibration of a numerical model for the transport of floating wooden debris

The paper describes the calibration of a numerical model to simulate the 2D motion of floating rigid bodies. The proposed model follows a one-way coupling Eulerian-Lagrangian approach, in which the solution of the Shallow Water Equations (SWE) is combined with the Discrete Element Method (DEM) to compute the displacement of rigid bodies. Floating bodies motion is computed by adapting the Maxey-Riley equation to the case of semi-submerged bodies at high Reynolds number. In order to account for the flow velocity distribution along the body axis, the elements are divided into shorter subsections. A specific formulation is proposed to calculate the rotation of wooden cylinders, by computing the angular momentum. The model includes also a term of added inertia, which accounts for the resistance to rotation and requires the calibration of a specific inertia coefficient. A series of flume experiments is performed to calibrate the model. The 2D trajectories of floating spheres and the linear and angular displacement of cylinders are recorded in stationary conditions. The comparison between the experimental data and the simulation shows that the numerical results are in agreement with the experimental ones, although less accuracy is observed in the reproduction of the angular displacement.


Introduction
During floods, wooden elements can be transported by the flow and may influence the effects of the inundation.Some authors [1-3] have performed post-event surveys highlighting that wood transport may increase the hazardousness of the flood, especially in presence of urban structures.Some attempts of physically-based design of safety structures, built to avoid, or reduce, the backwater effect, can be found in the literature [4][5][6].In order to predict the displacement of floating bodies during floods, numerical models that consider the wooden elements as cylinders were developed [7][8][9][10][11].Until now, no detailed studies have been performed to analyse the changes in orientation of transported cylinders, even if this datum is strongly significant to estimate the interaction with in-channel structures.The aim of this paper is (i) to introduce the new rotation formulation implemented in the 2D model and (ii) to present the outcome of the simulation of floating bodies transport.

Numerical model
To compute the hydrodynamic forces exerted on the body, the hydraulic variables of the flow are needed.Flow velocity and water level are computed through a full Eulerian 2D Shallow Water Equation (SWE) model, solved by the finite volume code ORSA2D [12,13].It implements a Roe's Riemann solver, 1 st order accurate in time and space [14], applying an upwind discretization to the bottom slope source term [15], while the friction slope is evaluated in a semi-implicit way [16].The time step is evaluated according to the Courant-Friedrichs-Lewy (CFL) condition [17].The motion and rotation of the rigid bodies are computed with a Lagrangian method, by considering each body as a single element (Discrete Element Method, DEM).The DEM method follows a dynamic approach and is one-way coupled with ORSA2D.The translation equation is adapted from the Maxey Riley equation [18] and includes the drag force, the side force, the pressure gradient and added mass forces.Each force is computed in four points along the body main dimension (Fig. 1), in order to take into account the flow velocity gradient.The translation equation reads: In Eq. 1, the subscript i refers to the four segments on which the forces are computed, while subscripts b and f refer to the rigid body and the flow, respectively.V is the velocity vector, A is the area on which forces are applied, C D , C S and C AM are the drag force coefficient, the side force coefficient and the added mass coefficient, m is the mass, ρ is the density and i z is the unit vector normal to the 2D flow plane.In order to take into account both the time and space variation of the flow velocity, the flow acceleration is computed with the total derivative DV f Dt � [18].For elongated bodies, the drag and side coefficients vary with the body orientation [10].
For the computation of the angular displacement an original formulation is implemented, involving two terms: the angular momentum and the added inertia.The first one is computed by taking into account the forces on the four segments of the body and calculating the angular momentum with reference to the centre of mass (CM, Fig. 1).The second term, which is a resistance term as suggested by [19], is computed as a torque due to the added moment of inertia, based on the difference among the flow and the body angular acceleration.The equation implemented is: The suffixes are tha same as for Eq. 1.The area A i and the volume Vol i , used to compute the forces, change according to the body shape; u and v are the streamwise and cross-sectional velocity components.To calculate the angular momentum, all the previously mentioned forces are considered.Equation 2 includes the complete added mass force term, which is made up of both the flow and the body linear acceleration (Eq.3).
The forces are multiplied by the distances between the centre of mass of the body and the four points of application of the forces (b CMxi and b CMyi ).The added inertia term is computed with the total derivative of the flow angular velocity, Dω Dt � , and C AI is the added inertia coefficient, which has to be calibrated.

Experimental campaign
In order to provide reliable data for the model calibration, an experimental campaign was carried out in the laboratory of the Fluid Mechanics group of the University of Zaragoza, in a prismatic flume with horizontal bottom.The tests were performed under steady nonuniform conditions, with a discharge equal to 15.3 m 3 h -1 , and different flume configurations to provide different flow fields.Figure 2 shows the scheme of the flume for the case with two side obstacles.It also displays the flow velocity (streamwise and crosssectional component) distribution.In order to check the accuracy of the hydraulic model, the flow velocity was measured at the points shown in Fig. 2a.
Two sets of experiments are here presented for the model calibration.In the first case, one curved side obstacle was placed in the flume and the experiments were performed with a wooden sphere (density 690 kg m -3 ).The second set of experiments was performed with two rectangular side obstacles and with 8 wooden cylinders (average density 770 kg m -3 ).The cylindrical samples were released one by one by a semi-automatic device, from 0.02 m above the water surface and about 1.25 m downstream from the inlet section.The sphere was released manually at the same distance from the inlet.The experiments were recorded, orthorectified and the information was extracted using Tracker [20].For the first set of experiments, with one smooth side obstacle in the channel, the outcome of 9 replications with the wooden sphere is shown in Fig. 3.The diagram corresponds to the planar view of the channel, with the right side along the x-axis.The releasing method (which is manual for the sphere) strongly affects the final trajectory.In all cases, the planar rotation of the sphere was not significant.
Focusing on the second experimental case, i.e. cylinders released in the flume with two rectangular side obstacles, Fig. 4 shows the trajectory of the centre of mass for 8 repetitions and corresponds to the planar view of the channel, with the right side along the x-axis.Fig. 5 displays the variation of the angle in relation with the x position of the body.Note that the experiments are named after the name of the sample (from T0 to T7) followed by the video number.As regards trajectories, two main paths appear, one reaching a maximum y position around 0.20 m and the other that reaches 0.16 m (for samples T3 and T6).Overall, the cylinder trajectories tend to follow the main stream in the flume and do not spread much.
As regards orientation, the cylinders are released perpendicularly to the flow, i.e. with an angle of 180° according to the selected convention, and tend to rotate to align with it.In the first part of the flume (up to 1.95 m from the inlet), the orientation varies smoothly as the flow velocity, with the exception of T5 which probably suffered some alteration during the release and is thus influenced by the initial conditions more than by the flow.When the cylinders encounter the first obstacle, some changes are shown, but when they are near the second obstacle (placed at x = 2.50 m) they abruptly vary their orientation and tend to align with the channel axis in the final part (90° or 270° indicate alignment, according to the selected convention).The final angles are markedly grouped around 50°-100° and 280°-300°, but are not clearly correlated to the initial angle.

Numerical simulation of laboratory experiments
The simulations are performed under steady conditions on an unstructured mesh of about 16000 triangular elements.The channel sides are represented by a solid wall condition, the upstream boundary condition is constant discharge and the downstream boundary condition is the critical state.The channel roughness is modelled by the Manning coefficient, n = 0.01 s m -1/3 .

Comparison of measured and simulated flow velocity
For the second experimental case (flume with two side rectangular obstacles), the streamwise and the cross-sectional components of velocity, measured at the points shown in Fig. 2a, are compared with the simulated results in Fig. 6.
The measurements are grouped longitudinally, and some differences can be observed.There is a better agreement of the streamwise component with the expected values, especially along the channel axis.The correlation coefficients are 0.939, 0.739 and 0.666 for the axis, the right side and the left side respectively.The correlation coefficients for the cross-sectional component are 0.846, 0.875 and 0.182, showing a lower correlation for the left side of the channel.Overall, the hydraulic simulation performed with ORSA2D provides results similar to the measured ones and confirms that the geometry, the discharge and the channel roughness are correctly modelled.

Modelling floating body transport
The simulation of floating bodies transport is performed by implementing in ORSA2D the translation and rotation equations, integrated in time with a Runge-Kutta method and with the same time step as the hydraulic simulation.In Fig. 7, 5 over 9 of the simulated trajectories of the sphere in the flume with one curved side obstacle (first experimental case) are compared with the experiments.To assess the outcome of the rotation formulation, the transport of elongated bodies with two side rectangular obstacles is simulated (second experimental case).Figures 8, 9 and 10 show the correlation among the measured and simulated x position, y position and orientation, respectively.Note that the added inertia coefficient is calibrated to obtain the best correlation for the considered experiments and is set equal to 1.8.
A confidence interval is outlined for each figure (solid lines), set as ±5% of the possible range of variation, which is 2.25 m in x, 0.24 m in y and 360° for ϑ.Figs. 8, 9 and 10 show that the best correspondence among experimental and simulated displacement for the second experimental case is obtained in x, i.e. in the direction of the channel axis.The lateral displacement appears to be less accurately represented, with simulated y positions being higher, i.e. nearer to the left side of the channel, than the experimental one.As regards rotation, the simulated angles are initially close to the experimental results but tends to differ strongly from the expected values in the final part of the trajectory downstream of the obstacles.Overall, the percentage of data included in the confidence interval is 92% for x displacement, 70% for y displacement and 77% for the angle.

Conclusions
The paper presents an Eulerian-Lagrangian model for the transport of floating rigid bodies.Their linear and angular displacements are computed with a dynamic approach.The forces exerted on the body are estimated from the flow velocity and acceleration, leading to the computation of the body motion.The translation equation is adapted from the literature, while the rotation equation is obtained by combining the angular momentum due to the distribution of the forces on the body and a new term, named added inertia.This term takes into account the relative angular acceleration and the added inertia coefficient, which needs to be calibrated.The method is first applied to the case of a floating wooden sphere, showing that it can reproduce the experimental trajectories with high accuracy.Since the sphere rotation was not appreciable in the experiment, eight tests with a wooden cylinder were also simulated to calibrate the rotation formulation.The results are in agreement with the experimental data, despite the simulation of the transversal and angular displacements are less accurate.
The simulation of the trajectory of floating bodies is an initial step towards the modelling of the transport of wood during floods.Further experiments and simulations are needed to assess the effect of rotation for the case of floating spheres, to verify that the rotation formulation provides good results for different flow fields and to calibrate the added inertia coefficient under different conditions.In addition, the interaction among bodies and the simulation of the backwater effect have to be implemented to provide an effective instrument for the prediction of the effect of floating bodies transport.

Fig. 1 .
Fig. 1.Scheme of the division of the transported body.Forces are evaluated at the centre of each segment (points 1 to 4); the corresponding volume and frontal areas are highlighted by dashed lines.

Fig. 2 .
Fig. 2. Flow field for the second experimental case.(a) Sketch of the channel with two rectangular obstacles and location of the points where the velocity is measured; (b) streamwise velocity field; (c) cross-sectional velocity field.

Fig. 3 .
Fig. 3. Trajectories of a wooden for the first experimental case (one curved side obstacle).Three main trends are shown: a lower trajectory, nearer to the x-axis, which presents the mostly curved trace, then an intermediate trajectory between y = 0.16 m and y 0.20 m and, for one experiment, a trajectory very close to the left side of the flume.The latter presents a straighter trend, probably due to the effect of the boundary.The trajectories tend to spread and to occupy nearly 70% of the channel width, despite being released in the same location.The releasing method (which is manual for the sphere) strongly affects the final trajectory.In all cases, the planar rotation of the sphere was not significant.Focusing on the second experimental case, i.e. cylinders released in the flume with two rectangular side obstacles, Fig.4shows the trajectory of the centre of mass for 8 repetitions and corresponds to the planar view of the channel, with the right side along the x-axis.Fig.5displays the variation of the angle in relation with the x position of the body.Note that the experiments are named after the name of the sample (from T0 to T7) followed by the video number.

Fig. 4 .
Fig. 4. Trajectories of the centre of mass of wooden cylinders for the second experimental case.

Fig. 5 .
Fig. 5. Experimental angle evolution of wooden cylinders for the second experimental case.

Fig. 6 .
Fig. 6.Comparison of the measured and simulated flow velocity for second experimental case.For the distinction among series refer to the points in Fig. 2a.(a) Streamwise component; (b) Crosssectional component.Beware different axes limits.

Fig. 7 .
Fig. 7. Numerical trajectories (symbols) and experimental trajectories (lines) for the first experimental case.The figure shows that the model can reproduce quite well the motion of spherical objects, although only two main paths are obtained, one nearer to the obstacle and the other which reaches y = 0.20 m.The intermediate path (between y = 0.20 and y = 0.16 m) is not reproduced.The simulated trajectories end between x = 2.40 and x = 2.50 m because the sphere stops floating and starts rolling on the bottom (a mode of transport not included in the model).The average correlation coefficients among all the experimental and simulated data are 0.998 and 0.780 in x and y, respectively.To assess the outcome of the rotation formulation, the transport of elongated bodies with two side rectangular obstacles is simulated (second experimental case).Figures8, 9and 10 show the correlation among the measured and simulated x position, y position and orientation, respectively.Note that the added inertia coefficient is calibrated to obtain the best correlation for the considered experiments and is set equal to 1.8.A confidence interval is outlined for each figure (solid lines), set as ±5% of the possible range of variation, which is 2.25 m in x, 0.24 m in y and 360° for ϑ.