Modelling of transport and collisions between rigid bodies to simulate the jam formation in urban flows

This study deals with the simulation of transport and interaction between bodies considered as a rectangular shape particles, in urban flow. We used an hydrodynamic two-dimensional finite elements model coupled to the particles model based on Maxey-Riley equations, and taking into account of contact between bodies. The finite element discretization is based on the velocity field richer than pressure field, and the particles displacements are computed by using a rigid body motion method. A collision strategy is also developed to handle cases in which bodies touch.


INTRODUCTION
Many observations of damages caused by flash floods or dam-break floods indicate that floating debris and debris jams are a substantial additional risk during flooding.The jams result from accumulations of materials (Brushes, uprooted trees, carried cars. . .), which temporarily causes a water reserve to the upstream.When these jams yield or breakup, the rupture causes the violent discharge of a great quantity of water and floating bodies.With the downstream, the characteristics of the flood change then brutally (acceleration of the current, abrupt rise of water, solid transport), and this in a not easily foreseeable way.Thus, they often to the origin of the principal damages generated by the quick swellings (small water courses to strong slope, pluvial urban streaming).
The jam generally occurs in the rivers and to proximity of the bridges, and often generates a flood (Fig. 1).Thus, the motion and the orientation of transported bodies, which form this jam, is a problems of practical interest to understand how these natural barriers are formed.In particular for shallow water, the orientations for example of tree's branches in river and their collisions can form stopping between bridge piers and then perturb the flow (Fig. 2).
In the past decade, many researchers developed numerical methods for direct simulation (DNS) of fluid-particle interaction (see [1], [2], [3], [4]).These problems are fluid and solid interaction problem and employ remeshing techniques or fictitious domain methods which involve a lot of computational storage and long calculation time at each time step.The purpose of this paper is to simulate the Debris-Jam formation by floating bodies.The method used is based on the computation of the pathway of bodies of rectangular shape, by including the translation and rotation movements, without the remeshing technique.Thus, a numerical simulation of incompressible viscous flow past moving rigid bodies, when the motion of the bodies is not known in advance but results from the hydrodynamical coupling and external forces such as gravity and collisions, is performed.A simple strategy to take into account body/body and body/wall collisions is then used.The outline of this paper is as follows.The hydrodynamic finite element model is presented in section 2. The modelling of the Fluid-Rigid body interaction is presented in section 3. Section 4 is devoted to Particle-Particle collisions for rigid bodies.The computational procedure is presented in section 5 and the numerical tests are presented in section 6. Results and discussions are given in 7.

HYDRODYNAMIC FLOW MODEL
The two-dimensional hydrodynamical model was based on the resolution of shallow-water equations, which were obtained by using hydrostatic and Boussinesq approximations and by integrating Navier-Stokes equations over total water depth.Notice, that for many practical surface-water flow applications, knowledge of the full three dimensional flow behaviors is not needed and it is sufficient to use the depth average flow quantities in two perpendicular horizontal directions.Thus, the depth-averaged momentum, and continuity equations lead to the following Saint-Venant set of equations ( 1)  where u f = (u f , υ f ) the depth average velocity with u f and υ f are the components in the x and y coordinate directions, defined by: ( H (x, y, t) is the water depth, z b (x, y, t) the bed elevation and h(x, y, t) = z b (x, y, t) + H (x, y, t) is the water surface elevation, v is the effective viscosity, which includes the dispersion and the turbulence contributions.F = (Fx, Fy) integrate volume forces (Coriolis), actions exerted on the bottom (friction) or on the free surface (wind) and also forces as radiation constraints.Thus: ( where k = (0, 0, 1) is the vertical unit vector; ω is the earth velocity in rd/sec, ϕ the latitude in degrees, r is a coefficient equal to 0.0026, W = (ω x , ω y ) are wind components in x and y directions, Q a and Q are densities of air and water respectively.β depends on the value of the following coefficient γ defined by: ( Thus: (5) For each finite element, the turbulent viscosity is calculated using the Prandtl-Smagorinsky formula: (6) For each finite element, is the averaged velocities, the averaged free surface, the friction velocity and l m the mixing length defined by: (7)

FINITE ELEMENT DISCRETIZATION AND VARIATIONAL FORM
The weak form associated to the set of equations ( 1) is: (8) with: (9) are the Galerkin test functions, [r ] the viscous constraints tensor given by: (10) the edge friction tensor and S is the Neumann stress boundary.
The discretised model associated to the set of equations ( 9) is obtained by using a P 2 -P 1 element for finite element approximation [10]- [13] (Fig 3).

MODELLING OF THE FLUID-RIGID BODY INTERACTION
A rigid body (Particle) in a non-uniform time-dependent fluid flow responds to the forces imposed by fluid velocity and body forces.A delineation of the particle motion and the interaction between the particle and fluid depends on the correct calculation of the forces.Thus, the Lagrangian equation of a rigid body moving in an arbitrary 2D flow are governed by the Newton's laws, that includes in the force balance Buoyancy, pressure gradient of ambient flow, added mass and drag forces.(11) m p is the particle mass, u p = (u x , u y ) is the translation velocity vector of the particle masscenter is the moment of inertia, I p , is the moment of inertia, ω p is the angular velocity and are total forces acting on the particle.G is the position of the center of mass and M the position of the center of pressure where are concentrate hydrodynamic forces.

FORCES ACTING ON THE RIGID BODY
The essence of the present simulations is to treat the total forces imposed on the fluid by particles of non spherical form, which can be obtained from the particle equation of motion-equation (11), as a body force.Since the effect of the particle on the fluid is to generate a disturbance in the fluid velocity field, it would also be appropriate to couple the particle effect into the fluid in terms of a velocity disturbance.Since there is usually a relative motion between the particle and the fluid, and also the particle can be rotating, the liquid flow around individual particle leads to local variations in pressure and shear stress.Here fives different contributions for the interaction force term can be taken into account, a friction (Stokes drag) force , an added mass force, the pressure gradient force, the Basset history force, and the Buoyancy force, leading to the following approximation for the force of interaction, governed by the Maxey-Riley equation [12]: where:

•
, is the Stokes drag term, where C d is the drag coefficient which depends on the particle form u p and m f are the velocity and the masse of the particle respectively, whereas, u f is the fluid velocity.
• , is the pressure gradient related to fluid acceleration, • , is the virtual mass contribution (or added mass effect).It arises because acceleration of the particle requires acceleration of the fluid surrounding the body, and the pressure gradient force accounts for the acceleration of the displaced fluid (i.e. the force a fluid particle of the same size would experience in absence of the particle).
is the Basset history term.It adjusts the particle acceleration by taking into account the past acceleration on the particle motion, including the effect of the conditions that prevailed during development of the flow.In the twodimensional study, where gravity is absent, this term may be negligible, especially because the particle is not released from rest but with an initial velocity corresponding to the ambient fluid velocity.
Elsewhere, in the above equation, the Buoyancy term: F Bo as well as the pressure gradient F pg because of the 2D horizontal space problem.We also neglect neglect the Faxen correction term, which becomes significant only in the event of large curvature in the velocity profile.We notice that the derivative d/dt is used here to denote a time derivative following the moving particle, thus: where S eff is the particle area normal to the direction of the drag force.It changes with the incidence angle α between relative velocity (u f -u p ) and particle major axis direction and is determined by L and l are the length and large of the particle respectively, β is the drag coefficient which is determined by [5]: where is the equal projected area circle diameter, Re is the particle Reynolds number, defined on the basis of the relative velocity between particle and the surrounding fluids and equal-volume sphere diameter; is the particle sphericity, where s is the surface of a sphere having the same volume as the particle and S is the actual surface area of the non spherical particle.The equations (15) calculate the translation velocity vector at the center of mass.The velocity vector for an arbitrary point M on a rotating and translating rigid body is then: G is the position of the center of mass.
As the center of pressure P does not coincide with the center of mass at nonzero incidence angles, hydrodynamics forces described above, which act at the center of pressure P rather than at the center of mass G, will give rise to a torque on the particle.Thus, the distance ||d || between G and P is given by (see [5] ): (20)

PARTICLE-PARTICLE COLLISIONS FOR RIGID BODIES
The inter-particle collision depends mainly on the body motion and on it's size, and on the fluid dynamic transport effects, that is, drag and pressure gradient forces.However, the calculation of the particle linear (u p ) and angular (ω p ) velocities change caused by an interparticle collision is related to the following assumptions:

•
Bodies are assumed to be rigid, and no shape deformation of the bodies during the collision process is considered.

•
Bodies move in a two-dimensional plane.There is no linear movement and rotation in a third-dimensional space.

•
Collisions are instantaneous, and during the collision calculation, only impulsive forces are considered.Thus, in what follows and for the sake of simplicity, we will present only the case of contact between two bodies.This approach can be extended to more complex cases by treating them as series of contacts between two bodies.Suppose a vertex on particle A is colliding into an edge of particle B at the point M, (see Fig. 4) GA and GB are the centers of mass of particles A and B respectively d A , and d B are the distances vectors from GA and GB to point P for particles A and B respectively.n is the normal vector to edge of particle B. Two solutions exist to model collision between particles.The first one is to introduce an impulsion force directly in equation (15) when collision occurs.That means to know details about the materials of the particles, their exact geometry and how they deform under stress, etc.This approach is recommended when one wants to know more details on collisions but require lot of CPU time cost.The second approach consists on modifying the directions and magnitudes of velocities vectors after the impact of the particles by introducing the concept of impulse [5].An impulse is the change in momentum of an object when a force is applied over a very brief period of time.We imagine that during the collision there is a force acting for a very brief period of time.If you integrate that force over that brief time, you get the impulse.This is resumed on the following equation: are the initial and final translational velocities of particles A and B respectively.This equation says that the velocity at which the particles fly apart is proportional to the velocity with which they were coming together.The proportionality factor is the parameter e.The elasticity contact parameter e is equal to 0 if the collision is inelastic and 1 if the collision is perfectly elastic.Our goal is to identify the collisions between particles and to know approximately how stopping are formed in river for example when tree's branch contact bridge piers and perturb the flow.Because of its relative simplicity, our choice is focussed for the second method.If we neglect the friction force (which is parallel to the edge) and assume that the only force during the collision is in direction perpendicular to the edge, which is given by the vector n, and using the equation ( 21), the pre-and postvelocities are related by the following equations: (22) are respectively the initial and final angular velocities of particles A and B. is I A is the moment of inertia od the particle A, α AB is the impulse parameter and it's expression is calculated by introducing equations (21) into equation ( 22) and assuming that n.n = 1.The expression of α AB is then: We can use this last equation ( 22) for calculating collisions with a wall by assuming the mass of the wall is infinite, hence it becomes: (24) Notice that all those equations calculate velocities at center of mass.For an arbitrary point C in particle, one uses the following equation to estimate the velocities: (25)

COMPUTATIONAL PROCEDURE
We discretize the system of equations (15) using a finite element procedure, with an additional constraint given by equation ( 19) and by taking into account of inter-particles collisions.

THE TIME-INTEGRATION SCHEME
The time-integration scheme consists on using the Runge-Kutta scheme with an adaptive step size control algorithm.The goal of this adaptive time discretization step control is to achieve some predetermined accuracy in the solution with minimum CPU-time.Thus, after writing the set of equations (15) in the form: (26) the basic idea of the Runge-Kutta with adaptive step size control algorithm, may be given as follow: • Combine a 4 th order Runge-Kutta method and a 5 th order Runge-Kutta method; • Estimate error for each local step combining solution of the 4th and 5th order; • If the local truncation error is larger than a given accuracy, reduce the step size; • If the local truncation error is smaller than a given accuracy, increase the step size For the 5 th order Range-Kutta, the solution is: For the 4 th order Range-Kutta, the solution is: (28 where ∆t is the time discretisation step, y is the function to be approximated, S i are computed as regular Range-Kutta methods: (29) f is the right hand side, the parameters introduced in this method are given in the table 1, where the error reads: (30) Thus, the adaptative time discretisation step Range-Kutta algorithm is as follow: 1. Input f, t 0, t last, ∆t, y 0, tol 1and go to 1 C S as y i as i a alculate error (30) err as Modelling of transport and collisions between rigid bodies to simulate the jam formation in urban flows Geometry of the bend tube.

SEARCH ALGORITHM FOR INTERPOLATION
To get information for each particle moving through the finite element mesh and solve numerically equations ( 15) and ( 21) -( 26) at each time step, one need to interpolate data associated with unstructured grids.This interpolation can take a non-negligible portion of total CPU-time, especially for large applications.To optimize this CPU-time, we use a fast neighbor-to-neighbor technique [6] that accelerates speedup by a factor 5 compared to usual techniques.
Consider an unstructured finite element mesh (see Fig. 5), as well as a point P with coordinates (xp, yp).A straightforward way to determine if the point M is inside a given element i is to determine the shape-function values of P (Ni) with respect to the coordinates of the points belonging to element i as follow: (31) where N i are the shape functions, that verify: • for a 3 nodes triangle: • for six nodes finite element, we split each element into four triangles and evaluate each of these sub-triangles.Then the point M is in element i, if and only if: Starting with "starting element", the algorithm jump from an element to its neighbor in the known grid.If the element into which x falls can be found in a few attempts (< 10), this procedure, will outperform all interpolation algorithms.The algorithm is as follows: 1. Form the list of elements adjacent to element for the given mesh 2. Loop over all particles 3. Obtain good starting element "starting element" 4. For "starting element", evaluate Ni If ( Min (N i, 1 -N i ) ≥ 0, i = 1, 3 Exit else Set "starting element" to neighbor associated with Min(N i ) Goto 4

FAST LINE-SEGMENT INTERSECTION
For each time step, one must detect if the particles, considered as rigid bodies, intersect each other.Especially for large number of particles, this intersection's search can take a lot of total CPU-time and must be optimized to reduce computational time.A simple way to calculate intersection between two segments is to construct the convex hull of the segments.Two segments intersect if and only if the convex hull is a quadrilateral whose vertices alternate between the two segments.Nevertheless, this technique could be used only when the number of particles is small.In this section, we propose the sweep line algorithm witch is more complicated to program than the first one, but accelerates drastically speedup especially for large number of particles.Two segments AB and CD intersect (Fig. 6) if and only if : • the endpoints A and B are on opposite sides of the line CD, and • the endpoints C and D are on opposite sides of the line AB.In order to detect whether there is an intersection in a set of more than just two segments, we use the sweep line algorithm (Fig. 7).The segments which the sweep line intersect form a "working set" of line segments, amongst which we test for intersections.If the working set is ordered from left to right, then only neighbors in the working set could intersect each other.As the sweep line moves down, segments enter and leave the working set, and the working set is maintained in left-to-right order.Every time the working set changes, we search again for possible intersections.

VALIDATION OF THE HYDRODYNAMICAL MODEL
The validation of the hydrodynamical model (Eq.1) is carried out by using a bend tube as a benchmark test.The geometry used in this study was a 180 o U-bend, with a diameter B, as shown in Fig. 8, where the mesh grid used in the present 8 calculations is given in Fig. 9.This geometry has been used in several experimental 9 and numerical investigations [7], [8], [9].The computed water elevation, the velocity magnitude contours and iso-vectors are presented respectively in Fig. 10, Fig. 11 and Fig. 12.As illustrated in Fig. 10, in the external part of the curve, the pressure force is in balance with the centrifugal force, however, the  fluid particle velocity decreases.Consequently, the centrifugal force acting on it declines.The maximum velocity tends to be found closer to the inner curve (Fig. 11 and Fig. 12).It is also worthwhile to mention that for square-sectioned ducts or for non-circular ducts, the velocity vectors may give rise to secondary motion, which may be strong enough to generate small vortexes near the boundary layer.

TRANSPORT AND JAM FORMATION: RESULTS AND DISCUSSIONS
We consider the simulation of the motion of tree's branches in river and their collisions with four bridge piers.The tree's branches have rectangular shape with 1m large and 20 m long.The channel has 1500m length and 170 m large and inclined at slope 0.01 (Fig. 13).The obstacles have circle shape and their radius equal 30 m.The computational domain is discretized using triangular mesh for a total of 4468 finite elements and 9357 nodes with 15000 degrees of freedom (Fig. 14).At each node, the unknowns of the problem are the water surface elevation h(x,t) and the velocities u and v.The physical parameters for the fluid are as follows: the density, Q = 1:Kg/m 3 ; the Reynolds number, Re = 1000 and the Manning coefficient's friction= 0.03m 1/3 /s -1 .Initial solution, u = v = h = 0; boundary conditions • on AB, h = 10m, • on CD, h = 9.80m, • on AC and BD, a reflection boundary condition is imposed.
The simulation takes about 300 second CPU, on a Windows based PC with 2:8Ghz Pentium.The computed water elevation, the velocity magnitude contours are presented respectively in Fig. 15   Fig. 16 shows that complex flow around the piers reflects eddies in the hydrodynamical flow field and Fig. 17 shows the collisional effects from the particles accumulation around the piers, and that particles tend to orient with long axes parallel to flow in the center of the channel.Furthermore, Fig. 17 shows that interactions among multiple particles have tendency to both increase and decrease particle transport: particle-to-particle collisions entrain previously deposited particles, while immobile particles can obstruct moving ones, causing deposition, and formation of log jams.
Nevertheless, in order to examine interactions among hydraulics, channel geometry, transport distance and deposition of floating bodies, and to assess the potential effects of floating bodies on the dynamic morphology of streams, we need to understand where these   bodies will deposit and how long they will remain there.Furthermore, the natural formation of debris jams, is a build-up of floating bodies of variable sizes and quantities.Thus it is necessary to simulate numerically the pathway of particles of different size and shapes, while taking into account channel geometries.

CONCLUSION
In this paper, a model is successfully derived to track the motion of non spherical rigidbodies in a non-uniform flow field.The model reproduces successively the bodies pathway by including translation and rotation movements, and by taking into account of the Particle-Particle collisions.By using the fast neighbor-to-neighbor technique [6], the results showed the capacity of the model to reproduce the collisions between bodies and walls without increasing the computing time.
This model can help to identify potential problems for flood management in urban areas as well as to improve the accuracy of flow modelling of these areas.It can also be useful to researchers during investigations of the interaction between the floating debris, the flood wave propagation and the jam formation.Indeed, the information available consists mainly of after-flood observations of the remains of debris jams.

Figure 1
Figure 1 Schematic presentation of potential Debris-jam sites.

Figure 2
Figure 2 Conceptual sketch of the formation and breaching of debris-jams at potential locations (Bridges with piers).

Figure 5
Figure 5 Sketch of the search line algorithm.

Figure 6
Figure 6 Schematic presentation of the sweep line algorithm.

Figure 9
Figure 9 Mesh of the bend tube.

Figure 10
Figure 10 Water surface elevation.

Figure 11
Figure 11 Isco-contour of the velocity magnitude.

Figure 12
Figure 12 isco-vectors of the velocity.
and Fig.16.Snapshots of the positions and orientations of the particles are shown in Fig.17

Figure 13
Figure 13 Geometry of the channel.

Figure 14
Figure 14 Mesh of the channel: two piers are located in the zoomed picture.

Figure 15
Figure 15 Water surface elevation.

Figure 17
Figure 17 Particle transport passing through piers, located in the channel.Particles consist of floating long sticks.
To test whether two points are on opposite sides of a line through two other points, we use a counter clockwise (CCW) test.A and B are on opposite sides of line CD if and only if exactly one of the two triples A, C, D and B, C, D is in counter clockwise order.To test the order of the triple A, B, C, we use the simple algorithm CCW: