Constitutive modelling of non-cohesive soils under high-strain rates: a consistency approach

Recent advancements in computational techniques to model large deformations, such as the arbitrary Lagrangian–Eulerian finite-element method (ALE) (Benson, 1989; DSS, 2012) or the material point method (MPM) (Sulsky et al., 1995; Al-Kafaji, 2013; Yerro Colom, 2015) have made possible the simulation of high-strain-rate (HSR) loading on large-deformation problems in geo-systems (Kim et al., 2015; Tran & Soowski, 2019). However, only a handful of constitutive models are available for granular materials that account for HSR effects (Liingaard et al., 2004; Lu & Fall, 2018).

One reason is that sand is generally hypothesised to be insensitive to rate effects. This is because sand's viscous effects are not significant at low strain rates and relatively large confinement pressures (e.g. more than one metre of overburden) (da Cruz et al., 2005; Hurley & Andrade, 2015). However, there is abundant evidence that for strain rates in the impact regime (e.g. free-fall or projectile impact), HSR effects are significant (Lade et al., 2009; Yamamuro et al., 2011; Omidvar et al., 2012; Suescun-Florez & Iskander, 2017). This is illustrated in Fig. 1, which presents the results of vacuum-dry triaxial tests conducted by Yamamuro et al. (2011) at different strain rates and void ratios. The effects of increasing the loading rate can be summarised as: (a) increased strength (Fig. 1(a)); (b) early manifestation of the peak stress ratio (⁠$η=q/p$⁠, with $p$ = mean stress and $q$ = deviatoric stress) (Fig. 1(a)); and (c) increased dilation and reduced contraction (Fig. 1(b)). These salient features of soil behaviour are caused by several HSR-induced micro-mechanical effects occurring at the particle level (e.g. modified trajectory of particles, inertial effects, particle breakage) (da Cruz et al., 2005; Yamamuro et al., 2011; Andrade et al., 2012; Omidvar et al., 2012; Hurley & Andrade, 2015; Suescun-Florez et al., 2015). Understanding the aforementioned phenomena is of critical importance, for example, for the analysis of impact and blasting problems (e.g. Higgins et al., 2013; Zambrano-Cruzatty & Yerro, 2020). Thus, the development of constitutive equations incorporating the effects of HSR is highly desirable, especially for the simulation of non-cohesive soils (e.g. sands).

Numerous researchers have calibrated viscoplastic constitutive models using the Yamamuro et al. (2011) data set (Higgins et al., 2013; Xu & Zhang, 2015; Mukherjee et al., 2020). The majority of these developments are based on phenomenological equations and are typically constructed using overstress methods such as those proposed by Perzyna (1966) and Duvaut & Lions (1976). In those methods, an explicit calculation of the viscoplastic strain (Perzyna, 1966) or the viscoplastic stress (Duvaut & Lions, 1976) is used to integrate the model. Thus, the approach eliminates the need to satisfy the consistency condition and allows stresses laid outside the yield surface; hence, the term ‘overstress’. For example, Perzyna viscoplasticity requires an overstress function to explicitly calculate the viscoplastic strains (⁠$dεvp$⁠), as shown in equation (1)

with $F$ = yield function, $ΦF$ = overstress function, $∂F/∂σ$ the normal vector to the yield function, and $x=x$ if $x>0$ or 0 otherwise (being $x$ a scalar variable). This approach adds an extra parameter that modifies the stress–strain response of the model. Overstress viscoplasticity can reproduce time-dependent behaviours, such as relaxation and creep (Katona, 1984; Wang et al., 1997; Heeres et al., 2002; Liingaard et al., 2004; An et al., 2011). However, overstress methods are not suitable for stress reversal paths (Heeres et al., 2002), rendering them unsuitable for dynamic problems like earthquake loading.

An alternative method to overstress is to enforce the consistency condition using strain-rate hardening/softening equations; hence, the name ‘consistency approach’. Wang et al. (1997) coined the term and used the consistency approach for a von Mises model with simple work-hardening and strain-rate hardening/softening laws which dynamically modify the yield loci. This dynamic change of the yield surface is the so-called ‘non-stationary’ yield surface, which was generalised by Olszak & Perzyna (1966). Consistent-viscoplasticity models are capable of reproducing relaxation and creep behaviour (Heeres et al., 2002; Liingaard et al., 2004). In addition, due to their mathematical structure, consistent-viscoplasticity models are more efficient when tracking stress reversals than overstress methods (Heeres et al., 2002). However, unlike overstress methods, consistency-viscoplasticity requires extra strain-rate hardening/softening equations to describe the evolution of the model's internal state variables.

In this study, a consistent-viscoplasticity approach is proposed and is implemented on a novel and simple non-associative modified Mohr–Coulomb model (NAMC) with coupled stress–dilation. The document is organised as follows. First, the reference laboratory data set (Yamamuro et al., 2011) and corresponding data-processing methods are presented. Then, the NAMC is established together with the equations governing the state variable of the model (i.e. plastic dilation $Dp$⁠), which are derived using a work–energy balance approach by way of stress dilatancy and inertial effects. Following that, the generalised consistency–viscoplasticity formulation is presented in such a way that it may be applied to any inviscid elastoplastic model. Two procedures for the integration of such a framework on a given constitutive equation are presented: the first one is based on the original ideas of Wang et al. (1997), and the second one is a particularisation of Olszak & Perzyna (1966) for strain rates termed the ‘dashpot method’. Subsequently, the experiments of Yamamuro et al. (2011) are replicated using the proposed viscoplastic NAMC. Furthermore, the viscoplastic NAMC is evaluated under relaxation and impact loading paths. Finally, a discussion and conclusion section, delineating the strengths, limitations and needed research, is provided.

As with previous studies, the data set of Yamuro et al. (YEA2011) (Yamamuro et al., 2011) was utilised in this work to evaluate the behaviour of sands under different strain rates. It consists of 14 triaxial tests performed on precrushed coral sand with specific gravity of 2·68 and void ratios ranging from 1·20 to 0·74. It is classified as poorly graded sand (SP) according to the Unified Soil Classification System (Howard, 1984). Yamamuro et al. (2011) presented results of vacuum-dry drained triaxial tests, which prevented any generation of air pore pressure. Strain rates ranging between 0·0022%/s and 1764%/s were reached with a free-fall weight that pushed an axial piston to shear the samples. Two confinement pressures of 98 kPa and 350 kPa were applied to sand samples having initial void ratios of 0·93 and 1·04.

For the purpose of this work, the data set was manually digitised and a Savitzky–Golay filter (Savitzky & Golay, 1964) was utilised to denoise the filtered data (Klotz & Coop, 2002; Jefferies et al., 2015). The filtering process permitted the rendering of consistent relationships between the deviator stress (⁠$σd$⁠), the volumetric strain (⁠$εv$⁠) and the axial strain (⁠$εa$⁠). A filtered data set displayed in Fig. 2 demonstrates the consistency of the original and filtered data. It is important to note that the YEA2011 data set has some limitations due to the nature of HSR triaxial tests. One concern is that global strain measurements contain significant oscillations due to wave propagation, as described in Abrantes & Yamamuro (2002). For this reason, the strains were measured locally, on a disc of 1·25 cm positioned at the specimen's mid-height, as it appeared to filter out the strain oscillations (Abrantes & Yamamuro, 2002). As will be shown in further sections, few oscillations persisted, particularly for HSR (e.g. test with $ε˙a=1495$%/s in Fig. 1(b)), which makes the dilation calculation (⁠$dεv/dεq$ with $εq$ the deviatoric strain) relatively challenging. Moreover, in local strains, dilation and contraction are enhanced compared to the global ones (Abrantes & Yamamuro, 2002; Klotz & Coop, 2002). Therefore, to partially correct these distortions, the local strains are adjusted according to the approach detailed in Klotz & Coop (2002), with the steps for this adjustment explained in detail in Zambrano-Cruzatty (2021).

This section proposes a non-associative viscoplastic model based on the Mohr–Coulomb formulation. The objective of this exercise is to develop a straightforward and robust model that fundamentally combines friction and dilatation. While the standard Mohr–Coulomb model has disassociated friction and dilation angles, the model given here couples friction and dilation by way of using a stress–dilatancy equation. This gives a valuable and straightforward approach for capturing the non-linear behaviour of granular materials.

The yield surface is defined as the following, which is equivalent to the non-cohesive Mohr–Coulomb failure criterion

where the state of plasticity is attained if $F(σ,ηy)=0$⁠, with $σ$ the stress tensor with invariants $q$ and $p′$⁠, and $ηy$ the mobilised stress ratio (Muir Wood, 2003). The mobilised stress ratio can be calculated using the stress–dilatancy relationship in equation (3), which is based on the seminal work of Nova & Muir Wood (1982).

with $M$ being the critical stress ratio for shearing at constant volume; $N$ is Nova's volumetric coupling coefficient; and $Dp=dεvp/dεqp$ is the plastic dilatancy with $dεvp$ and $dεqp$ the increment of plastic volumetric and deviatoric strain, respectively. In this work, equation (3) is directly connected to the deviatoric strain using a simple hardening/softening equation of the type (Andrade et al., 2012)

where $Dminp$ is the minimum dilatancy and $h$ is the hardening parameter.

Equations (2)–(4) set the basis for an inviscid elasto-plastic model that, up to this point, does not incorporate the effects of the loading rate. To overcome this limitation, the HSR effects on the model's state variable $λvp=(n⋅Dvel⋅dε+L⋅b⋅dε˙)/(n⋅Dvel⋅m+H)$ must be determined. First, the inertial coefficient $I$ (da Cruz et al., 2005) is introduced as a new state variable including the rate effects (i.e. $Xs=f(I)$⁠). da Cruz et al. (2005) showed that the critical friction angle and dilation were connected to this parameter for rheology problems through a series of discrete simulations. Numerous studies have reproduced and validated these observations (Andrade et al., 2012; Hurley & Andrade, 2015). The inertial coefficient is defined by

where $D$ is the particle diameter; $ε˙q$ is the deviatoric strain rate, which replaces the shear strain rate (⁠$γ˙$⁠) for stress path generalisation; and $ρs$ the particle or solid density.

The inertial coefficient takes into account the effect of mobilising the weight of the particle at the enforced strain rate $ε˙q$ normalised by the confinement pressure $p′$⁠. The evolution of the inertial coefficient plotted against the deviatoric strain is depicted in Fig. 3 for the YEA2011 data set, calculated with $D=0⋅32$ mm and $ρs=2860$ kg/m³.

Note that $I$ decreases during the test as the axial strain increases, resulting from an increase in $p′$⁠, until it reaches a constant value during the last part of the test. Smooth changes in $I$ are beneficial in ensuring the stability of the mathematical integration of the model. In the following subsection, the dependency of the model's elastic parameters and the state variable (⁠$ηy(Dp)$⁠) on $I$ is established.

According to the studies of Abrantes & Yamamuro (2002) and Yamamuro et al. (2011), HSR increases the Young's modulus up to 115%. Consistently, the implementation of a non-stationary yield surface requires expressions relating the state variables and the strain rate. These expressions work as strain-rate hardening/softening laws and they can be phenomenological or analytical.

Figure 4 presents the stress–dilatancy plots for the YEA2011 data set grouped by confinement stress and initial void ratio (Figs 4(a)–4(d)). The purpose of these graphs is to provide an overview of the phenomena associated with strain rate and inertial effects. Generally, the stress ratio increases in lockstep with the strain rate in all subplots. Moreover, the paths show a vertical line that shifts to the left, indicating less contraction. In contrast, enhanced dilation is visible only in Fig. 4(b), whereas apparent reduced dilation is observed in the remaining tests in Figs 4(a), 4(c) and 4(d). This is attributed to shear banding preventing samples from reaching minimum dilatancy (e.g. tests with $e0=1⋅03$ in Figs 4(a) and 4(c)), or severe dilation reversal due to oscillations in the volumetric strain response (e.g. tests with $ε˙a=1425$%/s in Fig. 4(d)). Therefore, the minimum dilation for tests in Figs 4(a), 4(c) and 4(d) is not fully reliable. Nonetheless, it is clear that HSR affects stress–dilatancy. As a result, the strain-rate hardening/softening equations can be derived using a work–energy dissipation procedure.

Consider the rate of elastic strain work (⁠$W˙$⁠) in a viscoplastic material using equation (6) (see Appendix 1)

with superscript $(⋅)v$ reserved for viscous quantities. By the superposition principle, equation (6) can be particularised for viscoelastic components and reworked to the form shown in equation (7) (see Appendix 1)

with $Mvel$ the elastic energy–work dissipation density. According to equation (7), the $G/K$ ratio (⁠$∂εvel/∂εqel=G/K$⁠) is independent of the stress ratio, which is consistent with the vertical lines observed at the beginning of the stress–dilation paths corresponding to elastic behaviour (Fig. 4). Increased strain rate results in increased $G/K$ ratio (vertical lines move toward the left in Fig. 4) and prolonged elastic behaviour (vertical lines are higher). Assuming that the initial soil fabric and properties are identical, this is only possible if the yield surface expands as the strain rate increases, which indicates a non-stationary yield surface.

From equation (7), the change in the $G/K$ ratio is expressed as

where the function $fI$ connects $G/K$ and $I$ (see Appendix 1). Integrating equation (8) with limits $G/K=G0/K0→I=I0$⁠, with $I0=$ reference inertial coefficient, results in a power function that relates the dilation and inertial coefficient ratios as described by equations (9) and (10), which satisfy equation (8).

with $κK>κG$ to produce reduction of $G/K$ with increasing strain rate. If $κK=κG$ the model will produce rate-independent Poisson's ratio and does not capture the elastic contraction behaviour depicted in Fig. 4. Note that the inertial coefficient ratio is equivalent to the strain rate ratio (⁠$ε˙q/ε˙ref$⁠) with $ε˙ref$ = the reference strain rate.

The Young's modulus is obtained from the YEA2011 data set. In particular, the Young's modulus values are calculated as the secant modulus at 1% axial strain and converted to shear and bulk modulus using $G=E/21+ν$ and $K=E/31−2ν$⁠, where Poisson's ratio (⁠$ν$⁠) is assumed to be equal to 0·2. It is important to note that a back-calculation of Poisson's ratio using Muir Wood (2003) and Duncan & Chang (1970) yields disparate negative values, probably related to how the strains were measured by Abrantes & Yamamuro (2002). Only for the test with $σ′3=98$ kPa and $e0=0⋅94$ is the Poisson's ratio equal to 0·2, and hence the assumed value. Elastic parameters for quasi-static strain rates serve as a baseline for comparison with those with HSR.

Figure 5(a) shows the evolution of the shear modulus with the inertial coefficient. It is seen that the YEA2011 data points fit the proposed model well with $κG=0⋅06$⁠. Uncertainty for high $I/I0$ ratios is attributed to pressure dependency not being considered and wave oscillations in the YEA2011 data set (as already discussed in the previous section) that could slightly alter the calculation of the secant modulus. In addition, the influence of the inertial coefficient on the $G/K$ ratio is depicted in Fig. 5(b). It is observed that there is a significant scatter, probably due to assuming that all samples have the same Poisson's ratio.

Utilising various Poisson's ratios for each test and adding pressure-dependent elasticity adds additional complexity to the constitutive model. Henceforth, a simple relation between $κK$ and $κG$ is proposed. In Fig. 5(b) it can be observed that $κK=2κG$ yields acceptable results for a subset of points comprising primarily tests with loose sand. However, $κK=3κG$ is more appropriate for dense configuration tests. The relation $κK=2⋅5κG$ is proposed since it returns an average value of the dilatancy ratio and henceforth is used in further simulations.

Following the same strategy as for the elastic parameters, equation (11) is derived from the dissipation of plastic work (see Appendix 1).

with $Mvp$ the viscoplastic energy dissipation rate.

The behaviour observed in Fig. 4 appears to be consistent with equation (11). It is seen that as the strain rate increases, the observed stress ratio increases as all paths move upward. When the critical state is attained $Dp=0$⁠, then equation (11) is transformed into

with the critical condition denoted by the subscript $(⋅)c$⁠. It is deduced from equation (12) that the critical stress ratio is masked by inertial effects, as numerous authors have demonstrated (da Cruz et al., 2005; Jop et al., 2006; Andrade et al., 2012; Hurley & Andrade, 2015). Under relatively large confinement stress, in contrast, inertial effects become less important (as $p↑$⁠, $I<0⋅001$ in equation (5)) and the critical stress ratio remains theoretically unchanged. The latest is the case with the YEA2011 data set, as previously reported by Zambrano-Cruzatty (2021).

Additional evidence for a rate-independent critical stress ratio is presented in Fig. 6(a), which contains the maximum stress ratio (⁠$ηmax$⁠) plotted against the minimum dilatancy (⁠$Dmin$⁠) for all experiments. The colour coding in Fig. 6(a) denotes different strain-rate levels, the white crosses indicate the confinement pressure, and the dashed line with arrows separates data points with different initial void ratios. As depicted in Fig. 6(a), a linear fitting model with $Mtc=1⋅4$ can be constructed for most points (line 1 in Fig. 6(a)), implying that the critical stress ratio is strain-rate independent. This is consistent with the widely accepted assumption that $Mtc$ is commensurate with the mineral-to-mineral friction coefficient. A smaller group of tests, mainly corresponding to $σ′3=350$ kPa, are slightly off trend. Overall, an average $Mtc=1⋅31$ is obtained when all points are considered, which is in good agreement with the critical stress ratio used in Higgins et al. (2013) and Xu & Zhang (2015) for the YEA2011 data set. The outlier points correspond to tests that suffer from artefacts attributed to shear banding that prevent samples from reaching minimum dilatancy or severe dilation reversal due to oscillations in the volumetric strain response, as discussed previously in this section. Therefore, the dilation measurements of these four points have been reasonably considered unreliable.

Theoretical (equation (11)) and experimental evidence (e.g. Fig. 1) indicate that strain-rate hardening must occur as a result of changes in the plastic dilation. This is observed in Fig. 6(a) where the majority of points with the same initial void ratio and confinement pressure shift to the right (i.e. $Dmin$ increases) as the strain rate increases. This indicates that a higher stress ratio is achieved precisely as a result of the enhanced dilation.

By analogy with the $G/K$ ratio, equation (8) serves to model the changes in minimum plastic dilation under HSR.

with $Dmin,0p$ being the inviscid minimum plastic dilation and $κD$ the ‘dilation viscosity’ coefficient. Fig. 6(b) shows a plot of dilation ratio (⁠$Dmin/Dmin,0$⁠) against the inertial coefficient ratio (⁠$I/I0$⁠) to validate equation (13). The points for the tests at which minimum dilation is achieved (Fig. 4(b)) are observed to follow a growing trend (i.e. $Dmin/Dmin,0>1$⁠) as indicated by the regression line with $κD=0⋅04$⁠. Consistently, the dilation ratio for the remaining data sets has significant scatter, since the comparison with the quasi-static dilation is not reliable (i.e. $Dmin/Dmin,0<1$⁠). Moreover, results from other studies (Omidvar et al., 2012; Suescun-Florez & Iskander, 2017) show that the increase in the maximum deviatoric strain, and therefore the dilation, is highly correlated with the initial confinement pressure and void ratio, which are not directly considered in equation (13). In the subsequent section, equation (13) will be revisited to show that its formulation is consistent with Perzyna's viscoplasticity to add a supporting proof for its suitability.

Finally, Table 1 summarises the model's equations and Table 2 (presented in the calibration section below) summarises the value of the parameters determined from the YEA2011 data set.

In this study, viscous effects are incorporated by updating the internal state variables due to changes in the strain-rate tensor (⁠$ε˙$⁠), which is defined as the time derivative of the strain tensor $ε˙=dε/dt$⁠. It is convenient to assume that the stress–strain (⁠$σ$ plotted against $ε$⁠) of soils follows an isotach behaviour. This means that the stress–strain response at different steady strain rates is unique (Fig. 7(a)). A stress path can move from one strain rate and inertial level, say at $ε˙0$ and $I0$⁠, to another strain rate $ε˙f$ and $If$⁠, as illustrated by the arrow in Fig. 7(a). Extrapolating to three dimensions, the increase in viscoplastic stress (⁠$dσvp$⁠) produced by this path is expressed by equation (14).

where $Dvel$ is the viscoelastic constitutive matrix; $dεvel$ is the viscoelastic portion of the strain increment; $dεvp$ is the viscoplastic strain increment, and $dε$ is the total strain increment.

The relationship between strain and strain rate is also required, particularly when substepping is needed for numerical integration. Fig. 7(b) shows an idealised relationship between the strain and the strain rate for a one-dimensional response. A linear relationship can be assumed with sufficient accuracy providing a small strain substepping approach as it is expressed in equation (15).

where $c$ is a constant scalar. In three-dimensional space, equation (15) implies that the increments of strain rate and strain are parallel, which can be assumed with reasonable accuracy if the increment of strain is sufficiently small.

Consider a material model with a yield function $Fσ,Xs$ depending on the stress tensor (⁠$σ$⁠) and a set of internal state variables $Xs=Xs,1,Xs,2,Xs,3,…,Xs,n$⁠, in which the internal state variables can harden/soften as a function of the plastic strain and the strain rate such that $Xs,i=fiεvp,ε˙$⁠. The consistency equation enforces that any set of stress and state variables remain on the yield surface once plasticity is attained (i.e. $dF=0$⁠).

The viscoplastic increment of strain (⁠$dεvp$⁠) is computed by the flow rule, which in the case of the non-associated flow rule is expressed by equation (16).

where $λvp$ is the so-called plastic multiplier and $m=∂P/∂σ$ is the vector normal to the plastic potential function $P$⁠; defined in terms of the stress and internal state variables (⁠$P=f(σ,Xs)$⁠). Recall that Perzyna's viscoplasticity uses equation (1) explicitly to find the plastic strain tensor. If $P=q+Dpp′$ is the plastic potential function, $∂P/∂p=Dp$⁠. For the consistency–viscoplasticity $∂P/∂p=D0p(I/I0)κD$ and henceforth $dε˙∝D0p(I/I0)κD$⁠. This shows that, for the proposed model, the consistency–viscoplasticity is equivalent to Perzyna's with $Φ=(I/I0)κD$⁠. Similar parallelism has been demonstrated by other researchers (Wang et al., 1997; Heeres et al., 2002; Liingaard et al., 2004).

Combining equations (15) and (16) with the consistency condition (⁠$dF=0$⁠), the consistency equation considering strain-rate effects is obtained as (see Appendix 1)

where $n=∂F/∂σ$ is a vector normal to the yield surface, $L=∂F/∂Xs$ is the direction of maximum change of the yield function with respect to the internal state variables, $a=∂Xs/∂εvp$ is the derivative of the internal state variables with respect to the plastic strain, and $b=∂Xs/∂ε˙$ is the derivative of the internal state variables with respect to the strain rate.

At this point, two integration strategies are considered. The first takes into account the original assumption of Wang et al. (1997), in which the elastic portion of the increment of strain (⁠$dεve$⁠) is insensitive to the strain rate (⁠$dε˙vel=0$⁠), hence the increment of strain rate can be expressed using equation (18).

Note that Wang et al. (1997) proposed this approach specifically for the von Mises model, which contrasts with the generalised mathematical implementation presented herein. In this work, this approach is called the ‘original’ consistency viscoplasticy. Alternatively, a variation of this method is proposed, where the increment of strain rate (⁠$dε˙$⁠) is used without further decomposition between plastic and elastic portions. For which reason, no assumptions are needed regarding the elastic strain rate. This variation is termed the ‘dashpot method’. The treatment of the strain rate as an internal parameter is a special case of the non-stationary yield surface theory proposed by Olszak & Perzyna (1966).

The increment of viscoplastic stress is found by substituting equation (18) in equation (17) and subsequently solving for the plastic multiplier (⁠$λvp$⁠). Equation (19) shows the result of this operation

where the viscoplastic hardening/softening term (⁠$Hvp$⁠) is given by equation (20).

where $H$ is a hardening/softening term expressed by

Finally, equations (19), (16) and (14) are combined to obtain equation (22)

which can be simplified as shown in equation (23).

where $Dvp$ is the viscoplastic constitutive matrix which is given by equation (24).

In which $⊗$ represents the tensor product operator.

Equation (25) is retrieved by solving for $λvp$ in equation (17).

Replacing equation (25) in equation (16) and then in equation (14), one can obtain the increment of viscoplastic stress.

Equation (26) can be rearranged and simplified as shown below

where $Dep$ is the inviscid elastoplastic constitutive matrix, and $Cvp$ is a viscoplastic damping matrix described by equation (28)

Equation (27) resembles a three-dimensional Kelvin model, in which a dashpot and a spring are connected in parallel. Hence, the name of this variant of the consistency approach.

The numerical implementation follows the refined explicit integration scheme (REIS) proposed by Sloan et al. (2001). This was developed for inviscid elastoplasticity. The REIS implementation is chosen because it provides a robust framework that is more stable for consistency approaches; however, the consistency method could also be used with other integration schemes. This section describes the numerical integration of constitutive equations. Detailed pseudo-algorithms are presented in Appendix 2.

In the following, the symbol $Δ⋅$ represents the change of a variable between two consecutive steps (⁠$(⋅)i+1−(⋅)i$⁠), the subscript $(⋅)0$ represents initial or inviscid values and the subscript $(⋅)u$ denotes the updated state.

The calculation starts with the initialisation of parameters including the initial stress tensor (⁠$σ0$⁠), state variables (⁠$Xs,0$⁠), strain-rate tensor (⁠$ε˙0$⁠), increment of strain (⁠$Δε$⁠), time increment (⁠$Δt$⁠) and elastic properties included in the elasticity matrix (⁠$D0vel$⁠). Once the model variables are initialised, the current strain rate and strain-rate invariants can be computed.

As mentioned in the previous section, strain-rate hardening/softening controls the yield function size dynamically; thus, negative strain rate increments might cause over-softening of internal state variables. As a result, it is convenient to establish a strain-rate threshold under which the material behaviour is considered to be quasi-static and no additional softening is permitted. Fig. 8(a) depicts the three strain-rate paths for which the internal state variables need to be updated. Paths 1 and 3 cross the threshold and path 2 indicates strain-rate change. Contrarily, Fig. 8(b) shows examples of paths where the internal state variables do not need updating. The algorithm 1 is designed to evaluate the condition above using the deviatoric strain-rate invariant (⁠$ε˙q$⁠) and can be found in Appendix 2.

Subsequently, the elastic-predictor stress (⁠$σvel$⁠) can be calculated with the updated state and elastic parameters (i.e. $Duvel$ and $Xs,u$⁠), and then the yield function can be evaluated on the updated stress (⁠$F(σ0+Δσvel,Xs,u)$⁠) to determine if plasticity is attained within the bounds of a numerical tolerance for the yield function (⁠$F(σ0+Δσvel,Xs,u)>FTOL$⁠, where FTOL denotes yield function tolerance).

After this, the yield function is evaluated on the initial stress (⁠$F(σ0,Xs,0)$⁠) to determine the proportion of elastic to plastic strains (⁠$α$⁠). Three scenarios are possible: (a) the stress state undergoes an elastoplastic transition (⁠$0<α<1$⁠); (b) the stress state experiences pure plasticity (⁠$α=0$⁠); or (c) the stress state transitions from plasticity to elasticity, and then goes back to plasticity. These states can be assessed using the value of the yield function in the initial state and elastic parameters and using either the Pegasus algorithm described in Sloan et al. (2001) or the Newton–Raphson method used here.

To check for unloading in plasticity, the angle between the normal to the yield surface (⁠$n$⁠) and the elastic predictor increment of stress are evaluated (⁠$Δσvel$⁠) as is show in equation (29).

where $||⋅||2$ is the L-two norm of the tensors. If $βF$ is larger than 90° the material experiences elastoplastic unloading, which will require the determination of the elastic unloading proportion (⁠$αΔε$⁠). More detail is presented in the main REIS algorithm (algorithm 2 in Appendix 2).

Calculation of the elastic proportion (⁠$αΔεvel$⁠) is needed when there is elastoplastic loading or elastoplastic unloading (algorithm 2) from an elastic state. An example of the former is illustrated in Fig. 9. When the stress is inside the elastic area, the yield surface (⁠$F$⁠) can shrink or expand dynamically because of the strain-rate softening/hardening processes. A situation that may arise is that the strain rate softens the yield surface while the stress increases. For instance, in Fig. 9(a), the initial yield surface (black dashed line) is dragged to its final position (grey dashed line). The stress path will satisfy $FT=Fσ0+Δσvel,Xs,u>FTOL$ and $F0=Fσ0,Xs,0<FTOL$⁠. Hence, there is a point at which the stress path and the shrinking yield surface intersect, as is illustrated by the continuous grey yield surface. Similarly, this situation also arises if the yield surface expands at a slower pace than the stress rate (Fig. 9(b)). To calculate the proportion $α$ at which this event occurs, a Newton–Raphson method is used.

Assume a positive $α$ such that the strain and strain rates can be decomposed into their plastic and elastic proportions (equations (30) and (31)).

Equation (31) is applicable if the strain and strain rates are linearly related as was established in equation (15).

If the strain-rate hardening/softening equations are bijective, there is a unique set of internal state variables that correspond to a certain level of strain rate, such that one can determine those variables using $Xs,α=fXs,0,ε˙q,0+αΔε˙$⁠, where $f$ represents the set of hardening/softening equations (e.g. equation (13)). Similarly, there exists a stress state $σα=σ0+Δσαvel$ that satisfies $Fσα,Xs,α=0$⁠, where $Δσαvel=αDαvel⋅Δε˙$⁠, and $Dαvel$ is the updated viscoelastic constitutive matrix at $Δεvel$ and $Δε˙vel$⁠. Hence, it is possible to find $α$ by iterating equation (32)

where $Fαi$ is the yield function evaluated at $αi$⁠, and $F′ααi$ is the partial derivative of $F$ with respect to $α$ evaluated at $αi$⁠, and $i$ represents the iteration step.

The derivative of the yield function with respect to $α$ can be computed using equation (33)

Note that the terms $∂Dvel/∂ε˙⋅Δε˙$ and $∂Xs/∂ε˙⋅Δε˙$ are equivalent to their change caused by strain rate $ΔDε˙vel$ and $ΔXs,ε˙$ respectively; which is easier to implement than computing the strain-rate derivatives each time step. The Newton–Raphson procedure is detailed in algorithm 3 in Appendix 2.

Once the proportions of elastic and plastic strains have been computed, the viscoplastic stress can be approximated using an explicit modified Euler algorithm with substepping (Sloan et al., 2001). In principle, integration follows the same algorithm described in Sloan et al. (2001) with minor changes to include the selected viscoplastic framework (i.e. the dashpot or the original consistency–viscoplasticity).

The substepping algorithm (algorithm 4) provides an automatic adequate size of the strain increase. The procedure to subdivide the strain is based on the calculation of the relative error between the two approximations of the increment of stress and state variables (⁠$Δσ1$⁠, $Δσ2$⁠, $ΔXs,1$⁠, $ΔXs,2$⁠) using the modified Euler's algorithm. The step stress and state variables are calculated using equations (34) and (35), respectively.

and the step relative error is calculated using equation (36).

If the relative error is larger than a user-defined tolerance (⁠$STOL$⁠), the step ‘failed’, and the increment of strain and strain rate is subdivided using a scaled pseudo-time (⁠$ΔT$⁠). Sloan et al. (2001) propose to calculate the scaling factor (⁠$qR$⁠) using equation (37).

and subsequently $ΔT$ is updated using $ΔT←qRΔT$⁠.

After the calculation of the updated stress and internal state variables, stress drift can occur (i.e. $Fσu,Xs,u>FTOL$⁠). This drift can be corrected using an algorithm that finds the error in the stress and internal state variables, updating the stress back to the tolerance threshold. The algorithm uses a first-order approximation of the yield function, as shown in equation (38).

The operator $δ⋅$ refers to the error in estimating the stress and state variables. In theory, $δε˙=0$ because the strain rate is imposed. Therefore, the stress drift correction coincides with the procedure proposed in Sloan et al. (2001). Nevertheless, for the original method, the strain-rate error can be computed using $δε˙=δλvpm/ΔtΔT$⁠.

With the considerations described above, the error in the viscoplastic multiplier can be expressed by setting $F=0$ and using the flow rule in equation (16), rendering

where $F0$ is the yield function evaluated at the point where the stresses and state variables are known. Finally, using the flow rule (equation (16)), and the stress–strain relationship (equation (14)) the corrected increments of plastic strain and stress are calculated (⁠$δεp$ and $δσ$⁠). The corrected stresses and state variables are then computed by

where $δXs$ can be computed using equation (53). Equations for the NAMC's derivatives $a$⁠, $b$⁠, $n$⁠, $m$ and $L$ are provided in Appendix 1.

Simulations are performed using the IncrementalDriver programme (Niemunis & Grandas-Tavera, 2017). IncrementalDriver is used for single-element tests, in which stress or strain increments can be prescribed and passed to the constitutive model using the Abaqus user-defined material syntax (UMAT). All simulations are performed with $ΔTmin=1×10−10$⁠, $STOL=1×10−3$ and $FTOL=1×10−9$⁠. The discretisation is achieved by subdividing the maximum strain into 1000 steps to achieve the convergence of the prescribed boundary stresses after 200 iterations per step.

First, the model's parameters $h$ and $N$ are calibrated under quasistatic conditions for each test. Subsequently, the results of the dashpot and original viscoplastic approaches are compared in terms of relative error, maximum uncorrected stress drift and computation cost. A comparison between physical and simulated triaxial tests under steady axial strain rate is presented. Finally, the model's capacity to reproduce relaxation and pulse testing under transient strain rates is investigated.

The hardening modulus $h$ and Nova's coefficient $N$ are calibrated using quasistatic triaxial tests. $h$ and $N$ are chosen by balancing the area below the numerical response and the physical tests using the definition of fitness given in Pal et al. (1996). The averaged $Mtc=1⋅31$ is selected with the generally accepted assumption that $Mtc=1⋅31$ should be unique for a specific material. Fig. 10 shows the comparison between physical tests and their numerical counterparts using the calibrated properties shown in Table 2. The experimental stress–strain response is observed to be well captured in all simulations. However, some volumetric strain responses show more dilation than physical tests. This might be due to the assumption of a unique Poisson's ratio, as mentioned in the viscoplastic NAMC section.

The proposed integration strategies are compared using a triaxial simulation with the same properties as the sample tested at $σ′3=98$ kPa (Table 2), with $κD=κG=0⋅04$ and $κK=2⋅5κG$ values. The test was simulated at $ε˙a=2⋅5×10−3$ %/s and $ε˙a=100$%/s strain rate. As illustrated in Fig. 11(a), both the original and dashpot integration techniques provide the same solution at low strain rates. However, when the rate and axial strain increase, the original and dashpot approaches diverge (Figs 11(a) and 11(c)). This is also evident in Fig. 11(b), which plots the maximum uncorrected stress drift (⁠$F0,max$⁠) against the axial strain. Not surprisingly, the assumption of discarding elastic strain rates unbalances the consistency condition, resulting in the observed drift. A more precise integration technique, such as Runge–Kutta, may be used to resolve the problem mentioned above, although with an increased computational cost.

Similarly, the relative residual error of the modified Euler's procedure (Fig. 11(d)) for the original viscoplasticity approach is higher than the one for the dashpot method. As a result, the dashpot approach is slightly faster than the original (12·1 s to 11·56 s for the simulation shown in Fig. 11). For comparison, the dashpot approach saves around 4·5% of the computation time per integration point.

Apart from the observations made previously, the simulated stiffness and strength also rise when the strain rate increases. Consistently, the material becomes more dative in Fig. 11(c), and the contraction is almost non-existent. This behaviour is further examined in the following section, which employs the model to replicate the YEA2011 data set.

A series of drained triaxial tests are simulated that correspond to the boundary conditions of the tests carried out in Yamamuro et al. (2011), with the properties shown in Table 2. The simulations are conducted only with the dashpot approach. Preliminary simulations using constant viscosity coefficients, $κD=κG=0⋅4$ with $κK=2⋅5κG$⁠, resulted in a poor HSR match between YEA2011 and the simulation, as shown in Fig. 12. It is observed that for tests with $σ′3=98$ kPa, the prediction is closer to the experimental results. In contrast, the maximum deviatoric stress is greatly underestimated for tests with $σ′3=350$ kPa. Similarly, the results are more underestimated when the material is loose than when it is dense. The variability of the viscosity coefficient $κD$ found to fit the YEA2011 data set has been determined and can be summarised in equation (42).

with $Dr$ the relative density and $patm$ the atmospheric pressure. According to equation (42), as the mean stress increases, the viscosity increases. This relationship indicates a stronger relationship between the grain contacts and the viscoplasticity features observed at HSR.

Figure 13 shows deviatoric stress (⁠$q$⁠) plotted against axial strain (⁠$εa$⁠) and volumetric strain (⁠$εv$⁠) plotted against axial strain (⁠$εa$⁠). It is observed that the model captures the increase in the peak deviatoric stress and stiffness as the strain rate increases (Figs 13(a-1)–(d-1)). The model captures an early manifestation of peak deviatoric stress, enhanced dilation and less compression, as observed in the subplots of volumetric strain plotted against axial strain (Figs 13(a-2)–(d-2)). The NAMC performs better for dense conditions than for loose conditions, where the model predicts a better fit to both the prediction of deviatoric stress and the prediction of volumetric strain (e.g. Fig. 13(a)). However, one limitation is that the experimental results show a stiffer response in the plastic regime (⁠$qmax$ occurs early in the strain domain), indicating a need to increase the hardening modulus with the strain rate. In addition, the $q$ plotted against $εa$ curves of simulations with high-strain-rate ratios (i.e. $ε˙a/ε˙ref>105$⁠) seem to coincide. This is because, at high strain rates of similar order of magnitude, equation (11) loses sensitivity (e.g. results of simulations for $ε˙a=897$%/s and $ε˙a=1761$ %/s in Fig. 13(a)).

Relaxation tests are conducted to investigate the numerical behaviour of the model for different viscosities $κD$⁠. Five simulations of drained triaxial test are considered with $σ′3=350$ kPa, the properties in Table 2 and viscosities $κD$ ranging from 0·02 to 0·10. The simulations consist of four loading stages (Fig. 14(a)). First, 8% of axial loading is imposed at 1 %/s, after which pure relaxation occurs for 6 s. Subsequently, the sand is unloaded elastically at 1 %/s for 1 s, and finally relaxation occurs for 4 s.

The results of stopping the loading are observed in Figs 14(b)–14(d). As expected, as the viscosity increases, the maximum deviatoric stress increases. However, during the second stage of the simulation, the deviatoric stresses relax and decrease to their quasistatic value independently of their viscosity (Fig. 14(b)). This relaxation is appreciated in the plot of deviatoric stress against axial strain (Fig. 14(c)), in which a vertical drop occurs at constant axial strain, similar to the one-dimensional viscoplastic model of Smith (1960), which is widely used in driven pile applications. During elastic unloading, the model predicts a lower plateau of deviatoric stress, since the modulus increases with the rate of loading. However, there is little effect on the relaxation of stresses in the elastic range during the last stage of the simulation. Finally, as viscosity increases, volumetric dilation also increases (Fig. 14(d)).

A set of drained triaxial tests with the properties presented in Table 2, with $σ′3=350$ kPa and $e=0⋅94$ is performed using half the period of a cosine function to simulate the pulse loading shown in Fig. 15(a). The deviatoric stress plotted against axial strain (Fig. 15(b)), and the volumetric strain plotted against axial strain (Fig. 15(c)) responses are similar to their counterparts at constant strain rate (Fig. 13(c)), with the inclusion of a decrease in deviatoric stress due to the strain rate softening process in the last part of the test (Fig. 15(b)). This is because in logarithmic space a sharp reduction of strain rates produces a relaxation effect, manifested as a rapid drop in deviatoric stress at the end of the pulse.

Overall, the proposed model captures several features of soil behaviour under HSR that are also reported in models developed with overstress methods (Katona, 1984; An et al., 2011; Higgins et al., 2013; Xu & Zhang, 2015; Mukherjee et al., 2020). A simple NAMC is used with semi-empirical strain-rate hardening/softening equations, which could be used in several applications in geomechanics with minimal calibration effort. In additional, setting $κK=2⋅5κG$ produced acceptable results. Because most geotechnical engineers are familiar with the Mohr–Coulomb model, the implementation of the proposed viscoplastic NAMC in boundary value problems can be expedited. However, the NAMC has limitations owing to its simplicity. For example, NAMC cannot capture plastic deformation under pure volumetric compression, requires individual calibration for different densities, and particle breakage is not considered in the derivation of the model equations.

The YEA2011 data set has limitations due to the difficulties in measuring strains in a highly dynamic triaxial test. Here, it is important to note that volumetric strains were measured locally, which caused distortions in the data that were not feasible to analyse. The reader is directed to Abrantes & Yamamuro (2002) for more information on data acquisition and to Klotz & Coop (2002) for more information on the effect of measuring local strains in defining dilation and the critical state line. Future efforts to extend the experimental data of HSR triaxial testing, whether by physical or discrete-element simulations, are critical in this regard.

A simple power law is selected based on a comparison analysis between the current and some reference values as proposed by equation (8). This functional form, in addition to being simple, captured the trends observed in the data. It provides physical meaning to the viscosity parameter, which can be seen as a damping coefficient. A similar functional form has been selected in other viscoplastic models (Mukherjee et al., 2020). It is also important to highlight that the model lacks sensitivity for strain rates of the same order of magnitude. Further development using the approach in Appendix 1 could be used to construct a more sophisticated model similar to those proposed by Nova & Muir Wood (1982), Jefferies (1993) or Muir Wood (2003).

Finally, the reference strain rate is selected on the basis of the minimum strain rate used in the experiments. More research is required to establish selection criteria. However, a reference strain rate on the order of $ε˙ref=2⋅5×10−5$%/s is suggested to start a calibration process, which is equivalent to a typical loading rate in conventional triaxial tests.

In this study, a simple viscoplastic model based on a non-associated Mohr–Coulomb model is presented. The model is constructed using Nova's stress–dilatancy equation and through an analysis of the effects of strain rates on the work–energy dissipation. Using the approach above, a strain-rate hardening equation is proposed based on a comparison between the current state variables and the reference or quasistatic state. Besides, it is proposed that the inertial coefficient (⁠$I$⁠) is used as strain-rate state variable, since it is proportional to $ε˙q/p$⁠. The model equations are calibrated using impact triaxial data sets published in Yamamuro et al. (2011) with a good fit between the proposed model and the data.

Two integration approaches based on the enforcement of the consistency equation are presented named: (a) the dashpot method (an adaptation of Olszak & Perzyna (1966)) and (b) the original consistency method (Wang et al., 1997). The theoretical framework is developed and expressed in generic terms so that any inviscid model can be adapted following the proposed numerical implementation contained in Appendix 2.

The following conclusions are drawn from the numerical simulations.

• (a)

  Consistency–viscoplasticity offers a reliable alternative to overstress methods. Moreover, since it does not vary significantly from inviscid elasto-plastic models, little effort is required to adapt an inviscid model to form a viscoplastic one.

• (b)

  The consistency approach is shown to capture several features of sand behaviour under HSR. For example, the proposed viscoplastic NAMC reproduces increased stiffness under HSR, enhanced maximum deviatoric stress, improved dilation and reduced compression under constant and transient axial strain rates.

• (c)

  The original and the dashpot viscoplastic variants yielded comparable results. However, the stress drift for the original method is more significant, necessitating more iterations to correct the drift. Computationally, the dashpot method is slightly superior (at least 4·5% more efficient); however, it may render significantly more rapid and stable solutions in boundary value problems.

• (d)

  The consistency–viscoplasticity proposed herein can reproduce relaxation.

• (e)

  The calibration against experimental data shows reasonable comparisons; however, more data are needed for various levels of strain rates, confinement, void ratios, mineralogy and particle morphology.

This material is based in part upon work supported by the National Science Foundation under grant number 1735139. Any opinions, findings and conclusions or recommendations expressed in this material are those of the authors and do not necessarily reflect the views of the National Science Foundation. The authors express their gratitude to the National Science Foundation.

The energy balance equation in a continuum medium will now be considered.

By superposition it is possible to subdivide the inviscid and viscous components of the stress tensor

where superscripts $(.)0$ and $(.)v$ represent inviscid and viscous components, respectively.

The equation above can be reworked using Roscoe's invariants as in

where the inviscid superscript has been suppressed.

By normalising equation (45) by the mean stress and the deviatoric stress one obtains

with $M$ = critical stress ratio, $η=p/q$⁠, $D=ε˙v/ε˙q$ = dilation, and $KvD2+3Gv$ an unknown viscoplastic constitutive relationship.

Recall that the inertial coefficient $I$ scales with $1/p$⁠, and henceforth the unknown terms can be replaced such that, generally

with $fI$ a linear term of $I$⁠. Moreover, the normalised change of $D$ is expressed as

For elastic conditions $Del=G/K$⁠. Integrating equation (48) with limits $G0/K0I→I0$⁠, an equation is obtained of the form

with $κ$ as a viscosity coefficient. From equation (49) is seen that for $G/K$ to change and produce the behaviour observed in Figs 4 and 5, Poisson's ratio must change. That is only possible if equation (49) is decomposed into

with $κ=κG−κK$ with $κK>κG$ to offset the vertical lines in the stress–dilatancy plot to the left.

By analogy, for the plastic range, the following is obtained

The differential of the yield function (⁠$dF$⁠) can be calculated using a first-order approximation shown in equation (52).

where $n=∂F/∂σ$ is a vector normal to the yield surface and $L=∂F/∂Xs$ is the direction of maximum change of the yield function with respect to the internal state variables. Similarly, the increment of internal state variables (⁠$dXs$⁠) can be calculated with a first-order approximation (equation (53)).

where $a=∂Xs/∂εvp$ is the derivative of the internal state variables with respect to the plastic strain and $b=∂Xs/∂ε˙$ is the derivative of the internal state variables with respect to the strain rate.

Combining equations (52) and (53), the following is obtained

Finally, using equations (14) and (16) as replacements in equation (54), the following is obtained

Define $Xs={ηy(Dp)}$⁠. Thus

Then, developing the derivatives of the model equations in Table 1 

with $ϵvp=εvp−1/3εvp1$ the viscoplastic deviatoric strain tensor.