\makeatletter
\@ifundefined{HCode}
{\documentclass[screen,CRMECA,Unicode,biblatex,published]{cedram}
\addbibresource{crmeca20240483.bib}
\newenvironment{noXML}{}{}%
\def\tsup#1{$^{{#1}}$}
\def\tsub#1{$_{{#1}}$}
\def\sfrac#1#2{{#1}/{#2}}
\def\sofrac#1#2{({#1}/{#2})}
\def\thead{\noalign{\relax}\hline}
\def\endthead{\noalign{\relax}\hline}
\def\tbody{\noalign{\relax}\hline}
\def\botline{\\\hline}
\usepackage{etoolbox}
\usepackage{upgreek}
\newcounter{runlevel}
\let\MakeYrStrItalic\relax
\def\jobid{crmeca20240483}
\def\morerows#1#2{#2}
%\graphicspath{{/tmp/\jobid_figs/web/}} 
\graphicspath{{./figures/}} 
\def\refinput#1{}
\def\back#1{}
\csdef{Seqnsplit}{\\}
\DOI{10.5802/crmeca.282}
\datereceived{2024-06-12}
\dateaccepted{2024-12-19}
\ItHasTeXPublished
}
{\documentclass[crmeca]{article}
\usepackage{upgreek}
\usepackage[T1]{fontenc} 
\def\CDRdoi{10.5802/crmeca.282}
\def\hyperlink#1#2{#2}
\def\hypertarget#1#2{#2}
\def\sfrac#1#2{{#1}/{#2}}
\def\sofrac#1#2{({#1}/{#2})}
}
\makeatother

\usepackage{upgreek}

\dateposted{2025-01-27}
\begin{document}

\begin{noXML}

\CDRsetmeta{articletype}{research-article}

\title{On the behavior of a granular soil deposit subjected to
horizontal vibration. A discrete element modeling}

\alttitle{Comportement d'un d\'{e}p\^{o}t de sol granulaire soumis
\`{a} une vibration horizontale : Mod\'{e}lisation par \'{e}l\'{e}ments
discrets}

\author{\firstname{Said} \lastname{Derbane}\CDRorcid{0009-0003-7773-5806}\IsCorresp}
\address{Department of Civil Engineering, LMGHU Laboratory, University of 20 August 1955 Skikda, BP 26, 21000, Skikda, Algeria}
\email[S. Derbane]{derbane\_said@yahoo.fr}

\author{\firstname{Mouloud} \lastname{Mansouri}\CDRorcid{0000-0002-4179-7488}}
\address{Department of Civil Engineering, University of Ferhat Abbas S\'{e}tif 1, Setif 19000, Algeria}

\author{\firstname{Salah} \lastname{Messast}\CDRorcid{0000-0002-0864-2257}}
\addressSameAs{1}{Department of Civil Engineering, LMGHU Laboratory, University of 20 August 1955 Skikda, BP 26, 21000, Skikda, Algeria}

\author{\firstname{Aboulghit} \lastname{~El~Malki~Alaoui}\CDRorcid{0000-0003-1357-8509}}
\address{ENSTA Bretagne, IRDL, UMR 6027, CNRS, 29806 Brest, France}

\begin{abstract}
This paper presents an analysis of the behavior of a non-cohesive
granular material deposit excited at its base by a horizontal harmonic
vibration. The analysis is carried out numerically by means of a 2D
discrete element model. The performed simulations highlighted some
aspects of vibration behavior in non-cohesive deposits, such as the
shape of the vertical profile of the displacement, notably in the case
of large displacements. The analysis particularly focused on the
amplification of the movement at the free surface of the deposit, as
well as its dependence on some parameters such as the excitation
frequency and the excitation amplitude of the deposit confinement. The
obtained results showed that the behavior of the deposit following the
change in the excitation frequency is similar to
the case of an elastic deposit excited by a
harmonic displacement at the base, i.e. the Dynamic Amplification
Factor (DAF) initially increases with the frequency increase,
it reaches a peak of resonance then it decreases. The resonance
frequency estimated from this analysis is close to the fundamental
frequency for low excitation amplitudes, but becomes smaller as the
excitation amplitude increases. On the other hand, for a fixed
frequency, increasing the amplitude of the excitation induces greater
amplification. It has been shown that this increase results from the
degradation of the shear modulus due to the increase in the level of
involved shear strain. Therefore, unlike elastic deposits, for
non-cohesive granular deposits, increasing strain leads to a
degradation of the shear modulus, resulting in a downward shift of the
resonance frequency and can induce a significant increase in
amplification. The confinement of the deposit is achieved by increasing
the gravitational acceleration; it has been shown that increased
confinement makes the deposit stiffer, and therefore reduces the
amplification of the introduced movement. 
\end{abstract}


\begin{altabstract}
Cet article pr\'{e}sente une analyse du comportement d'un
d\'{e}p\^{o}t de mat\'{e}riau granulaire non coh\'{e}sif soumis \`{a}
une vibration harmonique horizontale appliqu\'{e}e \`{a} sa base.
L'analyse est r\'{e}alis\'{e}e num\'{e}riquement \`{a} l'aide
d'un mod\`{e}le bidimensionnel bas\'{e} sur la m\'{e}thode des
\'{e}l\'{e}ments discrets. Les simulations effectu\'{e}es ont mis en
\'{e}vidence certains aspects du comportement vibratoire des
d\'{e}p\^{o}ts non coh\'{e}sifs, notamment la forme du profil vertical
du d\'{e}placement, en particulier dans le cas de grands
d\'{e}placements. L'analyse s'est particuli\`{e}rement
concentr\'{e}e sur l'amplification du mouvement \`{a} la surface
libre du d\'{e}p\^{o}t, ainsi que sur sa d\'{e}pendance \`{a} certains
param\`{e}tres tels que la fr\'{e}quence et l'amplitude de
l'excitation ou encore le confinement du d\'{e}p\^{o}t.

Les r\'{e}sultats obtenus montrent que le comportement du
d\'{e}p\^{o}t face \`{a} la variation de la fr\'{e}quence
d'excitation est similaire \`{a} celui d'un d\'{e}p\^{o}t
\'{e}lastique soumis \`{a} un d\'{e}placement harmonique \`{a} sa base.
Plus pr\'{e}cis\'{e}ment, le facteur d'amplification dynamique (DAF)
augmente initialement avec la fr\'{e}quence, atteint un pic de
r\'{e}sonance, puis diminue. La fr\'{e}quence de r\'{e}sonance
estim\'{e}e dans cette analyse est proche de la fr\'{e}quence
fondamentale pour de faibles amplitudes d'excitation, mais devient
plus basse \`{a} mesure que l'amplitude augmente. Par ailleurs, pour
une fr\'{e}quence donn\'{e}e, l'augmentation de l'amplitude de
l'excitation engendre une amplification plus importante. Il a
\'{e}t\'{e} d\'{e}montr\'{e} que cette augmentation r\'{e}sulte de la
d\'{e}gradation du module de cisaillement due \`{a} une augmentation du
niveau de d\'{e}formation en cisaillement. Contrairement aux
d\'{e}p\^{o}ts \'{e}lastiques, pour les d\'{e}p\^{o}ts granulaires non
coh\'{e}sifs, l'augmentation de la d\'{e}formation conduit \`{a} une
d\'{e}gradation du module de cisaillement, entra\^{i}nant un
d\'{e}calage vers le bas de la fr\'{e}quence de r\'{e}sonance et
pouvant induire une amplification significative.

Le confinement du d\'{e}p\^{o}t est simul\'{e} par une augmentation de
l'acc\'{e}l\'{e}ration gravitationnelle. Il a \'{e}t\'{e} montr\'{e}
qu'un confinement accru rend le d\'{e}p\^{o}t plus rigide, ce qui
r\'{e}duit l'amplification du mouvement introduit.
\end{altabstract}

\shortrunauthors

\keywords{\kwd{Sand deposit} 
\kwd{Behavior of a granular soil}
\kwd{Discrete element method} 
\kwd{Shear wave}
\kwd{Propagation}}

\altkeywords{\kwd{D\'{e}p\^{o}t de sable}
\kwd{Comportement d'un sol granulaire}
\kwd{M\'{e}thode des \'{e}l\'{e}ments discrets}
\kwd{Onde de cisaillement}
\kwd{Propagation}}

\maketitle

\end{noXML}

\section{Introduction}

Seismic wave propagation from the bedrock to the free ground surface
depends strongly on the nature of the traversed soil layers. In
particular, it is known that stiffness and damping of these layers have
significant effects on the transmission of the motion. For granular
soils these two properties i.e.\ stiffness and damping depend on the
applied movement characteristics such as amplitude and frequency, as
well as on certain soil state parameters, such as confining and density
states. The classical methods used in the analysis of movement
transmission and its amplification from the bedrock to the surface are
usually based on continuum mechanics. These methods work well for
cohesive soils, however for granular soils, given the discontinuous
nature of the material, discrete element modeling could provide more
understanding of the involved mechanisms.

In the context of this modeling approach, transient wave propagation in
granular materials was studied using a discrete element model by Martin
et~al.~\cite{1}, three particular cases were investigated including dry
cohesionless material, elastic cemented particulate media and fluid
saturated granular material, the results of this study indicated that
the wave velocity depends upon the stiffness of the interparticle
contacts and the distribution of branch vectors along the propagation
direction. Furthermore, it is shown that the attenuation of the wave
amplitude depends on the number of branch vectors in the propagation
direction. The influence of force chains on the wave scattering and
attenuation in granular soils was also studied by Peters and Muthuswamy
\cite{2}. As part of a study on the instabilities of force chains,
Campbell \cite{3} treated the destabilizing effect of wave propagation
using discrete element modeling. Zamani and El Shamy \cite{4} have
carried out an investigation of the vertical propagation of a shear
wave in a granular deposit through discrete element modeling, they
studied the change of the shear modulus and the damping ratio for low
frequencies. The properties extracted from the DEM simulations were
used to calculate the deposit response assumed as a continuum with
linear viscoelastic behavior using the SHAKE software. The obtained
results were then compared to those of the DEM modeling. Ning et al.
\cite{5}, after developing an appropriate method for measuring wave
velocity in a discrete element model, have studied the effects of
particle size and elastic properties on the wave propagation velocity.
O'Donovan et~al.~\cite{6} have conducted an experimental study on a soil
model in a cubical cell, they then compared the results with both
discrete element (DEM) simulations and continuum modeling. In their
study, they used point source transmitters and receivers in order to
assess the shear and compression wave velocities in the samples, from
which some of the elastic moduli can be deduced. As a result, a
satisfactory agreement between experimental observations and DEM
simulations is confirmed. Recently, Arran and co-authors \cite{7}
modeled the flow of spherical particles on a rough inclined base. By
comparison with their experimental results on physical models dealing
with seismic waves emitted by granular flows \cite{8,9} they
highlighted the effect of velocities of basal particle, which are
difficult to measure experimentally, on the flow properties. In order
to better understand how the dynamic load affects the evolution of the
microstructure and the whole granular assembly behavior, Longlong~and
Shunhua \cite{10} proposed a criterion to recognize the major
propagation path of dynamic load in 2D granular materials, called the
``dynamic force chain''. This analysis showed that the spatial
distribution of dynamic force chains in the indentation of granular
materials provides a direct measure of dynamic load diffusion. It is
demonstrated that the statistical evolution of dynamic force chains has
a strong correlation with the indentation behaviors. Besides, discrete
element modeling of wave propagation in saturated granular materials
has also received interest, particularly in the context of the study of
sand liquefaction phenomenon under seismic vibrations \cite{11,12,13}.

Despite all these studies, movement propagation in granular soils still
has many unclear behavioral aspects. Indeed, for these soils,
deformation is particularly related to intergranular slips, which
themselves depend on several parameters such as the strain amplitude,
the vibration frequency, the confining pressure, etc. Hence, these
complex aspects deserve to be subjected to more analysis. This work
consists of an analysis of the behavior of a non-cohesive granular
deposit subjected to dynamic excitation at its base using a 2D discrete
element model. The exciting motion is harmonic with a controlled
amplitude and frequency. We particularly aim to highlight and analyze
some involved features, such as the movement amplification at the free
surface and its dependence on the frequency and amplitude of excitation
as well as the deposit confinement. We also show the macroscopic
behavior of the deposit during large displacements or near the
resonance. 


\section{Discrete element modeling} \label{sec2}

The discrete element method used in this work was first initiated by
Cundall in 1992 \cite{14} and is based on molecular dynamics. It
models granular materials at the micromechanical scale using
independent elements. The grains of the material interact with each
other through contact forces that are calculated through simple models
based on a slight overlap of the grains \cite{15,16}. It is assumed
that the global deformation of the medium is mainly due to the relative
movements of the grains considered as rigid bodies \cite{17,18}.
Because of their simplicity and the savings in calculation time they
involve, circular element shapes are still the most widely used. 

The movement of each grain $i$ is governed by the following
Newton's second law, which allows to obtain the accelerations of the
grains. 
{\begin{eqnarray}\label{eqn1}
\begin{cases} 
\displaystyle
m_{i}\vec{\ddot{x}}_{i}=\sum _{j}{\vec {F}}_{ij}^{\mathrm{contact}}+m_{i}\vec {g}
\vspace*{6pt}\\
\displaystyle
I_{i}\vec {\ddot{\varphi }}_{i}=\sum _{j}{\vec {M}}_{ij}^{\mathrm{contact}} 
\end{cases} 
\end{eqnarray}}\unskip
with $\vec {\ddot{x}}_{i}$ and $\vec {\ddot{\varphi }}_{i}$ represent
respectively the translational and rotational accelerations of the
grain $i$, $m_{i}$ and $I_{i}$ are the mass and the moment of   
inertia respectively, ${\vec {F}}_{ij}^{\mathrm{contact}}$ is the
interaction force applied by a grain $j$ in contact with grain $i$,
$\vec {g}$ is the acceleration due to gravity and ${\vec
{M}}_{ij}^{\mathrm{contact}}$ is the torque with respect to the center
of gravity of grain $i$ applied by the contact forces ${\vec
{F}}_{ij}^{\mathrm{contact}}$.

Contact forces are decomposed into a normal component $\vec {F}_{n}$
and a tangential component $\vec {F}_{s}$ contained in the plane
tangent to the grain at the point of contact \cite{19}. In this work,
the normal contact force is computed through the viscoelastic model
(Figure~\ref{fig1}a) as:
{\begin{eqnarray}\label{eqn2}
\vec {F}_{n}=(-k_{n}D_{n}-\nu _{n}V_{n})\cdot\vec {n}
\end{eqnarray}}\unskip
where $k_{n}$ is the elastic constant, $v_{n}$ is the viscous damping
constant, $V_{n}$ is the normal velocity and $D_{n}$ is the overlap of
the two grains $i$ and $j$, defined geometrically by:
{\begin{eqnarray}\label{eqn3}
D_{n}=\Vert\vec {x}_{j}-\vec {x}_{i}\Vert -r_{i}-r_{j}
\end{eqnarray}}\unskip
in which $r_{i}$ and $r_{j}$ are the radii of grains $i$ and $j$.



\begin{figure}
\includegraphics{fig01}
\caption{\label{fig1}(a) Normal force model. (b) Shear force model.}
\end{figure}


The tangential contact force $(\vec {F}_{s})$ is computed through a
viscoelastic model with friction (Figure~\ref{fig1}b) as:
{\begin{eqnarray}\label{eqn4}
\vec {F}_{s}=\min (k_{s}D_{s}+v_{s}V_{s}, \mu _{d}F_{n})\vec {s}
\end{eqnarray}}\unskip
where $k_{s}$ is the tangential stiffness, $v_{s}$ the viscous damping
coefficient, $\mu _{d}$ is the inter-particle coefficient of friction,
$D_{s}$ is the grain deformation due to shear force, $V_{s}$ is the
tangential velocity of grain $j$ with respect to grain $i$ and $\vec
{s}$ is the tangential unit vector as indicated in Figure~\ref{fig1}a.

The duration of a contact between two grains ($t_{c}$) is equal to the
natural half-period of the equivalent oscillator. Using this contact
force model, this half-period is equal to:
{\begin{eqnarray}\label{eqn5}
t_{c}=\uppi \sqrt{\frac{m_{\mathrm{eff}}}{k_{n}}}
\end{eqnarray}}\unskip
It is crucial in practice to consider this contact period $t_{c}$. In
fact, the integration of the movement equations is only stable if the
integration time step $\Delta t$ is sufficiently small in comparison
to $t_{c}$, i.e.\ the evolution of the contact should be accurately
described. In practice, in order to prevent certain instabilities, the
integration time step $\Delta t$ is taken as \cite{19}:
{\begin{eqnarray}\label{eqn6}
\Delta t_{\max}\approx 0.1\uppi \sqrt{\frac{m_{\mathrm{eff}}}{k_{n}}}
\end{eqnarray}}\unskip
where $m_{\mathrm{eff}}$ is the smallest effective mass in the system.

In order to account for the damping of the rolling motion of natural
grains, a rolling resisting torque is incorporated in the model. This
torque is calculated by placing the normal force off-center of the
contact, such eccentricity is adjusted to reduce or increase rolling
resistance.

It should be noted that in order to obtain quantitative results
comparable to experimental ones, a calibration of the model parameters
is necessary. However, for a phenomenological analysis such as in our
case, these parameters are usually selected in such a way as to save
computation time without departing from the normal behavior of the
model.       


\section{Analytical solution of wave propagation in a viscoelastic deposit}

For some mechanical and/or movement conditions where the deposit can be
assumed as a continuum with elastic damped behavior, the equation
governing the propagation of a horizontal harmonic displacement ($u$)
in the vertical direction ($z$) is given by \cite{20} as:
{\begin{eqnarray}\label{eqn7}
\frac{1}{c^{2}}\frac{\partial ^{2}u}{\partial t^{2}}=\frac{\partial ^{2}u}{\partial z^{2}}+\frac{2\xi }{\omega }\frac{\partial ^{3}u}{\partial t\partial z^{2}}
\end{eqnarray}}\unskip
where $c$ is wave propagation velocity, $\xi $ is the hysteretic
damping ratio and $\omega $ is the circular frequency of the movement.
For a deposit of thickness $h$ resting on an oscillating bedrock, the
final solution of this equation is presented in \cite{20} as:
{\begin{eqnarray}
\frac{A \cdot u(z,t)}{u_{0}}&=&\cosh (ph)\cos (qh)\cosh (pz)\cos (qz)\sin (\omega t)  
\nonumber\\&&
+\,\sinh (ph)\sin (qh)\sinh (pz)\sin (qz)\sin (\omega t)  
\nonumber\\&&
+\,\cosh (ph)\cos (qh)\sinh (pz)\sin (qz)\cos (\omega t) 
\nonumber\\&&
-\,\sinh (ph)\sin (qh)\cosh (pz)\cos (qz)\cos (\omega t)
\label{eqn8}
\end{eqnarray}}\unskip
with $A$, $p$ and $q$ are calculated from equations:
{\begin{eqnarray}
& \displaystyle
A=\cosh^{2}(ph)- \sin^{2}(qh)
\label{eqn9}\Seqnsplit
& \displaystyle
(ph)^{2}+(qh)^{2}=\frac{(\uppi ^{2}/4)(\omega ^{2}/{\omega }_{1}^{2})}{1+4\xi ^{2}}
\quad \mbox{and}\quad
(ph) (qh)=\xi \frac{(\uppi ^{2}/4)(\omega ^{2}/{\omega }_{1}^{2})}{1+4\xi ^{2}}                                     
\label{eqn10}
\end{eqnarray}}\unskip

The fundamental frequency of the deposit being the first resonance
frequency for very low damping, it is given by:
{\begin{eqnarray}\label{eqn11}
\omega _{1}=\frac{\uppi c}{2h}
\end{eqnarray}}\unskip
This solution indicates that for a harmonic excitation, the
displacement of any point of the deposit is harmonic. In addition, it
states that the shape of the displacement, i.e.\ the plot of the
displacement versus of the depth at a given time, depends on the
excitation frequency. As examples, Figure~\ref{fig3} shows the displacement
shapes for a deposit excited at its base by a harmonic displacement of
amplitude $u_{0}$ and with excitation frequencies $\omega /\omega
_{1}=0.42$ and $\omega /\omega _{1}=1.12$. For these examples, the
damping ratio is set to $\xi =8$\%. Note that the system's damping
ratio can be estimated from the dynamic amplification factor at
resonance as shown below.

This figure shows that for an excitation frequency lower than the
fundamental frequency, the displacement profile has the shape of the
first natural oscillation mode (Figure~\ref{fig3}a). However, when the
excitation frequency exceeds the fundamental frequency, the second mode
shape appears by the formation of a zero-displacement node at a
specific height that rises from the base. This height increases with
the increase in the excitation frequency (see Figure~\ref{fig2}b) until
the second natural frequency is reached, then the shape of the third
mode appears and so on. 


\begin{figure}
\includegraphics{fig02}
\caption{\label{fig2}Normalized displacements of analytical solution for two
different excitation frequencies.}
\end{figure}


\begin{figure}
\includegraphics{fig03}
\caption{\label{fig3}Dynamic Amplification Factor at the deposit's free surface
(DAF) in terms of the frequency ratio $(\omega /\omega
_{1})$ damping ratios: 5\%, 10\% and 15\%.}
\end{figure}


In this study, we define the dynamic amplification factor at the top
(i.e.\ the free surface) as $\mathrm{DAF}=\sfrac{u_{s}}{u_{0}}$, where $u_{s}$ and
$u_{0}$ are the displacement amplitudes at the free surface and at the
base respectively. Using the displacement solution (Equation~(\ref{eqn8})),
the DAF is plotted in Figure~\ref{fig3} in terms of the
frequency ratio $(\omega /\omega _{1})$ ranging from 0 to 1.12 and for
damping ratios of 5\%, 10\% and 15\%.


This figure clearly shows the typical evolution of amplification
against the excitation frequency increasing. The DAF
increases starting from 1 for low frequencies, it reaches a maximum
value at the first resonance frequency then it decreases. The maximum
DAF value occurring at resonance is inversely proportional
to the damping ratio, it can be approximated for low damping ratios by
$\mathrm{DAF}_{\max}\approx \sfrac{2}{\uppi \xi }$~\cite{20}. It is worth noting
that for common dampings, the resonance occurs at frequencies very
close to the fundamental frequency.


\section{Details of the numerical model of wave propagation}

In the present work, an analysis of the vertical propagation of a shear
wave in a sand deposit is carried out using 2D discrete element
modeling, the used model is implemented in a C${{+}{+}}$ code. In a first step,
a dense deposit is built by pluviation under gravity on a horizontal
rigid plateau (considered as a bedrock), with introducing periodic
boundaries to model an infinite layer in the horizontal direction.
After stabilization, the positions of the grains of the deposit are
stored to become the reference positions for displacements. Next, the
deposit is subjected to a horizontal vibration though the supporting
bedrock, while the upper surface of the deposit is kept free. In order
to ensure good transmission of movement to the deposit, the bedrock is
roughened by bonding the lowest deposit grains to it. Therefore, these
grains undergo the same imposed bedrock displacement. In order to
describe the response, the deposit is divided into a number of layers
for which the average grain displacements are monitored during the
excitation.

The properties of the grains as well as the micromechanical parameters
used in the computation of the intergranular contact forces are
presented in Table~\ref{tab1}. 


\begin{table}
\caption{\label{tab1}Model parameters}
\begin{tabular}{cc}
\thead
Parameters & Value \\ 
\endthead
Number of grains  & 5000  \\ 
Grains radii $r_{i}$ & $1.00 \times 10^{-3}$ to $2.00 \times 10^{-3}~\mathrm{m}$  \\ 
Density $G_{s}$ & $2600~\mathrm{kg}/\mathrm{m}^{3}$ \\ 
Normal intergranular stiffness $k_{n}$ & $1.2\times 10^{6}$~N/m  \\ 
Shear intergranular stiffness $k_{s}$ & $9.6\times 10^{5}$~N/m  \\ 
Intergranular friction coefficient $\mu $ & 0.5   
\botline
\end{tabular}
\end{table}



The width of the deposit (i.e.\ the period) is initially set to 2~m,
the resulting height for the used number and size of the grains is
about 0.54~m. Figure~\ref{fig3} shows the obtained sand deposit model after
pluviation with the different boundary conditions as well as the bonded
chain of grains called ``excitation chain''. The subdivision of the
deposit into control layers is also illustrated in the same figure.



In this work, the applied excitation is a harmonic horizontal
displacement of specified amplitude and frequency          
$u (t)=u_{0}\sin (\omega t)$, where  $u_{0}$ is the amplitude
and $\omega $ is the circular frequency. In order to avoid the
transient response disturbances, the excitation is applied in three
stages according to the same idea presented in reference \cite{4}. In
the first step, the amplitude is gradually increased from zero to the
selected steady-state vibration amplitude. In the second stage, the
amplitude is maintained constant to produce a steady state oscillation.
Finally, in the third stage the amplitude is decreased from the
steady-state amplitude until zero. In the following simulations, the
durations of the three stages are respectively set to 4.5~s, 6~s
and 1.5~s. In this way, the transient response should become almost
attenuated in the steady-state oscillation phase. The results that will
be discussed below, essentially belong to this later phase. The deposit
is divided into 10 control layers and the recorded displacement is the
average displacement over all grains belonging to the layer.


\section{Simulation results and discussion}

\subsection{Displacement time history and movement amplification}

In this simulation, the parameters of the exciting displacement are set
to $u_{0}=4 \times 10^{-4}$~m  and $\omega =25$~rad/s. In the first
stage (i.e.\ from $0$ to 4.5~s), the displacement amplitude is increased
proportionally to time from zero until $4\times 10^{-4}$~m. Then,
maintained at the same level for 6~s and finally decreased gradually to
reach zero during 1.5~s. Figure~\ref{fig4} shows the displacement time
history at the bedrock (Figure~\ref{fig5}a) and at the top layer of the
deposit (Figure~\ref{fig5}b).


\begin{figure}
\includegraphics{fig04}
\vspace*{2pt}
\caption{\label{fig4}Illustration of the sand deposit's model, boundary conditions
and subdivision into a number of monitoring layers.}
\vspace*{4pt}
\end{figure}



\begin{figure}
\includegraphics{fig05}
\caption{\label{fig5}Displacement time history, (a) at the bedrock (excitation),
(b) at the top layer of the deposit.}
\end{figure}


These figures clearly show the different phases of the movement at the
bedrock and at the top layer (which represents approximately the
surface of the deposit). Figure~\ref{fig4}b shows that the amplitude at
the deposit's surface in the steady state phase is almost constant,
i.e.\ the effect of the transient response disappears in this phase. The
comparison of the two displacement plots (Figures~\ref{fig5}a, b)
in the steady state phase, shows that the
introduced displacement at the bedrock (excitation) having an amplitude
of 0.4~mm is amplified in the top layer where the amplitude becomes
about 0.85~ mm. Besides, Figure~\ref{fig4}b shows that the deposit
experiences a free elastic damped oscillation after ending the
excitation (i.e.\ after $t=12$~s). 

The amplification of the movement in the bulk of the deposit is
illustrated in Figures~\ref{fig5}a,~b. The first shows a
part of the steady state response $(7~\mathrm{s}\leq t \leq
8~\mathrm{s})$ at the base, the middle and the top layers, and the
second presents the vibration amplitude versus the layer's height from
the bedrock ($y$-coordinate). 






These figures indicate that the movement amplification increases from
the base up to the top layer. On the other hand, Figure~\ref{fig6}b shows that
the displacement profile has the characteristic shape of the first mode
of vibration; this reveals that the excitation frequency is lower than
the first resonance frequency of the deposit.


\begin{figure}
\includegraphics{fig06}
\caption{\label{fig6}(a) Displacement time history for steady state oscillation at
three different depths. (b) Displacement amplitude versus the layer's
height from the bedrock during steady state oscillation.}
\vspace*{6pt}
\end{figure}


\subsection{Fundamental frequency of the deposit}\label{sec5.2}

Soils are known for their non-linear mechanical behavior, their shear
modulus degrades when the level of strain involved increases, as shown in 
{Figure~\ref{fig7}a.} In this figure, $G_{\max}$ is the small strain
shear modulus and $G$ is the modulus corresponding to a specified
strain $\gamma $. Figure~\ref{fig7}b shows the degradation curve of $G$
in terms of the involved level of strain ($\gamma _{0}$) according to
reference \cite{21}. It is clear that $G_{\max}$  corresponds to
strains approximately smaller than $10^{-5}$.

\begin{figure}
\includegraphics{fig07}
\caption{\label{fig7}(a) Illustration of the degradation of the shear modulus with
increasing applied strain. (b) Shear modulus degradation curve versus the
strain level \cite{21}.}
\end{figure}



The fundamental frequency is defined as the first natural frequency of
a small strain (pseudoelastic) oscillation. To obtain this frequency,
we subjected the deposit to a harmonic vibration in the same way
described above, with the amplitude $u_{0}=4\times 10^{-4}~\mathrm{m}$
and the frequency of 25~rad/s. The fundamental frequency of the deposit
is taken as the frequency of free vibration remaining after the end of
the excitation. Figure~\ref{fig8} presents a zoom of the time history
of displacements corresponding to the free vibration phase (after the
end of excitation) at different depths of the deposit. This figure
shows that the displacements are harmonic and without a phase offset
for the different depths, which corresponds to the oscillation shape of
the fundamental mode. The fundamental period can be estimated from the
displacement plot, it is about $T=0.166~\mathrm{s}$. Hence, the
fundamental  frequency is  $\omega _{1}=\sfrac{2\uppi
}{T}=37.82~\mathrm{rad}/\mathrm{s}$.

\begin{figure}
\includegraphics{fig08}
\caption{\label{fig8}Zoom of time history response at different depths for free
vibration phase.}
\end{figure}


For an elastic oscillation of the deposit, the shear wave velocity can
be calculated from the fundamental period through the expression
$c=\sfrac{4h}{T}$~\cite{20}, hence for our deposit 
$c=12.91~\mathrm{m}/\mathrm{s}$. 

Besides, it could be noted that the wave propagation velocity for this
model is lower than the characteristic velocities of natural soils.
This is reasonable due to the low intergranular stiffness deliberately
adopted in the model in order to reduce the calculation time. However,
it appears that this wave velocity has a negligible effect on the
phenomenological aspects of the response. O'Donovan et~al.~\cite{6},
following their simulations on wave propagation in granular media,
indicated that the change in intergranular stiffness changes the
propagation velocity but has no noticeable effect on the nature of the
response.



\subsection{Influence of the excitation frequency on the response of the deposit}

In order to show the excitation frequency effect on the deposit's
response, several simulations are performed for frequencies ranging
from 10 to 40~rad/s, while the excitation amplitude is set to $u_{0}=4
\times 10^{-4}~\mathrm{m}$.  Figure~\ref{fig9} shows the time history
of displacements during the steady state excitation stage at the
bedrock and at the top layer of the deposit for frequencies of 10 and
20~rad/s. It can be noted that for the frequency of 10~rad/s, there is
no movement amplification at the top layer, thus the deposit exhibits a
behaviour of a rigid body that follows the bedrock's displacement.
However, for the frequency of 20~rad/s, the movement is amplified at
the top layer, therefore this frequency causes the deformation of the
deposit.

\begin{figure}
\includegraphics{fig09}
\caption{\label{fig9}Steady state displacements time history at the bedrock and at
the top layer of the deposit for excitation frequencies: (a) 10~rad/s,
(b) 20~rad/s. Excitation amplitude $u_{0}=4 \times
10^{-4}$~m.}
\vspace*{1.5pc}
\end{figure}

\begin{figure}
\includegraphics{fig10}
\caption{\label{fig10}Time history displacements of the top layer for the excitation
frequencies 20, 30, 32 and 35~rad/s. Excitation amplitude 
$u_{0}=4 \times 10^{-4}~\mathrm{m}$.}
\end{figure}

Figure~\ref{fig10} shows the displacement time history of the top layer
for the entire excitation duration. Even though the displacement time
histories were plotted for all frequencies, only the plots
corresponding to frequencies 20, 30, 32 and 35~rad/s are shown. In
fact, these frequencies seem to clearly trace the evolution of
behaviour. It can be observed that going from 20~rad/s to 30~rad/s
there is an increase in the amplification of the movement; the
amplitude of the top layer goes from 0.6~mm for 20~rad/s to
approximately 1.6~mm for 30~rad/s. At frequencies of 32~rad/s and
35~rad/s, and in particular for the top layer but not for the other
layers, the displacement plot shapes indicate the occurrence of
irreversible grain sliding. This sliding could be attributed to the
strong movement amplification and the non-confining of the considered
layer (top layer). Note that this sliding begins when the vibration
amplitude becomes sufficiently high, i.e.\ around the end of the
excitation increasing stage. Thereafter, the layer continues to
oscillate around a new position which may evolve over time. It may be
noted from the plots, that the vibration amplitude tends to undergo a
slight decrease with time; this behaviour may be due to the
densification of the deposit which results from the significant
relative movements of the grains. Note that densification leads to a
reduction in the deformability of the deposit.

Meanwhile, for frequencies where top layer sliding has occurred, the
steady-state vibration phase is considered to be the short phase
immediately after the amplitude increase phase (i.e.\ after 4.5~s), it
is shown in the time history displacement plots (red windows in
Figure~\ref{fig10}). This phase is shortened for these frequencies and
it should correspond to the case where the sample is close to its
initial state before densification. To account for the possible shift
from the initial position (zero position), the amplitude of the
displacement is calculated as half the difference between the maximum
and minimum peak values in the taken time interval.

The displacement profiles (i.e.\ the deformation shape of the deposit at
a given time) for the different excitation frequencies are represented
in Figure~\ref{fig11} for the steady state vibration at a positive peak
time and the next negative peak time (positive peak time $+$ a half
period). 

\begin{figure}
\includegraphics{fig11}
\caption{\label{fig11}Displacements profiles for different excitation frequencies.
Excitation amplitude $u_{0}=4 \times 10^{-4}~\mathrm{m}$.}
\vspace*{-3pt}
\end{figure}


This figure demonstrates that for low frequencies (up to 30~rad/s), the
displacement profiles have the shape of the first vibration mode. For
the frequency of 32~rad/s, we notice the appearance of the node with
zero displacement at the base then it rises for the frequency of
35~rad/s. According to the relative discussion in Section~\ref{sec2},
this indicates reaching and then exceeding the first resonant
frequency. These shapes of displacement profiles indicate that the
first resonance frequency for the excitation amplitude $u_{0}=4 \times
10^{-4}~\mathrm{m}$, is around 32~rad/s.

We define in the following the dynamic amplification factor
(DAF) at the top layer in the steady state vibration phase as:
{\begin{eqnarray}\label{eqn12}
\mathrm{DAF}=\frac{u_{s}}{u_{0}}
\end{eqnarray}}\unskip
where $u_{s}$ denotes the displacement amplitude at the top layer
(surface layer) and $u_{0}$ is the displacement amplitude at the
bedrock (excitation). 

In order to show the evolution of the amplification with frequency, the
DAF as defined here is represented in Figure~\ref{fig12}a in
terms of the excitation frequency. 

\begin{figure}
\includegraphics{fig12}
\caption{\label{fig12}Dynamic Amplification Factor versus excitation frequency for
excitation amplitude $u_{0}=4 \times 10^{-4}~\mathrm{m}$. (a) DAF
obtained from numerical modelling. (b) DAF comparison
between numerical modelling and analytical solution.}
\end{figure}


Figure~\ref{fig12}a shows that the DAF evolves in a
similar way to the case of an elastic deposit excited by a harmonic
displacement at the base; with the increase of the excitation
frequency, the DAF increases up to a maximum value that
corresponds to the resonance then it decreases. The resonance frequency
that can be estimated from this graph is around 32~rad/s, which agrees
well with the displacement profiles in Figure~\ref{fig11}.

For a comparison with the analytical solution, Figure~\ref{fig12}b
shows the evolution with frequency of the DAF obtained numerically and
that calculated from the analytical solution presented~above. The data
used in this solution are those corresponding to the numerical deposit,
namely $h=0.54~\mathrm{m}$ and $\omega
_{1}=37.82~\mathrm{rad}/\mathrm{s}$. The damping ratio is chosen so
that the maximum value (at resonance) is equal to that obtained from
the numerical model which is approximately $\mathrm{DAF}_{\max}=4.75$. Using the
equation $\mathrm{DAF}_{\max}\approx \sfrac{2}{\uppi \xi }$, we obtain the damping
ratio $\xi =13.4$\%. In this figure, the excitation frequencies are
normalized by the fundamental frequency of model $\omega _{1}$. This
figure shows that the evolution of the numerically computed
DAF is very similar to that of the analytical solution.
However, resonance occurs at a frequency about 13\% lower. We believe
that this difference can be attributed mainly to the degradation of the
shear modulus due to the increase in shear strains near resonance.

The analysis of the effect of excitation frequency on deposit response
is repeated with an excitation amplitude of $u_{0}=4 \times
10^{-6}~\mathrm{m}$ in the same way as for the amplitude $4 \times
10^{-4}~\mathrm{m}$. It is observed that for all frequencies, there is
no sliding of the surface layer, and that the movement remains harmonic
and regular for all layers and frequencies. The plot of the
displacement profiles showed that the resonance frequency increases to
around 37.5~rad/s. For this case we present only the evolution of the
DAF with excitation frequency (Figure~\ref{fig13}), in the same way as
Figure~\ref{fig12}b. 

\begin{figure}
\includegraphics{fig13}
\caption{\label{fig13}Dynamic Amplification Factor versus excitation frequency for
excitation amplitude $u_{0}=4 \times 10^{-6}~\mathrm{m}$.  DAF obtained
from numerical modelling and analytical solution.}
\end{figure}


This figure shows that the resonance frequency becomes very close to
the fundamental frequency estimated in Section~\ref{sec5.2} ($\omega
_{1}=37.82~\mathrm{rad}/\mathrm{s}$). In addition, the maximum dynamic
amplification factor increases significantly ($\mathrm{DAF}_{\max}=17.08$).
These results indicate that there is almost no degradation of the shear
modulus and that the damping ratio becomes very low for small
amplitudes. We believe that the damping rate depends on intergranular
slippage, which decreases with decreasing excitation amplitude. It can
therefore be concluded that the resonance frequency depends on the
excitation amplitude.


\subsection{Influence of the excitation amplitude on the movement
amplification}

It should be remembered that for an elastic deposit excited at the
bottom by a harmonic displacement, the \textit{FAD} is independent of
the amplitude of the excitation. However, as the granular deposit does
not behave elastically in all situations, it becomes interesting to
show how it behaves for different excitation amplitudes. For this
purpose, simulations of the deposit response under harmonic excitation
at fixed frequency (25~rad/s) and with amplitudes ranging from $2\times
10^{-6}~\mathrm{m}$ to $1\times 10^{-3}~\mathrm{m}$ are carried out.
For all simulations, the displacements in the steady-state phase are
plotted and the DAF as defined above is calculated.


Figure~\ref{fig14} shows the evolution of the DAF with the
excitation amplitude. In this plot, the logarithmic scale is used for
the $x$-axis (amplitudes) given the wide range of amplitude values.

\begin{figure}
\includegraphics{fig14}
\caption{\label{fig14}Evolution of the DAF at the top layer versus excitation
amplitude. Excitation frequency 25~rad/s.}
\end{figure}


It is clear that unlike the case of an elastic deposit where the
DAF is theoretically constant, for the present granular
deposit the DAF is almost constant for small amplitudes,
but it increases sharply for large excitation amplitudes. This
behaviour could be attributed to the degradation of the shear modulus
due to the increase of the shear strain level when the amplitude
increases. In order to check this presumption, we were interested in
evaluating the shear strain level involved for the different excitation
amplitudes. Observing that the displacement profile has an almost
linear shape for the used excitation frequency (25~rad/s) as shown in
Figure~\ref{fig11}, the shear strain level ($\gamma_{0}$) can be approximated
by the ratio of the difference between the horizontal displacements of
two layers to the vertical distance between them. In order to obtain an
average value over the whole deposit, $\gamma_{0}$ is calculated in
the following using the penultimate layer (layer 9) and the first
layer. Considering the possible phase shift (even very small) between
the displacements of the layers, $\gamma _{0}$ is expressed as:
{\begin{eqnarray}\label{eqn13}
\gamma _{0}=\frac{{u}_{9}^{p}-u_{1}}{y_{9}-y_{1}}
\end{eqnarray}}\unskip
where ${u}_{9}^{p}$ denotes a peak displacement of layer 9, $u_{1}$
denotes displacement of layer 1 at the same time (corresponding to
$u_{9}^{p}$) and $y_{i}$ is the vertical coordinate of the $i$ layer's
center. Average shear strain levels are calculated for all the
excitation amplitudes. Figure~\ref{fig15}a shows the evolution of the
shear strain level with the excitation amplitude, while
Figure~\ref{fig15}b shows the DAF versus the shear strain
level involved by the different excitation amplitudes.

\begin{figure}
\includegraphics{fig15}
\caption{\label{fig15}(a) Average shear strain level versus excitation amplitude. (b)
DAF versus the shear strain level involved by the different
excitation amplitudes.}
\end{figure}

Figure~\ref{fig15}a shows that the shear strain level increases
monotonically (almost linearly) with increasing excitation amplitude.
Taking into account the shear modulus degradation with the shear strain
level curve (Figure~\ref{fig7}b), this increase reveals that the shear
modulus of the deposit decreases with the increase of the excitation
amplitude, which leads to a change in the DAF. Such a change
is well demonstrated by Figure~\ref{fig15}b, this later shows that the
DAF increases very slightly for small deformations, while it
takes on an increasingly significant variation slope for large strains.
This aspect is very similar to the shear modulus degradation curve.
Furthermore, the strain threshold which separates the two phases of
variation of the DAF agrees well with that which can be drawn
from the degradation curve of the shear modulus (Figure~\ref{fig7}b),
it is about $10^{-4}$. This result indicates that for a granular
deposit, the movement can be largely amplified even when far from the
fundamental frequency. This is due to the degradation of the shear
modulus, which degrades with increasing excitation amplitude.

\subsection{Effect of the confining on the movement amplification}

When a layer of granular soil is deep and subjected to the weight of
the upper layers, it becomes confined and its motion transmission
properties change. This section aims to highlight the effect of
confinement on the amplification of movement at the free surface of a
granular deposit excited at its base. A simple way to simulate
confinement is to increase the gravity acceleration, so in this work we
analyzed the evolution of the DAF for gravity accelerations
ranging from g up to 20~g. The properties of the model are the same as
those presented above and the base excitation is harmonic with an
amplitude of 0.4~mm and a frequency of 25~rad/s. 

Figure~\ref{fig16} shows the evolution of the DAF as defined
above, with the increase of the gravity acceleration in 
$(\mathrm{m}/\mathrm{s}^{2})$.

\begin{figure}
\includegraphics{fig16}
\caption{\label{fig16}Evolution of the DAF with the increase of the gravity acceleration.}
\end{figure}

It is clear that the increase in the gravity acceleration decreases the
amplification of the displacement at the free surface. Indeed,
increasing gravity leads to an increase in effective stresses and
therefore more frictional forces and less sliding, which makes the
granular deposit more rigid and therefore reduces the amplification of
surface movement. It can be noted that the DAF tends towards
unity for large values of gravitational acceleration. This result
agrees with theoretical knowledge according to which very stiff
deposits (like rocks) do not amplify seismic movements.

\section{Conclusion}

This work consisted of a numerical analysis of the vertical propagation
of shear waves in a sand deposit through 2D discrete elements modeling.
The used discrete element model is based on the molecular dynamics
method with circular shaped elements. The intergranular normal forces
at contacts are calculated through a linear viscoelastic law while the
tangential forces are calculated through a viscoelastic perfectly
plastic model. Rolling friction is incorporated to account for the
damping of the grains rolling motion. The exciting motion is introduced
through the deposit's base (bedrock), the deposit's upper surface is
free and periodic boundaries are implemented on lateral boundaries. The
used excitation is harmonic with variable frequencies and amplitudes.
Qualitatively, the simulations carried out highlighted some behavioral
aspects of the vibration of non-cohesive deposits, such as the shape of
the vertical displacement profile and the large displacements near the
resonance. In the quantitative aspect, first the deposit's fundamental
frequency is obtained from the free vibration at low amplitudes.
Thereafter, the analysis of the effect of the excitation frequency on
the response of the deposit showed in the range of the studied
frequencies, that the Dynamic Amplification Factor (DAF)
increases starting from 1 for low frequencies, it reaches a maximum
value at the first resonance frequency then it decreases. This
analysis was used to estimate the resonance frequency. It appeared that
the latter is close to the fundamental frequency for low excitation
amplitudes, but becomes smaller as the excitation amplitude increases.
This finding may be mainly due to the degradation of the shear modulus
as shear strains increase with increasing excitation amplitude. Unlike
elastic deposits where the DAF is theoretically insensitive to
the amplitude, the analysis of the effect of excitation amplitude on
motion amplification confirmed that for a granular deposit, the
DAF increases with increasing excitation amplitude. It has
been shown that this behavior is due to the degradation of the shear
modulus following the increase in shear strain. Furthermore, it is
found that the strain threshold beyond which degradation becomes
significant is in the order of $10^{-4}$. Finally, the analysis of the
confinement effect on the response of the deposit, is carried out by
increasing the gravitational acceleration. It is shown that increasing
the confinement makes the granular deposit stiffer and reduces the
movement amplification at the free surface.

\section*{Ethical Approval} 
Not applicable. 

\section*{Declaration of interests} 
The authors do not work for, advise, own shares in, or receive funds
from any organization that could benefit from this article, and have
declared no affiliations other than their research organizations. 

\section*{Authors' contributions} 
All authors participated equally in all aspects of the study, including
modeling, writing, and reviewing the manuscript.

\section*{Availability of data and materials} 
The authors declare that the data supporting the findings of this study
will be made available on reasonable request. 

\back{}

\printbibliography
\refinput{crmeca20240483-reference.tex}

\end{document}
