We discuss the multiscale modeling of a granular material trapped between continuum elastic domains. The amorphous granular region, usually termed “gouge”, is under high confinement pressure, to represent the loading of faults at depth. We model the granularity of gouge using the Discrete Element Method (DEM), while the elastic regions surrounding it are represented with two continuum domains modeled with the Finite Element Method (FEM). We resort to a concurrent coupling of the discrete and continuum domains for a proper transmission of waves between the discrete and continuum domains. The confinement pressure results in the appearance of a new kind of ghost forces, which we address via two different overlapping coupling strategies. The first one is a generalization to granular materials of the Bridging Method , which was originally introduced to couple continuum domains to regular atomic lattices. This method imposes a strong formulation for the Lagrange constraints at the coupling interface. The second strategy considers a weak formulation. Different DEM samples sizes are tested in order to determine at which scale a convergence of the elastic properties is reached. This scale sets the minimal mesh element size in the DEM/FEM interface necessary to avoid undesirable effects due to an elastic properties mismatch. Then, the two DEM/FEM strategies are compared for a system initially at equilibrium. While the performance of both strategies is adequate, we show that the strong coupling is the most stable one as it generates the least spurious numerical noise. Finally, as a practical example for the strong coupling approach, we analyze the propagation of pressure and shear waves through the FEM/DEM interface and discuss dispersion as function of the incoming wave frequency. 1 Introduction The behavior of seismic faults depends on the response of a gouge region squeezed between solid rock masses. The micro-constituents of the gouge play an important role in the fault friction stability [22]. Numerical models representing granular materials subject to shearing provide useful insights to understand stability and emerging stick-slip behavior [1, 15, 16, 22, 29]. Such models generally utilize the Discrete Element Method (DEM) to represent the gouge, which is composed of an amorphous grain structure subjected to high confining pressure. While the gouge is made of grains, the surrounding regions can be modeled as an elastic continuum domain to represent efficiently and accurately the far-field boundary conditions. To obtain a proper transmission of waves between the discrete and continuum domains, the discrete and continuum domains need to be coupled together. Here, we choose to handle the interaction between the surrounding solids and the gouge with a FEM/DEM coupling [17, 20]. FEM/DEM couplings have the advantage of conserving the physical mechanisms of both domains, while drastically reducing the computational cost when compared to pure DEM simulations. During the last decades, two main classes of coupling methods have emerged: overlapping methods [28, 31] and edge-to-edge methods [11, 24, 26]. Edge-to-edge methods do not consider a reconciliation zone between the two domains. The coupling is implemented by imposing reciprocal boundary conditions. Contrariwise, overlapping methods consider a reconciliation zone where the coupled domains are mixed together. The Arlequin method [13] introduced by H. B. Dhia, considers an energy weighing in order to control the 1 influence of each domain, while a set of kinematic constraints enforces coherency thanks to Lagrange multipliers. Strong and weak Lagrange constraints lead to different formulations of kinetic constraints. The Bridging Method [31] is a classical overlapping method, based on the Arlequin method, which considers a strong constraint and, has been used to couple atomistic and continuum domains. Edge-to-edge and overlapping methods have been initially developed to couple perfect crystal structures (through molecular dynamics) with continuum domains (MD/FEM). DEM/FEM coupling differs from MD/FEM coupling on several aspects and brings unique challenges. First, the targeted particles obey different interaction laws: atoms interact through energy potentials (leading to both attractive and repulsive forces that are function of the distance between particle centers), while discrete particles interact according to macroscopic contact laws (short range interactions), which require the consideration of rotational degrees of freedom. Second, MD is more often than not used with a well-organized crystal structure, while DEM considers generally an amorphous structure, which leads to heterogeneous effective elastic properties. These differences make the extension of MD/FEM to DEM/FEM a far-from-trivial task. Several approaches have been put forward to couple continuum and discrete (DEM) domains, including overlapping [10, 12, 30, 32] and edge-to-edge methods [17, 20]. However, to the best of our knowledge, none of these methods have been tested for a three-dimensional amorphous arrangement of DEM particles (e.g. a gouge) subjected to a confinement pressure. The confinement pressure has a destabilizing effect, which challenges the enforcement of Lagrange multipliers and results in Ghost Forces of a new kind, which have not been previously reported in the literature. In this paper we demonstrate the emergence of these spurious effects, and provide a simple modification of the Bridging Method allowing to cancel all Ghost Forces both for a weak and a strong formulation of the constraints. As a by-product, this work also provides simple guidelines to the multiscale modeling of a gouge at depth, in particular in view of seismic applications. Section 2 presents two FEM/DEM coupling methods. The first one enforces a strong coupling and is an extension of the Bridging Method [31] developed for MD/FEM to DEM/FEM. The second one considers a weak formulation of the Lagrange constraint. In Section 3, the behavior of a granular material is studied to determine its effective elastic properties and help choose the finite-element mesh size at the DEM/FEM interface. Section 4 compares the coupling strategies in order to determine the most appropriate one for stability. Then, in Section 5, the passing through of pressure and shear waves is analyzed for different wave frequencies. 2 Coupling strategy for FEM-DEM Overlapping methods consider a reconciliation zone (also called the Bridging Zone ), where the continuum and discrete domains are mixed together. We consider two different overlapping method formulations. The first method imposes a strong coupling with a Lagrange constraint where each particle is constrained to the equivalent position coming from the FEM interpolation. The second is a weak coupling that enforces such a constraint in an average sense within the volume of each element. The first step of Arlequin couplings consists in defining an energy weighting, which permits to con- trol the relative influence of each model. Thus, the global energy in the system is defined using the Hamiltonian H=Z Ωα(X)EC(X) + (1 −α(X))ED(X)dX, (1) where ECandEDare the energy densities in the continuum and discrete regions. These are weighted in the bridging zone using a scaling function αwhich defines the influence of the continuum domain α(X)for a given spatial location X. 2 However, kinematic consistency between the coupled models cannot be enforced by the energy weighting. Therefore, the two models are tied together by means of Lagrange constraints. The energy weighting and the Lagrange constraint ( g=0) lead to a Lagrangian problem which can be summarized as: HL=H+λTg, (2) where His the original Hamiltonian, HLits expanded version including the coupling restrictions and λ the Lagrange multiplier matrix. This leads to new equations of motion:   M¨u=F+λ·∂g ∂u m¨d=f+λ·∂g ∂d(3) where u,FandMare displacements, forces and masses of the FEM domain, and d,fandmare the displacements, forces and masses of the particles. The energy weighting leads to the altered masses m andM, which are defined such that mi= (1−α(Xi))miandMI=α(XI)MIwithXiandXIthe initial positions of particles and continuum FEM nodes. The Verlet[5] scheme is employed to perform the time integration of the uncoupled domains, which produces test velocities ( ˙u⋆,˙d⋆) not respecting the coupling constraint. ˙u⋆,n+1/2=˙u⋆,n+∆t 2M−1Fn un+1=u⋆,n+ ∆t˙u⋆,n+1/2 Computation of Fn+1 ˙u⋆,n+1=˙u⋆,n+1/2+∆t 2M−1Fn+1 ˙d⋆,n+1/2=˙d⋆,n+∆t 2m−1fn dn+1=d⋆,n+ ∆t˙d⋆,n+1/2 Computation of fn+1 ˙d⋆,n+1=˙d⋆,n+1/2+∆t 2m−1fn+1(4) The correction of the test velocities needs to be done within the reconciliation zone thanks to the Lagrange multipliers, which gives at the n-th step the following form [31]: ˙un+1=˙u⋆,n+1−M−1∂g ∂uT·λ ˙dn+1=˙d⋆,n+1−m−1∂g ∂dT·λ(5) The velocities produced by Eq. 5 are chosen to satisfy the time derivative of the Lagrange constraint ˙ g=0. This choice imposes that the Lagrange multipliers λare obtained by solving: ˙ g=Aλ, (6) where Ais the constraint matrix, which depends on the constraint function gand can deliver distinct formulations of the coupling, such as the strong and weak coupling detailed in sections 2.1 and 2.2. 2.1 Strong coupling The strong coupling considers a Lagrange constraint, where each particle is constrained to the position interpolated from the FEM displacement field. Therefore, the constraint applied to the particles is expressed as: g=NTu−d=0, (7) where NTuis the interpolated displacement for all coupled particles, computed with the standard finite-element shape functions evaluated at each of the particles’ positions (Xi): 3 N= φ1(X1)···φ1(Xi) ... ···... φI(X1)···φI(Xi)  (8) where φIis the FEM shape function for node Iemployed to discretize the continuum domain. The expression of corrected velocities from Eq. 5 together with the time derivative of Eq. 7 allows to identify the particular form of Eq. 6 specific to the strong coupling: NT˙u⋆,n+1−˙d⋆,n+1 | {z } ˙g= NTM−1N−m−1 | {z } Aλ, (9) Additionally, to reduce the computational cost and provide a dissipative treatment of wave reflections [7], the constraint matrix Acan be lumped. Therefore, the Lagrange multipliers can be obtained by solving Eq. 9. Finally, the velocities can be corrected using the computed Lagrange multipliers, which leads to the following algorithm to be applied in the reconciliation zone: ˙ g=NT˙u⋆,n+1−˙d⋆,n+1 A=NTM−1N−m−1 Optional lumping of A λ=A−1˙ g ˙un+1=˙u⋆,n+1−M−1Nλ ˙dn+1=˙d⋆,n+1+m−1λ(10) 2.2 Weak coupling The idea behind the weak coupling is to relax the constraint so that it is respected in an integral sense. The intended effect is that the constraint is imposed in average over each element and not anymore point wise over the particles. For this, the standard collocation method taken from finite elements is employed, leading to the following weak constraint: gW I=Z ΩφI(X)G(X)dX=X iφI(Xi)gi=0 (11) where the integral becomes a discrete sum due to the discrete nature of the constraint, which can be formally written by defining the continuous constraint G(X) =P iδ(X−Xi)giwith Dirac δdistributions. Therefore, there are as many Lagrange constraints as the number of coupled FEM nodes, and it takes the following matricial form: gW=N NTu−d =0 (12) This leads to the weak formulation of Eq. 6: N NT˙u⋆,n+1−˙d⋆,n+1 | {z } ˙gW=NAz }| { NTM−1N+m−1 NT | {z } AWλW(13) In the case of weak coupling, the constraint matrix AWcannot be lumped as it would introduce too important approximations in the resolution of the Lagrange multipliers, a poor control of the influence of each coupled model and a spurious time evolution. Considering the previous developments, the algorithm employed to correct weakly for the velocities in the reconciliation zone becomes: 4 ˙ gW=N NT˙u⋆,n+1−˙d⋆,n+1 A=NTM−1N−m−1 AW=NANT λW= [AW]−1˙ g ˙un+1=˙u⋆,n+1−MNNTλW ˙dn+1=˙d⋆,n+1+mNλW(14) 3 Granular medium and constitutive behavior In order to test the presented coupling approaches, we consider a system where the continuum and discrete domains have a matching material behavior within the range of small (elastic) deformations. It is therefore necessary to determine the macroscopic elastic properties of the granular ensemble. The granular medium represents a dry gouge, with no adhesion between the particles, and subjected to a specific confinement pressure coming from the surrounding medium. Several DEM box sizes are generated to determine if the granular material can be considered as linear elastic, and to observe if the extracted material properties are converging when the sample size becomes large enough. The material properties of the FEM and the mesh size at the FEM/DEM interface will be chosen based on these single scale DEM results, with the ambition to avoid any mismatch in the material response between the coupled domains. 3.1 DEM sphere packing Nine different box sizes, with periodic boundary conditions in the three directions, are considered. In order to account for variability due to the random particle distribution, twenty samples have been generated for each box size. The DEM simulations are run using LAMMPS software [27], with DEM grains parameters corresponding to glass bead values as described in Table 1. Glass beads are widely used as a model material in simulations [4, 15] and in experiments aiming to represent the behavior of granular materials [14, 19]. Two grain sizes were considered, with an average diameter d0= 3mm. The dimensionless stiffness is defined as κ= (E∗/σn)2/3with E∗=E/(1−ν2)[2, 3] and is characterizing the equilibrium overlap between particles δ,i.e.δ/d0≃κ−1[3]. To be in the so-called rigid limit, δ≪d0, it is necessary to have δ/d0=κ−1≪1. In our case κ= 495 ≫1, for which the granular assembly can be considered fairly stiff. Note that a too high overburden pressure could lead to a non-elastic behavior. Agnolin and Roux [3] observed an elastic behavior until a pressure of 10 MPa. As we are considering an overburden pressure ofσn= 5MPa, it is safe to assume an elastic contact behavior. The time step is chosen according to the criterion ∆t <0.14p mmin/kmax[21], with mminthe minimum grain mass and kmaxthe maximum normal contact stiffness. We considered the Hertz contact law [18] with tangential forces between two grains (details on the contact law are given in the supplementary material). To generate the DEM samples, the initial boxes are first filled with particles such that they occupy a volume fraction of 60%. Then, the boxes are compressed in the three directions and relaxed, iteratively, until an overburden pressure of 5MPa is reached. During the process, the final equilibrium conditions are approached as the particles’ excess kinetic energy is dissipated via viscous damping forces. For a box size of 12d0, approximately three million time steps are necessary to reach equilibrium at the desired overburden pressure. 5 Young’s modulus, E 50 (GPa) Grain density, ρ 2500 (kg/m3) Poisson’s ratio, µ 0.3 Overburden pressure, σn 5 (MPa) Grain minimal diameter, dmin 2.5·10−3(m) Grain maximal diameter, dmax 3.5·10−3(m) Grain mean diameter, d0 3·10−3(m) Friction between grains, µg 0.5 Elastic constant for normal force, kn 3.663·1010(N/m2) Elastic constant for tangential force, kt 4.5249·1010(N/m2) Viscoelastic damping constant for normal force, γn 0.2 (1/(ms)) Viscoelastic damping constant for tangential force, γt 0.1 (1/(ms)) Time step, dt 10−8(s) Table 1: DEM grain parameters. See supplementary material for the definitions of kn, kt, γnandγt. 3.2 Homogeneous Elastic properties The boxes containing the particles are then subjected to positive and negative displacements equivalent to particular strain states, where each component ( εxx,εyy,εzz,γxy,γxz,γyz) is evaluated separately (see elastic examples of LAMMPS [27]). After each displacement, the boxes are relaxed, before measuring the overall stress. The 36 elastic tensor coefficients ( Cij) are computed with the corresponding applied strain direction (an average of positive and negative strains is performed to get a better accuracy). The Zener [33] anisotropy index 2C44/(C11−C12)can be extracted as a measure of anisotropy. Figure 1a shows that this anisotropy coefficient is very close to unity for a box size larger than 6d0. This implies that for a large enough box, our granular samples can be considered isotropic. Figure 1b shows the elastic tensor of a box of length 12d0for a specific realization. The structure of the elastic tensor confirms the isotropy of the granular medium. 02468101214161820222426 Lx/d00.000.250.500.751.251.501.752.00 1.00Zener ratio (a) Ci1Ci2Ci3Ci4Ci5Ci6C1j C2j C3j C4j C5j C6jλ+2µ 1.18 λ 0.4 λ 0.4 ∼0 ∼0 ∼0λ 0.4 λ+2µ 1.18 λ 0.4 ∼0 ∼0 ∼0λ 0.4 λ 0.4 λ+2µ 1.18 ∼0 ∼0 ∼0∼0 ∼0 ∼0 µ 0.39 ∼0 ∼0∼0 ∼0 ∼0 ∼0 µ 0.39 ∼0∼0 ∼0 ∼0 ∼0 ∼0 µ 0.39Elastic tensor 12 d0 0.00.20.40.60.81.0GPa (b) Figure 1: (a) represents the Zener anisotropy index in terms of the number of particles per box length Lx/d0. (b) represents the elastic tensor of the box of length 12d0at a specific realization. The vertical error bars represent the standard deviation due to the 20 box generations per box length. 6 The elastic properties (Young’s modulus, bulk modulus, shear modulus and Poisson’s ratio) are then determined. In Figure 2a, the Young’s modulus is shown as a function of the number of particles per box length, and reaches a plateau for box sizes larger than 12d0. The shear modulus and the bulk modulus are also following the same trend (see supplementary material Figures 9a and 9b), while the Poisson’s ratio (Figure 2b) and density (supplementary Figure 9c) appear practically constant with respect to the box size. The tests conducted suggest a convergence of the elastic properties for box sizes larger than 12d0. Therefore, any linear-elastic continuum domain discretized using a characteristic mesh size equal or larger than 12d0should display the same constitutive behavior (see Table 2 for numerical values) than the homogenized DEM samples. Young’s modulus, E 0.977 (GPa) Density, ρ 1543 (kg/m3) Poisson’s ratio, µ 0.255 Overburden pressure, σn 5 (MPa) Table 2: Effective constitutive parameters for a granular domain. Such values will be employed as material properties for the FE simulations. 02468101214161820222426 Lx/d00.700.750.800.850.900.951.001.051.101.151.20 0.977Young’s modulus E(GPa) (a) 02468101214161820222426 Lx/d00.00.10.20.30.40.5 0.255Poisson’s ratio ν (b) Figure 2: Convergence of elastic properties as the sample size becomes larger. (a) shows the Young’s modulus in terms of the number of particles per box length Lx/d0. (b) shows the Poisson’s ratio in terms of the number of particles per box length Lx/d0. The vertical error bars represent the standard deviation due to the 20 box generations per box length. 4 Simulations of coupled FEM/DEM 4.1 Geometry The geometry considered in all the following coupled simulations is that of an elongated 3D bar of constant section, made of a DEM region surrounded by two finite-element meshes, as presented in Figure 3. According to the results from the previous section, finite elements with a characteristic size h≈12d0have elastic properties matching with DEM granular packings. The DEM dimensions are chosen to ensure at least two elements, ≈24d0, in the section ( xandzdirections), a gouge thickness ≈20d0(the unconstrained 7 DEM), and two bridging regions covering the length of one element along the ydirection. The FEM dimensions are chosen in such a way that the dimensions in xandzmatch the DEM region, while the length of each FEM meshes are ≈400d0. Such a long bar will allow studying the propagation of large wave lengths. Periodic boundary conditions are imposed to the system in the lateral directions ( xandz). As previously stated, the DEM particles are packed with a hydrostatic confinement pressure, which needs to be balanced with the FEM subdomains. This is achieved with eigen-strains pre-stressing the FEM domains. The concurrent coupling will be in charge of exchanging forces between interface nodes and particles. As will be discussed below, the efficiency of the proposed coupling approaches differ notably. Figure 3: Schematic of the FEM/DEM coupling. 4.2 Strong/Weak coupling and ghost-force correction The strong and weak coupling strategies have been implemented in the open-source LibMultiScale software [6–9]. Note that a parallel version of the strong coupling is available in LibMultiScale . To determine whether the strong or weak coupling is more appropriate, this coupled geometry is first considered with no added perturbation, and run over long simulation times. Three different coupling set-ups, composed of three different DEM regions, have been tested to consider the effect of the amorphous structure. The displacement within a chosen cross-section in the FEM (indicated in red dashed line in Figure 3) is monitored in order to establish the performance of each coupling strategy. For the weak coupling, after a simulation time of 500,000timesteps, a displacement noise is observed to be in the order of 10−5m∼10−2·d0(Figure 4b), while it is in the order of 10−6m∼10−3·d0(Figure 4d) for the strong coupling. Such spurious displacements will keep growing as a consequence of a lack of equilibrium in the coupling region, due to the imposed pre-stress. The applied eigen-strains, which induce forces on the free surface of the mesh where nodes should be associated with a vanishing energy weight, are 8 at the origin of this issue. The approximations in the integration of the Lagrange constraint led to the equations presented in section 2, which precludes a null weight on the continuum interface. This becomes obvious when considering equations 10 and 14, where a division by the weighted mass ( ¯M) is performed, forbidding a perfectly null weight. In Anciaux [7], it was shown that an adequate, non-zero, choice of the weight for these limiting nodes could also diminish wave reflections substantially, in the case of crystalline structures evolving around their ground state configuration. With the presented amorphous and pre-stressed configuration, the same coupling approach does not compensate well the interface forces, giving rise to ghost forces . To balance these, a strategy similar to the one used by Shenoy et al. [25] has been employed: initial nodal forces were saved at the simulation onset, and then subtracted throughout the entire simulation. With this force correction, the weak coupling shows a diminished noise displacement in the order of 10−7m∼10−4·d0 (Figure 4a), and of the order 10−11m∼10−8·d0(Figure 4c) for the strong coupling. Therefore, it appears that the strong coupling with ghost force correction yields the most stable system, which will be used in the following section to study propagating waves through the coupling interface. 9 0 1052.1053.1054.1055.105 Timestep−1.0−0.50.00.51.0Displacement (10−7m)Weak coupling and substraction of initial force dx dy dz(a) 0 1052.1053.1054.1055.105 Timestep−0.50.00.5Displacement (10−5m)Weak coupling and no substraction of initial force dx dy dz (b) 0 1052.1053.1054.1055.105 Timestep−202Displacement (10−11m)Strong coupling and substraction of initial force dx dy dz (c) 0 1052.1053.1054.1055.105 Timestep−4−20Displacement (10−6m)Strong coupling and no substraction of initial force dx dy dz (d) Figure 4: Testing different coupling strategies. (a, b) show the average displacement at a specific cross- section considering the weak coupling. (c, d) show the average displacement at a specific cross-section considering the strong coupling. (a, c) consider the subtraction of the initial force to correct ghost forces, while (b, d) do not. The shaded regions represent the minimum and maximum displacement values obtained through the three different simulations. 5 Waves passing through the coupling interface 5.1 Wave train characterization Ricker [23] demonstrated that a seismogram can be decomposed into a succession of wavelets, which are nowadays called Ricker wavelets . To excite our coupled system and form a compressional wave propagating along the ydirection, we impose through time a Gaussian displacement corresponding to a ricker wavelet: 10 uy(t) =Iexp −1 22πc(t−t0) λ2! , (15) with cthe velocity of the wave, λits wavelength, Iits intensity and t0=λ/(2c). The Fourier transform is used to analyze the spectrum of the wave: F[uy(y)](k) =1√ 2πZ+∞ −∞uy(y)eikydy=1√ 2πZ+∞ −∞Iexp −1 22π(y−ct0) λ2! eikydy(16) The Fourier transform has a Gaussian shape localized around the fundamental wavelength λ. If the fundamental wavelength λis chosen to be much larger than the mean diameter of the particles d0, all the harmonics of significant intensity that are contained in the spectrum will be much larger than the diameter of a particle d0. Additionally, if a wavelength is larger than the mean diameter of the particles d0, it will be experienced by the DEM grains as a smooth transition from one grain to its neighbors, and it will not trigger rearrangements of the fabric. Therefore, the DEM portion will undergo small deformations and will behave as a linear-elastic medium. Also, a large wavelength reduces numerical dispersion in the FEM. Let us now derive a condition for a compressional wave to be considered low-amplitude. In what follows, let us note ξ=2π λ(y−c(t−t0)). To quantify the maximum value of the strain, we apply the chain rule: ε= ∂uy(y) ∂y = ∂uy(y) ∂ξ∂ξ ∂y ∼|uy,ξ(y)| λ∼I λ. (17) The condition for the continuum to remain linear and elastic is ε≪1(small strains), which consequently leads to I/λ≪1. While this result was derived for compressional waves, a similar approach can be pursued for shear waves with identical conclusions. 5.2 Wave propagation A wave propagation through the strong coupling system presented earlier (see Figure 3) will be monitored in order to characterize the artifacts produced by the coupling. A Ricker’s wavelet is generated at the extremity of the left FEM portion, by applying the displacement given by Eq. 15 as a boundary condition. The propagation of the wave is studied for two different mesh sizes h:12d0and3d0, for both compressional and shear Ricker waves. We consider two large wavelengths ( 200d0,400d0) and two small strains ( 8.3·10−4, 8.3·10−5) for compressional and shearing waves. Table 3 summarizes the wavelengths, intensities, and strains of the generated waves. For all the conducted simulations, the displacement of the wave will be monitored within the cross-section at y= 0.138m. This cross-section is highlighted with a dashed red line in the Figure 3 in the case of element size 12d0. Wavelength λ(m) 200d0 200d0 400d0 Intensity I(m) 5·10−45·10−51·10−4 Strain ε 8.3·10−48.3·10−58.3·10−5 Table 3: Waves characteristics. 11 −1 0 1 position in the bar (m)−0.50.00.51.0|uy/ umax|Wavelength λ=200·d0andh=12d0 ε= 8.3e-04 ε= 8.3e-05(a) −1 0 1 position in the bar (m)−0.50.00.51.0|uy/ umax|Wavelength λ=200·d0andh=3d0 ε= 8.3e-04 ε= 8.3e-05 (b) −1 0 1 position in the bar (m)−0.50.00.51.0|ux/ umax|Wavelength λ=200·d0andh=12d0 ε= 8.3e-04 ε= 8.3e-05 (c) −1 0 1 position in the bar (m)−0.50.00.51.0|ux/ umax|Wavelength λ=200·d0andh=3d0 ε= 8.3e-04 ε= 8.3e-05 (d) Figure 5: Waves passing through the coupling interface. (a, c) consider a characteristic mesh size of 12d0, while (b, d) consider a characteristic mesh size of 3d0. The top row figures (a, b) have been generated considering push waves, while the bottom row figures (c, d) have been generated considering shear waves. The red curves represent the waves before reaching the DEM, and the blue curves correspond to the waves after passing through the DEM. The dash vertical lines are visual guides for the location of the FEM-DEM interfaces. In the case of compressional and shear waves for a large mesh size of 12d0, we observe that the intensity of the wave after passing the DEM is close to 1 for every strain (see Figure 5a, Figure 5c). This means that we have an almost-complete transmission of the wave, regardless of its strain level. Contrariwise, for a mesh size of 3d0, we observe that the intensity of the transmitted wave is diminished by a small factor ( ∼ 0.88), independently of the strain (see Figure 5b, Figure 5d). The less-accurate transmission of the wave in the case of the smallest mesh size 3d0is due to the mismatch of material properties between FEM and DEM. In fact, it was shown in Section 3, that a mesh characteristic size greater than 12d0is necessary 12 to homogenize the elastic properties of the DEM. If this is not the case, the constitutive behavior varies among elements. It can be observed that a reflected wave has been generated in all cases (see Figures 5a, 5b, 5c, 5d). In order to determine the origin of the reflected wave we pursue an analogy with MD/FEM coupling. It has already been observed in MD/FEM coupling [7] that a reflected wave emerges at the coupling interface. The dispersion relation that relates the wave frequency to its wavelength can differ between the MD and the FEM depending on the FE mesh size (Figure 2.8 in [7]). If the mesh size is equal to the distance between atoms, the dispersion relation will be the same. However, when the mesh size increases, MD and FEM domains admit different wave velocities for the same harmonic. This difference of propagation velocities (dispersion) creates the impedance contrast that is at the origin of the reflected wave. The same mechanism is at the origin of the reflected wave that is observed in our DEM/FEM simulations, the main difference being that the dispersion relation of the granular medium is much more convoluted than the simple 1D atomic chain that is used in the illustrative example of [7]. The intensity of the reflected wave can be reduced by increasing the wavelength, as the dispersion relation in both domains will be similar in the long-wavelength limit if the material properties of the homogenized granular medium are assigned to the continuum. By increasing the wavelength compared to the mesh size, we reduce the difference of wave speeds between the two domains and therefore reduce the intensity of the reflected wave. In order to validate our discussion as to the origin of the reflected wave in the DEM/FEM coupled system, we study the propagation of two compressional waves, each with the small strain 8.3·10−5, but with different wavelengths ( 200d0,400d0). The considered discretization used characteristic element sizes 12d0. Figure 6 shows the displacement of these two compressional waves in the system. In both cases, a reflected wave is created during the propagation through the coupling and DEM regions. Though, a higher reflected intensity is observed for λ= 200 d0, than for a wavelength λ= 400 d0. This indicates that the reflected wave is due to dispersion happening at the FEM/DEM interface, which is strongly dependent on the wavelength. −3−2−1 0 1 2 3 position in the bar (m)0.00.20.40.60.81.0|uy/ uy,max|Strainε= 8.3·10−5 λ= 400d0 λ= 200d0 Figure 6: Transmission through the coupling interface of two compressional waves with λ= 200 d0(orange and dashed-dark blue for before and after the interface) and λ= 400 d0(red and dashed-light blue for before and after the interface). The dash vertical lines are visual guides for the location of the FEM-DEM interfaces. To strengthen our analysis, we plot the ctdiagrams for λ= 200 d0. Figure 7 shows the results for a large mesh size 12d0and Figure 8 for a small mesh size 3d0. These ctdiagrams show a better transmission of 13 the compressional and shear waves for a large mesh size. One can nonetheless observe reflected waves in all cases. Note that a residual displacement at the unconstrained DEM region is observed for a large mesh size 12d0, in the case of a compressional wave of strain 8.0·10−4. For the smallest mesh size 3d0, a residual displacement is observed independently of the strain, with higher values for compressional waves than for shear waves. Overall, the propagation of compressional and shear waves confirms that a mesh size of at least 12d0is necessary to avoid a mismatch of material properties between FEM and DEM. It also confirms the ability of our numerical framework to satisfactorily transmit waves of different kinds yet only for small strains, which is a standard assumption of concurrent coupling methods. -1.2 0 1.2 Position in the bar (m)1052.1053.1054.105Timestepε= 8.3·10−5 012345Displacement|uy|(10−5m) (a) -1.2 0 1.2 Position in the bar (m)1052.1053.1054.105Timestepε= 8.3·10−4 012345Displacement|uy|(10−4m) (b) -1.2 0 1.2 Position in the bar (m)1052.1053.1054.1055.105Timestepε= 8.3·10−5 012345Displacement|ux|(10−5m) (c) -1.2 0 1.2 Position in the bar (m)1052.1053.1054.1055.105Timestepε= 8.3·10−4 012345Displacement|ux|(10−4m) (d) Figure 7: Wave transmission through the coupling interface for a mesh size of 12d0andλ= 200 d0. (a, b) represent the ctdiagrams of compressional waves at different strains. (c, d) represent the ctdiagrams of shear waves at different strains. 14 -1.2 0 1.2 Position in the bar (m)1052.1053.1054.105Timestepε= 8.3·10−5 012345Displacement|uy|(10−5m)(a) -1.2 0 1.2 Position in the bar (m)1052.1053.1054.105Timestepε= 8.3·10−4 012345Displacement|uy|(10−4m) (b) -1.2 0 1.2 Position in the bar (m)1052.1053.1054.1055.105Timestepε= 8.3·10−5 012345Displacement|ux|(10−5m) (c) -1.2 0 1.2 Position in the bar (m)1052.1053.1054.1055.105Timestepε= 8.3·10−4 012345Displacement|ux|(10−4m) (d) Figure 8: Wave transmission through the coupling interface for a mesh size of 3d0andλ= 200 d0. (a, b) represent the ctdiagrams of compressive waves at different strains. (c, d) represent the ctdiagrams of shear waves at different strains. 6 Conclusions We have investigated the transmission of elastic compression and shear waves through a DEM/FEM interface. Two different coupling strategies have been proposed with regard to the formulation of the kinematic constraints, termed “strong coupling” and “weak coupling". The effect of subtracting the initial force at the FEM interface has also been studied. For verification purposes, we devised a simulation set-up in which both DEM and FEM are meant to represent the same material when deforming in small strain conditions. We first determined that our granular samples can be considered as homogeneous, linear-elastic and isotropic for sample sizes larger than 12times the representative grain size. This size constrains the minimum finite-element mesh size to be used in order to match elastic properties between the continuum and the discrete domains. Next, we developed a test to determine the most stable approach (in the sense of introducing the least noise in the system). We observed that the “strong” version of the coupling (e.g., constraints defined 15 exactly at each grain that overlaps with FEs) outperforms the “weak” one (which weights the constraints differently depending on the relative position of the grains with respect to the nodes of the element that engulfs it). Additionally, correcting the spurious forces induced by pre-stressed FEs was proved remarkably beneficial when it comes to stabilize the DEM/FEM overlapping region. We concluded that the strong coupling considering the subtraction of initial force is the proper choice for the problem at hand. Finally, longitudinal and shear waves were sent through coupled DEM/FEM model. The results reveal the presence of a reflected wave, due to the dispersion happening at the FEM/DEM interface, whose amplitude decreases as the wavelength increases. Additionally, the propagation of waves confirms that a mesh size of at least 120is necessary to avoid a mismatch of material properties between FEM and DEM. Overall, the results confirm the aptitude of our numerical framework in the LibMultiScale software to satisfactorily transmit small-strain waves of different kinds (compressional/shear) between the two domains. Acknowledgements J.-F.M., G.A., J.G.-S., and M.V.-L. acknowledge financial support from the Swiss National Science Foundation (Grant 200021_197152, “Wear across scales”).