On a hemi-variational formulation for a 2D elasto-plastic-damage strain gradient solid with granular microstructure

1.   Introduction

A large amount of scientific literature deals with non-conservative physical systems, where it is necessary to use methods that are able to handle the related dissipation phenomena ^{[3,16,28,38,54]}. In engineering applications ^[39,40,62] a wide range of materials, like steel and concrete, experience dissipative phenomena such as damage and plasticity. That is why they are especially interesting for the engineering community. But building accurate description of these phenomena can be very difficult, especially for complex material systems like granular ^[39] or lattice-type ^[25,58] microstructure. Continuum damage ^[1,8] and plasticity ^{[4,9,17,22,23,24,32,39,40,49,50]} modeling have been vigorously pursued in the literature, considering the improvements based on phase field models ^{[10,13,14,33,35,55]} for shear bands and fracture. Multi-scale approaches ^[6,27,34,42], in addition to phenomenological ones ^[56,57], have been proposed. Their purpose is to link low-scale descriptions with continuum ^[15,29,31] to include complex emerging behaviors in the continuum. Besides, the strain gradient regularization of the elastic response ^[26,30] should be considered also in this non conservative context. In this paper we recap, in a new and better way, recently developed continuum model for granular materials undergoing damage and plastic deformations. Damage and plasticity are irreversible phenomena ^[2,32,49,50]. In this model, the irreversibility is taken into account by assuming that damage and plastic variables are non-decreasing quantities in time ^[59]. It has been extensively discussed in the literature, that deformed shapes of granular materials may be described in a continuum model by using the relative displacements of the barycenters of the grains, regardless their deformations. Thus, in a coarse-grained description, for a material with a granular micro-structure, the deformation energy as well as the energy dissipated due to damage and plasticity, are expressed in terms of these relative movements. The volumetric energy of deformation, i.e., the deformation energy per unit volume, is assumed to be the sum of deformation energies associated with each intergranular interaction. Each grain-grain interaction is identified by the orientation of the grain-grain direction that, in the continuum limit, are infinite in number. This approach has shown its efficiency in describing granular systems, both at the discrete and at the continuum levels ^{[5,7,18,37,51,60,61]}. We also use a variational approach ^[20,21]. First of all, we define an objective and reversible kinematic vector variable to measure the relative displacements between grains. Thus, we decompose the objective relative displacement between the grains ^[39,48] into two components. The first one is directed along the vector connecting the grain centroids, and it is called the normal component. The second is directed along the orthogonal direction and is called the tangent component. These components are decomposed into elastic and plastic parts. The functionals of elastic deformation and of dissipation energy are defined in terms of these reversible components ^[52,53] and of irreversible damage and plastic variables. Damage is defined by two variables, i.e., the normal and the tangential damage variables, that are both a function of grain-grain orientation. On the one part, plastic displacement does not have to be non-decreasing in time. On the other part, it is characterized as the difference between two non-diminishing plastic variables, that are the accumulated plastic relative displacement in tension and in compression, respectively. The evolution of the form of the body that is obtained without external load defines in this approach the evolution of the plastic strain. The elastic evolution of the damaged material is defined by the total elastic strain energy. This type of energy is explained in terms of the elastic energies related to each intergranular orientation. Therefore, dependencies are obtained for the standard elastic modules (4th order stiffness tensor), of the second gradient elastic modules (6th order stiffness tensor) and of the chiral interaction modules (5th order stiffness tensor). They are functions of parameters describing the damage-elasto-plastic state, describing micro-mechano-morphology, of each orientation. Besides, because of plastic effects, we obtain an expression also of the pre-stress (2nd order tensor) and of the pre-hyperstress (3rd order tensor). The hemi-variational approach ^{[43,44,45,46]} has been used to obtain Karush-Kuhn-Tucker (KKT) kind conditions, that drive the evolution of damage and plastic irreversible kinematic variables. According to the same derivation, also the Euler-Lagrange equations for the progress of the reversible placement function is obtained. It is worth to be noted that, since grain-pairs are oriented in different directions, for a given loading-sequence, they experience different loads and hence different damage and plastic evolution. Thus, the macroscopic response will be complex and with an intrinsically dependence upon path.

The content of this paper is organized as follows. Section 2 provides a rational recap of the model that was introduced in papers ^{[12,41,44,47,59]}. Section 3 provides the scheme of a possible numerical (or analytical) implementation of the model. Section 4 is devoted to the representation of the results for the homogeneous case, where no Finite Element implementation is necessary for the illustration of the results. Section 5 addresses a few concluding comments and future viewpoints.

2.   A rational recap of the model

In this section we will recap the model that has been investigated by the authors more widely in different papers ^{[12,41,44,47,59]}.

2.1.   The elastic energy per unit area

Let $\mathcal{B}$ be the $2D$ reference configuration of a strain gradient elastic body. Its elastic energy $U$ per unit area is assumed to take the following form

where the elastic part $u_{\eta}^{el}$ of the normal displacement $u_{\eta}$ is postulated be equal to the difference of the total normal displacement $u_{\eta}$ and its plastic part $u_{\eta}^{pl}$

where $u_{\tau}$ is the tangential displacement and both are defined as follows,

where the unit vector $\hat{c}$ gives the direction of the considered grain-pair interaction and the domain $\mathcal{S}^{1}$ is the unit circle to which it belongs; the Green-Saint-Venant tensor $G$ and its gradient are tensors of a 2nd and 3rd order,

respectively, where the deformation gradient $F$ and its gradient are defined

in terms of the placement function $\chi\left(\boldsymbol{X}, t\right)$, that is a function of the position $X$ and of time $t$. In Eq (2.3) and (2.4), $L$ is the averaged grain-pair distance. Besides, the damaged tangent stiffness is $k_{\tau, D}$ and the damaged normal stiffness is $k_{\eta, D}$. The damaged normal stiffness is assumed to be asymmetric in tension and compression, i.e.,

where $k_{\eta, D}^{t}$ is the stiffness in tension that is assumed to be smaller than the stiffness in compression $k_{\eta, D}^{c}\gg k_{\eta, D}^{t}$. Besides, the dividing line between tension and compression is given by the sign of the elastic normal displacement $u_{\eta}^{el}$. Thus, the Heaviside function $\varTheta$ is here used. Damage is modeled, as we have already pointed out, with two variables, i.e., the normal damage $D_{\eta}$, and the tangent damage $D_{\tau}$. The damage variables $D_{\eta}$ and $D_{\tau}$ have the role to reduce the damaged normal stiffness $k_{\eta, D}$ (2.7) and the damaged tangent stiffness $k_{\tau, D}$, respectively,

where the non-damaged normal stiffness $k_{\eta}$ and the non-damaged tangent stiffness $k_{\tau}$ have been introduced. Definitions of non-damaged tension ($k_{\eta}^{t}$) and compression ($k_{\eta}^{c}$) normal stiffness through the following expressions

yield the analogous of (2.7) for the non-damaged normal stiffness, i.e.,

We therefore obtain, by insertion of (2.10) into (2.8)$_{1}$, the following expression for the damaged normal stiffness

Insertion of (2.2), (2.3), (2.4) and (2.8) into (2.1) yield the elastic energy per unit area in a more compact form as

where, accounting for the symmetrization induced by the symmetry of the strain tensor $G$, the elastic stiffnesses $\mathbb{C}$, $\mathbb{M}$, $\mathbb{D}$, $\mathbb{P}$ and $\mathbb{Q}$ are identified as follows

According to the legacy of strain gradient continua ^[11,19], a consequence of the expression (2.12) for the elastic energy per unit area is the form of both the stress tensor $S$ and the hyper stress tensor $T$, i.e.,

where $\mathbb{P}$ and $\mathbb{Q}$ take the roles of the pre-stress and the pre-hyper stress, respectively. We oversee that (i) the normal plastic displacement $u_{\eta}^{pl}$ has a direct influence, as expected, from (2.16)–(2.17) on the pre-stress and on the pre-hyper stress and (ii) damage variables $D_{\eta}$ and $D_{\tau}$ has a direct influence from (2.13)–(2.17) on all the stiffness tensors.

2.2.   The dissipation energy per unit area

Damage and plastic variables are dissipative in nature and their evolution are related to the form of the dissipation energy. The dissipation energy per unit area $W$ is the energy dissipated because of irreversible phenomena. An additive decomposition of the dissipation energy is assumed in terms of $W_{D} = W_{D}^{\eta}+W_{D}^{\tau}$, the energy dissipated because of damage phenomena (where $W_{D}^{\eta}$ is the part that is due to the normal phenomena and $W_{D}^{\tau}$ is the part that is due to tangential phenomena), and $W_{pl}$, the energy dissipated because of plasticity phenomena, i.e.,

where $B_{\eta}^{c}$ and $B_{\eta}^{t}$ are two characteristic displacements associated with normal damage dissipation in compression and in tension, respectively. The complicated forms of the assumed dissipated energy in (2.20)–(2.21) are devoted to obtain, in the subsection 2.5, an exponential and/or an arctangential damage evolution, as it will be proved in (2.42) and (2.43).

For cementitious materials it is intuitive that $B_{\eta}^{t}\ll B_{\eta}^{c}$. The reason is that, in tension, a smaller amount of elastic normal displacement is needed to activate damage mechanisms. Besides, $\widetilde{B}_{\tau}\left(u_{\eta}^{el}\right)$ is the characteristic displacement associated with tangent damage dissipation. It is assumed to depend on $u_{\eta}^{el}$, i.e., on the elastic part of the normal displacement as follows,

where $B_{\tau0}$, $\alpha_{1}$ and $\alpha_{2}$ are necessary constitutive parameters needed to express the function $\widetilde{B}_{\tau}\left(u_{\eta}^{el}\right)$. These parameters have the role to couple the two terms, namely the addends (2.20) and (2.21), of the damage dissipation energy $W_{D}$ per unit area. It is worth to be noted that usually, for cementitious materials and in elastic tension, the characteristic tangential displacement $B_{\tau} = B_{\tau0}$ is much lower than the one $B_{\tau} = \alpha_1 B_{\tau0}$ it is necessary in elastic compression ($\alpha_{1}\gg1$). Indeed, a smaller amount of elastic tangential displacement is needed in elastic extension to activate tangential damage mechanisms with respect to the tangential displacement that is needed in elastic compression. In other words, referring to Eq (2.23), this means both that $B_{\tau0} < B_{\tau0}-\alpha_{2}u_{\eta}^{el}$ (which implies $\alpha_{2} > 0$, as $u_{\eta}^{el} < 0$ in compression) and $B_{\tau0} < < \alpha_{1}B_{\tau0}$ (which implies $\alpha_{1}\gg1$).

The plastic dissipation energy function per unit area $W_{pl}$ in (2.22) is assumed to linearly depend on the plastic multipliers $\lambda_{\eta}^{t}$ and $\lambda_{\eta}^{c}$ that are the plastic accumulated displacement in tension and in compression, respectively. The plastic normal displacement is defined as the following difference,

We will show at the end of Subsection 2.5 that the scalars $\sigma_{\eta}^{t}$ and $\sigma_{\eta}^{c}$, defined in (2.22), dictate the yielding conditions (more specifically, they are the characteristic force that is necessary to apply to the grain-pair to activate plastic deformation for no-damage case) of the damage-elasto-plastic grain-pair interaction in tension and compression, respectively. It is worth to be noted that the identification of newly introduced constitutive parameters, i.e. of $L$, $B_{\eta}^t$, $B_{\eta}^c$, $B_{\tau0}$, $\alpha_1$ and $\alpha_2$, is necessary for the application of the present approach for modeling the mechanical behavior or real materials such as, e.g., concrete.

2.3.   The energy functional

The energy functional is defined as the sum of the elastic and dissipation energy,

integrated over the $2D$ reference configuration $\mathcal{B}$. It is a functional of the fundamental kinematical descriptors of the model, i.e. the placement

that is a function of position $\boldsymbol{X}$ and time $t$, and the $4$ irreversible descriptors

that are functions not only of position $X$ and time $t$ but also on the orientation $\hat{c}$.

Damage ($D_{\eta}$ and $D_{\tau}$) and plastic ($\lambda_{\eta}^{t}$ and $\lambda_{\eta}^{c}$) variables are defined both by two variables that are non-decreasing in time. These inequality assumptions,

imply a generalization of standard variational principle into a so-called hemivariational principle.

2.4.   The hemivariational inequality principle

Let us introduce a monotonously increasing time sequence $T_{i}\in\left\{ T_{i}\right\} _{i = 0, \ldots, M}$ with $T_{i}\in\mathbb{R}$ and $M\in\mathbb{N}$ and give initial datum on each of the fundamental kinematic quantities for $i = 0$, i.e., for time $T_{0}$. A family of placements $\chi$ defines the motion for each time $t = T_{0}, T_{1}, \ldots, T_{M}$. The set $AM_{t}$ of kinematically admissible placements is defined for a given time $t$ and the set $AV_{t}$ is defined as the corresponding space of kinematically admissible variations, i.e., $\upsilon = \delta\chi\in AV_{t}$. Admissible variations $\beta$ of the irreversible kinematic quantities $\left(D_{\eta}, D_{\tau}, \lambda_{\eta}^{t}, \lambda_{\eta}^{c}\right)$ must be positive, namely

By definition, the first variation $\delta\mathcal{E}$ of the energy functional (2.25) is calculated as

Besides, the increment of (2.26-2.27), i.e. of the fundamental kinematic quantities, at $t = T_{i}$ is given by the difference between these quantities as evaluated at times $t = T_{i}$ and $t = T_{i-1}$, namely

The same definition is utilized for the increment $\Delta\mathcal{E}$ of the energy functional

Finally, as a matter of facts, the hemi-variational principle is formulated as follows

It is worth to be noted here that introducing the three vectors

where $A$ is intended as the Frechet derivative of the energy functional, the first variation $\delta\mathcal{E}$ of the energy functional in (2.30) and its increment $\Delta\mathcal{E}$ in (2.32) are represented as linear functional of the variation $C$ and the increment $B$ as follows,

so that the hemi-variational principle (2.33) can be formulated

As remarked in ^[36], the inequality (2.33) states that the actual energy release rate is not smaller than any possible one. Thus, it constitutes a kind of principle of maximum energy release rate.

2.5.   The consequence of the hemivariational principle

First of all, the reversibility of the admissible placement variation $\upsilon = \delta\chi\in AV_{t}$ implies

that correspond to standard strain gradient elasticity equations for fixed values of irreversible kinematic quantities $\left(D_{\eta}, D_{\tau}, \lambda_{\eta}^{t}, \lambda_{\eta}^{c}\right)$. Equation (2.37) is derived simply evaluating the inequality (2.36) both for $\delta\chi = \Delta\chi+\hat{\delta}\chi$ and $\beta = \left(\Delta D_{\eta}, \Delta D_{\tau}, \Delta\lambda_{\eta}^{t}, \Delta\lambda_{\eta}^{c}\right)$ and for $\delta\chi = \Delta\chi-\hat{\delta}\chi$ and $\beta = \left(\Delta D_{\eta}, \Delta D_{\tau}, \Delta\lambda_{\eta}^{t}, \Delta\lambda_{\eta}^{c}\right)$, where $\hat{\delta}\chi$ is another arbitrary variation that in (2.37) takes the same symbol $\delta\chi$ just for the sake of simplicity. Secondly, following ^[47], the variational inequality (2.36) implies the following KKT conditions on the 4 irreversible kinematic descriptors $\left(D_{\eta}, D_{\tau}, \lambda_{\eta}^{t}, \lambda_{\eta}^{c}\right)$

where the derivation of (2.38) is done simply evaluating the inequality (2.36) both for $\delta\chi = \Delta\chi$ and $\beta = \left(2\Delta D_{\eta}, \Delta D_{\tau}, \Delta\lambda_{\eta}^{t}, \Delta\lambda_{\eta}^{c}\right)$ and for $\delta\chi = \Delta\chi$ and $\beta = \left(0, \Delta D_{\tau}, \Delta\lambda_{\eta}^{t}, \Delta\lambda_{\eta}^{c}\right)$, the derivation of (2.39) is done simply evaluating the inequality (2.36) both for $\delta\chi = \Delta\chi$ and $\beta = \left(\Delta D_{\eta}, 2\Delta D_{\tau}, \Delta\lambda_{\eta}^{t}, \Delta\lambda_{\eta}^{c}\right)$ and for $\delta\chi = \Delta\chi$ and $\beta = \left(\Delta D_{\eta}, 0, \Delta\lambda_{\eta}^{t}, \Delta\lambda_{\eta}^{c}\right)$, the derivation of (2.40) is done simply evaluating the inequality (2.36) both for $\delta\chi = \Delta\chi$ and $\beta = \left(\Delta D_{\eta}, \Delta D_{\tau}, 2\Delta\lambda_{\eta}^{t}, \Delta\lambda_{\eta}^{c}\right)$ and for $\delta\chi = \Delta\chi$ and $\beta = \left(\Delta D_{\eta}, \Delta D_{\tau}, 0, \Delta\lambda_{\eta}^{c}\right)$ and the derivation of (2.41) is done simply evaluating the inequality (2.36) both for $\delta\chi = \Delta\chi$ and $\beta = \left(\Delta D_{\eta}, \Delta D_{\tau}, \Delta\lambda_{\eta}^{t}, 2\Delta\lambda_{\eta}^{c}\right)$ and for $\delta\chi = \Delta\chi$ and $\beta = \left(\Delta D_{\eta}, \Delta D_{\tau}, \Delta\lambda_{\eta}^{t}, 0\right)$. In (2.38)–(2.41) the auxiliary threshold functions $\tilde{D}_{\eta}(u_{\eta}, \lambda_{\eta}^{t}, \lambda_{\eta}^{c})$, $\tilde{D}_{\tau}(u_{\tau})$, $\tilde{\lambda}_{\eta}^{t}(u_{\eta}, \lambda_{\eta}^{c}, D_{\eta}, D_{\tau})$ and $\tilde{\lambda}_{\eta}^{c}(u_{\eta}, \lambda_{\eta}^{t}, D_{\eta}, D_{\tau})$ have been defined as follows,

From Eqs (2.42) and (2.43) the meaning of $B_{\eta}^{t}$, $B_{\eta}^{c}$ and $B_{\tau}$ as characteristic displacements for the activation of the damage phenomena is evident at least for the loading case. Besides, we observe from (2.44) and (2.45) that with no damage, the meaning of the scalars $\sigma_{\eta}^{t}$ and $\sigma_{\eta}^{c}$ as those characteristic forces that dictate the yielding conditions in tension and compression, is also explained. However, the presence of the normal damage $D_{\eta}$ makes higher such effective characteristic displacement that, in the failure case (with $D_{\eta}\rightarrow1$), becomes infinite.

3.   Implementation of the model

In this Section, the implementation of the model previously presented is divided in the following 5 steps.

$1)$ Null initial, i.e., at time $t = 0$, conditions is assumed on the displacement field for all the points of the body

and and the same for damage and plastic irreversible descriptors both for all the points of the body and for all the directions,

$2)$ Initial isotropy is assumed, that means that non-damaged stiffnesses $k_{\eta}^{c}$, $k_{\eta}^{t}$ and $k_{\tau}$ are assumed to be initially constant $\forall\hat{c}\in\mathcal{S}^{1}$ and $\forall\boldsymbol{X}\in\mathcal{B}$

where $\bar{k}_{\eta}^{c}$, $\bar{k}_{\eta}^{t}$ and $\bar{k}_{\tau}$ are the averaged non-damaged initial stiffnesses.

$3)$ Numerical values of the parameters of the model are assumed and here are reported in Tables 1 and 2.

$4)$ The elastic stiffness tensors ($\mathbb{C}, \mathbb{M}, \mathbb{D}$), as well as the pre-stress and pre hyperstress tensors ($\mathbb{P}, \mathbb{Q}$) are calculated according to Eqs (2.13)–(2.17) at time $t = 0$ with the initial input (3.1)–(3.3). These ingredients with proper boundary conditions are the input for a standard variational principle in (2.37), where the dissipation energy becomes simply an additive constant that does not influence the minimization process. Such a minimization can be performed analytically (as in the homogeneous cases of Section 4) or, more generally, with the use of a Finite Element Method (FEM). Thus, we obtain the displacement field at $i = 1$,

$5)$ With (3.4) we compute the new irreversible descriptors $\left(D_{\eta}, D_{\tau}, \lambda_{\eta}^{t}, \lambda_{\eta}^{c}\right)$ via the KKT conditions (2.38)–(2.41) at $i = 1$,

Thus, we iterate the points 4 and 5 for all those time sequence $T_{i}\in\left\{ T_{i}\right\} _{i = 0, \ldots, M}$ with $T_{i}\in\mathbb{R}$ and $M\in\mathbb{N}$ of the researched time history defined at the beginning of Subsection 2.4.

4.   Homogeneous response

Constitutive parameters are depicted in Table 1. Thus, from ^[12,59] we have an equivalent initial Young modulus in compression $E_{c}$ or in tension $E_{t}$ and Poisson ratio in compression $\nu_{c}$ or in tension $\nu_{t}$, that yields in compression

or in tension

Let us solve the problem in Figure 1. The imposed displacement is

and, consequently, the displacement is trivially deduced in all the body,

Thus, the strain and strain gradient are

and the relative displacements from (2.3) and (2.4) are

where a standard parameterization of the unit vector $\hat{c}$ has been used in terms of an angle $\theta$, i.e.,

The stress response is given in terms of the components of the stress tensors in (2.18),

that implies

4.1.   Strain tension control

In tension

the normal displacement $u_{\eta}\geq0$ is positive in any direction (i.e., $\forall\theta\in\left[0, 2\pi\right]$) and therefore from (2.45) we have

From (2.42) and (2.44) we have

At the beginning of the time history (i.e., with $0 < u_{\eta}\ll L$ from (4.6)) we have from (4.18)$_{2}$ that $\tilde{\lambda}_{\eta}^{t}(u_{\eta}, \lambda_{\eta}^{c}, D_{\eta}, D_{\tau}) < 0$ that means an analytical solution for the accumulation in tension and for the normal damage,

The solution (4.19) is valid only before the threshold condition (2.40) with (2.44)

is satisfied. Thus, by insertion of (4.19)$_{2}$ into (4.20) we have,

The nonlinear algebraic equation (4.21) defines a Lambert function, does not have an analytical solution and can be solved only numerically or graphically, e.g., in Figure 2, from which it is clear, with the constitutive parameters that we have chosen in Table 1, that the condition (4.21) is satisfied for no values of the normal displacement $u_{\eta}$, that implies that (4.19)$_{1}$ is always valid in this investigated case and no plastic behavior takes place. Besides, a different choice of the normal damage characteristic displacement $B_{\eta}^{t}$ would give a different result. From Figure 2 it is in fact also shown that by assigning a value of the normal damage characteristic displacement $B_{\eta}^{t}$ equal to 20 times that assigned in Table 1, there exist a value for the normal displacement $u_{\eta} = \check{u}_{\eta}$ that satisfies (4.21),

that implies that (4.19)$_{1}$ is not valid and plasticity takes a role in this new investigated case. The new analytical solution comes from (4.18), from which we derive the nonlinear algebraic equation for normal damage

and for accumulation of tension

We observe from (4.23) that the normal damage $D_{\eta} = \check{D}_{\eta}$ at which the condition $u_{\eta} = \check{u}_{\eta} > 0$ holds is constant with respect to the normal displacement $u_{\eta}$ and therefore with respect to time. This implies that the accumulation in tension in (4.24) evolves linearly with the normal displacement.

The response is calculated by (4.13)–(4.15) and it is graphically represented in Figure 3. We observe ($S_{11}$ in the left-hand side of Figure 3) that, after a first part where we have an increasing function of time, in a second part of the response the softening induced by damage is evident and a peak reaction is observed as well as a descending curve. A positive reaction is observed also in the orthogonal direction ($S_{22}$ in the right-hand side of Figure 3). The behavior in the orthogonal direction is almost equivalent. However, a second hardening stage is observed because of a non trivial evolution of the equivalent (the response is not anymore isotropic) Poisson ratio.

The softening behavior is due the non homogeneous evolution of the damage variables with respect to the grain-pair orientation $\hat{c}\in\mathcal{S}^{1}$ or, equivalently, to the angle $\theta\in\left[0, 2\pi\right]$, as it is shown in Figure 4.

The anisotropic behavior is also explicated in the evolution of the ratio $\mathbb{C}_{2222}/\mathbb{C}_{1111}$ between the vertical and the horizontal stiffness in Figure 5.

In the left-hand side picture of Figure 4 normal damage evolution is shown and it is evident that the horizontal grain-pair orientation, i.e., around $\theta = k\pi$, with $k\in\mathbb{Z}$, are the most affected by the damage effect. In the right-hand side picture of Figure 4 tangential damage evolution is shown and it is evident that the oblique grain-pair orientation, i.e., around $\theta = \pi/4+k\pi/2$, with $k\in\mathbb{Z}$, are the most affected by the damage effect. It is also observed that the velocity of the damage evolution is constitutively driven by the damage characteristic displacements $B_{\eta}^{t}$, $B_{\eta}^{c}$ and $B_{\tau}$. It is also worth to be noted from Figure 5 that the anisotropy between vertical and horizontal stiffness on the one hand in an elastic isotropic simulation should be maintained at $1$ for the entire history of deformation. Here, on the other hand, it goes from $1$ at $t = 0$ (i.e., initial isotropic behavior) to $60$ at $t = 20s$ (i.e., we have anisotropic behavior induced by deformation).

As it is shown in Figure 2 on the one hand no plastic behavior occurs in this case. On the other hand, increasing $20$ times the normal damage characteristic displacement ($B_{\eta}^{t}\rightarrow20B_{\eta}^{t}$), the response changes dramatically. In Figure 6 the stress response is shown in this second case and the softening behavior is very much attenuated.

In Figure 7 normal and tangent damage evolution are also different. In fact, e.g., maximum normal damage is not any more equal to the admissible value $D_{\eta}\simeq1$ but to the constant value $D_{\eta} = \check{D}_{\eta}\simeq0.4$ that was analytically calculated in (4.23).

The non trivial plastic displacement evolution $u_{\eta}^{pl} = \lambda_{\eta}^{t}-\lambda_{\eta}^{c}$ is therefore shown in Figure 8. It is evident that the plastic displacement occurs only around the horizontal direction.

4.2.   Strain compression control

In compression

the relative displacement $u_{\eta}\leq0$ is negative in any direction (i.e., $\forall\theta\in\left[0, 2\pi\right]$) and therefore from (2.44) we have

From (2.42) and (2.44) we have

At the beginning of the time history (i.e., with $0\geq u_{\eta}\gg-L$ from (4.6)) we have from (4.27)$_{2}$ that $\tilde{\lambda}_{\eta}^{c}(u_{\eta}, \lambda_{\eta}^{c}, D_{\eta}, D_{\tau}) < 0$ that means an analytical solution for the accumulation in compression and for the normal damage as follows,

The solution (4.28) is valid only before the threshold condition (2.41) with (2.45)

is satisfied. Thus, by insertion of (4.28)$_{2}$ into (4.29) we have,

With the constitutive parameters that we have chosen in Table 1 this condition happens at $u_{\eta} = \bar{u}_{\eta} < 0$ and therefore we have (4.28) for $u_{\eta} < \bar{u}_{\eta}$

By insertion of (4.31) into (4.30) we also obtain a nonlinear algebraic equation

the solution of which gives the normal damage $D_{\eta} = \bar{D}_{\eta}$ at which the condition $u_{\eta} = \bar{u}_{\eta} < 0$ holds and that is constant with respect to the normal displacement $u_{\eta}$ and therefore with respect to time. This implies that the accumulation in compression evolves linearly with the normal displacement by insertion of (4.32) into (4.31)$_{2}$,

This response is calculated by (4.13)–(4.15) and it is graphically represented in Figure 9. We observe ($S_{11} < 0$ in the left-hand side of Figure 9) that, after a first part where we have an decreasing function of time, in a second part of the response the softening induced by damage is evident and a peak reaction is observed as well as a slightly increasing curve. A negative reaction is observed also in the orthogonal direction ($S_{22} < 0$ in the right-hand side of Figure 9). The behavior in the orthogonal direction is almost equivalent except for the slightly increasing part again because of a non trivial evolution of the equivalent (the response is not anymore isotropic) Poisson ratio.

The softening behavior is due the non homogeneous evolution of the damage variables with respect to the grain-pair orientation $\hat{c}\in\mathcal{S}^{1}$ or, equivalently, to the angle $\theta\in\left[0, 2\pi\right]$, as it is shown in Figure 10.

The anisotropic behavior is also explicated in the evolution of the ratio $\mathbb{C}_{2222}/\mathbb{C}_{1111}$ between the vertical and the horizontal stiffness in Figure 11.

In the left-hand side picture of Figure 10 normal damage evolution is shown and it is evident that the horizontal grain-pair orientation, i.e., around $\theta = k\pi$, with $k\in\mathbb{Z}$, are the most affected by the damage effect. In the right-hand side picture of Figure 10 tangential damage evolution is shown and it is evident that the oblique grain-pair orientation, i.e., around $\theta = \pi/4+k\pi/2$, with $k\in\mathbb{Z}$, are the most affected by the damage effect. This behavior is nevertheless obscured by the non constant dependence of the damage tangential characteristic displacement $B_{\tau}$ with respect to the normal displacement in compression that was made explicit in (2.23). We have that the higher is the compression the lower is the damage velocity evolution. Thus, the orientation that are more in compression, i.e., $\theta = k\pi$, with $k\in\mathbb{Z}$, have lower damage velocity evolution with respect to the orientation that are less in compression, i.e., $\theta = \pi/2+k\pi$, with $k\in\mathbb{Z}$, where damage have new peaks.

It is also worth to be noted from Figure 11 that the anisotropy between vertical and horizontal stiffness on the one hand in an elastic isotropic simulation should be maintained at $1$ for the entire history of deformation. Here, on the other hand, it goes from $1$ at $t = 0$ (i.e., initial isotropic behavior) to 13 at $t = 20s$ (i.e., anisotropic behavior induced by deformation).

The non trivial plastic displacement evolution $u_{\eta}^{pl} = \lambda_{\eta}^{t}-\lambda_{\eta}^{c}$ is therefore shown in Figure 12. It is evident that the plastic displacement occurs only around the horizontal direction.

5.   Conclusions

This paper recaps and updates a recently developed continuum model for granular materials in order to handle with those important dissipative phenomena as damage and plasticity. The novelty of such a recap is original. The reason is that here for the first time the dissipation energy and the variational inequality are defined directly integrated over the orientation space. This makes the Euler Lagrange equations related to the displacement field (i.e., the elastic Partial Differential Equations PDEs and Boundary Conditions BCs) already integrated over that space. The advantage is that we derive directly the equations that we use for the numerical integration and we do not need to make an artificial integration after their derivation. Its application to an analytical homogeneous case is considered. Plasticity is given by two separate kinematic descriptors (i.e., the accumulation in tension and the accumulation in compression), that are the two plastic multipliers, for every position, time, and grain pair orientation. A hemi-variational principle was adopted to derive the governing equations, from which we obtain Karush-Kuhn-Tucker (KKT)-kind conditions that specify the progression of damage and plasticity relating to each pair of grains interaction. For the case of homogeneous deformation, an analytical solution for the displacement field is assumed and damage and plastic evolution have been derived. It is worth to remark that for non homogeneous deformation, the computation of non homogeneous strain can be reached, e.g., with a Finite Element method according to the scheme developed in ^[47,59], where the presence of strain gradient terms in the PDEs related to the elastic evolution guarantees the overlook of the problem of the mesh-dependence results. Different loading patterns are experienced by different grain-pairs that are oriented in different directions, resulting in complex anisotropic behavior due to damage and plastic evolution. Competition between damage and plasticity dissipative phenomena is demonstrated in these simulations. We show that, for specific parameters, the evolution of plasticity may stop the growth of damage and vice versa. Besides, the presence of newly conceived constitutive parameters that are present in the dissipation energy functional, imposes a fundamental outlook related to their identification. In the presented model we have the inclusion of simple local plastic interactions that contribute to a complex plastic response of the material as a whole. Finally, no additional assumptions, out of the form of the dissipation energy and such as flow rules, are required to describe the plastic behavior and it is worth to point out that the plastic strain is compatible with the existence of a placement function.

Acknowledgments

This project has received funding from the European Union's Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No 899546.

Conflict of interest

All authors declare no conflicts of interest in this paper