# Hyperelastic or Hypoelastic Granular Circular Chain Instability in a Geometrically Exact Framework

Noël Challamel, François Nicot, Antoine Wautier, Félix Darve, Jean Lerbet

► To cite this version:

Noël Challamel, François Nicot, Antoine Wautier, Félix Darve, Jean Lerbet. Hyperelastic or Hypoelastic Granular Circular Chain Instability in a Geometrically Exact Framework. Journal of Engineering Mechanics - ASCE, 2022, 148 (9), pp.04022053. 10.1061/(ASCE)EM.1943-7889.0002139. hal-04012052

HAL Id: hal-04012052

https://hal.inrae.fr/hal-04012052

Submitted on 2 Mar 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Hyperelastic or Hypoelastic Granular Circular Chain Instability in a Geometrically Exact Framework

Noël Challamel, M.ASCE<fnref n="1" />; François Nicot<fnref n="2" />; Antoine Wautier<fnref n="3" />; Félix Darve<fnref n="4" />; and Jean Lerbet<fnref n="5" />

**Abstract:** This paper investigates several granular interaction laws used in the modeling of discrete granular media. In the considered model, each grain interacts with its neighbors with a coupled shear-normal interaction law. The analysis is performed in a geometrically exact framework allowing large rotation and displacement evolutions, without any geometrical approximations. It is shown that most of the granular interaction laws available in the literature are classified as hypoelastic interaction laws, and we precise the requirements to build some hyperelastic interaction laws that avoid artificial dissipation. We also show that the uncoupled granular interaction law is hyperelastic for all the studied models. The analysis is applied to a paradigmatic elementary system of a granular loop with a diamond pattern (a four-grain cyclic granular chain) loaded by concentrated forces. Instabilities are observed for large displacement of the diamond chain for all the classified models. It is observed that the discrepancies between each model may grow during the deformation process. The instability phenomenon is associated with the appearance of a limit load for this granular structural problem due to large nonlinear geometrical effects. Blocking phenomena may also appear for such granular structural systems due to secondary granular contacts. DOI: 10.1061/(ASCE)EM.1943-7889.0002139. © 2022 American Society of Civil Engineers.

**Author keywords:** Discrete element method; Granular interactions; Dissipative phenomena; Hyperelastic interactions; Instability; Circular pattern; Granular chain; Geometrical nonlinearity.

## Introduction

The modeling of granular media using the so-called Discrete Element Method (DEM) emerged in the 1970s with the development of computational facilities (Cundall 1971). DEM applied to discrete granular media is based on a modeling of each grain as a rigid body (each grain can be idealized by a rigid disk in a two-dimensional framework, or a rigid sphere in three dimensions), which interacts with its neighbor through normal and shear

fundamental interaction laws (Serrano and Rodriguez-Ortiz 1973; Cundall and Strack 1979). The principle of this physically based numerical method is to capture the complex macroscopic response of granular media through elementary granular interactions. In most of the available codes used currently, the interaction law between two adjacent grains is defined in a rate form, generally based in its simplest formulation, on a linear relation between the rate of elementary forces at the grain interface (rate of normal force and shear force) to the normal component and tangential components of the relative velocities between two points belonging to the interfacial region between the two grains. This linear relation involves some stiffness parameters, namely the normal and tangential stiffness parameters. Inelastic granular interaction laws are commonly used to upscale the complex features of grain interactions to the macroscopic scale. However, even in the linear range for the granular interaction, some unexpected results are found, and still poorly documented. The pioneer formulations of Serrano and Rodriguez-Ortiz (1973) and Cundall and Strack (1979) are only based on rate equations, and are not introduced from energetic arguments. In fact, it can be shown on simple examples (this will be done in this paper), that this rate formulation is not a total differential of the state variables, so that these models can be considered only as hypoelastic. In such a case, the model can be also classified as strongly hypoelastic (see Lerbet et al. 2018 for this terminology). It is not possible with such models to express the internal force variables in terms of direct granular kinematic parameters by integration. This may lead to artificial energy dissipation phenomenon through ratcheting effects, even if the granular interaction law is linear. The so-called paradox of Green (Green and Naghdi 1971; Knops 1982) with infinite dissipation in closed kinematic cycles may be observed for such discrete granular models. McNamara et al. (2008) specifically highlighted this paradox for the hypoelastic model of Cundall and Strack (1979) on closed kinematic cycles. McNamara et al. (2008) proposed a correction for the initial shear-normal interaction law by relating the shear force multiplying by the distance between the two grains, to the relative rotations

<footnote n="1">¹Professor, Univ. Bretagne Sud, Institut de Recherche Dupuy de Lôme—Université Bretagne Sud—Unité Mixte de Recherche—Centre National de la Recherche Scientifique 6027, Centre de Recherche, Rue de Saint Maudé-BP 92116, Lorient cedex 56321, France (corresponding author). ORCID: https://orcid.org/0000-0002-7122-0700. Email: noel.challamel@univ-ubs.fr</footnote>

<footnote n="2">²Research Director, Université Savoie Mont Blanc, Environnements Dynamiques et Territoires de la Montagne, Unité Mixte de Recherche—Centre National de la Recherche Scientifique 5204, Le Bourget du Lac 73376, France. Email: francois.nicot@irstea.fr</footnote>

<footnote n="3">³Researcher, Aix-Marseille Université, Institut National de Recherche pour l'Agriculture, l'Alimentation et l'Environnement, Risques, ECOsystèmes, Vulnérabilité, Environnement, Résilience, Aix-en-Provence 13182, France. ORCID: https://orcid.org/0000-0002-2551-103X. Email: antoine.wautier@inrae.fr</footnote>

<footnote n="4">⁴Professor, Univ. Grenoble Alpes, Grenoble Institut National Polytechnique, Centre National de la Recherche Scientifique, lab 3SR, Grenoble 38000, France. ORCID: https://orcid.org/0000-0002-1276-1929. Email: felix.darve@3sr-grenoble.fr</footnote>

<footnote n="5">⁵Professor, Laboratoire de Mathématiques et Modélisation d'Evry, Unité Mixte de Recherche—Centre National de la Recherche Scientifique 8071, Univ. Evry, Université Paris-Saclay, Evry 91037, France. Email: jlerbet@gmail.com</footnote>

<footnote>Note. This manuscript was submitted on February 22, 2022; approved on May 2, 2022; published online on July 12, 2022. Discussion period open until December 12, 2022; separate discussions must be submitted for individual papers. This paper is part of the *Journal of Engineering Mechanics*, © ASCE, ISSN 0733-9399.</footnote>
between each grain. With such a correction, the granular interaction becomes hyperelastic, in the sense that it can be derived from a potential. This model is equivalent to consider that the grain interactions is composed of normal translational springs with shear rotational springs (see also Challamel 2015). This granular interaction model can be also classified as a Green elasticity model (see again Green and Naghdi 1971; Knops 1982 or more generally on elasticity concepts, Rajagopal 2011). The model of McNamara et al. (2008) is variationally founded and does not present the paradox of possible artificial dissipation in close cycles. Furthermore, it gives the same results as the hypoelastic DEM model of Cundall and Strack (1979) for uncoupled granular motion (pure normal or pure shear interaction modes).

Recently, Turco (2018) and Turco et al. (2019) developed another hyperelastic formulation based on an energy expression constituted of normal, shear, and bending kinematic interactions. Turco (2018) studied the shear behavior of packing of nonuniform circular grains, and Turco et al. (2019) mostly focused on uniform circular grains with and without geometrical defects (based on the same granular interaction model). When restricted to pure shear and normal interactions, this granular interaction model is equivalent to consider normal translational springs with shear translational springs (see also Challamel 2015 who also presents such hyperelastic model compared to the one based on shear rotational springs). The model can be also classified as a Green elasticity model (see again Green and Naghdi 1971; Knops 1982 or more generally on elasticity concepts, Rajagopal 2011). The model of Turco et al. (2019), or equivalently the model of Turco (2018), is also variationally founded and does not present the paradox of possible artificial energy dissipation in close cycles. However, this hyperelastic model (model of Turco et al. 2019) presents a surprising behavior where the normal force is also affected by the rotation variables of each grain. Moreover, the shear force with such a model depends on both the shear stiffness and the normal stiffness parameters. Turco (2018) and Turco et al. (2019) restricted their analysis to two-dimensional packing of grains, whereas three-dimensional packing of grains has been recently considered by Turco (2022).

In the present paper, we will present the response of a granular loop (circular granular chain) under horizontal and vertical forces, with hypoelastic and hyperelastic granular interaction models. Starting with the incremental formulation of Serrano and Rodriguez-Ortiz (1973) or Cundall and Strack (1979), which may be also labeled as classical DEM model, we show how these strongly hypoelastic models can be modified into a weakly hypoelastic model and a hyperelastic one. Eventually, the granular interaction models may coincide with the initial DEM model in case of uncoupled interaction laws (for pure normal interaction, or pure shear interaction). By embedding the three interaction models into the four-grain cyclic chain studied in detail by Nicot et al. (2016), we analyze the relative contribution of the contact law and the nonlinear geometry effects. This evolution problem has been formulated in rate form using the initial DEM hypoelastic model, and integrated by a time-discretization process. The four-grain structural problem has been also reconsidered by Lerbet et al. (2018) with a discussion related to stability of this so-called hypoelastic model. Nicot et al. (2017) numerically studied the behavior of a predeformed granular column with periodic initial pattern, and based on the hypoelastic granular interaction law with shear and normal force coupling. The behavior of a cyclic granular chain composed of six grains has been also extensively studied by Nicot and Darve (2011) in both the elastic and the elastoplastic regimes. More recently, Challamel and Kocsis (2021) studied the stability and the bifurcation behavior of a granular column, based on a shear-bending granular interaction law, in a geometrically exact framework. The instabilities of this granular column

was previously numerically studied by Hunt et al. (2010) in presence of elastic confinement. Challamel and Kocsis (2021), following the works of Challamel et al. (2020) and Lerbet et al. (2020) also showed the link between the discrete granular problem and the finite-difference formulation of a continuous Engesser-Timoshenko column (also labeled Bresse-Timoshenko beam; see Challamel and Elishakoff 2019). In the present paper, only a few grains are considered in order to investigate analytically the nonlinear evolution problem and point out the relative effect of the type of contact law (hypoelastic or hyperelastic) on the mechanical response of an archetypal grain assembly. The classification of several granular interaction laws available in the literature in terms of hyper or hypoelasticity has not been investigated in detail, to the authors' knowledge. Furthermore, new, exact solutions for the static response of a circular granular chain composed of a few rigid grains connected by some hyperelastic or hypoelastic interaction law are presented and commented on as a paradigmatic granular system representative of some more complex granular structures.

## Description of the Granular Model

The differences between each granular interaction model will be highlighted more easily on a simple two-degree-of-freedom granular system (see Figs. 1 and 2). A regular assembly of four rigid disks initially in contact is considered (Fig. 1). This structural problem has been already studied in detail by Nicot et al. (2016) for this so-called diamond pattern.

The elementary assembly of these four identical grains with the same radius is assumed to be loaded by a system of two perpendicular forces, aligned along the symmetry axis of the diamond pattern, as illustrated in Figs. 1 and 2. For symmetrical reasons, this problem is reduced to a two-degree-of-freedom system as highlighted by Fig. 2. We consider the biaxial compression of a two-degree-of-freedom granular system, with normal and shear interactions between the two grains. Each grain has a radius denoted by $R$. The two displacements ($\Delta_1, \Delta_2$) are associated with the external forces ($F_1, F_2$) in the same direction. It is possible to express

<figure><img src="image_2.png" /><figcaption>**Fig. 1.** A granular symmetrical cyclic chain composed of four identical rigid grains: definition of kinematic parameters.</figcaption></figure>
<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Interaction between the symmetrical four-grain system reduced to a two-grain problem: initial configuration (each grain is constrained in translation).</figcaption></figure>

this two-degree-of-freedom system using some alternative kinematic variables, for instance the distance $l$ between the center of each grain and the orientation $\alpha$ of the axis defined by the line intersecting the center of each grain with respect to the vertical axis. We assume that in the initial configuration, the two grains are in contact with no overlap, and the initial orientation angle is denoted by $\alpha_0$.

The two equivalent kinematic variables $(\alpha, l)$ are expressed as a function of the global displacements $\Delta_1$ and $\Delta_2$, and the geometrical characteristics $R$ and $\alpha_0$

$$ \tan \alpha = \frac{2R \sin \alpha_0 - \Delta_2}{2R \cos \alpha_0 - \Delta_1} \quad \text{and} \quad l = \sqrt{(2R \sin \alpha_0 - \Delta_2)^2 + (2R \cos \alpha_0 - \Delta_1)^2} \quad (1) $$

In the particular case for which grains are assumed rigid without normal interaction, the problem is reduced to a single-degree-of-freedom system, with the following kinematic constraint

$$ l = \sqrt{(2R \sin \alpha_0 - \Delta_2)^2 + (2R \cos \alpha_0 - \Delta_1)^2} = 2R \quad (2) $$

But, in the more general case considered in this paper, assuming an extension or a penetration condition between the two grains, leads to a true two-degree-of-freedom system. The following kinematic relationship may be considered as well, for relating each set of kinematic parameters:

$$ \cos \alpha = \frac{2R \cos \alpha_0 - \Delta_1}{l} \quad \text{and} \quad \sin \alpha = \frac{2R \sin \alpha_0 - \Delta_2}{l} \quad (3) $$

The equilibrium of this granular system can be directly written applying Newton's second law of motion, relating the external forces $(F_1, F_2)$ to the internal forces $(N, T)$ at contact

$$ \begin{aligned} F_1 &= 2N \cos \alpha + 2T \sin \alpha \\ F_2 &= 2N \sin \alpha - 2T \cos \alpha \end{aligned} \quad (4) $$

$N$ = normal force, which acts between the two grains (directed toward the line connecting the center of each grain); and $T$ = tangential force, which is orthogonal by definition, to the normal force. Note that the forces $N$ and $T$ as displayed in Fig. 1 correspond to the forces applied by the left grain on the top grain. $F_1$ is the vertical force and $F_2$ is the horizontal force. These equilibrium equations are expressed in a geometrically exact framework, which accounts for the geometrical change of the granular

evolution problem. Note that all four contacts in the diamond cell have the same $(N, T)$ for symmetry reasons. These expressions of direct equilibrium forces may be equivalently obtained from the principle of virtual work

$$ \begin{aligned} \delta W_{ext} = \delta W_{int} \quad \text{with} \quad \delta W_{ext} &= \frac{F_1}{2} \delta \Delta_1 + \frac{F_2}{2} \delta \Delta_2 \\ \text{and} \quad \delta W_{int} &= -N \delta u_n + T \delta u_t \end{aligned} \quad (5) $$

where the variations of normal displacement and tangential displacements are related to the kinematic variables through the differential relation

$$ \delta u_n = \delta l \quad \text{et} \quad \delta u_t = l \delta \alpha \quad (6) $$

We note that it is not possible to choose $(u_n, u_t)$ as the state variables of this system if $(l, \alpha)$ are kinematic variables describing the state of the system. In particular, we have

$$ \begin{aligned} \delta u_t &= A(\alpha, l) \delta l + B(\alpha, l) \delta \alpha \quad \text{with} \quad A(\alpha, l) = 0 \quad \text{and} \quad B(\alpha, l) = l \\ \text{It is clear that} \end{aligned} \quad (7) $$

$$ \frac{\partial A(\alpha, l)}{\partial \alpha} \neq \frac{\partial B(\alpha, l)}{\partial l} \quad (8) $$

And then $\dot{u}_t = l\dot{\alpha}$ is not a total differential of $(l, \alpha)$. As a conclusion, $(u_n, u_t)$ cannot be considered as the state variables of this two-degree-of-freedom system.

That being underlined, the principle of virtual work is finally expressed as

$$ -N \delta l + T l \delta \alpha = \frac{F_1}{2} \delta \Delta_1 + \frac{F_2}{2} \delta \Delta_2 \quad (9) $$

From the definition of the kinematic variables in Eq. (1), one calculates the rate form of the kinematic variables

$$ \begin{aligned} \dot{l} &= \frac{\partial l}{\partial \Delta_1} \dot{\Delta}_1 + \frac{\partial l}{\partial \Delta_2} \dot{\Delta}_2 = \frac{\Delta_1 - 2R \cos \alpha_0}{l} \dot{\Delta}_1 + \frac{\Delta_2 - 2R \sin \alpha_0}{l} \dot{\Delta}_2 \\ &= -\cos \alpha \dot{\Delta}_1 - \sin \alpha \dot{\Delta}_2 \end{aligned} \quad (10) $$

The rate of the orientation variable is calculated as well as

$$ \dot{\alpha} = \frac{\partial \alpha}{\partial \Delta_1} \dot{\Delta}_1 + \frac{\partial \alpha}{\partial \Delta_2} \dot{\Delta}_2 \Rightarrow l\dot{\alpha} = \sin \alpha \dot{\Delta}_1 - \cos \alpha \dot{\Delta}_2 \quad (11) $$

Injecting both Eqs. (10) and (11) in the weak formulation of equilibrium gives

$$ \begin{aligned} N(\cos \alpha \delta \Delta_1 + \sin \alpha \delta \Delta_2) + T(\sin \alpha \delta \Delta_1 - \cos \alpha \delta \Delta_2) \\ = \frac{F_1}{2} \delta \Delta_1 + \frac{F_2}{2} \delta \Delta_2 \end{aligned} \quad (12) $$

which of course is equivalent to the strong form of the equilibrium equations Eq. (4).

In this paper, we investigate three classes of granular interaction laws. The first class, which is widely used in the literature, is a model introduced in a rate form, initially by the pioneer researchers in the computation of granular material through DEM applied to granular material (Serrano and Rodriguez-Ortiz 1973; Cundall and Strack 1979). The simplified model based on the following incremental relationship between the internal forces (normal force and tangential force) and the kinematic variables:
$$ \begin{aligned} \dot{N}(\alpha, l) &= -k_n \dot{l} \\ \dot{T}(\alpha, l) &= k_t l \dot{\alpha} \end{aligned} \qquad (13) $$

$k_n$ and $k_t$ = normal and tangential stiffnesses, which characterizes this coupled granular interaction. Serrano and Rodriguez-Ortiz (1973) also mentioned that the stiffness parameters may be dependent on the internal force variables, which make the model even more complex. In the following, we will assume that the stiffness parameters $k_n$ and $k_t$ are constant. Eq. (13) expresses that the rate of normal force is proportional to the normal component of the relative velocities of each grain initially in contact, whereas the rate of tangential force is assumed to be proportional to the tangential component of this relative velocity.

It is worth mentioning that Eq. (13) can be partially integrated, for constant values of the stiffness parameters $k_n$ and $k_t$, so that

$$ \begin{aligned} N(\alpha, l) &= -k_n(l - l_0) \\ \dot{T}(\alpha, l) &= k_t \left( -\frac{N}{k_n} + l_0 \right) \dot{\alpha} \end{aligned} \qquad (14) $$

Eq. (14) can be reformulated in a classical rate form

$$ \begin{cases} \dot{N} = K_n(N, T) \dot{l} \\ \dot{T} = K_T(N, T) \dot{\alpha} \end{cases} \quad \text{with} \quad \begin{cases} K_n(N, T) = -k_n \\ K_T(N, T) = k_t \left( -\frac{N}{k_n} + l_0 \right) \end{cases} \qquad (15) $$

This definition of hypoelastic interaction law is in agreement with the one of Truesdell in the restricted assumption of small strain (Truesdell 1955; Ericksen 1958; Bernstein 1960). Usual DEM models based on such constitutive laws are then undoubtedly of a hypoelastic nature because $l\dot{\alpha}$ is not a total differential of $(l, \alpha)$. In the following, it is referred to as strongly hypoelastic. It is also possible to propose some corrections to this contact law to build other classes of granular interaction laws, such as hyperelastic granular interaction laws.

Starting back from Eqs. (5) and (6), the variation of internal work is presented in the following form:

$$ \delta W_{\text{int}} = -N\delta l + Tl\delta\alpha \qquad (16) $$

So that the hyperelasticity condition for the elastic granular interaction is

$$ -\frac{\partial N(\alpha, l)}{\partial \alpha} = \frac{\partial [lT(\alpha, l)]}{\partial l} \qquad (17) $$

This condition is fulfilled with the model of McNamara et al. (2008)

$$ \begin{cases} N(\alpha, l) = -k_n(l - l_0) \\ lT(\alpha, l) = 4R^2 k_t (\alpha - \alpha_0) \end{cases} \quad \text{with} \quad l_0 = 2R \qquad (18) $$

It is possible for the model of McNamara et al. (2008) to build explicitly the internal energy

$$ W_{\text{int}} = \frac{1}{2} k_n (l - l_0)^2 + \frac{1}{2} k_t 4R^2 (\alpha - \alpha_0)^2 \qquad (19) $$

It is worth mentioning that Eq. (18) can be also written with the notation of McNamara et al. (2008) as

$$ \left\{ \begin{array}{l} N(\alpha, l) = k_n D_n \\ T(\alpha, l) = \beta k_t D_t \end{array} \right. \quad \text{with} \quad \left\{ \begin{array}{l} D_n = 2R - l \\ D_t = 2R(\alpha - \alpha_0) \end{array} \right. \quad \text{and} \quad \beta = \frac{2R}{l} \qquad (20) $$

In McNamara et al. (2008), the factor $\beta = 2R/l$ is denoted by $\alpha$ [see Eq. (30) of McNamara et al. 2008]. It is worth mentioning that the hyperelastic granular interaction law can be also presented in a rate form as

$$ \begin{aligned} \dot{N}(\alpha, l) &= -k_n \dot{l} \\ \dot{l}T(\alpha, l) + l\dot{T}(\alpha, l) &= 4R^2 k_t \dot{\alpha} \end{aligned} \qquad (21) $$

or equivalently

$$ \begin{aligned} \dot{N}(\alpha, l) &= -k_n \dot{l} \\ \dot{T}(\alpha, l) &= 4R^2 k_t \frac{\dot{\alpha}l + (\alpha_0 - \alpha)\dot{l}}{l^2} \end{aligned} \qquad (22) $$

The model of McNamara et al. (2008) coincides with the model of Cundall and Strack (1979) for uncoupled granular interaction motion, for a pure tangential motion

$$ k_n \rightarrow \infty \Rightarrow \begin{cases} l = 2R \\ \dot{T}(\alpha) = k_t 2R \dot{\alpha} \end{cases} \qquad (23) $$

or for a pure normal motion

$$ k_t \rightarrow \infty \Rightarrow \begin{cases} \alpha = \alpha_0 \\ \dot{N}(l) = -k_n \dot{l} \end{cases} \qquad (24) $$

Some other possibilities can be built. In the specific case where $N(\alpha, l)$ only depends on $l$, we have the hyperelastic constraint on the shear force constitutive law

$$ \frac{\partial N(\alpha, l)}{\partial \alpha} = 0 \Rightarrow T(\alpha, l) + l \frac{\partial T(\alpha, l)}{\partial l} = 0 \qquad (25) $$

In this last case, the shear force $T(\alpha, l)$ cannot depend solely on the relative rotations of each grain. For instance, it is tempting to replace the intergranular distance $l$ by the initial intergranular distance $l_0 = 2R$ in the tangential rate equation of the initial DEM model, so that $\dot{u}_t = 2R\dot{\alpha}$ and the tangential force may be integrated

$$ \begin{aligned} \dot{N}(\alpha, l) &= -k_n \dot{l} \\ \dot{T}(\alpha, l) &= k_t 2R \dot{\alpha} \end{aligned} \qquad (26) $$

which may be equivalently written as

$$ \begin{aligned} N(\alpha, l) &= -k_n(l - l_0) \\ T(\alpha, l) &= -2Rk_t(\alpha_0 - \alpha) \end{aligned} \qquad (27) $$

However, such a granular interaction law Eq. (27) which violates the hyperelastic criterion Eq. (25) is not hyperelastic. In other words, Eq. (27) is not a Green elasticity granular interaction law (or equivalently not a hyperelastic constitutive granular law). It is a Cauchy elasticity granular interaction law, which cannot be derived from a potential (also classified as a nonconservative elastic law). In this sense, it is referred to as weakly hypoelastic. Even if the contact forces can be expressed in a total form, energy dissipation can occur in closed cycles of loading.

An alternative hyperelastic granular interaction law is the translational model of Turco (2018) or Turco et al. (2019; see also Challamel 2015 for the introduction of this translational granular model), which can be presented for the granular problem studied in this paper as
$$
\begin{align}
N(\alpha, l) &= -k_n [l - 2R \cos(\alpha - \alpha_0)] \nonumber \\
lT(\alpha, l) &= k_n 2R \sin(\alpha - \alpha_0) [l - 2R \cos(\alpha - \alpha_0)] \nonumber \\
&\quad + k_t 2R \sin(\alpha - \alpha_0) 2R \cos(\alpha - \alpha_0) \tag{28}
\end{align}
$$

It can be checked that the hyperelastic criterion Eq. (17) is verified for the hyperelastic model of Turco et al. (2019). Indeed, the model of Turco et al. (2019) is based on the writing of the following internal energy, related to translational normal and shear interactions:

$$
W_{\text{int}} = \frac{1}{2} k_n (l - 2R \cos(\alpha - \alpha_0))^2 + \frac{1}{2} k_t 4R^2 \sin^2(\alpha - \alpha_0) \quad (29)
$$

Nonlinear Evolution for Hypo- and Hyperelastic
Granular Models

Hypoelastic Model

The hypoelastic model is not integrable. The rate form of the equi-
librium equations are obtained from Eq. (4)

$$
\begin{align}
\dot{F}_1 &= 2\dot{N} \cos \alpha + 2\dot{T} \sin \alpha - F_2 \dot{\alpha} \nonumber \\
\dot{F}_2 &= 2\dot{N} \sin \alpha - 2\dot{T} \cos \alpha + F_1 \dot{\alpha} \tag{30}
\end{align}
$$

which can be equivalently written in a rate form, for the hypoelastic
model

$$
\begin{align}
\dot{F}_1 &= -2k_n \cos \alpha \dot{u}_n - 2k_t \sin \alpha \dot{u}_t - F_2 \dot{\alpha} \nonumber \\
\dot{F}_2 &= -2k_n \sin \alpha \dot{u}_n + 2k_t \cos \alpha \dot{u}_t + F_1 \dot{\alpha} \tag{31}
\end{align}
$$

It is then possible to express this rate form of the load–displace-
ment relationship in case of hypoelastic granular interactions

$$
\begin{pmatrix} \dot{F}_1 \\ \dot{F}_2 \end{pmatrix} = 2 \begin{pmatrix} k_n \cos^2 \alpha + k_t \sin^2 \alpha & \cos \alpha \sin \alpha (k_n - k_t) \\ \cos \alpha \sin \alpha (k_n - k_t) & k_n \sin^2 \alpha + k_t \cos^2 \alpha \end{pmatrix} \begin{pmatrix} \dot{\Delta}_1 \\ \dot{\Delta}_2 \end{pmatrix} \\
+ \dot{\alpha} \begin{pmatrix} 0 & -1 \\ 1 & 0 \end{pmatrix} \begin{pmatrix} F_1 \\ F_2 \end{pmatrix} \tag{32}
$$

which can be equivalently presented, using Eq. (11), in rate form
thanks to the following tangent stiffness matrix

$$
\begin{pmatrix} \dot{F}_1 \\ \dot{F}_2 \end{pmatrix} = \begin{pmatrix} K_{11} & K_{12} \\ K_{21} & K_{22} \end{pmatrix} \begin{pmatrix} \dot{\Delta}_1 \\ \dot{\Delta}_2 \end{pmatrix} \quad \text{with} \quad K_{21} \neq K_{12} \quad (33)
$$

The terms of the tangent stiffness matrix of the hypoelastic
model are given by

$$
\begin{align}
K_{11} &= 2k_n \cos^2 \alpha + 2k_t \sin^2 \alpha - \frac{\sin \alpha}{l} F_2 \nonumber \\
K_{12} &= 2(k_n - k_t) \cos \alpha \sin \alpha + \frac{\cos \alpha}{l} F_2 \nonumber \\
K_{21} &= 2(k_n - k_t) \cos \alpha \sin \alpha + \frac{\sin \alpha}{l} F_1 \nonumber \\
K_{22} &= 2k_n \sin^2 \alpha + 2k_t \cos^2 \alpha - \frac{\cos \alpha}{l} F_1 \tag{34}
\end{align}
$$

which is consistent with the results presented by Nicot et al. (2016;
note that the notation of the present paper slightly differs from the
one used by Nicot et al. 2016), as Δ₁ and Δ₂ stand for the displace-
ment of each grain in the direction of application of forces F₁ and
F₂, in the present work.

It is remarkable that the tangent stiffness matrix of this incre-
mental equation is not symmetric due to the hypoelastic nature of
the evolution problem. In general, no analytical solutions are avail-
able for this nonlinear rate equation. The evolution problem can be
numerically integrated with a discrete incremental approach, for
various loading paths. We mainly investigated the response of
the diamond pattern for fixed values of a lateral force and displace-
ment-controlled system in the other direction.

For a one-dimensional loading system, characterized by a fixed
value of F₂, for instance, the rate equation relating the rate of force
to the rate of displacement is obtained from

$$
\dot{F}_2 = 0 \Rightarrow \dot{F}_1 = \frac{K_{11}K_{22} - K_{12}K_{21}}{K_{22}} \dot{\Delta}_1 \quad (35)
$$

This incremental relation is controlled by the sign of the deter-
minant $\det \underline{K}$. A limit load is obtained when

$$
\dot{F}_1 = 0 \Rightarrow \det \underline{K} = K_{11}K_{22} - K_{12}K_{21} = 0 \quad (36)
$$

This is the stability criterion of the associated granular structural
system along the biaxial loading path with a constant lateral force.
It is possible to consider the dimensionless variables

$$
\begin{align}
F_1^* &= \frac{F_1}{k_n 2R}; & F_2^* &= \frac{F_2}{k_n 2R}; & l^* &= \frac{l}{2R}; \nonumber \\
\Delta_1^* &= \frac{\Delta_1}{2R}; & \Delta_2^* &= \frac{\Delta_2}{2R} \quad \text{and} \quad \gamma &= \frac{k_t}{k_n} \tag{37}
\end{align}
$$

The rate equations of the hypoelastic granular model expressed
in dimensionless form are written

$$
\begin{pmatrix} \dot{F}_1^* \\ \dot{F}_2^* \end{pmatrix} = \begin{pmatrix} K_{11}^* & K_{12}^* \\ K_{21}^* & K_{22}^* \end{pmatrix} \begin{pmatrix} \dot{\Delta}_1^* \\ \dot{\Delta}_2^* \end{pmatrix}
$$

with

$$
\begin{align}
K_{11}^* &= 2 \cos^2 \alpha + 2\gamma \sin^2 \alpha - \frac{\sin \alpha}{l^*} F_2^* \nonumber \\
K_{12}^* &= 2(1 - \gamma) \cos \alpha \sin \alpha + \frac{\cos \alpha}{l^*} F_2^* \nonumber \\
K_{21}^* &= 2(1 - \gamma) \cos \alpha \sin \alpha + \frac{\sin \alpha}{l^*} F_1^* \nonumber \\
K_{22}^* &= 2 \sin^2 \alpha + 2\gamma \cos^2 \alpha - \frac{\cos \alpha}{l^*} F_1^* \tag{38}
\end{align}
$$

Hyperelastic Model of McNamara et al. (2008)

Using both Eqs. (22) and (30), it is possible to express the incre-
mental equations for the hyperelastic granular interaction model of
McNamara et al. (2008)

$$
\begin{pmatrix} \dot{F}_1 \\ \dot{F}_2 \end{pmatrix} = \begin{pmatrix} K_{11} & K_{12} \\ K_{21} & K_{22} \end{pmatrix} \begin{pmatrix} \dot{\Delta}_1 \\ \dot{\Delta}_2 \end{pmatrix} \quad \text{with} \quad K_{21} = K_{12} \quad (39)
$$

As this problem is conservative, the tangent stiffness matrix is
now symmetric for this (hyperelastic) structural problem. The terms
of the tangent stiffness matrix are given by
$$K_{11} = 2k_n \cos^2 \alpha + k_t \frac{8R^2 \sin^2 \alpha}{l^2} - \frac{4Rk_t 2R(\alpha_0 - \alpha)}{l^2} \\ \times \cos \alpha \sin \alpha - F_2 \frac{\sin \alpha}{l}$$

$$K_{12} = 2k_n \cos \alpha \sin \alpha - k_t \frac{8R^2 \sin \alpha \cos \alpha}{l^2} - \frac{4Rk_t 2R(\alpha_0 - \alpha)}{l^2} \\ \times \sin^2 \alpha + F_2 \frac{\cos \alpha}{l}$$

$$K_{21} = 2k_n \cos \alpha \sin \alpha - k_t \frac{8R^2 \sin \alpha \cos \alpha}{l^2} + \frac{4Rk_t 2R(\alpha_0 - \alpha)}{l^2} \\ \times \cos^2 \alpha + F_1 \frac{\sin \alpha}{l}$$

$$K_{22} = 2k_n \sin^2 \alpha + k_t \frac{8R^2 \cos^2 \alpha}{l^2} + \frac{4Rk_t 2R(\alpha_0 - \alpha)}{l^2} \\ \times \cos \alpha \sin \alpha - F_1 \frac{\cos \alpha}{l} \quad (40)$$

It is easy to check that

$$K_{12} - K_{21} = - \frac{4Rk_t 2R(\alpha_0 - \alpha)}{l^2} + \frac{F_2 \cos \alpha}{l} - \frac{F_1 \sin \alpha}{l} = 0 \quad (41)$$

due to the equilibrium equation

$$T = \frac{F_1 \sin \alpha - F_2 \cos \alpha}{2} = - \frac{2R}{l} k_t 2R(\alpha_0 - \alpha) \quad (42)$$

Because the model is integrable, the force-displacement relationship can be expressed directly using the following nonlinear equations:

$$\begin{cases} F_1 = 2N \cos \alpha + 2T \sin \alpha \\ F_2 = 2N \sin \alpha - 2T \cos \alpha \end{cases} \quad \text{and} \quad \begin{cases} N(\alpha, l) = -k_n(l - 2R) \\ lT(\alpha, l) = 4R^2 k_t (\alpha - \alpha_0) \end{cases} \quad (43)$$

It is possible to extract $l$ as a function of $\alpha$ as

$$l = 2R + \frac{F_1 \cos \alpha + F_2 \sin \alpha}{-2k_n} \quad (44)$$

Multiplying the first equilibrium equation Eq. (43) by $l$ leads to

$$F_1 l = -2k_n l(l - 2R) \cos \alpha + 2 \sin \alpha 4R^2 k_t (\alpha - \alpha_0) \quad (45)$$

The nonlinear equations are then obtained in a dimensionless form [dimensionless variables defined in Eq. (37)] as

$$\begin{cases} l^* = 1 - \frac{1}{2}(F_1^* \cos \alpha + F_2^* \sin \alpha) \\ F_1^* l^* = -2l^*(l^* - 1) \cos \alpha - 2\gamma(\alpha_0 - \alpha) \sin \alpha \\ \Delta_1^* = \cos \alpha_0 - l^* \cos \alpha \end{cases} \quad \text{and}$$

In the specific case $F_2^* = 0$ (uniaxial compression of the diamond granular chain), Eq. (46) is reduced to the following form

$$-(F_1^*)^2 \frac{\sin \alpha \cos \alpha}{2} + F_1^* \sin \alpha + 2(\alpha_0 - \alpha)\gamma = 0 \quad (47)$$

We finally find the nonlinear load-rotation constitutive law

$$F_1^*(\alpha) = \frac{1 \pm \sqrt{1 + \frac{4\gamma(\alpha_0 - \alpha)}{\tan \alpha}}}{\cos \alpha} \quad (48)$$

We then obtain for the initial load-rotation variable, the branch with the minus sign

$$F_1^*(\alpha) = \frac{1 - \sqrt{1 + \frac{4\gamma(\alpha_0 - \alpha)}{\tan \alpha}}}{\cos \alpha} \quad (49)$$

It is easy to check from Eq. (49) to check that the initial load vanishes for the initial slope, i.e., $F_1^*(\alpha_0) = 0$. Furthermore, the linearized load-rotation relationship can be derived from Eq. (49) using an asymptotic expansion

$$F_1^*(\alpha) = \frac{2\gamma}{\sin \alpha_0} (\alpha - \alpha_0) + \dots \quad (50)$$

It is also possible to express equivalently the evolution problem in rate form. The rate equations of the hyperelastic granular model expressed in dimensionless form are written

$$\begin{pmatrix} \dot{F}_1^* \\ \dot{F}_2^* \end{pmatrix} = \begin{pmatrix} K_{11}^* & K_{12}^* \\ K_{21}^* & K_{22}^* \end{pmatrix} \begin{pmatrix} \dot{\Delta}_1^* \\ \dot{\Delta}_2^* \end{pmatrix}$$

with

$$\begin{align*} K_{11}^* &= 2\cos^2\alpha + 2\gamma \frac{\sin^2\alpha}{(l^*)^2} - \frac{2\gamma(\alpha_0 - \alpha)}{(l^*)^2} \cos\alpha \sin\alpha - F_2^* \frac{\sin\alpha}{l^*} \\ K_{12}^* &= 2\cos\alpha \sin\alpha - 2\gamma \frac{\sin\alpha \cos\alpha}{(l^*)^2} - \frac{2\gamma(\alpha_0 - \alpha)}{(l^*)^2} \sin^2\alpha + F_2^* \frac{\cos\alpha}{l^*} \\ K_{21}^* &= 2\cos\alpha \sin\alpha - 2\gamma \frac{\sin\alpha \cos\alpha}{(l^*)^2} + \frac{2\gamma(\alpha_0 - \alpha)}{(l^*)^2} \cos^2\alpha + F_1^* \frac{\sin\alpha}{l^*} \\ K_{22}^* &= 2\sin^2\alpha + 2\gamma \frac{\cos^2\alpha}{(l^*)^2} + \frac{2\gamma(\alpha_0 - \alpha)}{(l^*)^2} \cos\alpha \sin\alpha - F_1^* \frac{\cos\alpha}{l^*} \end{align*} \quad (51)$$

The uniaxial compression behavior of this cyclic granular chain is first illustrated for $F_2 = 0$, $\gamma = 2$, and $\alpha_0 = \pi/4$. The case of the initial configuration $\alpha_0 = \pi/4$ will be mainly investigated in this paper as a paradigmatic configuration. It is shown in Fig. 3 that the response of the hyperelastic granular chain following Eq. (49) presents a limit point with a maximum rotation $\alpha$ value, whereas the load parameter is still continuing to increase. Fig. 4 represents the same numerical test in the load-displacement ($F_1^*, \Delta_1^*$) space (for the hyperelastic granular chain). Each branch of Eq. (48) is highlighted by the two colors in Figs. 4 and 5. Fig. 4 gives the impression that there exists a bifurcation point, as the two solution components apparently intersect. However, in Fig. 5, a three-dimensional representation of the equilibrium paths is shown, which clearly indicates here that there is no bifurcation, even if multiple solutions may be obtained for a given load parameter. The responses of both models, the hypoelastic and the hyperelastic granular model, are compared in Fig. 6 in the load-displacement space. The two models have the same initial load-displacement
<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Load-rotation curve for the hyperelastic granular model: $\gamma = 2$; $F_2 = 0$; and $\alpha_0 = \pi/4$.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Fig. 4.** Load-displacement curve for the hyperelastic granular model: $\gamma = 2$; $F_2 = 0$; and $\alpha_0 = \pi/4$.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Fig. 5.** Three-dimensional representation of the load-displacement-rotation curve for the hyperelastic granular model: $\gamma = 2$; $F_2 = 0$; and $\alpha_0 = \pi/4$.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Fig. 6.** Comparison of the load-displacement curves for both the hyperelastic and the hypoelastic granular model: $\gamma = 2$; $F_2 = 0$; and $\alpha_0 = \pi/4$.</figcaption></figure>

slope, which confirms that both models coincide when the nonlinear effects of geometrical changes can be neglected. However, the response for large displacements of both models may significantly differ. In Fig. 6, a typical case of $\gamma = 2$ was chosen, which means that the shear stiffness $k_t$ was assumed to be larger than the normal stiffness $k_n$. This is probably unrealistic for most granular interactions and a more realistic value of $\gamma = 0.5$ (where the normal stiffness is larger than the shear stiffness) was investigated in Fig. 7. It is seen that the two hypoelastic and hyperelastic models for this loading configuration give very close results. Fig. 8 represents a parametric study with respect to the shear stiffness ratio $\gamma$. It can be observed that a limit state exists in terms of rotation $\alpha$ for large values of the shear stiffness ratio larger than unity. The asymptotic case
<figure><img src="image_8.png" /><figcaption>Fig. 7. Comparison of the load-displacement curves for both the hyperelastic and the hypoelastic granular model: $\gamma = 0.5$; $F_2 = 0$; and $\alpha_0 = \pi/4$.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Fig. 8. Comparison of the load-rotation curves for the hyperelastic granular model: $\gamma \in \{0.25; 0.5; 0.75; 1; 1.25; 1.5; 1.75; 2\}$; $F_2 = 0$; and $\alpha_0 = \pi/4$.</figcaption></figure>

$\gamma = 0$ corresponds to a pure shear interaction. The existence of two solutions for $F_1^*$ for some values of $\alpha$ corresponds indeed to two levels of grain interpenetration $l^*$. In the general case where $F_2^* \neq 0$, in presence of lateral confinement, considering both equations Eq. (46) leads to the second-order polynomial expression of $F_1^*$, written as

$$ -(F_1^*)^2 \frac{\sin \alpha \cos \alpha}{2} + F_1^* \left[ \sin \alpha - \frac{F_2^*}{2} + F_2^* \cos^2 \alpha \right] - F_2^* \cos \alpha \left( 1 - \frac{F_2^*}{2} \sin \alpha \right) + 2\gamma(\alpha_0 - \alpha) = 0 \quad (52) $$

which generalizes Eq. (47). The load-rotation function can be explicitly found from

$$ F_1^*(\alpha) = \frac{\sin \alpha - \frac{F_2^*}{2} + F_2^* \cos^2 \alpha \pm \sqrt{\delta}}{\sin \alpha \cos \alpha} \quad \text{with} \\ \delta = \left[ \sin \alpha - \frac{F_2^*}{2} + F_2^* \cos^2 \alpha \right]^2 \\ + 2 \cos \alpha \sin \alpha \left[ -F_2^* \cos \alpha \left( 1 - \frac{F_2^*}{2} \sin \alpha \right) + 2\gamma(\alpha_0 - \alpha) \right] \quad (53) $$

<figure><img src="image_10.png" /><figcaption>Fig. 9. Comparison of the load-rotation curves for the hyperelastic granular model with and without lateral load $F_2$: $F_2^* \in \{0; 0.75\}$, $\gamma = 0.5$; and $\alpha_0 = \pi/4$.</figcaption></figure>

As shown in Figs. 9–11, the response of the granular structural model is very sensitive to the presence of the lateral load $F_2$, which may induce some instabilities with the appearance of a limit load, both in the load-rotation diagram (Fig. 9) and the load-displacement diagram (Fig. 10). The loss of stability of the initial path is checked from the change in the sign of the determinant of the stiffness matrix. This phenomenon is more specifically highlighted in a three-dimensional space, where the load is shown to evolve with both the rotation and the displacement in the direction of application of the vertical force (Fig. 11). It is shown that, even in presence of hyperelastic law with quadratic potential with respect to both the intergranular distance and the shear rotation variable, strongly nonlinear geometrical effects may be responsible for the macroscopic nonlinear response including possible instability such as limit point behavior. Structural instabilities may arise even if the granular interaction is linearly elastic or hyperelastic (in the sense that it is associated to a quadratic potential of the kinematic variables).

## Hyperelastic Model of Turco et al. (2019)

The system of two nonlinear equations of $(l, \alpha)$ of this model (hyperelastic model based on normal and shear translational springs see also Challamel (2015) for the introduction of this translational-based interaction model) are summarized as

$$ F_1 = 2N \cos \alpha + 2T \sin \alpha \\ F_2 = 2N \sin \alpha - 2T \cos \alpha $$

and

$$ N(\alpha, l) = -k_n [l - 2R \cos(\alpha - \alpha_0)] \\ lT(\alpha, l) = k_n 2R \sin(\alpha - \alpha_0) [l - 2R \cos(\alpha - \alpha_0)] \\ + k_t 2R \sin(\alpha - \alpha_0) 2R \cos(\alpha - \alpha_0) \quad (54) $$

The nonlinear equations of this model can be equivalently rewritten as:
<figure><img src="image_11.png" /><figcaption>Fig. 10. Comparison of the load-displacement curves for the hypere-lastic granular model in presence of additional lateral load $F_2$: $F_2^* = 0.75$; $\gamma = 0.5$; and $\alpha_0 = \pi/4$.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Fig. 11. Three-dimensional representation of the load-displacement-rotation curve for the hyperelastic granular model: $F_2^* = 0.75$; $\gamma = 0.5$; and $\alpha_0 = \pi/4$.</figcaption></figure>

$$l^* = \cos(\alpha - \alpha_0) - \frac{1}{2}(F_1^* \cos \alpha + F_2^* \sin \alpha) \quad \text{and}$$

$$\frac{F_1^* l^* \sin \alpha - F_2^* l^* \cos \alpha}{2} = \sin(\alpha - \alpha_0)[l^* - \cos(\alpha - \alpha_0)] + \gamma \sin(\alpha - \alpha_0) \cos(\alpha - \alpha_0) \quad (55)$$

In absence of lateral load $F_2^* = 0$, these equations reduce to

$$F_2^* = 0 \Rightarrow -\frac{1}{4}\sin(2\alpha)(F_1^*)^2 + F_1^* \sin(2\alpha - \alpha_0) + \gamma \sin(2\alpha_0 - \alpha) = 0 \quad (56)$$

or equivalently, for the initial nonlinear branch

$$F_2^* = 0 \Rightarrow F_1^*(\alpha) = 2 \left[ \frac{\sin(2\alpha - \alpha_0) \pm \sqrt{\sin^2(2\alpha - \alpha_0) + \gamma \sin(2\alpha) \sin(2\alpha_0 - 2\alpha)}}{\sin(2\alpha)} \right] \quad (57)$$

Imposing the initial constraint $F_1^*(\alpha_0) = 0$, gives the branch with the minus sign

$$F_2^* = 0 \Rightarrow F_1^*(\alpha) = 2 \left[ \frac{\sin(2\alpha - \alpha_0) - \sqrt{\sin^2(2\alpha - \alpha_0) + \gamma \sin(2\alpha) \sin(2\alpha_0 - 2\alpha)}}{\sin(2\alpha)} \right] \quad (58)$$

The differences between the two hyperelastic granular models, namely the model based on shear rotational springs (model of McNamara et al. 2008) and the model based on shear translational springs (model of Turco et al. 2019; see also Challamel 2015) are highlighted on Figs. 12 and 13. As shown from Figs. 12 and 13 (and as confirmed by an asymptotic analysis; Eq. (50) is valid for

<figure><img src="image_13.png" /><figcaption>Fig. 12. Representation of the load-rotation curve for the two hyper-elastic granular models (model of McNamara et al. 2008 based on rotational shear springs, and model of Turco et al. 2019 based on translational shear springs): $F_2^* = 0$; $\gamma = 0.5$; and $\alpha_0 = \pi/4$.</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>Fig. 13. Representation of the load-displacement curve for the two hyperelastic granular models (model of McNamara et al. 2008 based on rotational shear springs, and model of Turco et al. 2019 based on translational shear springs): $F_2^* = 0$; $\gamma = 0.5$; and $\alpha_0 = \pi/4$.</figcaption></figure>

both hyperelastic models), the initial slope of both hyperelastic models also coincide, which means that each model can be distinguished from each other, only for large displacement values of the granular evolution process. While the two granular models behave qualitatively the same with respect to the rotation parameter $\alpha$, the behaviors of both models diverge when expressed as a function of the dimensionless displacement $\Delta_1^*$.

## Unilateral and Blocking Phenomenon

The present granular model may account for the so-called unilateral effect due to the possible asymmetry of the granular interaction in tension and in compression. The generalization of Eq. (19) with the so-called no-tension interaction law could be written from

$$W_{\text{int}} = \frac{1}{2} k_n \langle l_0 - l \rangle^2 + \frac{1}{2} k_t 4R^2 (\alpha - \alpha_0)^2 \quad \text{with} \quad \langle x \rangle = \frac{x + |x|}{2} \qquad (59)$$

where $\langle x \rangle = \text{Macaulay bracket.}$ Moreover, it is possible to account for the appearance of new contact conditions for the granular system. For the granular diamond structure considered in the paper, and in presence of pure shear interaction (which is asymptotically recovered for $\gamma = 0$), the hypoelastic model and the hyperelastic model of McNamara et al. (2008) coincides

$$\gamma = 0 \Rightarrow l = 2R \quad \text{and} \quad T = 2Rk_t(\alpha - \alpha_0) \qquad (60)$$

In this case, there is a specific phenomenon of blocking that appears at the emergence of the contact with the second neighboring grain. In case of pure shear interaction, the coupling of both the equilibrium equation and the elastic shear law gives the load-rotation relationship

<figure><img src="image_15.png" /><figcaption>Fig. 14. Load-rotation response of the cyclic granular chain under pure shear interaction; First and second contact emergence: $F_2^* = 0$; $\alpha_{\text{max}} = \pi/3$ (for the rotational shear spring model); and $\alpha_0 = \pi/4$.</figcaption></figure>

$$\begin{gathered} F_1 \sin \alpha - F_2 \cos \alpha = 2T = 4Rk_t(\alpha - \alpha_0) \quad \text{and} \\ \Delta_1^* = \cos \alpha_0 - \cos \alpha \qquad (61) \end{gathered}$$

which may be expressed in dimensionless form by

$$\begin{gathered} \hat{F}_1 = \frac{2(\alpha - \alpha_0) + \hat{F}_2 \cos \alpha}{\sin \alpha} \quad \text{where} \\ \hat{F}_1 = \frac{F_1}{2Rk_t} \quad \text{and} \quad \hat{F}_2 = \frac{F_2}{2Rk_t} \qquad (62) \end{gathered}$$

Here, the dimensionless force variables have been normalized with respect to the tangential stiffness parameters.

The blocking phenomenon also arises with the translational shear spring model of Turco et al. (2019) in case of pure shear interaction ($k_n \to \infty$)

$$\gamma = 0 \Rightarrow l = 2R \cos(\alpha - \alpha_0) \qquad (63)$$

In this case, the evolution problem in case of pure shear interaction is governed by the following system:

$$\begin{gathered} F_1 \sin(2\alpha - \alpha_0) - F_2 \cos(2\alpha - \alpha_0) = 2Rk_t \sin(2\alpha - 2\alpha_0) \quad \text{and} \\ \Delta_1^* = \cos \alpha_0 - \cos \alpha \cos(\alpha - \alpha_0) \qquad (64) \end{gathered}$$

which may be expressed in dimensionless form by

$$\hat{F}_1 = \frac{\sin(2\alpha - 2\alpha_0) + \hat{F}_2 \cos(2\alpha - \alpha_0)}{\sin(2\alpha - \alpha_0)} \qquad (65)$$

Figs. 14 and 15 illustrates the nonlinear evolution of the compression of the cyclic granular chain with pure shear interaction and in absence of additional lateral load $F_2 = 0$. For the three granular interaction models, the response is constituted of two branches, the first (nonlinear) elastic branch and the emergence of the second granular contact with a displacement-limited constraint. In this pure shear interaction mode, the hypoelastic model coincides with the hyperelastic granular model based on rotational springs. Each branch is analytically given by

$$\begin{gathered} \text{First branch: } \hat{F}_2 = 0 \Rightarrow \hat{F}_1 = \frac{2(\alpha - \alpha_0)}{\sin \alpha} \quad \text{and} \\ \text{Second branch: } \alpha = \alpha_{\text{max}} = \pi/3 \qquad (66) \end{gathered}$$

This rotation limited condition is equivalent to a displacement-limited condition, expressed by
<figure><img src="image_16.png" /><figcaption>**Fig. 15.** Load-displacement response of the cyclic granular chain under pure shear interaction—first and second contact emergence: $F_2^* = 0$; $\Delta_{I,\max}^* = \cos \alpha_0 - \cos \pi/3 = [(\sqrt{2}-1)/2]$; and $\alpha_0 = \pi/4$.</figcaption></figure>

$$ \alpha = \alpha_{\max} = \pi/3 \quad \text{or equivalently} $$

$$ \Delta_I^* = \Delta_{I,\max}^* = \cos \alpha_0 - \cos \pi/3 = \frac{\sqrt{2}-1}{2} \quad \text{for } \alpha_0 = \pi/4 \quad (67) $$

The granular system behaves like a displacement-limited material (blocking phenomenon), as investigated at the material scale by Challamel and Rajagopal (2016).

For the translational shear spring model of Turco et al. (2019), the response is very close with the same displacement-limited condition associated to a close rotation limited value

$$ \cos(\alpha_{\max}) \cos(\alpha_{\max} - \alpha_0) = \frac{1}{2} \quad \text{or equivalently} $$

$$ \Delta_I^* = \Delta_{I,\max}^* = \cos \alpha_0 - \cos \pi/3 = \frac{\sqrt{2}-1}{2} \quad \text{for } \alpha_0 = \pi/4 \quad (68) $$

Furthermore, the (pure shear) response is almost piecewise linear. We show that such macroscopic nonlinearities (and nonsmooth behavior) from a material point of view, may be induced by local geometric nonlinearities at the microscopic scale. Blocking and clogging phenomena are common phenomena observed in the behavior of granular materials. For the coupled system with both shear and normal interactions, the secondary contact is also responsible of a nonsmooth behavior with a piecewise nonlinear response.

## Conclusions

In this paper, we have investigated the nonlinear response of a cyclic granular chain in presence of strongly nonlinear geometrical effects. The general granular interaction laws considered here are classified as hypoelastic and hyperelastic granular interaction laws. Granular structural systems may present different kind of instabilities including symmetrical pitchfork bifurcations with initial stable branches or initial instable branches (as shown numerically by Hunt et al. (2010) in presence of elastic surrounding medium). Challamel and Kocsis (2021) have recently analytically shown that in absence of external elastic confinement, initial stable pitchfork bifurcation may also arise in the buckling and postbuckling of granular columns. In the present paper, we show another instability phenomenon related to a limit load with an exchange of stability. The granular system may present a limit point behavior induced by the geometrical nonlinear effects, as observed by Nicot and Darve (2011) for granular circular chains (Nicot and Darve 2011; Nicot et al. 2017) for a longitudinal granular column.

Nonlinear geometrical effects may be responsible of macroscopic granular instabilities and possible strongly nonlinear behaviors, even if the local granular interaction remains elastic and eventually linear. We also showed, with this elementary granular system, the appearance of some possible displacement-limited phenomena associated to the emergence of secondary granular contacts (and possible kinematics blocking constraints).

## Data Availability Statement

All data, models, and code generated or used during the study appear in the published article.

## References

Bernstein, B. 1960. "Hypoelasticity and elasticity." *Arch. Ration. Mech. Anal.* 6 (429): 89-104.

Challamel, N. 2015. *Letter to F. Nicot, F. Darve and J. Lerbet, on non-conservative and conservative interaction granular laws*, 27 September 2015. Internal Report. Lorient, France: Université Bretagne Sud.

Challamel, N., and I. Elishakoff. 2019. "A brief history of first-order shear-deformable beam and plate models." *Mech. Res. Commun.* 97 (103389): 1-7. https://doi.org/10.1016/j.mechrescom.2019.04.002.

Challamel, N., and A. Kocsis. 2021. "Geometrically exact bifurcation and post-buckling analysis of the granular elastic." *Int. J. Non Linear Mech.* 136 (103772): 1-15.

Challamel, N., J. Lerbet, F. Darve, and F. Nicot. 2020. "Buckling of granular systems with discrete and gradient elasticity Cosserat continua." *Ann. Solid Struct. Mech.* 12 (1): 7-22. https://doi.org/10.1007/s12356-020-00065-5.

Challamel, N., and K. Rajagopal. 2016. "On stress-based piecewise elasticity for limited strain extensibility materials." *Int. J. Non Linear Mech.* 81 (May): 303-309. https://doi.org/10.1016/j.ijnonlinmec.2016.01.017.

Cundall, P. A. 1971. "A computer model for simulating progressive large scale movements in blocky rock systems." In Vol. 1 of *Proc., Symp. of the Int. Society for Rock Mechanics*, 132-150. Lisbon, Portugal: International Society for Rock Mechanics and Rock Engineering.

Cundall, P. A., and O. D. L. Strack. 1979. "A discrete numerical model for granular assemblies." *Géotechnique* 29 (1): 47-65. https://doi.org/10.1680/geot.1979.29.1.47.

Ericksen, J. L. 1958. "Hypoelastic potentials." *Q. J. Mech. Appl. Math.* 11 (1): 67-72. https://doi.org/10.1093/qjimam/11.1.67.

Green, A. E., and P. M. Naghdi. 1971. "On thermodynamics, rate of work and energy." *Arch. Ration. Mech. Anal.* 40 (1): 37-49. https://doi.org/10.1007/BF00281529.

Hunt, G. W., A. Tordesillas, S. C. Green, and J. Shi. 2010. "Force-chain buckling in granular media: A structural mechanics perspective." *Philos. Trans. R. Soc. London, Ser. A* 368 (1910): 249-262. https://doi.org/10.1098/rsta.2009.0180.

Knops, R. J. 1982. "Instability and the ill-posed Cauchy problem in elasticity." In *Mechanics of solids, the Rodney hill 60th anniversary*, edited by H. G. Hopkins and M. J. Sewell, 357-382. New York: Pergamon Press.

Lerbet, J., N. Challamel, F. Nicot, and F. Darve. 2018. "Coordinate free nonlinear incremental discrete mechanics." *Z. Angew. Math. Mech.* 98 (10): 1813-1833. https://doi.org/10.1002/zamm.201700133.

Lerbet, J., N. Challamel, F. Nicot, and F. Darve. 2020. *Stability of discrete non-conservative systems*. Amsterdam, Netherlands: Elsevier.

McNamara, S., R. Garcia-Rojo, and H. J. Hermann. 2008. "Microscopic origin of granular ratcheting." *Phys. Rev. E* 77 (031304): 1-12.

Nicot, F., and F. Darve. 2011. "The H-microdirectional model: Accounting for a mesoscopic scale." *Mech. Mater.* 43 (12): 918-929. https://doi.org/10.1016/j.mechmat.2011.07.006.

Nicot, F., G. Veylon, H. Zhu, J. Lerbet, and F. Darve. 2016. "Mesoscopic scale instability in particulate materials." *J. Eng. Mech.* 142 (8): 04016047. https://doi.org/10.1061/(ASCE)EM.1943-7889.0001100.
Nicot, F., H. Xiong, A. Wautier, J. Lerbet, and F. Darve. 2017. "Force chain collapse as grain column buckling in granular materials." *Granular Matter* 19 (2): 18. https://doi.org/10.1007/s10035-017-0702-0.

Rajagopal, K. R. 2011. "Conspectus of concepts of elasticity." *Math. Mech. Solids* 16 (5): 536-562. https://doi.org/10.1177/1081286510387856.

Serrano, A. A., and J. M. Rodriguez-Ortiz. 1973. "A contribution to the mechanics of heterogeneous granular media." In *Proc., Symp. Plasticity and Soil Mechanics*, edited by A. C. Palmer, 215-227. Cambridge, UK: Cambridge Univ.

Truesdell, C. 1955. "Hypoelasticity." *J. Ration. Mech. Anal.* 4: 83-133.

Turco, E. 2018. "In-plane shear loading of granular membranes modeled as a Lagrangian assembly of rotating elastic particles." *Mech. Res. Commun.* 92 (Sep): 61-66. https://doi.org/10.1016/j.mechrescom.2018.07.007.

Turco, E. 2022. "Forecasting nonlinear vibrations of patches of granular materials by elastic interactions between spheres." *Mech. Res. Commun.* 122 (Jun): 103879. https://doi.org/10.1016/j.mechrescom.2022.103879.

Turco, E., F. dell'Isola, and A. Misra. 2019. "A nonlinear Lagrangian particle model for grains assemblies including grain relative rotations." *Int. J. Anal. Numer. Methods Geomech.* 43 (5): 1051-1079. https://doi.org/10.1002/nag.2915.