Coupling multibody system and granular dynamics application to a 2D benchmark

. We present a benchmark to experimentally validate the method we have developed to simulate mechanical devices interacting with granular media. Consequently the method we develop is able to solve problems involving bilateral and unilateral constraints. The mechanical device is modelled with the multibody system approach and is described by the use of relative coordinates. With these coordinates, the bilateral constraints are mainly due to loop-closure conditions (or speciﬁc bilateral constraints such as imposed screw-joint motion). The bilateral constraints are solved thanks to the coordinate partitioning method. The problem of contact dynamics, which introduces the unilateral constraints, is solved by a non-linear Gauss-Seidel procedure applied in the contacts coordinates. The procedure accounts for the constrained motion of the mechanism because the coordinate partitioning is also applied to the mapping between the contacts coordinates and the generalized coordinates. Consequently the solution of the unilateral constraints problem is locally compatible with the set of bilateral constraints coming from the mechanical system. The proposed planar benchmark consists of a modiﬁed four-bar mechanism blending a set of disks contained in a vibrating box. The numerical results reveals the inﬂuence of the vibration frequency on the granular compactness and thus the mechanism motion.


Introduction
Our research project deals with the tamping process of railway tracks. This process is part of the maintenance operation done on ballasted railway tracks. This particular operation restores the track geometry after some degradation. During the process the ballast grains are shacked at certain frequency range and then squeezed under the sleepers. Initially the ballast is a compact and rigid assembly of stones, when it is vibrated the compactness goes down which allows an easier squeezing. This application is an example of an articulated mechanism interacting with a granular media which takes advantage of the different state of behaviour of the ballast. In our work we want to simulate such process. To this purpose we suggest to couple the multibody systems dynamics (MBS) approach with the discrete element method (DEM) using non-smooth contact dynamics.
The objective of this paper is firstly to explain how we coupled MBS and DEM and secondly to present the experimental benchmark we are building to validate our methodology. The contact can be modelled either by the molecular dynamics (regularized contact law) or the nonsmooth contact dynamics (non-regularized contact law).
The contact dynamics has been chosen due to the high contact stiffness of the ballast and because it has proven effective to deal with railway ballast application. Previous works have been conducted on the coupling between MBS and non-smooth contact dynamic. For example Flores et al. [1] extend the approach of Jean [2] to unconstrained MBS. In Ref [3], Brüls and al. present the non-smooth generalized-α scheme which enables to enforce the unilateral constraints both at position and velocity levels for frictionless contact. This approach allows to solve the drift-off effect which is observed when the constraints are solved at velocity level only. The particularity of our work is the capacity of modelling granular media involving a lot of contacts and the ability of solving three-dimensional problem thanks to the NLGS procedure.
To validate our methodology, we are developing a planar benchmark that can be experimentally constructed and analyzed. The article will firstly describe this benchmark. Then the benchmark will be used to present the equation of motion of the MBS that will be coupled with the DEM equations of motion and contact equations. We will quickly describe the procedure used to solve the contact dynamics problem. Finally we will introduce some of our numerical results that were obtained in order to design the experimental benchmark. In near future, the benchmark will be build to compare experimental data's with numerical results.

Benchmark
The goal of the experimental benchmark is the consolidation of our development by experimental measurements. On one hand the benchmark has to represent a typical MBS that means containing at least one kinematic loop. On the other hand it has to be representative of the specific behaviour of granular media such as the fluidization or the jamming phenomena. To this purpose we suggest a benchmark (figure 1), composed of three parts: 1. The mechanical device is a four-bar mechanism i.e. a system composed of a floating rod connected to an input and output crank. Some compliance is introduced by using a telescopic bar instead of one arm. The position of the telescopic bar in the system (at the input, middle or output crank) is a design parameter to be investigated. The movement of the device is induced by a motor linked to the input crank. The last crank of the device is connected to a blender that will interact with the granular media.
2. The granular media is composed of disks of different radii. The disks are located in a box, with transparent lateral panels. This allows us to track the motion of the disks and observe the force network if they are made of a photoelastic material. The granular media is excited by the movement of the blender and also by a vertical movement of the box.
3. The box, containing the disks, can move vertically in order to introduce a vibration in the granular media. The upper part of the box is open in order to apply a confinement load on the disks.
Depending on the length of the different bars and on the position of the telescopic bar we may have a simple oscillation of the blender (as illustrated in figure 1) or a more complex movement including a full rotation. The motion of the mechanism, such as the compression of the telescopic bar or the movement of output crank, will be impacted by the compactness of the granular media. This compactness will be affected by the vibration generated by the box.

Equations of motion
Modelling this benchmark requires to couple two different techniques. On one hand the four-bar mechanism is well On the other hand the DEM is well suited to describe the granular dynamics. It is the contact interaction that causes the coupling between the two systems.

Multibody system
The MBS formalism is used to model the mechanical device of the experiment, i.e. the four-bar mechanism. The kinematics of the system is characterized thanks to the relative coordinates (q mbs ). This means that the position and orientation of a body is defined with respect to the previous body, the first body of the chain being the ground. The relative coordinates system implies that the assembly of tree-like mechanisms at the articulations is automatically satisfied. When a mechanism contains kinematic loops the generalized coordinates are completed by the following equations of constraints: = J (q mbs )q mbs h (q mbs ,q mbs ,q mbs ) = J (q mbs )q mbs +J (q mbs ,q mbs )q mbs (1) with J (q mbs ) = ∂h(q mbs ) ∂q T mbs . The kinematics of our four-bar mechanism is fully described by a set of 4 coordinates (figure 2) completed by a set of two constraints. In this example the set of constraints (h (q mbs )) imposes that the point D on the output bar is located at the position of the fixed point E on the ground. Such a constraint will only allow a rotation of the output bar about the point E. The dynamic equations of constrained MBS are: with λ being the Lagrange multipliers representing the forces and torques due to the constraints, Q is the joint forces, M mbs is the mass matrix of the system and c is the non linear dynamical vector that includes the effect of external forces frc and torques trq. The set of equations composed by equations 2 and 1 form a set of differential algebraic equation. This set of equations is reduced to a set of ordinary differential equation thanks to the coordinate partitioning method. Full details about the method are given in [4]. The main idea is to split the set of coordinates in two subsets. The first one contains the independent coordinates (q u ) that will be time integrated while the second  set contains the dependent coordinates (q v ) that will be determined thanks to the equation of constraints. The equation of motion for independant coordinates is then written in the reduced form:

Discrete Element Method
The DEM consists in analysing the granular media behaviour by computing the evolution of each particle of the media and its interaction with other particles. The position of each particles is described using the absolute coordinates (i.e. with respect to the ground fixed frame). Expressing and integrating the angular velocities in each body fixed frame leads to a constant mass matrix M: with F dem being the non linear dynamical vector including the external forces and torques. This equation must be completed to account for the interaction laws between grains and with the MBS, as explained in the next section.

Contact dynamics
To model the contact between bodies we choose the approach of non-smooth contact law which prevent the bodies from interpenetration. The equation in the normal direction is: if g ≤ 0 then Vn ≥ 0, Rn ≥ 0, Vn Rn = 0 where g is the gap, Vn and Rn are the normal component of the velocities and reaction force between the bodies. In the tangential direction a coulomb friction law is used: where μ is the friction coefficient, Vˆt and Rˆt are the tangential vector of the velocities and reaction force between the bodies. The equations 5 and 6 illustrated in figure 3 are referred as the "Signorini-Coulomb" condition. We use the Newton impact law V + n = −eV − n , to describe the dynamics of contact. In this equation e ∈ [0, 1] is the coefficient of restitution and the superscript − and + indicate the pre-and post-impact velocity. This law introduces discontinuities in the velocity field of the bodies with consequence that the equations of motion are no more integrable using standard methods based on the evaluation of the accelerations.

Coupled equations of motion
The coupling between MBS and DEM sets of equations is achieved by the contact conditions. This is taken into account by adding the unilateral constraints in the equations of motion of both system (eq. 3 & 4). These coupled equations: completed by the contacts equations (5 and 6) are written in terms of differential measures of the velocities dq u and dq dem . dI U are the differential measures of the contact reactions over the time step. The matrices H rmbs and H dem are used to project the quantities expressed in the contacts frame to the coordinate system of each domain: V = H Tq where H is the concatenation of H rmbs with H dem ,q the concatenation ofq u withq dem and V the velocities in the contacts (concatenation of Vn with Vˆt for each contact).

Resolution of the dynamics
We applied to our coupled equations of motion a time discretization which is close to the time stepping scheme of the Non-Smooth Contact Dynamics submitted by Moreau and Jean [2]. In this scheme the problem of the contact dynamics is written in the contact frame: where M and F groups the mass matrix, respectively the force vector, of MBS and DEM. h is the time step and θ a parameter of the time integrator. This problem is solved in agreement with the contact law (eqs. 5 and 6) with a Non-Linear Gauss-Seidel (NLGS) procedure. Full resolution of the contact dynamics is available in [2] and is no more detailed here. The NLGS procedure used to solve the contact problem enable us to deal with 3D problem.

Numerical analysis of the benchmark
We have applied the methodology described above to simulate some possible configuration of our four-bar mechanism. The results for the two configurations are summarized in the table 1. The radius of disks were randomly chosen between 5 equally spaced values between 6 and 10 mm. For physical data's (mass, inertia) the density of aluminium (2710 kg/m 3 ) and a width of 10 mm were chosen. The motor have an imposed operating speed of 15 rpm. The stiffness of the telescopic bar was chosen at 500 N/m. The amplitude of the box vertical motion is set at 4 mm. The confinement load is applied to increase the resistance of the media toward the blender motion and increase the amplitude of telescopic bar motion.
On the first configuration the complete rotation of the input crank will produce an oscillation on the output crank.  The figure 4 illustrates the length variation of the telescopic bar for different vibration frequencies of the box. We can observe at each rotation of the input crank that the solicitation of the telescopic bar is higher in absence of vibrations. Moreover a low frequency vibration (5 Hz) leads to similar results as the case without vibration. If we apply a higher vibration than 12 Hz, we have observed (for example with 25 Hz, not traced for readability reason) that the solicitation of the telescopic bar is not reduced any more.
On the second configuration a complete rotation of the input crank, which is the telescopic bar, will produce a complete rotation of the output bar. However if the length variation exceeds a certain value an oscillation will be produced at the output crank. The figure 5 illustrates the length variation of the telescopic bar for different vibration of the box. With this configuration the difference of solicitations for different vibration frequencies is not systematically observed at each turn of the input crank. On average we can still observe a higher solicitation in the absence of vibration or with a low frequency in comparison with the vibration of 12 Hz.
The solicitation of the telescopic bar is due to the resistance of the granular media. This resistance is an image of the media compactness. We expect a sudden decrease of the compactness when the vibration have induced a fluidized state of the granular media. As first approximation it has been considered that this state was reached when the maximal acceleration due to the box vibration matches the gravity acceleration. This corresponds to a frequency of 8 Hz with the amplitude of 4 mm. We can observe, on figures 4 and 5, that the solicitation of the telescopic bar is reduced when the frequency is over 8 Hz in comparison of a lower frequency. Moreover once the state of fluidization is reached, the solicitation does not decrease any more when the vibration frequency is increased.

Conclusion
We developed and applied a new methodology to solve problems involving constrained mechanical devices interacting with granular media. We modelled the fluidization of a granular media and observe the impact on the mechanical device. The results of the simulation were coherent with our expectations. We are now designing the experimental version of the benchmark in order to compare in a near future experimental data with numerical results.
In the future, we will challenge the time-stepping method for a better consideration of the non-linearity in the mechanical devices. A possible approach is to use a specific integrator for the smooth motion and then solving the non smooth part at each step as proposed in [5] and [3].