Mohamed Jebahi
email: mohamed.jebahi.1@ulaval.ca

Frédéric Dau

Jean-Luc Charles

Ivan Iordanoff

J.-L Charles

Multiscale modeling of complex dynamic problems: An overview and recent developments

Keywords: Multiscale modeling, dynamics, continuum methods, discrete methods, coupling methods, bridging domain decomposition

Multiscale modeling aims to solve problems at the engineering (macro) scale while considering the complexity of the microstructure with minimum cost. Generally, two scales are considered in multiscale modeling: small scale, which is designed to capture the mechanical phenomena at the atomistic, molecular or molecular cluster level, and large scale which is connected to continuous description. For each scale, well-established numerical methods have been developed over the years to handle the relevant phenomena. As a first part of this paper, the most popular numerical methods, used at different scales, as well as the coupling approaches between them are classified, according to their features and applications, so that the place of those used in multiscale modeling can be distinguished. Subsequently, the class of concurrent discrete-continuum coupling approaches, which is well adapted for dynamic studies of complex multiscale problems, is reviewed. Several techniques used in this class are also detailed. Among them, the bridging domain (BD) technique is used to develop a discrete-continuum coupling approach, adapted for dynamic simulations, between the Discrete Element Method (DEM) and the Constrained Natural Element Method (CNEM). This approach is applied to study the BD coupling parameters in dynamics. Several results giving more light on the setting of these parameters in practice are obtained.

Introduction

Numerical simulation has increasingly become a very important approach for solving complex practical problems in engineering and science. With the help of the increasing computers power, this modern numerical approach attacks the original problem in all its detail without making too many assumptions, rather than adopting the traditional theoretical practice using assumptions and approximations to simplify the studied problem. It plays a valuable role in providing validation for theories, offers insights to the experimental results and assists in the interpretation or even the discovery of new phenomena, which can be of different scales. Typically, four scales can be distinguished in the context of numerical simulation.

-The nanoscopic (or atomic) scale (∼ 10 -9 m), where phenomena related to the behavior of electrons become significant. At this scale, the interaction between particles (electrons, atoms, etc) is directly dictated by their quantum-mechanical state.

-The microscopic scale (∼ 10 -6 m), where phenomena related to the behavior of atoms are considered. The interaction between atoms is governed by empirical interatomic potentials, which are generally derived from Quantum Mechanical computations. The classic Newtonian mechanics is used to compute the atoms displacements and rotations.

-The mesoscopic scale (∼ 10 -4 m), where phenomena related to lattice defects are considered. At this scale, the atomic degrees of freedom are not explicitly treated, and only larger-scale entities (clusters of atoms, clusters of molecules, etc) are considered. The interaction between particles are also described by the classic Newtonian mechanics.

-The macroscopic scale (∼ 10 -2 m), where macroscopic phenomena which can be described by Continuum Mechanics are considered. At this scale, the studied physical systems are regarded as continua, whose the behavior is described by constitutive laws.

At each scale, several numerical methods have been developed to deal with the relevant phenomena.

These methods can be classified into: discrete methods (DMs) and continuum methods (CMs). The first class is based on discrete mechanics and typically covers the first three scales: quantum mechanical methods (QMs) for nanoscopic scale analyses, atomistic methods (AMs) for microscopic scale analyses and mesoscopic discrete methods (MDMs) for mesoscopic scale analyses. This class generally provides very accurate results, but it is very time-consuming and can only be used for small physical systems.

The class of continuum methods (CMs) is based on Continuum Mechanics and is rather adapted for problems at the macroscopic scale. Otherwise, additional handling is generally required to treat fine scale phenomena. Although less accurate than the discrete methods, this class is relatively inexpensive and can be applied on large systems.

The challenge in modern material science and engineering simulations is that real materials usually exhibit multiscale phenomena, which require on one scale a very accurate and computationally expensive description and on another scale coarse description to avoid prohibitively large computations.

Since none of the DMs and CMs methods alone is able to describe the entire system, it would be beneficial to develop coupling approaches which combine different methods at different scales, such that the computation effort can be distributed as needed. This has led to the development of various multiscale coupling approaches, able to predict the mechanical behavior of materials across all relevant scales, while keeping the accuracy of the individual approaches in their respective scales. These approaches can conceptually be classified into three categories: hierarchical approaches, concurrent approaches and hybrid hierarchical-concurrent approaches. The first category attempts to piece together a hierarchy of computational approaches, such that the large scale models use information obtained from the fine scale models. This approach can be applied for problems in which the different scales are decoupled or weakly coupled. The second category aims to couple methods at different scales in a combined model using some coupling procedure, such that the different scales communicate instantly during computation. This category is adapted for problems in which the different scales are strongly coupled, such as complex dynamic problems. The third category combines the features of the first two classes. A hierarchical modeling is used as long as the fine scale information is available, when this is no longer possible, a concurrent modeling takes place. From these categories, the concurrent coupling approaches seem more adapted for complex dynamic mechanical problems, which are the general framework of this review. Indeed, such problems generally involve highly dependent multiscale phenomena. The concurrent coupling approaches can be of different types: discrete-discrete coupling approaches, continuum-continuum coupling approaches and discrete-continuum coupling approaches. The latter present the advantage of combining dissimilar methods, which allows to obtain a good compromise between accuracy, efficiency and realistic description. This is why the concurrent discrete-continuum coupling approaches have become increasingly used in the literature [START_REF] Ben Dhia | The Arlequin method as a flexible engineering design tool[END_REF][START_REF] Guidault | On the L2 and the H1 couplings for an overlapping domain decomposition method using Lagrange multipliers[END_REF][START_REF] Guidault | Bridging domain methods for coupled atomistic-continuum models with L2 or H1 couplings[END_REF]112].

After reviewing the different classes of numerical methods, the present paper gives a birdseye view of the techniques commonly used in the concurrent coupling of discrete and continuum approaches.

From these techniques, the bridging domain (BD) is of particular interest for dynamic studies. As will be seen later, this technique can avoid spurious wave reflections without any additional filtering or damping, provided that its parameters are well adjusted. Several studies of these parameters, using 1D static models, can be found in the literature [START_REF] Bauman | On the application of the Arlequin Method to the Coupling of Particle and Continuum Models[END_REF][START_REF] Guidault | On the L2 and the H1 couplings for an overlapping domain decomposition method using Lagrange multipliers[END_REF][START_REF] Guidault | Bridging domain methods for coupled atomistic-continuum models with L2 or H1 couplings[END_REF]. Recently, these studies have been extended, in a previous paper [START_REF] Jebahi | 3D coupling approach between discrete and continuum models for dynamic simulations (DEM-CNEM)[END_REF], to dynamics using 3D models. The present paper attempts to give more light on the setting of these parameters for dynamic studies.

This paper is organized as follows. In Section 2, the numerical methods most commonly used at different scales are classified according to their underlying physics and applications. This is followed by a classification of the coupling approaches between these methods so that the place of the concurrent coupling approaches can be seen. In Section 3, the techniques commonly used for concurrent coupling of discrete and continuum methods are reviewed. Among these techniques, the bridging domain (BD)

technique, which presents a good candidate to study dynamic problems, is studied in Section 4, so as to corroborate the previous studies on the BD parameters and to give more light on their practical setting. This study is performed using DEM-CNEM coupling approach recently developed between Discrete Element Method (DEM) and Constrained Natural Element Method (CNEM) [START_REF] Jebahi | 3D coupling approach between discrete and continuum models for dynamic simulations (DEM-CNEM)[END_REF]. Section 5 presents some conclusions.

Classification of numerical methods

As mentioned in the introduction, several numerical methods used in Computational Mechanics can be found in the literature. These methods can mainly be classified into discrete methods (DMs), continuum methods (CMs). The main features and specificities of these two classes are reviewed in this section. Subsequently, the different types of coupling approaches between these methods are discussed so as to adequately position the notion of multiscale modeling, which is the scope of this review.

Discrete methods

The discrete methods (DMs), such as Molecular Dynamics (MD) [START_REF] Alder | Phase Transition for a Hard Sphere System[END_REF][START_REF] Alder | Studies in Molecular Dynamics. I. General Method[END_REF] and Discrete (or Distinct)

Element Method (DEM) [START_REF] Cundall | A discrete numerical model for granular assemblies[END_REF], are a family of numerical methods consisting in simulation of a set of rigid (or pseudo-deformable) bodies in interaction with each other by contact laws. Depending on the analysis scale, one can distinguish quantum mechanical (or ab initio) methods (QMs), atomistic methods (AMs) and mesoscopic discrete methods (MDMs) (Figs. 1 and2).

Quantum mechanical methods

The quantum mechanical (QMs) methods are used for material simulations at the atomic scale (∼ 10 -9 m), in which the electrons are the players (Fig. 2). They consist in solving the electronic Schrödinger equation for atoms and molecules, taking into account the electronic structure of each atom. In other words, the interaction laws between atoms are directly dictated by their quantum mechanical state. Examples of QMs methods are Quantum Monte Carlo (QMC) [START_REF] Foulkes | Quantum Monte Carlo simulations of solids[END_REF] and Quantum Chemistry (QC) [101], which are very accurate but computationally too demanding to handle more than a few tens of atoms. Other QMs methods are Density-Functional Theory (DFT) and Local Fig. 1: Classification of discrete methods Fig. 2: Characteristic length and time scales for numerical methods Density Approximation (LDA) [START_REF] Hohenberg | Inhomogeneous Electron Gas[END_REF][START_REF] Payne | Iterative minimization techniques for ab initio totalenergy calculations: molecular dynamics and conjugate gradients[END_REF]. These methods are less accurate than QMC or QC, but they can be readily applied to systems containing several hundred atoms for static properties. Dynamic simulations with DFT and LDA are usually limited to timescales of a few picoseconds. Overall, the QMs methods are generally very accurate since they hold out the possibility of performing simulations without prior need for tuning. However, they are extremely expensive and can only be applied on very small domains of a few nanometers size.

Atomistic methods

The atomistic methods (AMs) are used for material simulations at the microscopic scale (∼ 10 -6 m),

where atoms are the players (Fig. 2). The interaction laws between atoms can be described by empirical interatomic potentials that encapsulate the effects of bonding (mediated by electrons) between them.

These potentials can be obtained from fitting material properties from experimental data or quantum mechanical calculations. They may depend on the distance between the atoms, angles between bonds, angles between planes, etc. Equation [START_REF] Abraham | Simulating materials failure by using up to one billion atoms and the world's fastest computer: Brittle fracture[END_REF] gives the general form of the interatomic potential:

U (r 1 , r 2 , ..., r N ) = i U 1 (r i ) + i,j, j>i U 2 (r i , r j ) + i,j,k, k>j>i U 3 (r i , r j , r k ) + ... ( 1
)
where N is the total number of particles (atoms), U i is the i-particle part of U , r i the vector position of a particle i. Examples of AMs methods are Molecular Mechanics (Statics) (MM) [START_REF] Hehre | A Guide to Molecular Mechanics and Quantum Chemical Calculations[END_REF], Molecular Dynamics (MD) [START_REF] Alder | Phase Transition for a Hard Sphere System[END_REF][START_REF] Alder | Studies in Molecular Dynamics. I. General Method[END_REF] and Monte Carlo (MC) [START_REF] Metropolis | The Monte Carlo Method[END_REF]. Although they are less accurate than the QMs methods, the AMs methods are relatively inexpensive (compared to QMs methods) and are able to provide insight into atomic processes involving considerably large systems of up to 10 9 atoms [START_REF] Abraham | Simulating materials failure by using up to one billion atoms and the world's fastest computer: Brittle fracture[END_REF].

Dynamic simulations with AMs methods are generally limited to timescales of a few nanoseconds, which can be crippling for simulation of realistic mechanical problems.

Mesoscopic discrete methods

To overcome the timescale limitations of the quantum mechanical and atomistic methods, another generation of discrete methods has been developed: mesoscopic discrete methods (MDMs). These methods can be used for material simulations at the mesoscopic scale (∼ 10 -4 m), where lattice defects such as dislocations, crack propagation, and other microstructural elements are the players.

At this scale, the system is too small to be regarded as a continuum and too large to be simulated effectively using quantum mechanical or atomistic methods. More precisely, the mesoscopic scale can be defined as an intermediate scale at which the microscopic phenomena (e.g. particle motion) can be assumed in mechanical equilibrium. The MDMs can broadly be regarded as a generalization of the atomistic methods (AMs), where more complex interaction laws are used. These laws are usually derived by calibration or from phenomenological theories that encompass the effects of interactions between the atoms. In the MDMs, the atomic degrees of freedom are not explicitly treated and only larger-scale particles are modeled. The fundamental principle of dynamics (FPD) is applied to compute the displacements and rotations of these particles. Originally, this class of methods has been developed for problems of granular media in rocks mechanics [START_REF] Cundall | Computer Model for Simulating Progressive Large Scale Movements in Blocky Rock Systems[END_REF]. More recently, it has been applied

to study the damage of heterogeneous media such as concrete [START_REF] Hentz | Discrete element modelling of concrete submitted to dynamic loading at high strain rates[END_REF] and rocks [START_REF] Bobet | Numerical models in discontinuous media: review of advances for rock mechanics applications[END_REF], but also to study the damage of homogeneous media such as ceramics [104] and glasses [START_REF] André | Discrete element method to simulate continuous material by using the cohesive beam model[END_REF][START_REF] André | Using the discrete element method to simulate brittle fracture in the indentation of a silica glass with a blunt indenter[END_REF][START_REF] Jebahi | Simulation of Vickers indentation of silica glass[END_REF][START_REF] Jebahi | Discrete-continuum coupling method for simulation of laser-inducced damage in silica glass[END_REF]. In these studies, the material is modeled by an agglomerate of discrete elements which interact via bilateral cohesive links to ensure the cohesion of the medium. The cohesive links can be of different nature according to the physical properties of the material being modeled (Fig. 3). Figure 4 illustrates an example of a continuum modeled with MDM. Application of MDMs in continuum media simulations faces two major challenges:

choice of the cohesive links and identification of their microscopic parameters to ensure a correct macroscopic behavior;

construction of discrete domains able to take into account the structural properties of the problem domain, e.g. homogeneity and isotropy, such that the macroscopic mechanical properties are independent of the number of discrete elements.

Several works can be found in the literature which propose solutions of these points. Among them, one can cite the recent works of André et al. [START_REF] André | Discrete element method to simulate continuous material by using the cohesive beam model[END_REF].

Nowadays, the MDMs present an alternative way to study physical phenomena requiring a small scale analysis, for which the continuity assumption is no longer valid, or problems with discontinuities that cannot be easily treated by continuum methods (CMs). They have become widely used to study realistic complex mechanical problems, such as cracking behavior of silica glass [START_REF] André | Discrete element method to simulate continuous material by using the cohesive beam model[END_REF][START_REF] André | Using the discrete element method to simulate brittle fracture in the indentation of a silica glass with a blunt indenter[END_REF][START_REF] Jebahi | Simulation of Vickers indentation of silica glass[END_REF][START_REF] Jebahi | Discrete-continuum coupling method for simulation of laser-inducced damage in silica glass[END_REF]. The benefits of these methods have attracted several researchers, and consequently several variations of MDMs have been developed. These variations can be divided into four categories as shown in Figure 1. The fundamental concepts of each one is briefly recalled hereafter. Lattice models In lattice models, a solid is modeled by a set of nodes connected with truss or beam elements [START_REF] Schlangen | Experimental and numerical analysis of micromechanisms of fracture of cementbased composites[END_REF][START_REF] Schlangen | Simple lattice model for numerical simulation of fracture of concrete materials and structures[END_REF]. Typically, nodes have neither masses nor volumes (they do not occupy volumes).

Solving a mechanical problem with this class of discrete methods is based on the construction of a global stiffness matrix from the local connection properties. Both regular and irregular lattices were studied. Originally, the lattice models were used to represent elastic continuum; the equivalence was established for both truss [START_REF] Hrennikoff | Solution of problems of elasticity by the frame-work method[END_REF] and beam [START_REF] Schlangen | New method for simulating fracture using an elastically uniform random geometry lattice[END_REF] elements. Later on, obvious enhancements, such as brittle beam failure, were introduced. Lattice models nicely show the emergence of relatively complex structural behaviors, although fairly simple formulas are used to describe the governing local processes.

The major drawback of these models is that the nodes do not have volumes, which can cause numerical problems related to crack closure in post-fracture stage. Solutions are given to circumvent this problem. These solutions consist in associating to each element an equivalent volume, based on the spatial Voronoï decomposition [START_REF] Ragueneau | Inelastic behavior modelling of concrete in low and high strain rate dynamics[END_REF] for example.

Particle models This class of methods is very close to the first discrete approach proposed in the literature by Cundall and Strack [29,[START_REF] Cundall | A discrete numerical model for granular assemblies[END_REF]: Distinct (Discrete) Element Method (DEM). Contrary to lattice models, particle models consider elements with masses and volumes. These elements often have disc shape (in 2D) or spherical shape (in 3D): only a single parameter (the radius) is required to determine the elements geometry and there is only one possible contact easily detectable between them. Consequently, computer memory requirements and processing time are minimized with these element shapes, even when a relatively large number of elements is used. Nevertheless, discs and spheres can roll or rotate easily. This does not reflect the expected behavior for several materials, in the case of large shear processes for example. To solve this problem, more complex shapes such as ellipses [105], ellipsoids [START_REF] Lin | A three-dimensional discrete element model using arrays of ellipsoids[END_REF], polygons [START_REF] Issa | Numerical analysis of micromechanical behaviour of granular materials[END_REF] and polyhedra [START_REF] Cundall | Formulation of a three-dimensional distinct element model -Part I, A scheme to detect and represent contacts in a system composed of many polyhedral blocks[END_REF] were proposed in the literature to provide more flexibility for element characterization in particle models.

Basically, the associated algorithm involves two stages. In the first stage, interaction forces are computed when elements slightly interpenetrate each other. This force-interpenetration formulation is generally referred to as a "Smooth contact" method or "force-displacement" method. Actually, the interpenetration between discrete elements, which makes no mechanical sense, represents the relative deformation of the elements surface layers. In the second stage, Newton's second law is applied to determine the acceleration of each element, which is then integrated, using "dynamic explicit" schemes, to find the new elements velocities and positions. This process is repeated until the simulation is achieved.

Hybrid models

The features of the lattice and particle models are largely complementary. Indeed, in lattice models, additional treatments must be made to deal with fracture problems (crack closure).

This problem is not encountered in particle models in which elements have their own volumes. However, particle methods cannot correctly model a continuum using simply disc or spherical elements, especially when shear effects are significant. This problem can be solved using cohesive beams between elements, such as in lattice models. The complementary advantages of these models has given rise to a family of "hybrid" methods combining particle and lattice models, i.e. by considering sphere elements connected with cohesive beams [START_REF] Griffiths | Modelling of elastic continua using a grillage of structural elements based on discrete element concepts[END_REF]. Interesting recent works on the Hybrid models can be found in references [START_REF] André | Discrete element method to simulate continuous material by using the cohesive beam model[END_REF][START_REF] André | Using the discrete element method to simulate brittle fracture in the indentation of a silica glass with a blunt indenter[END_REF][START_REF] Jebahi | Simulation of Vickers indentation of silica glass[END_REF][START_REF] Jebahi | Discrete-continuum coupling method for simulation of laser-inducced damage in silica glass[END_REF].

Contact dynamics This class of methods provides an alternative approach based on a "non-smooth" formulation of mutual exclusion and dry friction between elements [START_REF] Jean | The non smooth contact dynamics method[END_REF][START_REF] Luding | Simulations of pattern formation in vibrated granular media[END_REF][START_REF] Moreau | Some numerical methods in multibody dynamics: application to granular materials[END_REF]. It introduces the notion of non-smooth (irregular) contact between elements which is, at present, subject of several studies.

Interpenetration between elements is prevented: no elastic contact law is used between them. There are two main numerical integrators for these methods, which are of "dynamic implicit" type: the event-driven integrators, also referred to as Even-Driven Method (EDM) [START_REF] Luding | Simulations of pattern formation in vibrated granular media[END_REF], and the so-called timestepping integrators, also referred to as Contact Dynamics Method (CDM) [START_REF] Jean | The non smooth contact dynamics method[END_REF][START_REF] Moreau | Some numerical methods in multibody dynamics: application to granular materials[END_REF]. In EDM, a collision or "event" occurs when two rigid elements touch each other and the post-collisional and angular velocities are prescribed by a collision operator [START_REF] Rapaport | The event scheduling problem in molecular dynamic simulation[END_REF]. Despite being very accurate, the even-driven integrators treat only one force at a time. Therefore, they are not well adapted for problems with many simultaneous contacts, as often encountered in mechanics. To overcome this limitation, Jean and Moreau [START_REF] Jean | The non smooth contact dynamics method[END_REF][START_REF] Moreau | Nonsmooth mechanics and applications[END_REF] have developed the Contact Dynamic Method (CDM) which has a specialized numerical scheme for problems with many contacts. The governing equations are expressed as differential inclusions and the accelerations are replaced by velocity jumps. In the generic CDM algorithm, an iterative process is used to compute forces and velocities. This process consists in solving a single contact problem with all other contacts kept constant, and iteratively updating the forces until a convergence criterion is fulfilled. Two basic kinematic constraints are used between elements in the CDM formulation:

-The Signorini conditions which state that the normal force f n is repulsive when the elements are in contact (distance between them is zero), and f n = 0 otherwise. To deal with persistent contacts, f n is reset to zero when no relative velocity exists between elements in contact.

-The Coulomb's friction law, which relates the sliding velocity and the friction force f t .

These kinematic constraints can also be complemented by a "rolling friction" constraint which introduces a moment resistance [START_REF] Bratberg | Dynamic rearrangements and packing regimes in randomly deposited twodimensional granular beds[END_REF]. Within the CDM, the time resolution is much larger than the collision characteristic time (unlike in the case of particle approaches). Instead, the time step represents a unit of time during which collisions can occur, causing velocity jumps. Although CDM has been successfully used for several geomechanical problems [START_REF] Donzé | Advances in Discrete Element Method applied to Soil, Rock and Concrete Mechanics[END_REF], it is much more difficult to implement than lattice and particle models. Also, the prediction of the contact forces and particle velocities in the following step from the current configuration is very problematic and is currently the subject of several studies.

The Contact Dynamics models are generally used to study quasi-static problems or problems with relatively low dynamic effects. This class of methods is perfectly suitable to study mechanical problems of granular mechanics. However, in the case of continuum media, the use of models based on regular or "smooth" interaction laws seems to be advantageous since the elasticity is naturally taken into account by these interaction laws.

Continuum methods

The continuum methods (CMs), such as Finite Element Method (FEM) [120, 121, 122] and Smoothed Particle Hydrodynamics (SPH) [START_REF] Liu | Smoothed particle hydrodynamics : a meshfree particle method[END_REF][START_REF] Lucy | Numerical approach to testing the fission hypothesis[END_REF], are used for material simulations at the macroscopic scale (∼ 10 -2 m), where continuum fields (e.g. density, temperature, displacement, stress, etc) are the players. The behavior of the studied domain, which is viewed as a continuum and completely filling the space it occupies, is governed by constitutive laws stemming from Continuum Mechanics. These constitutive laws are usually formulated so that they can capture the effects of the microstructure and lattice defects on materials properties. Solving mechanical problems using the CMs methods is typically a two-step process. In the first step, the continuum problem domain is approximated by a finite number of discrete components made up of reference points or nodes. Then, the continuous governing equations are changed into a system of algebraic equations. Solving this system gives the solutions at the reference points. For other points in the occupied space, the field variables can be approximated by interpolation or averaging (in the case of SPH) from the solutions at the reference points.

The CMs methods are well adapted to simulate macroscopic physical phenomena (having length scales much greater than the interatomic distance), in which the continuity assumptions are valid and remain valid during calculation. However, they are less suitable to study problems with occurrence of discontinuities, such as wear, fracture and abrasion problems. Indeed, additional treatments must be considered to deal with the new discontinuity surfaces. Compared to the discrete methods, the continuum methods are generally much less expensive. Consequently, they have become the most commonly used in Computational Mechanics and have been widely applied to various area of Computational Solid Mechanics (CSM) and Computational Fluid Dynamics (CFD). This resulted in development of a wide variety of continuum methods which can be divided into two classes: grid-based methods and meshless methods (Fig. 5). The main features of each class are reviewed hereafter.

Grid-based methods

In the grid-based methods the spatial domain is often represented by discretized elements. These elements are called "meshes" in the case of Finite Element Method (FEM) [122], grids in the case of Finite Difference Method (FDM) [START_REF] Smith | Numerical Solution of Partial Differential Equations: Finite Difference Methods[END_REF] and volumes or cells in the case of Finite Volume Method (FVM)

[108]. The terminologies of grids, volumes, cells, and elements carry different physical meanings linked to physical problems. However, all of them can be termed meshes according to the following definition.

A mesh, in the broadest sense of the word, is defined as an open space or interstice between the strands of a network which is formed by connecting nodes in a predefined manner. The key here is that, in the grid-based methods, the mesh must be predefined to provide a certain relationship between the nodes. This forms the formulation basis of this class of methods. Based on a properly predefined mesh, complex ordinary differential or partial differential governing equations can be approximated by a set of algebraic equations for each discretized element (mesh). The system of algebraic equations for the whole problem domain can be formed by assembling the elementary algebraic equations for all the meshes.

There are two main frames to describe the governing equations of mechanical problems: Lagrangian description and Eulerian description. The former, also called material description, is a way to look at material motion where the observer follows an individual material point as it moves through the space and time. The latter, also called spatial description, is another way to look at material motion where the observer focuses on a specific space position through which the material passes over time. The use of one or other of these descriptions leads to quite different governing equations. The difference between them arises from the definition of the total time derivative which is defined, in the case of

+ V .∇ convective derivative (2)
where V represents a velocity vector and ∇ denotes the vector differential operator (gradient). According to the used description, the grid-based methods can be divided into three subclasses: Lagrangian methods, Eulerian methods and Combined Lagrangian-Eulerian methods.

Lagrangian methods Lagrangian methods which are typically represented by the FEM [122] are based on the Lagrangian description. In this class of methods, the grid or mesh is attached to the material in the entire computation process, and moves with it (Fig. 6).

Fig. 6: Lagrangian mesh for a rabbit simulation Each node in the mesh follows the material path at the grid point. As a result, the mesh elements can be deformed during the simulation process. When the domain deforms, the mesh elements deform accordingly. Mass, momentum and energy are transported with the motion of the mesh. Because the mass within each cell (or element) remains constant, no mass flux is allowed through the mesh cells boundaries. The Lagrangian methods present several advantages:

-It is very easy to treat the boundary conditions at free surfaces, moving boundaries and material interfaces. Indeed, they are automatically imposed, tracked and determined simply by placing some grid nodes along the boundaries and the material interfaces.

-Complex geometries can be conveniently studied using an irregular mesh.

-Meshing is only required within the problem domain (no additional grids are required beyond it, as will be seen for Eulerian methods). Therefore, Lagrangian methods are computationally efficient.

-Since the mesh is attached to the material, the entire time history of all field variables at a specific material point can easily be obtained during the simulation.

-In the Lagrangian description, the total time derivative is the same that the local time derivative. Therefore, no convective terms exist in the associated governing equations. This makes Lagrangian simulations simpler and faster as no computational efforts are required to treat the convective terms.

Due to these advantages, Lagrangian methods have become very popular to simulate problems of Computational Solid Mechanics (CSM), where deformations are relatively small. However, application of these method to solve problems with large deformations presents a huge challenge. Indeed, in such problems, the mesh can become extremely distorted, and stability as well as convergence difficulties can arise, leading to computation breakdown. The accuracy of the formulation and consequently the solution can be severely affected. Also, the time step, which is controlled by the size of the smallest element grid, can become too small to be efficient for computation.

A solution to overcome these problems is to re-mesh (rezone) the problem domain or simply the regions where the initial mesh is hugely distorted. This technique involves overlying of a new (undistorted) mesh on the old one. The computation is then resumed on the new mesh. The field variables can be approximated at the new grid nodes via mass, momentum and energy transport using an Eulerian description. Despite the popularity of the rezoning techniques to simulate large deformation problems, they suffer from several difficulties. Indeed, rezoning procedure can be tedious and time-consuming. Besides, the transport of the field variables from old to new mesh is generally accompanied by material diffusion which can lead to a loss of the material history [START_REF] Benson | Computational methods in Lagrangian and Eulerian hydrocodes[END_REF].

Eulerian methods In Eulerian methods, the mesh is fixed on the space occupied by the simulated object which can move across the fixed mesh cells (grids) as shown in Figure 7. All the mesh cells as well as the associated nodes remain fixed in the space and do not move with the material moving across the mesh. The flux of mass, momentum and energy across mesh cell boundaries is simulated to compute the distribution of different field variables in the problem domain. As a result, the shape and the volumes of the grid cells remain constant during the entire computation process. Therefore, simulation of problems with large deformations is possible with Eulerian methods. Since they do not induce any cells deformation and thus do not cause the same kind of numerical problems as in the Lagrangian methods. For these reasons, Eulerian methods are dominant in Computational Fluid Dynamics (CFD) problems and problems with large deformations such as explosion, high velocity impacts (HVI), etc.

Despite the great success of Eulerian methods, there are many disadvantages which limit their application, from which one can cite:

-Since the Eulerian methods simulate the flux of mass, momentum and energy across mesh cells boundaries, it is very difficult to get out the time history of field variables at a fixed material point. Fig. 7: Eulerian mesh for a rabbit simulation -Eulerian methods require meshing not only within the problem domain, but also beyond it (Fig. 7). The mesh must be large enough to cover the area through which the simulated object can move. Therefore, coarse grid is sometimes recommended for computational efficiency, which can affect the solution accuracy.

-The free surfaces, moving and deformable boundaries and moving material interfaces are difficult to be determined accurately.

-It is difficult to treat irregular and complicated geometries of materials. Generally, tedious mesh generation procedure to convert the irregular geometry of the problem domain into a regular computational domain is necessary. This can become sufficiently expensive to be efficient. Efforts in this direction are still ongoing.

Combined Lagrangian-Eulerian methods The features of the Lagrangian and Eulerian methods are largely complementary. This has led to the development of several approaches involving in the same time the Lagrangian and Eulerian descriptions. The Coupled Eulerian Lagrangian (CEL) [START_REF] Felici | A coupled Eulerian/Lagrangian method for the solution of three-dimensinal vortical flows[END_REF][START_REF] Mair | Hydrocode methodologies for underwater explosion structure medium/interaction[END_REF]113] and the Arbitrary Lagrangian Eulerian (ALE) [START_REF] Benson | Computational methods in Lagrangian and Eulerian hydrocodes[END_REF][START_REF] Stein | [END_REF] can particularly be cited. The CEL approach consists in employing simultaneously both the descriptions in separate (or with overlap) regions of the problem domain. The Lagrangian region continuously interacts with the Eulerian one through a coupling module allowing for exchange of computational information. The ALE approach is conceptually close to the rezoning techniques for Lagrangian methods. In this case, the mesh moves independently of the material so as to minimize the mesh distortion. In the ALE approach, Lagrangian motion is firstly computed at each time step, followed by possible rezoning step in which the mesh is either rezoned to the original shape or rezoned to some more advantageous shape (between the previous shape and the current one). without any connectivity between them. Contrary to the grid-based methods, this class of methods normally do not require any predefined mesh to establish the system of algebraic equations for the whole problem domain. Instead, only a set of nodes scattered within the problem domain and along the boundaries are used to represent the simulated geometry (Fig. 8). In such a representation, no connectivity between nodes is required. The ideal requirements of the meshless methods are:

non-necessity of a predefined mesh, at least in field variable interpolation;

non-necessity of a predefined mesh at all throughout the computation process.

However, the meshless methods developed so far are not really ideal considering these two points. For example, some meshless methods require background cells (mesh) to derive the system of algebraic equations. The Element Free Galerkin (EFG) method [START_REF] Belytschko | Element-Free Galerkin methods[END_REF] and the Meshless Local Petrov-Galerkin (MLPG) method [START_REF] Atluri | A new meshless local Petrov-Galerkin (MPLG) approach in computational mechanics[END_REF], which use background cells for integration of system matrices derived from the weak formulation, belong to this category.

In the literature, one can find numerous meshless methods, some of which share several common features. It should be noted that most of the meshless methods are based on the Lagrangian description. Therefore, construction of shape functions is necessary before or during the computation process.

These shape functions will next be used to approximate the field variables at any point x within the problem domain as follows (3):

u(x) = n i=0 φ i (x)u i ( 3
)
where u is the field variable to be approximated (e.g. displacement), n is the number of all the nodes, u i is the nodal field variable at the node i, and φ i is the shape function associated to the node i. According to the type of shape functions, the meshless methods can be classified into two main categories: approximation methods and interpolation methods. The main features of each category are detailed hereafter.

Fig. 8: Meshless discretization for rabbit simulation

Approximation methods Contrary to the Lagrangian grid-based methods, the shape functions in meshless methods are constructed for the nodes (not for element grids). In the case of approximation methods, these functions generally do not satisfy the Kronecker delta function property (4) at the nodes.

φ i (x j ) = 1 i = j 0 i = j (4)
In other words, an approximated field variable based on these functions, called approximation shape functions, does not necessary pass through the nodal values (u

(x i ) = u i ).
The first method developed in this category is the famous Smoothed Particle Hydrodynamics (SPH) method [START_REF] Liu | Smoothed particle hydrodynamics : a meshfree particle method[END_REF][START_REF] Liu | Smoothed Particle Hydrodynamics (SPH): an Overview and Recent Developments[END_REF][START_REF] Lucy | Numerical approach to testing the fission hypothesis[END_REF][START_REF] Randles | Smoothed particle hydrodynamics: Some recent improvements and applications[END_REF]. This method was originally developed to model astrophysical phenomena, and later it has been widely extended for application to problems of continuum solid and fluid mechanics [START_REF] Lucy | Numerical approach to testing the fission hypothesis[END_REF]. The SPH method and its different variants, such as Corrective Smoothed Particle Method (CSPM) [START_REF] Liu | Smoothed particle hydrodynamics : a meshfree particle method[END_REF] and Discontinuous Smoothed Particle Hydrodynamics (DSPH) [START_REF] Liu | Smoothed particle hydrodynamics : a meshfree particle method[END_REF], are the major types of approximation methods, and have been incorporated into many commercial codes.

Several other methods belonging to this category can be found in the literature. Among them, one can cite the Generalized Finite Difference Method (GFDM) [START_REF] Liszka | The finite difference method at arbitrary irregular grids and its applications in applied mechanics[END_REF], Element Free Galerkin Method (EFGM) [START_REF] Belytschko | Element-Free Galerkin methods[END_REF], Reproducing Kernel Particle Method (RKPM) [START_REF] Li | Meshfree and particle methods and their applications[END_REF], etc. Several advantages have made the approximation methods very attractive:

-Since no predefined grid is required, approximation methods are able to deal with problems involving extremely large deformations.

-Approximation methods are more adapted than the Eulerian grid-based methods to study problems with free surfaces, deformable boundary and moving interfaces.

-The approximation methods are more adapted than the grid-based methods to study fracture of materials. Indeed, they do not use any kind of mesh which can affect the cracking mechanisms.

-Approximation methods are more adapted than grid-based methods to treat complex geometries which are represented only by a set of scattered nodes with no connectivity between them.

-Applications of this class of methods are very wide, ranging from microscale to macroscale and even to astronomical scale, and from discrete systems to continuum systems.

-Due to intensive research efforts, several approximation methods, such as SPH, are quickly approaching their mature stage.

-The accuracy of solutions is very weakly affected by the relative positions of the nodes. Hence, it is very easy to add or remove nodes. This point is particularly attractive to define an adaptive refinement strategy for problems with localization.

However, the approximation methods have several drawbacks related to the imposition of boundary conditions and implementation requirements:

-To impose essential boundary conditions in a straightforward way as in FEM, it is necessary that: (i) the approximated field variables pass through the nodal value (interpolation shape functions) and (ii) the influence of the internal nodes vanishes at the border of the problem domain (linear shape functions at the border). Since, for approximation methods, generally any of these conditions is satisfied, it is not easy to impose essential boundary conditions. Therefore, special techniques are required to overcome this problem. In the literature, one can find several works on this subject and several techniques have been proposed. Among them, one can cite the Lagrange multiplier method [START_REF] Belytschko | Element-Free Galerkin methods[END_REF], Penalty method [118] and the introduction of singular weight functions [START_REF] Kaljevic | An improved element free Galerkin formulation[END_REF]. But, in most cases, these techniques are tedious and involve additional computational efforts which can limit the application of these methods to problems with complex/moving interfaces and boundaries.

-The second drawback of the approximation methods is related to the numerical integration. Actually, almost all the shape functions used in meshless methods are rational (not polynomial).

This makes the Gauss integration schemes not well adapted for this class of methods. Dolbow et al. [START_REF] Dolbow | Numerical integration of the Galerkin weak form in meshfree methods[END_REF] have demonstrated that, using these integration schemes, a great number of Gauss points is required which can induce an extensive additional cost of computation. Concerning the other popular integration technique, Direct Nodal Integration, Chen et al. [START_REF] Chen | A stabilized conforming nodal integration for Galerkin mesh-free methods[END_REF] have demonstrated that the application of this technique can involve numerical instability problems. To overcome these problems, several solutions have been proposed in the literature, such as the Stabilized Conforming Nodal Integration scheme [START_REF] González | Numerical integration in Natural Neighbour Galarkin methods[END_REF][START_REF] Chen | A stabilized conforming nodal integration for Galerkin mesh-free methods[END_REF]. However, they are generally less accurate than the Gauss integration schemes for Lagrangian grid-based methods.

-Another drawback of this class of methods is related to the definition of the influence area of the particles (nodes), called support domain or also influence domain. In most cases, the support domain of a particle is taken as a sphere or parallelepiped centered on the particle. As discussed in several previous works [START_REF] Liu | Smoothed particle hydrodynamics : a meshfree particle method[END_REF][START_REF] Liu | Smoothed Particle Hydrodynamics (SPH): an Overview and Recent Developments[END_REF], this support should contain a sufficient number of neighbor particles for computation stability and then be sufficiently wide. However, a very wide support domain can affect the accuracy of the solutions. The choice of the size of the support domain is not trivial and represents, at present, an open question [START_REF] Lee | On error estimation and adaptive refinement for element free galerkin method : Part I : stress recovery and a posteriori error estimation[END_REF].

-As reported in several previous works [START_REF] Belytschko | Meshless methods : an overview and recent developments[END_REF][START_REF] Lee | On error estimation and adaptive refinement for element free galerkin method : Part I : stress recovery and a posteriori error estimation[END_REF], the application of approximation methods to nonconvex geometries, e.g. the presence of concave boundaries, cavity or cracks, is very challenging. This is mainly due to difficulties to define the support domain in the vicinity of the concave borders. In this case, the support domain can contain nodes which are not really neighbors, e.g.

nodes on opposite sides of a crack.

-Finally, compared to Lagrangian grid-based methods, approximation methods generally leads to poor and inaccurate solutions. The accuracy problem is much more amplified in the vicinity of the border of the problem domain.

Despite all the advantages of the approximation methods (meshless methods), the different drawbacks of these methods have limited their application to several mechanical problems.

Interpolation methods To overcome the limitations of the approximation methods, another generation of meshless methods have been developed: interpolation methods. In this generation, the field variables are approximated using interpolation shape functions, verifying the Kronecker delta function property (4) at the nodes. In other words, the approximated field variables pass through the nodal values (u(x i ) = u i ). The first interpolation meshless method that has been developed is the Natural Element Method (NEM) [START_REF] Braun | A numerical method for solving partial differential equations on highly irregular evolving grids[END_REF]107]. In this method, the interpolation shape functions are constructed based on the notion of Natural Neighbor (NN) interpolation. The NN interpolation is defined as a multivariate data interpolation scheme [START_REF] Sibson | A brief description of natural neighbour interpolation[END_REF], which has primarily been used in data interpolation and modeling of geophysical phenomena. It relies on the concepts of Voronoï diagrams [109] and Delaunay tessellations [START_REF] Delaunay | Sur la sphère vide. A la mémoire de Georges Voronoï[END_REF], which are widely used in computational geometry, to construct the shape functions (called NN shape functions). With this type of interpolation, the NEM method has overcome most of the difficulties encountered with the earlier generation of meshless methods, such as definition of the support domains of the nodes. However, application of this method for non-convex geometries remains difficult. In effect, the NN shape functions are not strictly linear over non-convex boundaries, which makes difficult the imposition of Dirichlet boundary conditions on such boundaries (the test functions do not vanish over non-convex boundaries). This is because NEM is closely related to the convex hull of a given set of points.

Recently, Cueto et al. [START_REF] Cueto | Imposing essential boundary conditions in the natural element method by means of density-scaled -shapes[END_REF][START_REF] Cueto | Overview and recent advances in Natural Neighbour Galerkin Methods[END_REF] have proposed an extension of NEM, called α-NEM, based on the concept of α-shapes or α-complexes, which are widely used in science visualization to give a shape to a cloud of points. The use of α-shapes in the context of NEM allows the construction of models entirely in terms of nodes and also ensures the linear precision of the shape function even over nonconvex boundaries. This is particularly true in the case of non-strongly non-convex domains. Indeed as mentioned by Cueto et al. [START_REF] Cueto | Imposing essential boundary conditions in the natural element method by means of density-scaled -shapes[END_REF][START_REF] Cueto | Overview and recent advances in Natural Neighbour Galerkin Methods[END_REF], the α-shapes used in such a method fail to ensure the linearity of the shape functions over highly non-convex boundaries (Fig. 9), e.g. around a crack tip. α-NEM is therefore limited to non-strongly concave geometries. Otherwise, the same problem related to the imposition of Dirichlet boundary conditions in the NEM, is raised. . This method can be regarded as another extension of the NEM which retains its attractive properties and overcomes the difficulties related to its application for strongly non-convex domains. The main originality of this approach is the introduction of a visibility criterion [115,116] in the NEM to select the natural neighbors of the nodes and construct the interpolation. A modified so-called constrained Voronoï diagram is introduced for the computation of the shape functions. As a meshless method, CNEM has practically all the advantages of this class of methods. Also, with the introduction of the visibility criterion, it overcomes almost all its drawbacks.

Coupling methods

As seen above, numerous numerical methods can be found in the literature. These methods can be classified into discrete methods and continuum methods. The first class is used to study phenomena ranging from nanoscale to mesoscale. The second class is generally used to study macroscopic phenomena, where the mechanical behavior can be described by Continuum Mechanics. The features of the different methods of these classes are largely complementary. This suggests that it would be computationally beneficial to combine two or more methods from these classes so as to strengthen their advantages and to avoid their drawbacks. This idea has attracted several researchers and several coupling methods have been developed. According to the analysis scale of the involved methods, two types of coupling approaches can be distinguished: monoscale coupling methods and multiscale coupling methods.

Monoscale coupling methods

The monoscale coupling methods are generally used to couple two or more continuum methods to study multiphysics problems, in which different continuum physical phenomena interact with each other. Among the multiphysics problems that can be treated by this type of coupling methods, one can cite the problems of dynamic fluid-structure interaction and high velocity impacts. Several works on the multiphysics coupling methods can be found in the literature. An example is the work of Fourey et al. [START_REF] Fourey | Violent Fluid-Structure Interaction simulations using a coupled SPH/FEM method[END_REF], who have developed a coupling method between the Smoothed Particle Hydrodynamics (SPH) and the Finite Element Method (FEM) to simulate violent fluid-structure interaction. The SPH method is applied to solve the Navier-Stokes equations in the fluid region, whereas, the FEM method is used to describe the mechanical behavior of the studied structure. Another example is the FEM-SPH coupling approach developed by Zang et al.

[117] to simulate high velocity impact of windshield. The SPH approach was only applied in the impact region to better characterize the fragmentation mechanisms. The monoscale coupling methods can also be applied to overcome the limitations of the approximation methods in the vicinity of the border of the studied domain. In such a coupling, an approximation method (e.g. SPH) is generally coupled to a Lagrangian grid-based method (e.g. FEM). This last method (Lagrangian grid-based method) is only used in the regions of application of the essential (displacement) boundary conditions, which cannot correctly be treated by approximation approaches. The remaining estate is modeled by the approximation method to better characterize the complex phenomena (e.g. large deformations) that can occur.

The monoscale coupling methods are beyond the scope of this review, which attempts to give light on the multicale modeling of mechanical problems. Therefore, it will not be examined here.

Multiscale coupling methods

The multiscale coupling methods are used to couple approaches at different scales. This type of coupling methods is applied to study problems involving on one scale phenomena that require very accurate and computationally expensive description, and on another scale phenomena for which coarser description is sufficient and, in fact, necessary to reduce the computation time. According to the nature of the coupled approaches, two cases can be distinguished: discrete-discrete coupling approaches (e.g. MD-DEM) and discrete-continuum coupling approaches (e.g. DEM-FEM). The discrete-discrete coupling approaches, in which only discrete methods at different scales are involved, are generally very expensive and can only be used to study small systems. The discrete-continuum coupling approaches, in which a combination of discrete and continuum methods are connected together, are less expensive and then more suitable to study realistic engineering applications.

In many mechanical problems, the notion of multiscale modeling arises quite naturally. In effect, the behavior of materials at the macroscopic scale, which is the scale of interest for technological applications, is determined by the microscopic interactions between atoms. This is why such a notion has become a special area of interest for many scientists, and then several multiscale coupling methods have been developed. These methods can mainly be divided into three classes: hierarchical, concurrent and hybrid hierarchical-concurrent multiscale coupling methods. It should be noted that there already exist a few review papers on multiscale modeling [START_REF] Lu | Handbook of Theoretical and Computational Nanotechnology, rieth, m. and schommers[END_REF]114], which use other classifications. In this paper, a three-class classification is proposed because it seems necessary and sufficient to include the main multiscale coupling approaches.

Hierarchical couping methods

The hierarchical coupling methods are the most widely used and computationally the most efficient. They are also referred in the literature to as sequential, serial, implicit, or message passing methods. These methods aim to piece together a hierarchy of numerical approaches in which the coarse scale models use information obtained by more detailed fine scale models. The homogenization methods for multi-phase elastic heterogeneous media are typical examples of hierarchical coupling methods. In these methods, the response of a representative volume element at the fine scale is first computed over a range of expected inputs, and from this, a stress-strain law is extracted which can be used at a coarser scale to describe the mechanical behavior of the homogenized material.

Other examples of hierarchical coupling approaches include almost all the Molecular Dynamics (MD) methods whose underlying potentials are derived from Quantum mechanical (QM) calculations.

The hierarchical technique is generally well suited for problems in which the different analysis scales are decoupled or weakly coupled. Also, it can be used when the large scale variations appear homogeneous and quasi-static from the fine scale point of view. However, application of this technique on nonlinear mechanical problems presents major challenges, and in some instances is invalid (i.e.

failure problems [START_REF] Belytschko | Multiscale aggregating discontinuities: A method for circumventing loss of material stability[END_REF]).

Concurrent coupling methods

The concurrent coupling methods consist in linking numerical approaches having different scales together in a single combined model, such that the fine scale model communicates directly and instantly with the coarse scale model though some coupling procedure.

These methods are also referred to as parallel or explicit coupling methods. Both compatibility and momentum balance are enforced across the interface between the coupled models. This type of coupling methods is well suited to study multiscale problems in which the behavior at each scale depends strongly on what happens at the other scale. Several works on the concurrent coupling methods can be found in the literature [START_REF] Aubertin | Energy conservation of atomistic/continuum coupling[END_REF][START_REF] Bauman | On the application of the Arlequin Method to the Coupling of Particle and Continuum Models[END_REF][START_REF] Belytschko | Coupling methods for continuum model with molecular model[END_REF][START_REF] Ben Dhia | The Arlequin method as a flexible engineering design tool[END_REF][START_REF] Broughton | Concurrent coupling of length scales: methodology and application[END_REF][START_REF] Fish | Concurrent AtC coupling based on a blend of the continuum stress and the atomistic force[END_REF][START_REF] Guidault | On the L2 and the H1 couplings for an overlapping domain decomposition method using Lagrange multipliers[END_REF][START_REF] Guidault | Bridging domain methods for coupled atomistic-continuum models with L2 or H1 couplings[END_REF][START_REF] Jebahi | On the H1 discrete -Continuum coupling based on the Arlequin method (DEM-CNEM)[END_REF][START_REF] Jebahi | 3D coupling approach between discrete and continuum models for dynamic simulations (DEM-CNEM)[END_REF]112]. A common feature of these works is that the problem domain is often partitioned into several subdomains characterized by different scales and physics. The challenge that arises here is how to ensure a smooth coupling between these subdomains.

To address this issue, several techniques have been developed. The main of these techniques will be reviewed in more detail in section 3.

A variation of concurrent coupling methods, generally referred to as semiconcurrent method, is that in which the coupled models run together and communicate instantly with each other but not intimately coupled. Compatibility and momentum balance are only satisfied approximately. The advantages of this approach lie in the fact that the coupled models can be computed by separate software, which can avoid numerical problems related to rapid change in element size. The FE 2 multiscale approach of Feyel and Chaboche [START_REF] Feyel | Multiscale approach for modelling the elastoviscoplasitic behavior of long fibre SiC/Ti composite materials[END_REF] is an example of this variation.

Hybrid hierarchical-concurrent coupling methods As seen above, the hierarchical methods are used for multiscale problems in which the coupled scales can be separated. The concurrent methods are used for strongly coupled multiscale problems. However, in some multiscale problems, the involved scales can be weakly coupled at the beginning of the computation up to a certain response limit, and subsequently become highly dependent. Therefore, it would be computationally beneficial to combine the above two coupling approaches to study these problems. A hierarchical approach is used as long as the requested fine scale information is available. When this information is no longer accessible, due to strain localization for example, a concurrent approach is invoked. If the fine scale response is history dependent, it is necessary to reconstruct at least an approximate history. This class of coupling methods is known in the literature as hybrid hierarchical-concurrent methods or hybrid hierarchical-semiconcurrent coupling methods (when the involved scales are not intimately coupled in the concurrent approach). An example of these methods is the adaptive multiscale approach developed by Akbari et al. [START_REF] Akbari | An adaptive multiscale method for fracture based on concurrent -hierarchical hybrid modelling[END_REF] to study the quasi-brittle crack propagation in metals. In this approach, the FE 2 technique is used in the safe regions of the studied domain. When strain localization appears, a domain decomposition (concurrent) coupling approach is used to solve the problem exactly at the material heterogeneities scale.

Compared to the first two coupling approaches, the hybrid hierarchical-concurrent coupling methods are relatively recent, and are, at present, the subject of several studies [START_REF] Akbari | An adaptive multiscale method for fracture based on concurrent -hierarchical hybrid modelling[END_REF][START_REF] Nuggehally | Adaptive Model Selection Procedure for Concurrent Multiscale Problems[END_REF] to alleviates the difficulties related to the transition from hierarchical to concurrent formulation.

3 Multiscale modeling of complex dynamic problems: concurrent discrete-continuum coupling methods Some of the most fascinating dynamic problems in all fields of science involve multiple spatial or temporal scales: processes that occur at a certain scale govern the behavior of the system across several (usually coarser) scales. Therefore, to simulate accurately such problems with minimum cost, it is often necessary to have recourse to multiscale modeling. Particularly, the concurrent multiscale coupling approaches between discrete and continuum methods seem best adapted to study these problems. Indeed, they offer several potential benefits, compared to other ones. On one hand, these approaches are of concurrent type. Therefore, as seen in the previous section, they are well adapted to study highly dependent multiscale phenomena, which are frequently encountered in complex dynamic problems.

On the other hand, both discrete and continuum methods are involved in these approaches, which allows to exploit their complementary advantages, while avoiding their drawbacks. The small scale phenomena can easily be treated by application of the discrete method in the associated regions.

Furthermore, the use of a continuum method, as a part of the coupling approach, allows to reduce considerably the computation time, and avoids prohibitively large computations. Consequently, realistic simulations of material science and engineering can be performed. Thanks to their advantages, the concurrent discrete-continuum coupling approaches have become very fashionable in the past decade.

After recalling the major coupling challenges usually encountered in these approaches, the techniques most commonly used in them will be listed hereafter. Although this list is not all-inclusive, it gives a clear vision of the topic. It should be noted that these techniques can also be applied to couple concurrently discrete methods (discrete-discrete concurrent coupling).

Coupling challenges

The major difficulty encountered in application of the discrete-continuum coupling approaches in dynamics is that, when the cutoff frequency of the continuum model is below that of the discrete one, spurious wave reflections can occur at the interface. This difficulty is almost always the case in this type of coupling approaches, since there is no benefit in using continuum model with discretization characteristic length of the same order of magnitude as for the discrete model. To illustrate simply the problem of spurious wave reflections, the behavior of waves in the discrete and continuum domains will be explored using simple 1D models. A chain of particles connected with equivalent springs is used for the discrete model, and its corresponding continuum counterpart is used for the continuum model (Fig. 10).

m dn = C (d n-1 -d n ) + C (d n+1 -d n ) = C (d n+1 + d n-1 -2 d n ) ( 5
)
where C is the stiffness of the springs, m and d n are respectively the mass and displacement of the particle n. Equation ( 5) can be solved by assuming a wavelike solution which is only defined on the particle positions:

d n = D e i (k n r-ω t) (6)
where D is the amplitude of the oscillations, k is the one dimensional wave vector (its amplitude gives the wavenumber κ = |k| = 2π λ ), λ is the wave length, r is the interparticle distance (x = n r is the position of the particle n along the chain) and ω is the angular frequency. Substituting (6) into [START_REF] Alder | Studies in Molecular Dynamics. I. General Method[END_REF] leads, after simplification, to:

ω 2 D = 4 C m sin 2 k r 2 (7)
which is known as dispersion relation. It characterizes the dependence of the wave frequency on the wavenumber κ = |k|. It can be shown from [START_REF] André | Using the discrete element method to simulate brittle fracture in the indentation of a silica glass with a blunt indenter[END_REF] that ω D (k) = ω D (k + 2 j π r ) for any j ∈ Z. Therefore, only the case of k ∈ [-π r , π r ] will be considered hereafter. Moreover, for symmetry reason, it is sufficient to restrict to k = κ > 0.

In the continuum model, the 1D equation of motion, in the absence of body forces, can be written as:

ρ ü = ∂σ ∂x ( 8
)
where ρ is the density, σ is the longitudinal Cauchy stress in x direction. Under the assumption of small deformations, equation ( 8) can be rewritten as:

ρ ü = E ∂ 2 u ∂x 2 (9)
where E is the Young's modulus. Assuming a harmonic wave propagating along x direction (u = U e i (k x-ω t) ), it is easy to find the continuum dispersion relation:

ω 2 C = E ρ k 2 (10)
To ensure equivalence between the discrete and continuum models, E and ρ must respectively verify the following conditions (the cross-sectional area S is assumed to be equal to the unit, for simplicity):

E = C r and ρ = m r (11)
Using these conditions 11, the continuum dispersion relation can be reformulated as follows:

ω 2 C = C r 2 m k 2 (12)
Influence of the spatial discretization on the dispersion relation As seen above, different dispersion relations are obtained for the discrete and continuum models. In this paragraph, the influence of the spatial discretization of the continuum model on the dispersion relation will be examined. To do this, the continuum domain is discretized into uniform segments of length h. Using linear interpolation functions, the discretized equation of motion of a node n can then be expressed:

ρ h ün = E h (u n-1 -u n ) + E h (u n+1 -u n ) = E h (u n+1 + u n-1 -2 u n ) (13)
where u n designates the nodal displacement of the node n. By analogy with the discrete model and using [START_REF] Belytschko | Meshless methods : an overview and recent developments[END_REF], the dispersion relation of the discretized form of the continuum model can be obtained:

ω2 C = 4 C r 2 m h 2 sin 2 k h 2 (14)
Using ( 7) and ( 14) the cutoff frequencies of the discrete and the discretized continuum models can be obtained:

f c D = ωmax D 2π = 1 π C m (15)
f c C = ωmax C 2π = 1 π C m r h (16)
The cutoff frequency of the discretized continuum model depends on h. In the case of coarse discretization (h ≫ r), the associated cutoff frequency ( f c C ) becomes smaller than that of the discrete one ( f c D ). As a result, if these models are coupled, the high frequency waves (f > f c C ) coming from the discrete model are not supported by the continuum model, and will be spuriously reflected at the interface (Fig. 11). This phenomenon has already been addressed using the finite element model with different element sizes [START_REF] Celep | Spurious reflection of elastic waves due to gradually changing finite element size[END_REF]. Fig. 11: Reflection of high frequency waves at the interface between the discrete and continuum models Taking into account the assumption k = κ > 0, the phase and group velocities in both the discrete model and discretized continuum model can be obtained form ( 7) and ( 14) as:

v ph D = ω D k = 2 κ C m sin κ r 2 and v gr D = ∂ω D ∂k = r C m cos κ r 2 ( 17
)
ṽph C = ω C k = 2 κ r h C m sin κ h 2 and ṽgr C = ∂ω C ∂k = r C m cos κ h 2 (18)
The use of different discretization characteristic lengths in the discrete and continuum models (r and h) leads to different wave velocities. Moreover, for a fixed κ, the wave velocities in the discretized continuum model decreases as h increases. This explains another mechanism of wave reflections, for which a portion of the main propagating wave is reflected at the interface (not only the high frequency waves). Figure 12 presents examples of a wave traveling between two models having different wave propagation velocities, due to different spatial discretization. influence of the temporal discretization on the dispersion relations will be examined. For a given time step ∆t, the second derivative of a field variable f with respect to time can be approximated as:

f ≈ f t+∆t -2 f t + f t-∆t (∆t) 2 (19)
Applying [START_REF] Benson | Computational methods in Lagrangian and Eulerian hydrocodes[END_REF] on the particles and nodes accelerations ( dn and ün , respectively) and using the equivalence conditions [START_REF] Belytschko | Meshless methods : an overview and recent developments[END_REF], Equations ( 5) and ( 13) can be approximated as:

m ∆t 2 D (d p+1 n + d p-1 n -2 d p n ) = C (d p n+1 + d p n-1 -2 d p n ) (20)
m ∆t 2 C (u p+1 n + u p-1 n -2 u p n ) = C r 2 h 2 (u p n+1 + u p n-1 -2 u p n ) (21)
where ∆t D and d p n are respectively the discrete time step and the displacement of a particle n at the p -th time step, and ∆t C and u p n are respectively the continuum time step and the displacement of a node n at the p -th time step. By analogy with the derivation of the previous dispersion relations and using the general form of the harmonic solutions, expressed at the p -th time step on the element (particle or node) n as d p n = D e i (k n r-ω p ∆t D ) and u p n = U e i (k n h-ω p ∆t C ) , the following equations can be obtained:

sin 2 ω D ∆t D 2 = C m ∆t 2 D sin 2 k r 2 (22)
sin 2 ωC ∆t C 2 = C m r 2 h 2 ∆t 2 C sin 2 k h 2 (23)
Using k = κ > 0 , equations ( 22) and ( 23) can be rewritten as:

sin ω D ∆t D 2 = C m ∆t D sin κ r 2 ( 24
)
sin ωC ∆t C 2 = C m r h ∆t C sin κ h 2 (25)
With the help of the first order and the third order Taylor expansions, an approximation of ω D and ωC can be obtained from ( 24) and ( 25) as follows:

ω D = C m k r 1 + k 2 r 2 24 C m ∆t 2 D -1 (26) ωC = C m k r 1 + k 2 h 2 24 C m r 2 h 2 ∆t 2 C -1 (27)
It can be concluded from these last equations that, for different discretization characteristic lengths (r and h) and for a given discrete time step ∆t D , it is possible to choose a suitable ∆t C such that the dispersion relations of both models are equivalent:

∆t C = ∆t 2 D + m C h 2 -r 2 r 2 (28)
Therefore, the inaccurate continuum dispersion relation caused by large discretization characteristic length h can theoretically be corrected by choosing larger continuum time step ∆t C . However, this solution is not always applicable in practice. For example in the case of "dynamic explicit" simulations, ∆t C must satisfy the CFL stability criterion [START_REF] Cundall | Computer Model for Simulating Progressive Large Scale Movements in Blocky Rock Systems[END_REF] to ensure convergence of the numerical approximation.

∆t C ≤ 2 π ωmax C (29)

Coupling techniques

Several techniques used to bridge discrete and continuum methods in a concurrent manner can be found in the literature. Some of them are detailed hereafter.

Edge-to-Edge coupling methods

A common technique for coupling discrete and continuum models is the edge-to-edge approach. This technique is also widely used in the finite element community for parallel implementation of the FEM method [119]. As shown in Figure 13, two types of discrete particles are defined in this technique. Aside from the real ones, virtual particles are defined to model the bond angle-bending for bonds (virtual)

between the discrete model and the continuum model. The virtual particles are connected on one side to the discrete model by virtual bonds, and on the other side are attached to the continuum model and move with it. The motion of these particles can be determined by interpolation of the continuum nodes motion. The most simple variation of this technique is that in which the particles (discrete elements)

and the continuum nodes on the interface are coincident and constrained to move together. This variation is known in the literature as direct coupling approach. In the general case, the coincidence of particles and nodes on the interface is not necessary. In this case, additional handling is required to ensure correct communication between the coupled models. This last point will be discussed in more detail below.

A domain Ω G is considered with boundary ∂Ω = ∂Ω u ∪ ∂Ω t (∂Ω u ∩ ∂Ω t = φ), such that essential (displacement) and natural (traction) boundary conditions are respectively prescribed on ∂Ω u and ∂Ω t . This domain is divided into two adjacent subdomains Ω C and Ω D , which are respectively modeled using the continuum approach and the discrete approach (Fig. 13). An isotropic linear elastic behavior and small deformations are assumed for simplicity. In the present development, an energetic (Hamiltonian) approach will be used to derive the governing equations of the coupling system. u g t g Fig. 13: Edge-to-edge coupling

The total energy in the discrete subdomain Ω D , assumed to be isolated from the continuum subdomain at this stage, is known as the Hamiltonian and is given by:

H D = E kinetic D + E int D -W ext D = n D i=1 1 2 m i p i D 2 - n D i=1 f int i d i + E virtual D - n D i=1 f ext i d i ( 30
)
where E kinetic D is the kinetic energy, E int D is the internal energy, W ext D is the external work, n D is the total number of particles, m i and p i D = m i ḋi are the mass and momentum of the particle i, d i and ḋi are respectively the displacement and velocity of the particle i, f int i is the total internal force exerted on the particle i by its neighbors, f ext i is the total external force acting on the particle i, and E virtual D is the bending potential of the virtual bonds. In some atomistic methods, E virtual D is simply a function of the bond angle change between the virtual bonds and their adjacent ones (Fig. 14). It should be noted that the stretching energy of the virtual bonds is automatically included in the continuum internal energy, since the virtual particles move with the continuum material. At this stage, f ext i does not include the coupling force which will be introduced later.

Fig. 14: Example of bending potential computation

Under the assumption of small deformations, the total energy in the continuum subdomain, supposed to be isolated at this stage, can be written:

H C = E kinetic C + E int C -W ext C = ˆΩC 1 2 ρ u2 dΩ + 1 2 ˆΩC σ : ε dΩ -´ΩC ρ b u dΩ -´∂Ω t C t g u dΓ (31)
where ρ is the density and u is the velocity field, σ and ε are respectively the Cauchy stress tensor and the strain tensor, b is the body force per unit mass, t g is the prescribed traction vector on ∂Ω t C . In the above development, the discrete and continuum energies are derived without taking into account the coupling conditions. In the following, it will be shown how these models are coupled.

Generally, the compatibility at the interface between the coupled models can be enforced using either velocity constraints or displacement constraints. However, in the case of highly dynamic problems, where dynamic effects become significant, it would be preferable to use velocity constraints to ensure a correct estimate of the kinetic energy at the interface. For the sake of consistency with the energetic (Hamiltonian) approach followed here to develop the edge-to-edge coupling method, displacement constraints will be used in the rest of this derivation. These constraints can be applied in different ways, as shown in [START_REF] Delaunay | Sur la sphère vide. A la mémoire de Georges Voronoï[END_REF] and [START_REF] Dobson | Analysis of a force-based quasicontinuum approximation[END_REF].

g weak i = u(x i ) -d i 2 = 0, i ∈ [1, n DI ] (32)
g strong i = u(x i ) -d i = 0, i ∈ [1, n DI ] (33)
u is the continuum displacement field, d i is the displacement of a particle i located on the interface at x i coordinates, and n DI is the number of particles on the interface. In [START_REF] Delaunay | Sur la sphère vide. A la mémoire de Georges Voronoï[END_REF], a single constraint is applied per particle belonging to the interface. This results in a system of equations with reduced number of degrees of freedom. However, using this equation, displacement compatibility is only satisfied approximately, which can lead to unacceptable large errors. In contrast, three constraints (one for each component) per particle are applied in [START_REF] Dobson | Analysis of a force-based quasicontinuum approximation[END_REF], which allows the coupled models to better communicate. Although it can lead to large system of equations, ( 33) is generally recommended and is retained here to enforce compatibility between the discrete and continuum approaches. In practice, these constraints can be introduced in the global system by using rigid models, such as the Lagrangian multipliers model (LM), elastic models such as the penalty model (PM), or a combination of these models such as the Augmented Lagrange Multipliers model (ALM). In the rest of this derivation, the ALM is used, because it is more general and includes the two other models. The associated coupling energy, in a continuous form, is given by:

H ALM I = ˆΓI λ (u -d) dΓ + p 2 ˆΓI (u -d) 2 dΓ ( 34
)
where Γ I is the interface between the two models, λ is the Lagrange multipliers field which can be regarded as a generalized coupling force field and p is the penalty parameter. Since the displacements d i in the discrete subdomain Ω D are only defined at the particle positions, to evaluate H ALM D (34), a continuous field d must be inferred, at least on the interface, from the particle displacements d i . This can be achieved using a kind of interpolation:

d(x i ) = n DI j=1 φ j D (x i ) d j (35)
where φ j D is the interpolation function associated to the particle j. Taking into account [START_REF] Dobson | Accuracy of quasicontinuum approximations near instabilities[END_REF], the total energy of the coupling system can be obtained:

H ALM G = H D + H C + H ALM I (36)
In [START_REF] Donzé | Advances in Discrete Element Method applied to Soil, Rock and Concrete Mechanics[END_REF], the Hamiltonian is given in a continuous form. To develop the corresponding discrete form HALM G , the continuum displacement and the Lagrange multipliers fields (u and λ, respectively) can be approximated, using for example a finite element method in the corresponding domains, as follows:

ũ(x) = n C i=1 φ i C (x) u i (37) λ(x) = n I i=1 φ i I (x) λ i ( 38
)
where φ i C and φ i I are the shape functions constructed respectively on ΩC and ΓI (the discretized forms of respectively Ω C and Γ I ), n C and n I are respectively the total number of continuum nodes and the total number of Lagrange multipliers nodes. In the general case, the Lagrange multipliers nodes do not necessary coincide with the discrete particles or the continuum nodes at the interface. Substituting [START_REF] Huanga | A dynamic atomistic-continuum method for the simulation of crystalline materials[END_REF] and ( 38) into (36) and using [START_REF] Dolbow | Numerical integration of the Galerkin weak form in meshfree methods[END_REF] on Γ I , the discrete Hamiltonian HALM G can be obtained, and then the discretized Hamiltonian equations of the coupling system can be derived:

ṗi D = - ∂ HALM G ∂d i , and ḋi = ∂ HALM G ∂ ṗi D for i ∈ [1, n D ] (39)
ṗi C = - ∂ HALM G ∂u i and ui = ∂ HALM G ∂ ṗi C for i ∈ [1, n C ] (40)
g i = ∂ HALM G ∂λ i = 0 for i ∈ [1, n I ] (41)
Replacing ḋi and ui with their expressions in respectively ṗi D and ṗi C , the global system of equations can be expressed, in matrix form, as:

m d = f int + f ext + f λ (42)
M ü = F int + F ext -F λ (43)
C I u -c I d = 0 ( 44
)
where m is the diagonal mass matrix of the particles, f int , f ext and f λ are the vectors of respectively the internal forces, the external forces and the coupling forces in the discrete model, M is the lumped mass matrix of the nodes, F int , F ext and F λ are the vectors of respectively the internal forces, the external forces and the coupling forces in the continuum model, C I and c I are the continuum and discrete coupling matrices.

The matrix systems (42), ( 43) and ( 44) are in the form of time-dependent ordinary differential equations (ODEs). These equations can be solved numerically using a time integration scheme. Many types of time integration schemes can be found in the literature. They can be classified into three headings: explicit, implicit and semi-implicit (or also semi-explicit). The explicit schemes calculate the state of a system at the current time step from its state at the previous time step. The implicit schemes use both the previous and current states to find the solution at the current time step. The semi-implicit schemes combine some features of the explicit and implicit schemes to find the solution. In dynamic studies, the explicit schemes are generally more recommended to capture the fine scale effects and ensure a proper momentum transfer. Many explicit schemes can be found in the literature, such as the Runge-Kutta, position and velocity Verlet schemes. A comparison between these schemes can be found in reference [START_REF] Rougier | Numerical comparison of some explicit time integration schemes used in DEM, FEM/DEM and molecular dynamics[END_REF]. According to this reference, the velocity Verlet scheme provides good results.

Moreover, it is easier to implement. Therefore, this scheme is a good candidate to solve the above matrix systems. This scheme gives an O ∆t 3 approximation for both velocities and displacements.

Thus, the choice of a displacement coupling (control of displacements at the interface) does not affect the coupling approach accuracy compared with a velocity coupling (control of velocities at the interface).

The edge-to-edge coupling methods are rather dedicated for quasi-static problems or problems with relatively low dynamic effects. Otherwise, fine discretization of the continuum subdomain, at the same level as the discrete subdomain, is required to deal with the spurious wave reflections at the interface.

Moreover, in most cases, this (use of fine continuum discretization) is not sufficient to overcome the reflections problem. Mostly, the particles belonging to the interface must be coupled with dampers, which are generally difficult to adjust. These difficulties make the edge-to-edge coupling approaches useless for dynamic studies, since the reduction in the computation effort from a fully discrete approach is not significant.

Bridging domain coupling methods

In the following, the bridging (overlapping) domain (BD) technique for coupling continuum models with discrete models is described. Contrary to the edge-to-edge coupling technique, in the present one, a bridging (overlapping) domain is considered between the coupled models. It should be noted that this technique is quite similar to the Arlequin approach developed by Ben Dhia et al. [START_REF] Ben Dhia | Problèmes mécanique multi-échelles: la méthode Arlequin[END_REF][START_REF] Ben Dhia | Analyse mathématique de la méthode Arlequin mixte[END_REF][START_REF] Ben Dhia | The Arlequin method as a flexible engineering design tool[END_REF].

They are based on practically the same mechanical ingredients which are:

-Decomposition of the global domain Ω G into two subdomains Ω D and Ω C , modeled respectively by the discrete model and the continuum model, with a bridging region Ω B .

-Weak coupling (based on the weak formulation): this second point states that, to couple the superposed models, the discrepancy between the mechanical states, e.g. displacement, deformation, strain, etc, must be controlled using some kinds of fictive forces. However, to allow each model to express its own wealth, the discrepancy should preferably be controlled in a weak manner using averaging operators. This point consists of:

-Definition of a gluing zone Ω GL (Ω GL ⊆ Ω B ), in which the coupling control will be performed.

In the general case, this zone can be different from the bridging region (Ω GL Ω B ). In the present development, Ω GL is chosen the same that Ω B . Hereafter, the term "bridging zone" will be used to designate both the bridging zone and the gluing zone.

-Definition of a mediator space M which can be defined as a space of the field variables, restricted to Ω B , that must be controlled in the bridging region. To ensure a correct dialogue between the coupled models, the control quantities must be carefully chosen. Generally, either the displacements or velocities are controlled in Ω B .

-Definition of a projection operator Π, which will be used to project the discrete and continuum field variables to be controlled, onto the mediator space M.

-Definition of a junction model, which will be used to ensure the compatibility of the controlled field variables in the bridging region. Besides the Lagrange, penalty and Augmented Lagrange Multipliers Methods, which are introduced earlier, another junction model has been proposed in the literature [START_REF] Bauman | On the application of the Arlequin Method to the Coupling of Particle and Continuum Models[END_REF][START_REF] Ben Dhia | Analyse mathématique de la méthode Arlequin mixte[END_REF][START_REF] Ben Dhia | Problèmes mécanique multi-échelles: la méthode Arlequin[END_REF][START_REF] Ben Dhia | The Arlequin method as a flexible engineering design tool[END_REF][START_REF] Ben Dhia | Further Insights by Theoretical Investigations of the Multiscale Arlequin Method[END_REF]:

< λ, f C -f D >= ˆΩB η 1 λ • (Πf C -Πf D ) + η 2 l 2 ε(λ) : ε(Πf C -Πf D ) dΩ (45)
where (Πf C -Πf D ) is the difference between the projected continuum and discrete control quantities on the mediator space M, λ is the Lagrange multiplier field and l is a parameter which is named "junction parameter" in this paper. This parameter which has the dimension of a length is added to ensure the homogeneity of the integral terms in [START_REF] Guidault | On the L2 and the H1 couplings for an overlapping domain decomposition method using Lagrange multipliers[END_REF]. η 1 and η 2 are non-negative weight parameters. These parameters can be chosen so as to scale the two integral terms in [START_REF] Guidault | On the L2 and the H1 couplings for an overlapping domain decomposition method using Lagrange multipliers[END_REF]. (η 1 , η 2 ) = (1, 0) refers to L 2 (Ω B ) inner (scalar) product which is the same that the Lagrange multiplier model, (η 1 , η 2 ) = (1, 1) refers to H 1 (Ω B ) inner (scalar) product and

(η 1 , η 2 ) = (0, 1) refers to H 1 (Ω B ) semi-inner product.
-Partition of energy between the discrete and continuum models in the bridging zone. The two models coexist in Ω B . Therefore, the discrete and continuum energies in this region must be weighted using a kind of unity partition functions (to avoid counting twice the associated energy).

Different weight functions can be used for the kinetic energy, the internal energy and the external work in Ω B . These functions must verify

f C : Ω G → [0, 1] x →      1 in Ω C \Ω B [0, 1] in Ω B 0 in Ω D \Ω B ( 46
)
in the continuum subdomain, and f D = fC = 1 -f C in the discrete subdomain. Figure 15 presents examples of weight functions.

Fig. 15: Examples of weight functions

These different ingredients will be detailed hereafter. As for the edge-to-edge technique, a domain Ω G is considered with boundary ∂Ω = ∂Ω u ∪ ∂Ω t (∂Ω u ∩ ∂Ω t = φ), such that essential (displacement)

and natural (traction) boundary conditions are respectively prescribed on ∂Ω u and ∂Ω t . This domain is divided into two subdomains Ω D , treated by the discrete model, and Ω C , treated by the continuum model, with a bridging region Ω B (Fig. 16). An isotropic linear elastic behavior and small deformations are assumed for simplicity.

To weight the energies in the continuum subdomain, three weight functions are assumed in this development: α for the internal energy, β for the kinetic energy and γ for the external work. In a complementary manner, the energies in the discrete subdomain are weighted using ᾱ, β and γ, for respectively the internal energy, the kinetic energy and the external work. Using these functions, the u g t g Fig. 16: Bridging domain coupling weighted discrete and continuum Hamiltonians can respectively be written:

H w D = E kinetic, w D + E int, w D -W ext, w D = n D i=1 1 2 m i βi p i D 2 - n D i=1 ᾱi f int i d i - n D i=1 γi f ext i d i ( 47
)
H w C = E kinetic, w C + E int, w C -W ext, w C = ˆΩC 1 2 β ρ u2 dΩ + 1 2 ˆΩC α σ : ε dΩ -´ΩC γ ρ b u dΩ -´∂Ω t C γ t g u dΓ (48)
To enforce compatibility between the coupled models in the bridging zone Ω B , the

H 1 (Ω B ) junction model ((45), with (η 1 , η 2 ) = (1, 1
)) is used. The associated coupling energy can be written:

H H1 B = ˆΩB λ • (Πu -Πd) + l 2 ε(λ) : ε(Πu -Πd) dΩ (49)
As explained above, to evaluate H H1 B , a continuous displacement field d must be approximated from the particle displacements d i in the bridging zone, using a kind of interpolation [START_REF] Dolbow | Numerical integration of the Galerkin weak form in meshfree methods[END_REF]. The global Hamiltonian of the coupling system can be obtained by summing [START_REF] Hall | The Boundary Element Method (Solid Mechanics and Its Applications)[END_REF], ( 48) and ( 49):

H H1 G = H w D + H w C + H H1 B ( 50
)
To obtain the corresponding discrete form HH1 G , the continuum displacement field u and the Lagrange multipliers field λ are approximated as given by ( 37) and [START_REF] Felici | A coupled Eulerian/Lagrangian method for the solution of three-dimensinal vortical flows[END_REF]. Using HH1 G , the global system of equations can be derived:

m β d = f int α + f ext γ + f λ (51) M β ü = F int α + F ext γ -F λ ( 52
)
C B u -c B d = 0 ( 53
)
where m β is the weighted diagonal mass matrix of the particles, f int α , f ext

β : Ω G → [0, 1] x →      1 in Ω C \Ω B [ε, 1 -ε] in Ω C|Ω B 0 in Ω D \Ω B ( 54
)
The global system of equations consisting of ( 51), ( 52) and ( 53) can finally be solved using a temporal integration scheme.

In the bridging domain (BD) coupling methods, the fine scale solution (of the discrete model) is continuously projected onto the coarse scale solution (of the continuum model) in the bridging zone.

This can cancel the high frequency waves (HFWs), and then avoids spurious wave reflections at the discrete/continuum interface, without any additional filtering or damping [START_REF] Ben Dhia | The Arlequin method as a flexible engineering design tool[END_REF][START_REF] Jebahi | 3D coupling approach between discrete and continuum models for dynamic simulations (DEM-CNEM)[END_REF]. However, this is conditional upon choosing correctly the coupling parameters, which is not a straightforward issue due to the large number of these parameters. Several studies aiming to alleviate this issue can be found in the literature. These studies have been extended in a previous paper [START_REF] Jebahi | 3D coupling approach between discrete and continuum models for dynamic simulations (DEM-CNEM)[END_REF] to dynamics using generalized 3D models, in which no coincidence conditions are imposed between the continuum nodes, the Lagrange multipliers nodes and the particles in the bridging zone. This last study [START_REF] Jebahi | 3D coupling approach between discrete and continuum models for dynamic simulations (DEM-CNEM)[END_REF] will be resumed and extended in section 4 to better investigate the effects of the main coupling parameters and to give recommendations and trends simplifying their setting in practice. Another interesting aspect of the BD coupling methods is that they can also be applied to obtain solution for nonlinear problems, in which devising an energy functional for the entire system is not possible, e.g. as in the presence of irreversible processes. In this context, Jacob et al. [START_REF] Fish | Concurrent AtC coupling based on a blend of the continuum stress and the atomistic force[END_REF] formulated an atomistic-continuum coupling method based on a blend of the continuum stress and the atomistic force.

Bridging scale coupling methods

The bridging scale coupling method was recently proposed by Wagner and Liu [110] to concurrently couple Molecular Dynamics (MD) and Finite Element Method (FEM). The fundamentals of this approach are briefly reviewed in this subsection, the reader can refer to [START_REF] Park | The bridging scale for two-dimensional atomistic/continuum coupling[END_REF][START_REF] Parka | Three-dimensional bridging scale analysis of dynamic fracture[END_REF]110] for more detail.

Two subdomains Ω D and Ω C , respectively modeled by discrete and continuum models, are considered.

To better understand the concepts of the bridging scale approach, these subdomains are first assumed to be completely superimposed in the global domain Ω G .

The key idea of such an approach is to decompose the total displacement field u at the particle positions into coarse and fine scales:

u T ot (x i ) = ū(x i ) + u ′ (x i ) (55)
where x i represents the position of a particle i, ū and u ′ are respectively the coarse scale and fine scale displacements (Fig. 17).

u u u Tot

Fig. 17: Scale decomposition of total displacement at the particle positions After discretization of the continuum subdomain Ω C , the coarse scale displacement at the particle positions x j is evaluated form the continuum nodal displacements as:

ū(x j ) = n C i=1 φ i C (x j ) u i (56)
where u i and φ i C are the continuum nodal displacement and the shape function associated to the node i, respectively. For the sake of clarity, it is convenient to employ a matrix representation for subsequent developments. In matrix form, (56) can be rewritten as:

ū = N u ( 57
)
where ū = {ū(x 1 ) ū(x 2 ) ... ū(x n D )} t is the vector of coarse displacements of the particles, n D is the total number of particles, u = {u 1 u 2 ... u n C } t is the vector of continuum nodal displacements, n C is the total number of nodes and N is the interpolation matrix.

The fine scale displacement is defined as the part of the total displacement that cannot be represented by the coarse scale. This quantity is computed from the vector d of particle displacements, which include the coarse scale parts, by subtracting (from d) the projection of the discrete solution onto the continuum solution.

u ′ = d -P d = Q d ( 58
)
where P is the projection matrix and Q = I -P is the complementary projector [110]. The term P d, which is called "bridging scale", is the part of the particle displacement that must be subtracted from the total displacement to completely separate the scales (i.e. the coarse and fine scales are orthogonal or linearly independent of each other). The total displacement u can finally be written as the sum of the coarse and fine scales as:

u T ot = N u + Q d (59)
To derive the coupled discrete and continuum equations, it is most convenient to adopt the Lagrangian approach, in which the multiscale Lagrangian is defined by:

L(u T ot , uT ot ) = E kinetic ( uT ot ) -E int (u T ot ) + W ext (u T ot ) = 1 2 ut M u + 1 2 ḋt m f ḋ -E int (N u + Q d) + f ext t (N u + Q d) (60)
where M and m f are respectively the coarse scale and fine scale diagonal mass matrices (defined in terms of the diagonal mass matrix of the particles m as:

M = N t m N and m f = Q t m Q = Q t m = m Q)
, and f ext is the vector of external forces acting on the particles. An important feature of the Lagrangian ( 60) is the absence of quadratic terms in the kinetic energy, which are canceled due to the presence of the "bridging scale". The coupled multi-scale equations of motion can obtained from L as follows:

d dt ∂L ∂ u - ∂L ∂u = 0 and d dt ∂L ∂ ḋ - ∂L ∂d = 0 (61)
Using ( 60), these equations can be rewritten after simplification as:

m d = f int + f ext (62)
M ü = N t f int + f ext (63)
The first (fine scale) equation is the equation of motion in the discrete model, m is the diagonal mass matrix of the particles. The second (coarse scale) equation is simply the equation of motion in the continuum model, while noting that the mass matrix M is consistent. The coupling between these equations is through the vector of the continuum internal forces F int = N t f int which is a function of the vector of the discrete internal forces (f int ).

In the above development, Ω D and Ω C are assumed to be completely superimposed. Therefore, the continuum equation of motion ( 63) is only an approximation of the first one [START_REF] Li | Perfectly matched multiscale simulations for discrete lattice systems: Extension to multiple dimensions[END_REF] and then is redundant. The total displacement and discrete displacements (u T ot and d, respectively) satisfy the same equation of motion [START_REF] Li | Perfectly matched multiscale simulations for discrete lattice systems: Extension to multiple dimensions[END_REF] and have the same initial conditions. Therefore, these quantities are identical for all time, and consequently, the continuum displacement u can be determined by [START_REF] Kaljevic | An improved element free Galerkin formulation[END_REF], which imply: N u = P d.

Assuming that the fine scale analysis is only required in a small region Ω B , the particles outside this zone will now be removed. In other words, the unnecessary fine scale degrees of freedom will be eliminated. By doing so, the discrete model is reduced to Ω B , where the two models coexist, and the remaining estate is only modeled by the continuum model (Fig. 18). The process of eliminating the unnecessary degrees of freedom results in a modified equation of motion in the discrete model, including an external force, called impedance force. This force is a function of a damping matrix or equivalently its time derivative, known as the time history kernel [START_REF] Adelman | Generalized Langevin equation approach for atom/solid-surface scattering: General formulation for classical scattering off harmonic solids[END_REF]110]. The damping kernel was first derived analytically by Adelman and Doll [START_REF] Adelman | Generalized Langevin equation approach for atom/solid-surface scattering: General formulation for classical scattering off harmonic solids[END_REF] for a harmonic 1D lattice. However, their analytical approach is very costly and is intractable above one dimension [START_REF] Parka | Three-dimensional bridging scale analysis of dynamic fracture[END_REF]. To overcome this limitation, other works [START_REF] Cai | Minimizing boundary reflections in coupled-domain simulations[END_REF][START_REF] Huanga | A dynamic atomistic-continuum method for the simulation of crystalline materials[END_REF] have proposed numerical approaches to evaluate this quantity, whose the effect is to dissipate the high frequency waves and to avoid wave reflections at the interface between the discrete and continuum models. In the first implementation of Wagner and Liu [110], the projection of the fine scale solution onto the coarse scale one is apparently only performed at the interface using the approach of Adelman and Doll [START_REF] Adelman | Generalized Langevin equation approach for atom/solid-surface scattering: General formulation for classical scattering off harmonic solids[END_REF] (without any Lagrange multipliers). This approach has been improved by computing the interface nodes using Green's functions for lattices [START_REF] Park | The bridging scale for two-dimensional atomistic/continuum coupling[END_REF][START_REF] Parka | Three-dimensional bridging scale analysis of dynamic fracture[END_REF]. More recently, other enhancements have been proposed to join the coarse and fine scales in a bridging domain, using the perfectly matched layer (PML) technique [START_REF] Li | Perfectly matched multiscale simulations for discrete lattice systems: Extension to multiple dimensions[END_REF]106]. As shown in reference [114], the performance of the bridging scale (BS) methods based on PML is comparable to that of bridging domain (BD) methods. However, the BD methods would be less costly, since they allow to avoid spurious wave reflections by only a correct setting of their parameters (no additional treatments are required).

Other coupling techniques

Quasicontinuum method The quasicontinuum (QC) method was originally developed in the context of lattice statics at zero temperature (Molecular Mechanics) using empirical interatomic potentials [102,103]. Later on, it was extended to dynamic problems using the coarse grained energy and the Hamilton's principle [START_REF] Li | On the Effect of Ghost Force in the Quasicontinuum Method: Dynamic Problems in One Dimension[END_REF][START_REF] Rodney | Mixed atomistic/continuum methods: static and dynamic quasicontinuum methods[END_REF][START_REF] Shenoy | Quasicontinuum models of atomic-scale mechanics[END_REF]. The chief objective of the theory is to systematically coarsen a particle (atomistic) description by introducing kinematic constraints. These constraints are selected and designed so that the fully atomistic model is preserved in the regions of fine scale effects and large numbers of particles (atoms) are collectively treated in the coarse scale regions. The fully discrete model is then represented by a set of representative particles (mainly located in the fine scale region),

where the equation of motion has to be solved. The displacement of the non-representative particles are determined by kinematic constraints based on finite element meshing, whose the nodes are coincident with the representative particles (Fig. 19). Although the imposed kinematic constraints significantly reduce the number of degrees of freedom (particles), the computational complexity of evaluating the generalized forces corresponding to the coarse-grained degrees of freedom (representative particles) still scales with the total number of particles in the system making computation on large systems intractable. Various approximations have been suggested to further reduce the complexity of forces computation [START_REF] Iyer | A field theoretical approach to the quasi-continuum method[END_REF][START_REF] Li | On the Effect of Ghost Force in the Quasicontinuum Method: Dynamic Problems in One Dimension[END_REF][START_REF] Ortiz | An analysis of the quasicontinuum method[END_REF][START_REF] Rodney | Mixed atomistic/continuum methods: static and dynamic quasicontinuum methods[END_REF][START_REF] Shenoy | Quasicontinuum models of atomic-scale mechanics[END_REF]. These include the mixed discrete (atomistic) and continuum formulations, or introduction of cluster summation rules on lattice sums. However, these approximations induce spurious forces, which can affect the solution accuracy. As shown in the literature [START_REF] Li | On the Effect of Ghost Force in the Quasicontinuum Method: Dynamic Problems in One Dimension[END_REF], the effects of theses forces are more significant in the case of dynamic studies. These effects may even be more severe than the spurious reflections at the interface.

Many strategies have been suggested to correct the errors incurred in these approximations [START_REF] Shimokawa | Matching conditions in the quasicontinuum method: Removal of the error introduced at the interface between the coarse-grained and fully atomistic region[END_REF]111],

but they introduce undesirable seams in the process. Further, recent numerical analysis suggests that the approximations introduced may not be consistent and stable, and can result in uncontrolled errors for rapid coarse-graining [START_REF] Dobson | Analysis of a force-based quasicontinuum approximation[END_REF][START_REF] Dobson | Accuracy of quasicontinuum approximations near instabilities[END_REF]. More recently, other works have proposed seamless quasicontinuum formulations. However, they generally suffer from a lack of systematic convergence [START_REF] Iyer | A field theoretical approach to the quasi-continuum method[END_REF].

ONIOM method

The ONIOM (our own N-layered integrated molecular orbital+molecular mechanics) approach was first developed by Svensson et al.

[100] to couple quantum mechanical methods and atomistic methods. Using this approach, the fine scale model is superimposed on the coarse scale model in the regions where high accuracy is required. Communication of the coupled models is ensured by another model (fictitious) added to cancel the coarse scale energy in the overlapping region. The ONIOM approach can also be used to couple discrete and continuum methods. Figure 20 gives a 1D illustration for a discrete-continuum coupling based on this approach. As can be seen, three models are involved in the zone of fine scale effects Ω B : the discrete model (DM), the continuum model (CM) and a fictitious model (DC) that is used to counteract the continuum energy. In Ω B , both the CM and DC models are fictitious and are only retained to simplify the coupling. The total internal energy as treated by ONIOM is given by:

E int = E int D + E int C -E int DC (64)
Fig. 20: ONIOM approach for coupling of discrete and continuum models

The discrete and continuum models are linked at the interface through dangling bonds (Fig. 20), which are connected on one side to ordinary particles at the boundary of the discrete model, and on the other side to fictitious link particles attached to the continuum model. These bonds represent a major source of error in the ONIOM methods, since their behavior can be polluted by incorrect placement of the link particles. This problem would even be more pronounced in the case of dynamic studies.

Coupling of discrete and continuum meshless methods In this paragraph, techniques used to couple discrete and continuum meshless methods are briefly reviewed, the reader is referred to [START_REF] Liu | Coupling meshfree particle method with molecular dynamics -a novel approach for multi-scale simulations[END_REF] for more detail. These techniques were first developed to couple Smoothed Particle Hydrodynamics (SPH)

and Molecular Dynamics (MD). As seen in section 2, in meshless methods, the studied domain is represented by a set of scattered particles (continuum particles) or nodes without any connectivity between them. Each continuum particle has its corresponding smoothing length R, representing the influence domain, and length scale h, characterizing locally the domain discretization. The idea is that the continuum particles close to the interface (transitional particles) interact on one side with their neighboring continuum and discrete particles (belonging the their influence domains), and act on the other side as virtual discrete particles which interact with the neighboring ordinary and virtual discrete particles (Fig. 21). To avoid interface problems, the length scale of the transitional particles must be graded down to the order of the discretization characteristic length of the discrete model.

The communication between the coupled models can further be enhanced by using large transitional zone, in which the particles act both as continuum and discrete particles (Fig. 21). As seen in section 3, the bridging domain (BD) coupling methods are well adapted to study dynamic problems. However, the effectiveness of these approaches remains the hostage of a correct setting of their coupling parameters. Several works studying these parameters can be found in the literature. Ben Dhia and Rateau [START_REF] Ben Dhia | Problèmes mécanique multi-échelles: la méthode Arlequin[END_REF][START_REF] Ben Dhia | Analyse mathématique de la méthode Arlequin mixte[END_REF][START_REF] Ben Dhia | The Arlequin method as a flexible engineering design tool[END_REF] were the first to study mathematically these ingredients for continuumcontinuum coupling. Their study is still conceptually valid for continuum-discrete coupling. The main results concerning the well-posedness of the coupling problem are recalled hereafter. The weight function α, applied to weight the internal energy, must be strictly positive in the bridging zone. Without this condition the coercivity of the internal energy cannot be verified. Another significant result concerning the junction models is that for the discretized problem, contrary to the H1 coupling ((45), with (η 1 , η 2 ) = (1, 1)) which yields a well-posed problem, the L2 coupling (l = 0) can lead to an ill-conditioned system of equations, especially in the case of highly refined discretization. In the same context, Bauman et al. [START_REF] Bauman | On the application of the Arlequin Method to the Coupling of Particle and Continuum Models[END_REF] have studied another coupling model, H1 semi-norm (associated to the semi-inner product given by [START_REF] Guidault | On the L2 and the H1 couplings for an overlapping domain decomposition method using Lagrange multipliers[END_REF], with (η 1 , η 2 ) = (0, 1)). This model leads to a well-posed problem, but it does not constrain enough the continuum and discrete displacements in the bridging zone.

Other works [START_REF] Bauman | On the application of the Arlequin Method to the Coupling of Particle and Continuum Models[END_REF][START_REF] Guidault | On the L2 and the H1 couplings for an overlapping domain decomposition method using Lagrange multipliers[END_REF][START_REF] Guidault | Bridging domain methods for coupled atomistic-continuum models with L2 or H1 couplings[END_REF] have studied numerically these parameters using 1D static models. Guidault et al. [START_REF] Guidault | On the L2 and the H1 couplings for an overlapping domain decomposition method using Lagrange multipliers[END_REF][START_REF] Guidault | Bridging domain methods for coupled atomistic-continuum models with L2 or H1 couplings[END_REF] have noted that, for the L2 coupling, the weight function α must be continuous at the boundary of the bridging (gluing) zone ∂Ω B . Indeed, the use of a discontinuous weight function can cause undesirable free conditions at ∂Ω B . Concerning the choice of the mediator space, Ben Dhia et al. [START_REF] Ben Dhia | Further Insights by Theoretical Investigations of the Multiscale Arlequin Method[END_REF][START_REF] Ben Dhia | Analyse mathématique de la méthode Arlequin mixte[END_REF] mentioned that in the case of continuous domains, it is convenient to choose M = H 1 (Ω B ); however, it is very difficult to choose the finite approximation space M. To address this difficulty, several works [START_REF] Bauman | On the application of the Arlequin Method to the Coupling of Particle and Continuum Models[END_REF][START_REF] Guidault | On the L2 and the H1 couplings for an overlapping domain decomposition method using Lagrange multipliers[END_REF][START_REF] Guidault | Bridging domain methods for coupled atomistic-continuum models with L2 or H1 couplings[END_REF] proposed static numerical studies of M, using 1D models. The different configurations that were studied are similar to those presented in Figures 22a, 22b and22c. The static studies of Guidault et al. [START_REF] Guidault | On the L2 and the H1 couplings for an overlapping domain decomposition method using Lagrange multipliers[END_REF] show that: (i) in the case of a fine mediator space M (Fig. 22c), the response of the structure does not depend on the weight functions and a locking phenomenon takes place, i.e. the fine solution exactly conforms to the coarse solution in the bridging zone; (ii) in the case of a coarse mediator space M (Fig. 22b), the weight functions have an influence on the solutions such that the larger the weight function on the fine scale model, the smaller the maximum solution jump between the two models becomes. These works have been recently extended to 3D dynamic studies [START_REF] Jebahi | 3D coupling approach between discrete and continuum models for dynamic simulations (DEM-CNEM)[END_REF].

In the following, it is proposed to corroborate these dynamics studies to better investigate the influence of the main coupling parameters and to give a methodology simplifying their setting in practice. It will be demonstrated that some of the results proven in statics using 1D models are not valid in 3D dynamic simulations. To do this, the BD coupling method recently developed between the Discrete Element Method (DEM) and the Constrained Natural Element Method (CNEM) is used.

The derivation of this approach is briefly recalled hereafter, the reader is referred to [START_REF] Jebahi | 3D coupling approach between discrete and continuum models for dynamic simulations (DEM-CNEM)[END_REF] for more detail.

Concurrent DEM-CNEM coupling method

The coupling ingredients used in this approach are as follows. Velocity constraints are applied in the bridging zone to enforce the DEM-CNEM compatibility. These constraints are introduced in the coupling system using the H 1 (Ω B ) inner product ((45), with (η 1 , η 2 ) = (1, 1)). The same weight function is used to weight the different energies (α = β = γ). Contrary to what is presented in subsection 3.2, where only particle displacements are considered in the discrete model, the present development is performed using another variation of discrete methods, where both displacement and rotation of particles are considered [START_REF] André | Discrete element method to simulate continuous material by using the cohesive beam model[END_REF]. The variational approach of virtual power is used to derive the coupling equations. This approach explicitly introduces velocities in the coupling formulation, which makes easier the application of the velocity constraints.

In an isolated system of the discrete subdomain, the governing equations can be written:

for i = 1..n D and t ∈ [0, t end ], find d i , θ i , f int i , c int i ∈ R 3 × R 3 × R 3 × R 3
such that, given the initial conditions:

f ext i + f int i = m i di c ext i + c int i = I i θi (65)
where d i , θ i , m i and I i represent respectively the displacement, rotation, mass and mass moment of inertia of the particle i. f int Applying the principle of the virtual power in Ω D , (65) can be reformulated as follows:

find (d, θ, f , c) ∈ D ad × O ad × F ad × C ad such that ∀ δ ḋ, δ θ ∈ Ḋad,0 × Ȯad,0 : n D i=1 f ext i δ ḋi + n D i=1 f int i δ ḋi + n D i=1 c ext i δ θi + n D i=1 c int i δ θi = n D i=1 m i di δ ḋi + n D i=1 I i θi δ θi (66)
The admissible spaces are defined as:

-D ad = {d = d i (t), i = [1..n D ] ∀ t ∈ [0, t end ]} -O ad = {θ = θ i (t), i = [1..n D ] ∀ t ∈ [0, t end ]} -F ad = f = f int i (t), i = [1..n D ] ∀ t ∈ [0, t end ] -C ad = c = c int i (t), i = [1..n D ] ∀ t ∈ [0, t end ] -Ḋad,0 = ḋ = ḋi (t), i = [1..n D ] ; ḋ = 0 on ∂Ω u D ; ∀ t ∈ [0, t end ] -Ȯad,0 = θ = θi (t), i = [1..n D ] ; θ = 0 on ∂Ω u D ; ∀ t ∈ [0, t end ]
Under the assumption of isotropic linear elastic behavior and small deformations, the governing equations in the continuum subdomain Ω C , which is assumed to be isolated at this stage, can be written as:

∀ x ∈ Ω C (t) and t ∈ [0, t end ], find (u, σ) ∈ H 1 (Ω C ) 3 × L 2 (Ω C ) 6 such that, given the initial conditions:                div(σ) + ρb = ρü in Ω C σ = A : ε(u) ε(u) = 1 2 (∇u + ∇ t u) u = u g on ∂Ω u C σ.n = t g on ∂Ω t C ( 67
)
where ρ is the density, u is the continuum displacement field, σ is the Cauchy stress tensor, ε is the stain tensor, A is the stiffness tensor, b is the body force per unit mass, u g and t g are respectively the given (prescribed) displacement and traction on ∂Ω u C and ∂Ω t C , and n is the outward-pointing unit normal vector. Using the principle of the virtual power, the continuum weak formulation associated with (67) can be derived: find u ∈ U ad such that ∀ δ u ∈ U ad,0 , given the initial conditions:

ˆ∂Ω t C δ u • t g dΓ - ˆΩC ε(δ u) : A : ε(u) dΩ + ˆΩC ρ δ u • b dΩ = ˆΩC ρ δ u • ü dΩ (68)
where δ u represent the test functions. The admissible solution spaces, U ad and U ad,0 , are defined as follows:

-

U ad = u = u(x, t) ∈ [H 1 (Ω C )] 3 ; u = u g on ∂Ω u C ; ∀ t ∈ [0, t end ] -U ad,0 = u = u(x, t) ∈ [H 1 (Ω C )] 3 ; u = 0 on ∂Ω u C ; ∀ t ∈ [0, t end ]
Based on the previous weak formulations of the isolated discrete and continuum subdomains ((66) and ( 68)), the coupling conditions can be introduced to obtain the global weak formulation of the whole problem. By weighting (66) and ( 68) and using the H 1 (Ω B ) junction model, the global weighted weak formulation can be derived: find (u, d, θ, λ) ∈ U ad × D ad × O ad × M such that ∀ (δ u, δ ḋ, δ θ, δλ) ∈ U ad,0 × Ḋad,0 × Ȯad,0 × M, given the initial conditions:

ˆΩC α ρ δ u • ü dΩ + ˆΩC α ε(δ u) : A : ε(u) dΩ - ˆ∂Ω t C α δ u • tg dΓ - ˆΩC α ρ δ u • b dΩ + n D i=1 ᾱi m i di δ ḋi + n D i=1 ᾱi I i θi δ θi - n D i=1 (ᾱ i f ext i + ᾱi f int i ) • δ ḋi - n D i=1 (ᾱ i c ext i + ᾱi c int i ) δ θi + δ ˆΩB λ • (Π u -Π ḋ) + l 2 ε(λ) : ε(Π u -Π ḋ) dΩ = 0 (69)
To obtain the corresponding discrete equations, the admissible spaces U ad and M as well as their dual geometric spaces must be discretized. In this development, the dual geometric space of M (mediator space) is taken independent from Ω D|Ω B and Ω C|Ω B . Therefore, a third subdomain specific to the bridging region is introduced which will be designated hereafter by Ω B (bridging subdomain). Consequently, the discretization of this subdomain can be made independently of the continuum and discrete subdomains. This makes easier the application of the coupling method for complex problems. Moreover, as will be seen in subsection 4. numerically in a previous paper [START_REF] Jebahi | 3D coupling approach between discrete and continuum models for dynamic simulations (DEM-CNEM)[END_REF] and will be resumed in subsection 4.2.

After discretization of the different geometries, the discretized form of the associated (dual) admissible spaces D ad , U ad and M can be derived. Their discretized forms are denoted by Dad , Ũ ad and M, respectively. Using the CNEM interpolation on the different discretized admissible spaces, Dad , Ũ ad and M, the displacement fields d and u and the Lagrange multipliers λ can be approximated as:

d(x) = n D i=1 φ i D (x) d i (70) ũ(x) = n C i=1 φ i C (x) u i (71) λ(x) = n B i=1 φ i B (x) λ i (72)
M α ü = F int α + F ext α -F λ (74)
with:

F λ = C L2 C + l 2 C H1 C λ = C C λ. C L2
C and C H1 C are the continuum L2 and H1 coupling matrices.

-Bridging domain

C B . u -c B . d = 0 (75)
with:

C B = C L2 C t + l 2 C H1 C t and c B = c L2 D t + l 2 c H1 D t
The global system of equations consisting of ( 73), ( 74) and ( 75) is solved using the velocity Verlet scheme. At each time step, ignoring the coupling conditions, predictive accelerations and velocities are determined for the two models (DEM and CNEM). These quantities are then sent to a DEM-CNEM Interface code to be corrected taking into account the coupling forces. Finally, the corrected quantities are returned to the discrete and continuum codes to continue computation. This process is repeated until the simulation is achieved.

Parametric study of the coupling parameters

The main coupling parameters, which are studied in this subsection, are (Fig. 23):

-The junction parameter l (45).

-The weight function α.

-The width of the bridging region L B .

-The discretization of the approximated mediator space M (or its associated geometric discretization characteristic length h B ). ). To generate a state of high frequency waves in propagation within the beam and investigate the wave reflexion at the DEM-CNEM interface, the free end (x = L) is submitted to a tensile loading with a very steep slope (Fig. 25). As shown in Figure 25b, the Fourier spectrum contains powerful high frequency waves (greater than f c C ). The material of the beam is the silica, whose the mechanical properties are given in Figure 24. The corresponding microscopic properties of the cohesive beam bonds in the DEM approach are given in Table (3). To control the wave propagation within the model, four check points are placed along this beam (Fig. 24) as follows:

-CnemCheckPoint: at the middle of the CNEM subdomain where the controlled quantities are computed using the CNEM nodes in this zone -OverlapCnemCheckPoint: at the middle of the bridging zone where the controlled quantities are computed using only the CNEM nodes in this zone.

-OverlapDemCheckPoint: at the middle of the bridging zone where the controlled quantities are computed using only the DEM particles in this zone.

-DemCheckPoint: at the middle of the DEM subdomain where the controlled quantities are computed using the DEM particles in this zone. (1, 1)), is mainly to the global coupling matrix A defined as:

A = ∆t 2 (C B [M α ] -1 C C + c B [m α ] -1 c D ) ( 76
)
where ∆t is the computation time step. The influence of this parameter l on the coupling results will therefore be investigated by studying its effect on the conditioning of A:

Cond(A) = A . A -1 (77)
which is an important indicator of the results accuracy. A coarse mediator space M (coarse ΩB , i.e.

at the same order as ΩC ) and continuous weight functions with ε = 0.05 (ε is defined in [START_REF] Jean | The non smooth contact dynamics method[END_REF]) are used in this study. Figure 27 shows the conditioning of A with respect to l. The A conditioning decreases with l and reaches a minimum for l = l opt . Beyond this value, the conditioning increases exponentially as l increases. Therefore, the H1 coupling ((45), with (η 1 , η 2 ) = (1, 1)) for a small value of l is better than the L2 coupling (l = 0). However, contrary to what is often presented in the literature, the H1 coupling becomes worse if l exceeds some optimal value. A Fig. 27: Conditioning of A with respect to l L B = 6 mm, coarse mediator space and continuous weight functions with ε = 0.05 Table 5 presents the optimal values of the junction parameter (l opt ) determined numerically for different L B and h B (width and discretization characteristic length of ΩB , respectively). l opt practically does not change with the width of the bridging region, except for very fine mediator space (case of h B = 32 µm). However, this parameter depends significantly on h B . In the case of relatively coarse mediator space (case of h B = 320 µm), l opt is roughly of the same order of magnitude as h B . In practice, the junction parameter l can be chosen as: l = h B . This choice generally results in acceptable conditioning of A. In the remainder of this paper, l is set equal to l opt which is determined numerically for each coupling parameters setting. Table 5: Variation of l opt with the width L B and discretization characteristic length h B of the bridging subdomain.

Influence of the weight functions α

The influence of the weight function α is studied in this subsection. The other coupling parameters are chosen as follow: a fine mediator space M (fine discretization of ΩB , i.e. at the same order as the discretization of ΩD ) and L B = 2 mm.

Constant weight functions α = ᾱ = 0.5

The DEM and CNEM models are supposed equivalent in the bridging region, i.e. weighted in the same manner (Fig. 15a). Figure 28 presents the free-end displacement obtained in this case. The mean displacement obtained using the coupling method is U mean = 0.081 mm. This is in agreement with the reference results (Tab. 4). However, the temporal coupling displacement (Fig. 28) presents several deviations with regard to the reference curves. To analyze this coupling result and to determine the causes of these deviations, the velocities at the different check points (Fig. 24) have been plotted

for the first round trip of the wave propagation (Fig. 29). This figure shows that the high frequency waves (HFWs) initially captured at the "DemCheckPoint" do not appear in "OverlapDemCheckPoint"

or "OverlapCnemCheckPoint". In other words, the major part of the HFWs are reflected without entering the bridging region. This explains the deviation in the temporal displacement each time the global wave crosses the bridging region. Thus, constant weight functions are not a good choice for dynamic simulations. Indeed, the projection mechanism, which occurs in Ω B , cannot dampen the HFWs, and additional filtering is required in this case. In contrast, in the case of static studies, Guidault et al. [START_REF] Guidault | On the L2 and the H1 couplings for an overlapping domain decomposition method using Lagrange multipliers[END_REF][START_REF] Guidault | Bridging domain methods for coupled atomistic-continuum models with L2 or H1 couplings[END_REF] have shown that constant weight functions can be used with H1 coupling. is worse with regard to the reflection phenomenon. In this case, a part of the global incident wave is reflected at the interface between the models (not only the HFWs which are not supported by the coarse scale model). This result proves that the 1D static studies available in the literature cannot be used to perform dynamic coupling. Guidault et al. [START_REF] Guidault | On the L2 and the H1 couplings for an overlapping domain decomposition method using Lagrange multipliers[END_REF] have shown that, in statics and using a fine mediator space, the solutions do not depend on the weight functions.

Continuous weight functions

As explained earlier, the weight functions must not vanish at the boundary of the bridging zone ∂Ω B (for the weighted lumped mass matrices to be invertible), and a small value ε must be applied rather than zero at ∂Ω B . Before studying the influence of the continuous weight functions, the influence of ε is investigated. Figure 32 presents the free-end displacement using continuous weight functions for different ε. This parameter ε, when less than 0.05, has no practical influence on the results, but a very small ε can lead to instability problems. Indeed, as shown in Table 6, the smaller the ε, the greater the conditioning of the coupling matrix A becomes. In the remainder of this section, ε = 0.05 is used each time a continuous weight function is applied. Figure 33 shows the free-end displacement for the case of a linear continuous weight function (Fig. 15b). The associated velocities measured at the different check points are presented in Figure 34. No high frequency waves (HFWs) are reflected at the interface between the two models. Using continuous weight functions, the HFWs enter the bridging zone, and then, they are dampened by the projection of the fine space solution (from the DEM model) onto the coarse space solution (from the CNEM model). Figure 34b evidences that with a fine mediator space, a small bridging zone is sufficient to cancel out all of the HFWs. The use of a continuous weight function significantly improves the results (Fig. 33). However, a small deviation from the reference results still persists and becomes greater each time the wave travels back (CNEM-DEM direction). Because of the very fine discretization of the DEM subdomain, the weight of the particles in the bridging zone decreases smoothly when approaching the CNEM domain. By examining the first round-trip in Figure 33, it can be observed that practically no deviation from the reference is noted when the wave travels from the DEM subdomain to the CNEM subdomain. Therefore, the incident wave correctly crosses the interface between the pure DEM (Ω D \Ω B ) and the bridging zone (Ω B ). However, in the CNEM subdomain, a coarse discretization is used. The jump between the weight values of two adjacent nodes is relatively large. Thus, the same reflection mechanism, observed using constant weight functions with α = 0.5, occurs when the wave travels back (CNEM-DEM direction). To enhance the results, one can increase the width of the bridging zone to reduce the slope of the weight functions and so the weight jump between adjacent nodes. Another solution consists in using continuous differentiable weight functions (Fig. 15c) to reduce the weight jump in the vicinity of ∂Ω B (boundary of the bridging zone). Figures 42,[START_REF] González | Numerical integration in Natural Neighbour Galarkin methods[END_REF] and 44 present the results using the two solutions. The wave correctly crosses Ω B without any significant deviation. In conclusion, the continuous differentiable weight functions (Fig. 15c) are best adapted to weight the different energies in the bridging zone. These functions ensure a better regularity (smoothing) in the vicinity of the ∂Ω B , and then the HFWs cross the bridging zone where they are dampened, provided that L B and M are well-chosen. Otherwise, the same reflection mechanism can take place but for other reasons as will be explained in the next subsections. In the previous subsection, a fixed fine mediator space has been used. In such a case, the fine scale model velocities are practically locked on those of the coarse scale model (Fig. 35a). Indeed, the velocity curve at the "OverlapDemCheckPoint" is superimposed on that at the "OverlapCnemCheckPoint". The same locking phenomenon is noted in the literature when the third configuration (Fig. 22c) is used [START_REF] Guidault | On the L2 and the H1 couplings for an overlapping domain decomposition method using Lagrange multipliers[END_REF]. Now, to study the influence of M on the coupling results, a coarse mediator space (coarse discretization of ΩB , i.e. at the same order as the discretization of ΩC ) is used. As shown in Figure 35b, the equality of velocities in the bridging region is satisfied only in a weak sense and not in each node of the mediator space. This allows the fine model (DEM model) to correctly act in Ω B .

However, in this case a small bridging zone is insufficient to correctly transmit the global incident wave and cancel the HFWs. A solution to reduce the wave reflections at the interface between the two models is to use a large bridging zone (large L B ). As a consequence of this solution, a large discrete subdomain has to be used which is costly and cumbersome. However, contrary to previous works on coupling approaches, the discretization of the mediator space M can be made independently from the discrete and continuum models in this work. Therefore, an intermediate discretization of M, between the discretization of the continuum and discrete models, can be used. As it will be seen in the next subsections, for any choice of L B , there is an optimal discretization of the mediator space h B that cancels the HFWs in the bridging region. As seen in the previous subsection, the influence of the bridging width (L B ) on the coupling results depends on the approximated mediator space M. Figures 36 and 37 provide the coupling results for different L B using fine and coarse mediator spaces, respectively. In the case of fine M, the width of the bridging region has no significant influence on the results because of the locking phenomenon.

In such a case, the DEM particles are strongly constrained in the bridging region which in turn will represent a rigid wall for the HFWs, regardless of L B . However, in the case of coarse M, the discrete particles can somewhat freely act in the bridging region. Therefore, using a small L B (narrow bridging region), the interface between the two models will act like a free surface for the HFWs. In this case, a large bridging region is required to sufficiently constrain the DEM particles in the bridging region and then to reduce the HFWs reflections. Figure 37 shows that, in the case of coarse mediator space, the coupling results are better as L B increases. (Fig. 39). However, in this case, the optimal value of h B is h B = 3 h D . This value depends, inter alia, on the ratio between h D and h C . The larger this ratio, the larger the optimal value of h B becomes.

To analyze the influence of L B on the optimal value of h B , L B = 6 mm is then used when h C = 10 h D is held fixed. Figure 40 presents the associated results. The same conclusions and remarks, as for L B = 2 mm, are pointed. However, in the present case, the optimal value is h B = 9 h D and the deviation mechanisms are much less sensitive to h B . In conclusion, the optimal value of h B depends also on L B , but the choice of this parameter is less sensitive when large bridging region is used. In a general case there is not an obvious method to determine, in a single way, the various coupling parameters to avoid wave reflections. This subsection gives several recommendations and trends to choose correctly these parameters. To minimize the conditioning of the coupling matrix A, the H1 coupling with l = l opt is recommended. l opt can be chosen at the some order of magnitude as the discretization characteristic length of the bridging subdomain (l opt ≈ h B ). Although, in the case of very fine mediator space M, the optimal value of l is not at the same order of magnitude as h B , the choice of l ≈ h B generally results in acceptable conditioning of A. The weight functions must have at least C 0 continuity in the bridging region, but more regular functions (C 1 continuity or more) are recommended. Indeed, with constant weight functions, the high frequency waves (HFWs) are reflected without entering the bridging region and cannot be dampened by the projection mechanism which occurs in this region. The choice of L B and M is more challenging. Indeed, these parameters depend on each other, and on the discretization of the discrete and continuum models. In practice, the choice of these parameters can be made as follows: the discrete and continuum models are discretized according to the physical effects that have to be captured; then, the width of the bridging region can be fixed as narrow as possible, containing sufficient number of CNEM nodes and DEM particles; finally, M can be chosen such that h D ≤ h B ≤ h C . Using an arbitrary choice of h B banded between h D and h C , a dynamic simulation can be performed to control the wave propagation between the models; if an inward deviation is found when the wave crossing the interface, h B must be increased; if an outward deviation is found, h B must be decreased; Otherwise, the optimal value of h B is luckily used.

Applications

The previous parametric study using 3D beam submitted to a dynamic tensile loading has highlighted several recommendations to perform a correct coupling. In this subsection, the results of this study are used to apply the DEM-CNEM coupling approach on several elementary tests: 3D beams submitted to dynamic tension, dynamic bending, dynamic torsion and initial velocity loading. The length and diameter of these beams are respectively L = 100 mm and D = 20 mm (L/D = 5). The DEM method is applied only for the portion located 20 mm from the right end (the section located at x = L) and the remainder of the domain is modeled using the CNEM method (Fig. 41). The different coupling parameters are chosen as follows: L B = 10 mm, continuous differentiable weight functions, ε = 0.05, coarse mediator space, l = h B . Figure 44 presents the temporal free-end rotation about the x -axis using the torsional loading. The coupling results are in good agreement with the DEM results (Tab. 7). Finally, the coupling method has been tested using an initial velocity loading (Test 4 in Figure 41). Figure 45 presents the associated free-end displacement with respect to x -axis. The coupling result is comparable to the reference one, obtained from fully DEM computation. In the context of numerical simulation, four scales can be distinguished: nanoscale, microscale, mesoscale and macroscale. For each scale, efficient computational methods have been developed to handle the relevant phenomena. An overview of these methods is given in section 2. These methods can be divided into two classes: discrete methods (DMs) and continuum methods (CMs). The discrete methods are used to study phenomena ranging from nanoscale to mesoscale. These methods give generally very accurate results, but are very time-consuming and can only be applied for small systems. The continuum methods are used to study macroscopic phenomena. Although less accurate than the DMs, these methods are relatively inexpensive and can be applied for large systems. The features and advantages of the different methods of these classes are largely complementary. Therefore, several coupling approaches involving numerical methods of different nature have been proposed so as to strengthen their advantages and to avoid their drawbacks. As shown in section 2, these approaches can be divided into three categories: discrete-discrete coupling methods, continuum-continuum coupling methods and discrete-continuum coupling methods. The last category is of particular interest in multiscale modeling of complex dynamic problems. Indeed, most of these problems involve multiscale phenomena that require, on one scale, a very accurate and computationally expensive description (by discrete method), and on another scale, a coarser description (by continuum method) to avoid prohibitively large computation. The discrete-continuum coupling methods can be of different nature: hierarchical, concurrent or hybrid hierarchical-concurrent. The hierarchical coupling approaches are adapted for problems in which decoupled or weakly coupled multiscale phenomena are involved. The concurrent approaches are well adapted for problems involving highly dependent multiscale phenomena. The hybrid hierarchical-concurrent approaches, which combine the features of the hierarachical and concurrent approaches, are used for problems in which the multiscale phenomena are weakly coupled up to a certain response limit, and then highly dependent. The concurrent coupling approaches are best suited to study complex dynamic problems, which generally involve highly dependent phenomena.

The most used techniques to couple concurrently discrete and continuum methods are reviewed in section 3. Among them, the bridging domain (BD) technique is a good candidate to study properly complex problems in dynamics. Indeed, this technique can avoid high frequency wave reflections at the interface between the coupled models without any additional filtering or damping. This is particularly true when the different BD parameters are well set up. To corroborate the existent works aiming to simplify the setting of these parameters in practice, a 3D dynamic study of these parameters is proposed in this paper. Based on this study, several recommendations have been drawn to address this issue.

Fig. 3 :Fig. 4 :

 34 Fig. 3: Example of cohesive links used in MDMs

Fig. 5 :

 5 Fig. 5: Classification of continuum methods

Fig. 9 :

 9 Fig. 9: Different geometries of problem domains

Fig. 10 :

 10 Fig. 10: 1D models for investigation of the spurious wave reflections

( a )Fig. 12 :

 a12 Fig. 12: Wave traveling between two models having different wave propagation velocities

γ

  and f λ are respectively the weighted vector of internal forces, the weighted vector of external forces and the vector of coupling forces in the discrete model, M β is the weighted continuum lumped mass matrix, F int α , F ext γ and F λ are respectively the weighted vector of internal forces, the weighted vector of the external forces and the vector of coupling forces in the continuum model. C B and c B are the continuum and discrete coupling matrices. For M β and m β to be invertible, the weight functions β and β must be strictly positive in Ω B and at the border ∂Ω B . Therefore, a small ε must be used instead of zero at the nodes assigned to ∂Ω B . The definition of the weight function β is slightly amended from[START_REF] Guidault | Bridging domain methods for coupled atomistic-continuum models with L2 or H1 couplings[END_REF] as:

Fig. 18 :

 18 Fig. 18: Removal of unnecessary degrees of freedom in bridging scale modeling

  Fig. 19: Quasicontinuum model

Fig. 21 :

 21 Fig. 21: Techniques for coupling discrete and continuum meshless methods

i and c int i are respectively the internal force and torque applied on the particle i. f ext i and c ext i are respectively the

  external force and torque applied on the particle i.

  2, this provides more flexibility to the choice of the discretization characteristic lengths in the continuum and discrete subdomains. The different discretized subdomains of Ω D , Ω C and Ω B are denoted by ΩD , ΩC and ΩB , respectively.Since Ω D is naturally discretized, ΩD is the same that Ω D ( ΩD = Ω D ). According to the configuration of the discretized subdomains in the bridging zone, four cases can be distinguished (Fig.22). The first three configurations (Figs.22a, 22b and 22c) present some regularity in the bridging zone and are widely studied in the literature. The present work considers the fourth configuration which is the most general and includes the first three configurations: no coincidence conditions are imposed in the bridging region. In fact, using this configuration, it is very difficult to prove mathematically the existence and uniqueness of the solutions. The well-posedness of the global problem has been analyzed

( a )Fig. 22 :

 a22 Fig. 22: Different configurations of the discretized subdomains in the bridging zone

Fig. 23 :

 23 Fig. 23: Main parameters of the bridging domain coupling approach

Fig. 24 :Fig. 25 :

 2425 Fig. 24: Beam model for the parametric study

Fig. 26 :

 26 Fig. 26: Free-end displacement obtained by fully DEM computation and fully CNEM computation

  h B (µm) l opt (µm)

Fig. 28 : 5 Fig. 29 :

 28529 Fig.28: Free-end displacement obtained using DEM, CNEM and the coupling method; L B = 2 mm, fine mediator space and constant weight functions α = 0.5

Fig. 30 :

 30 Fig. 30: Free-end displacement obtained using the coupling method with different constant weight functions: α = 0.3, α = 0.5 and α = 0.8; L B = 2 mm and fine mediator space

Fig. 32 :Table 6 :

 326 Fig.32: Free-end displacement obtained using the coupling method for different values of ε; L B = 2 mm, fine mediator space and continuous weight functions

Fig. 33 :

 33 Fig.33: Free-end displacement obtained using DEM, CNEM and the coupling method; L B = 2 mm, fine mediator space and continuous weight functions

Fig. 34 :

 34 Fig. 34: (a) Linear velocities at the different check points for the case of continuous α; (b) comparison of linear velocities at the "DemCheckPoint" obtained using continuous and constant (α = 0.5) weight functions; L B = 2 mm and fine mediator space

Fig. 35 :

 35 Fig. 35: Velocities comparison in the bridging zone using fine and coarse mediator spaces; L B = 2 mm, continuous weight function and ε = 0.05

Fig. 36 :Fig. 37 :

 3637 Fig. 36: Free-end displacements obtained using the coupling method for L B = 2 mm, L B = 4 mm and L B = 6 mm; Fine mediator space , continuous weight functions, ε = 0.05

Fig. 38 :Fig. 39 :Fig. 40 :

 383940 Fig. 38: Influence of the discretization characteristic length of the bridging subdomain h B ; L B = 2 mm and h C = 5 h D

Fig. 41 :

 41 Fig. 41: Validation model

Fig. 42 :Fig. 43 :Fig. 44 :Fig. 45 :

 42434445 Fig. 42: Tensile loading: Free-end displacements U x mean obtained by DEM, CNEM and the coupling method; L B = 10 mm, coarse mediator space , continuous differentiable weighting and ε = 0.05

  The combined Lagrangian-Eulerian methods are very promising and have received much research interest. However, they are generally difficult to implement. Besides, unexpected termination of the computation process is often encountered, due to a highly distorted Lagrangian element or a very small Eulerian cell. Table1presents a comparison between this category and the two first ones.

		Lagrangian methods	Eulerian methods	Combined L-E
	Discretization	Mesh attached to material	Mesh fixed in the space	Updated mesh
		FEM [120, 122, 121]	FVM [108]	CEL [38, 73, 113]
	Examples	X-FEM [76]	FDM [98]	ALE [99]
		BEM [47]
	Track	Moving of material points	Mass, momentum and	Combined solution
			energy flux accross mesh
			cells boundarries
	Implementation	⋆⋆⋆	⋆⋆⋆	⋆
	Time history at	⋆⋆⋆	⋆	⋆⋆
	material points
	Boundaries and	⋆⋆⋆	⋆	⋆⋆
	interfaces
	Complex	⋆⋆⋆	⋆	⋆⋆
	geometries
	Large	⋆⋆	⋆⋆⋆	⋆⋆
	deformation
	Discontinuities	⋆⋆	⋆	⋆
	Computation	⋆⋆⋆	⋆⋆	⋆⋆
	robustness
	⋆⋆⋆: Good/Easy	⋆⋆: Intermediate	⋆: Bad/Difficult
		Table 1: Comparison of grid-based methods
	2.2.2 Meshless methods
	Despite the great success of the grid-based methods in both CSM and CFD, their disadvantages have
	pushed researchers to seek for new numerical methods. One important goal of the initial research is to
	modify the internal structure of the grid-based methods to become more adaptive and more robust.
	Additional effort was given to problems which cannot easily be treated by this class of methods, such
	as problems with moving boundaries (for Eulerian methods) and large deformation (for Lagrangian
	methods). Such an effort has given rise to the next generation of computational methods: the meshless
	methods. The main key of the meshless methods is to provide accurate and stable solutions for integral
	equations or partial differential equations (PDE) with a set of arbitrary distributed nodes (or particles)

Table 2

 2 presents a comparison between the interpolation and approximation methods.

	Subclasses	Approximation methods	Interpolation methods
	Discretization	Scattered nodes attached to	Scattered nodes attached to
		the material	the material
		SPH [67, 71]	NEM [22, 107]
		PUFEM [74]	α-NEM [27, 28]
		EFGM [13]	CNEM [115, 116]
	Examples	MLPG [8]
		RKPM [61]
		PIM [66]
		FPM [80, 81, 82]
	Track	Movement of material points	Movement of material points
	Implementation	⋆⋆⋆	⋆⋆
	Time history at	⋆⋆	⋆⋆⋆
	material points
	Boundaries and	⋆⋆	⋆⋆⋆
	interfaces
	Complex	⋆⋆	⋆⋆⋆
	geometries
	Large	⋆⋆⋆	⋆⋆⋆
	deformation
	Discontinuities	⋆⋆	⋆⋆
	Computation	⋆⋆	⋆⋆⋆
	robustness
	⋆⋆⋆: Good/Easy	⋆⋆: Intermediate	⋆: Bad/Difficult

Table 2 :

 2 Comparison of meshless methods

Table 3 :

 3 The microscopic properties of the cohesive beam bonds in the discrete subdomain

		U mean (mm)	f 0 (Hz)	f 1 (Hz)	f 2 (Hz)
	Theory	0.087	71 757	215 272	358 787
	DEM	0.083	72 408	217 246	362 072
	CNEM	0.088	71 359	214 023	356 491

Table 4 :

 4 Comparison of DEM, CNEM and analytical resultsfrequencies. It can be seen that the results are in good agreement, and they are also in agreement with the beam theory results. This ensures the equivalence of the two models.

Table 7 :

 7 Comparison of results; L B = 10 mm, coarse mediator space , continuous differentiable weight functions, ε = 0.05

	Beam theory	DEM	CNEM	Coupling

A large difference between the results is observed. In the first case (α = 0.3), the magnitude of the free-end displacement is greater than the one obtained using α = 0.5. However, it is smaller in the case of α = 0.8. To provide an explanation for these results, the temporal velocities at the check points are presented in Figure 31.

For α = 0.8, a portion of the global incident wave is positively reflected at the interface between the two models without entering the bridging zone. Furthermore, only a complementary part is transmitted in the CNEM model. Quantitatively, the transmission and reflection coefficients are t num = 0.44 and r num = 0.56, respectively. By analogy with the wave propagation between media with different acoustic impedances, we define the transmission and reflection coefficients in terms of α and ᾱ as follows: t th = 2 ᾱ α+ ᾱ and r th = α-ᾱ α+ ᾱ . It can be verified that t num and r num are of the same order of magnitude as t th = 0.4 and r th = 0.6, respectively. Fig. 31: Linear velocities at the check points for different constant weight functions: α = 0.3, α = 0.5 and α = 0.8; L B = 2 mm and fine mediator space For α = 0.3, the same reflection mechanism takes place but with a negative coefficient. Indeed, the velocity magnitude of the transmitted wave (measured at "CnemCheckPoint") is greater than the velocity magnitude of the incident wave (initially measured at "DemCheckPoint"). Also, the numerical transmission and reflection coefficients, t num = 1.44 and r num = -0.44, are of the same order of