A basic hypoplastic constitutive model for sand

Hypoplastic constitutive models are based on nonlinear tensor functions and are characterized by simple formulation and few parameters. In its early stage, mainly basic hypoplastic constitutive equations were concerned, where the stress tensor is assumed as the only state variable. There followed some enhanced models based on the basic constitutive equation by including void ratio as an additional state variable. In this paper, we first show that the widely used hypoplastic model by Wolffersdorff is seriously flawed because the underlying basic equation does not perform properly. We proceed to develop a basic hypoplastic constitutive equation by introducing a new tensorial term, which preserves the critical state at large strain. The model performance is demonstrated by parameter study for some element tests. This simple and robust basic equation is well suited to build more sophisticated models.


Introduction
Recently, there is growing interest in hypoplastic models [2, 3, 6-8, 10-12, 15-17, 31]. Hypoplastic constitutive equations are based on nonlinear tensor functions and are characterized by simple formulation and few parameters. The early hypoplastic models contain four tensor polynomial terms with four material parameters as coefficients, usually with two linear terms and two nonlinear terms in strain rate. A major advantage of the basic model is the fact that it requires only four parameters, which can be easily identified with a single triaxial compression test. The stress tensor is considered as the only state variable in such basic models. As a consequence, the basic models cannot account for the complex history dependence of soil. Moreover, the constitutive model needs to be re-calibrated for the same material but with different initial densities.
The hypoplastic constitutive model with critical state presents a major achievement by introducing the void ratio as an additional state variable into the basic model [24,26]. For a given soil, the hypoplastic model with critical state requires a single set of parameters for the entire range of densities. Since the critical model is built on the basic model, the performance of the critical state model depends on the quality of the basic model.
The first critical state model [24,26] is based on a basic model by Wu [24] and Wu and Bauer [25]. However, this basic model has the shortcoming that the parameters calibrated at critical state in triaxial compression do not necessarily give rise to critical state in triaxial extension. In fact, the basic model shows excessive contractancy in triaxial extension at large strain. An improvement was made by Bauer [1] and Wolffersdorff [20] by requiring that the parameters (coefficients) of the two nonlinear terms are equal with opposite signs. The improved model shows vanishing rate of volumetric strain in triaxial compression and extension at large strain. Note that the number of terms (and of parameters) in the model are reduced from four to three. In this paper, we will show that the model by Wolffersdorff may show pathological behavior if parameters other than those in his paper are considered. We proceed to develop a robust basic model by introducing a new term, which is selected from the representation theorem for tensor polynomials.

Basic model: from Wu to Wolffersdorff
To gain perspectives, let us consider the following framework for basic hypoplastic constitutive equations by Wu and Kolymbas [22]: where L and N are tensorial functions and L is assumed to be linear in _ , r is the Cauchy stress tensor, _ is the strain rate (stretching) tensor, jj _ jj stands for the norm of the strain rate, r is the Jaumann stress rate defined by where _ r is the material time derivative of the Cauchy stress, and _ x is the spin tensor. The strain rate and the spin tensors are the symmetric and anti-symmetric parts of the velocity gradient.
The functions L and N must be isotropic to remain invariant under rigid body rotations (objectivity requirement). Since L is linear in _ , constitutive equations within the framework (1) are necessarily rate independent. This can be easily ascertained by the fact that (1) is homogeneous of the first degree in strain rate _ . To this end, we can write (1) in the following form by making use of Euler's theorem for homogeneous functions where L ¼ oL=o _ is a fourth-order tensor, the symbols and : denote an outer (dyadic) product and inner product between two tensors, respectively. The terms in the brackets of (3) represent the tangential stiffness, which depends on the direction of strain rate. At this stage, it is interesting to have a comparison with elastic and elastoplastic constitutive equations. For an elastic model, the same stiffness tensor is used for loading and unloading. For an elastoplastic model, two different stiffness tensors are used, one for loading and one for unloading. Loading and unloading are differentiated by the loading criteria. In elastic and elastoplastic models, the stiffness tensors may depend on stress but are independent of strain rate. In the hypoplastic model (3), loading and unloading are not explicitly stated but implicitly contained in the dependence on the direction of strain rate.
Based on the basic model (1), the critical state model can be obtained by multiplying the nonlinear term in (1) by a density function f(e) [26] with the following density function f ðeÞ ¼ e crt À e e crt À e min ða À 1Þ þ 1 ð5Þ The density function f(e) represents a linear interpolation between f ðe min Þ ¼ a and f ðe crt Þ ¼ 1, with e crt being the void ratio at critical state and e min the minimum void ratio. Both e crt and e min may depend on stress. Several relationships between void ratio and stress are available, e.g., logarithmic function in Cam Clay model [14], exponential function [26] and potential function [4]. Obviously, the critical state model reduces to the basic model at critical state with f ðe crt Þ ¼ 1. The void ratio is updated during deformation, and the density function serves as feedback to adjust the model response. To this end, the basic model is calibrated at critical state with vanishing rate of stress and void ratio. A specimen with an initial density looser or denser than critical is treated as an initial perturbation from the critical state. For continuing deformation, the model response is adjusted by the density function so that the critical state is approached at large strain. We proceed to consider the following basic constitutive equation proposed by [24], where c i (i = 1, 2, 3, 4) are dimensionless parameters. The deviatoric stress tensor r Ã in the above equation is defined by r Ã ¼ r À 1=3ðtrrÞI with I being the second-order unit tensor. The four parameters can be identified with a single triaxial compression test. The performance of the model is shown in great detail by Wu and Bauer [25]. Note that the basic model (6) is homogeneous of the first degree in stress. As a consequence, the tangential stiffness is proportional to stress with vanishing stiffness at null stress, which is reasonable for cohesionless soils. When the basic model (6) is used in the critical state model (4), the parameters c i (i = 1,2,3,4) are calibrated with a triaxial compression test with an initial density slightly looser than critical. In such a test, the critical state can be attained without apparent shear band formation [18,28]. Usually, the four parameters are determined by considering the stress state, the stress rate and strain rate at two distinguished points of the stress-strain curve in triaxial compression test, usually an isotropic stress state after consolidation and critical state, where the stress rate and the volumetric strain rate vanish simultaneously. The material behavior in an isotropic stress state is characterized by the initial tangential modulus E i and the Poisson ratio m i . The critical state is characterized by the friction angle and dilatancy angle u crt and w crt . Obviously, we have w crt ¼ 0.
It turned out, however, that the basic model, with the four parameters calibrated at critical state in triaxial compression, does not lead to critical state for other stress paths, e.g., triaxial extension. Figure 1 shows the stressstrain and volumetric strain curves for triaxial compression and extension (dotted lines) starting from an isotropic stress state. The dotted curves in Fig. 1 are obtained with the parameters E i ¼ 100r c , m i ¼ 0, u crt ¼ 30 and w crt ¼ 0 . r c stands for the consolidation stress in triaxial tests. These parameters are typical for loose sand except the initial Poisson ratio m i ¼ 0. However, let us leave this aside for a moment. Obviously, critical state is not attained in triaxial extension. The model shows excessive contractancy in extension at large strain. Moreover, the stress ratio at extension is higher than compression. The Mohr-Coulomb failure surface would give rise to equal stress ratios at compression and extension.
It was discovered quite fortuitously that critical state is reached for all stress paths if the parameters of the two nonlinear terms in (6) equal with opposite signs [1], i.e.
This is twofold fortuitous because the discovery was by chance and because the condition (7) is only valid for the basic model (1). To this end, the two nonlinear terms can be merged into one to give Note that we use c 3 in both (6) and (8) for simplicity. In general, c 3 assumes different values in (6) and (8). The performance of (8) in triaxial extension is shown by solid lines in Fig. 1. Apparently, critical state is reached at extension. The critical state model by Wolffersdorff [20] is based on the above basic equation, which is, however, seriously flawed. The detailed formulations of hypoplastic model given by Wolffersdorff are shown in Appendix. The critical state model by Wolffersdorff [20] used a different set of material parameters, which can be related to the three material parameters E i , m i and / crt . Note that we have w crt ¼ 0 at critical state. Moreover, the constitutive equation is homogeneous in stress, which means that the parameters h s in the model by Wolffersdorff can be scaled by E i . We are then left with two parameters m i and / crt . As a consequence, there exists a relationship between these two parameters. Figure 2 shows the relationship between m i and / crt according to the model by Wolffersdorff. A perusal of Fig. 2 shows that the initial Poisson ratio is close to zero for a friction angle / crt of about 30 . However, an initial Poisson ratio of m i ¼ 0 is at odds with most experimental observations. In fact, the Poisson ratio for sand usually lies between 0.2 and 0.3. One might ask what happens if the friction angle remains 30 and an initial Poisson ratio other than zero, say 0.2 is used? In this case, we need to modify another material parameter f d in the model, which depends on the stress state, initial density, critical density, minimal density and three other material parameters in the model (h s , n and a). The relationship between m i and f d is shown in Fig. 3. It can be seen from both figures that the initial Poisson ratio can even be negative for / [ 30 or f d [ 1. If the material parameters in the paper Wolffersdorff [20] are taken, the initial Poisson ratio m i ¼ 0:2 will lead to f d ¼ 0:4, which corresponds to an initial pore ratio of about 0.36, which is not realistic. It seems that the number of terms should not be less than four in order not to compromise the model performance.
We proceed to remedy the deficiency of (8) by adding a new term to the constitutive equation. The new term can be chosen based on the representation theorem for tensor  functions. To this end, let us have a look at the representation theorem in its simplest form, i.e., bilinear functions of r and _ [19].
where a i ; ði ¼ 1. . .5Þ are coefficients. The new term should not compromise the model at critical state, i.e., it should vanish identically for tr _ ¼ 0. Such a term will improve the model performance while leaving the critical state unchanged. As can be easily seen, the above theorem offers two such terms, namely rðtr _ Þ and ðtrrÞðtr _ ÞI. We have studied both terms by numerical simulation of laboratory tests and found both terms well suited to be included in the basic model. In this paper, however, we will focus on the former term. In fact, the representation theorems often offer more than one possible choices [29].
We proceed to remedy equation (8) by including the term rðtr _ Þ to obtain the following basic model: Again, the same notations for the four parameters c i , (i = 1 ... 4) are retained in the above equation without inducing confusion. The above equation has four terms and four parameters with the difference to (6) in that we have three linear terms and one nonlinear term.

Model performance
The performance of the basic model (10) is demonstrated by simulating some triaxial tests with different parameter combinations. The material parameters in (10) can be determined with a single triaxial compression test according to the procedure described by Wu and Bauer [25].

Parameter study with basic model
Several numerical simulations of triaxial compression and extension tests with different dilatancy angles are carried out. The other three material parameters are kept unchanged with E i =r c ¼ 170, m i ¼ 0:2, / crt ¼ 30 .
The stress-strain and volumetric strain curves are shown in Fig. 4. The following observations can be made. The stress-strain curves remain un-effected by the dilatancy angle. The volumetric strain curves show typical initial contractancy followed by dilatancy. The dilatancy is slightly larger in compression than in extension, which agrees well with experimental observations [23]. The next simulation series are carried out by varying the initial Poisson ratio from 0.0 to 0.3. The same initial stiffness and friction angle are used, i.e., E i =r c ¼ 170, / crt ¼ 30 . The numerical results with different Poisson ratios for two dilatancy angles w ¼ 0 (critical state) and w ¼ 10 are shown in Figs. 5 and 6. Again, the stress-strain curves are dictated by the friction angle and initial modulus and remain un-effected by the other parameters. The initial Poisson ratio has large influence on the volumetric strain, both the initial contractancy and the subsequent dilatancy. This offers more flexibility to fit model prediction to experimental data.

Simulation of stress path tests
The model performance along some stress paths other than compression and extension is shown by following the four stress paths in Fig. 7. Stress paths 1 and 4 represent triaxial compression and extension, respectively. In stress paths 2 and 3, both r 1 and r 3 increase with the ratio of 2 : 1 and 1 : 2, respectively. The following material parameters are used: E i =r c ¼ 170, m i ¼ 0:2, / ¼ 30 and w ¼ 10 . The numerical simulations are shown in Fig. 8. Stress path 1 has the highest strength and the lowest dilatancy. Stress path 4 has the lowest strength and the highest dilatancy. This is because the mean pressure in path 1 is higher than path 4. Dilatancy is suppressed by elevating mean pressure. With Eq. (8), we have a simple and robust basic model at hand to develop more sophisticated models. The four parameters can be easily determined with a single triaxial compression test. As will be shown later, the failure surface likes that of Drucker-Prager. We are tempted to compare our basic model with an ideally plastic model with the failure surface of Drucker-Prager, which comes also with four material parameters, namely elastic modulus, Poisson ratio, friction angle and dilatancy angle. However, our model is superior in many respects, e.g., single equation for loading and unloading, nonlinearity before failure and smooth transition from contractancy to dilatancy.

Model with critical state
The simplicity and robustness of the basic model (10) make it an interesting stand-alone model for many applications. However, the true strength of basic model can be best appreciated in connection with the critical state model. The focus of this paper is on the basic model. Therefore, we   (4) and (5) work. The parameters e crit and e min are assumed to be constant, i.e., independent of stress. Let us consider the following typical parameters for sand with a ¼ 0:8, e crit ¼ 0:7 and e min ¼ 0:5. The expression of f(e) can be obtained by inserting these parameters into (5), which can be further inserted into the constitutive equation (4) to give It can be easily seen in the above equation that void ratio e is included a state variable in addition to stress r. Comparison between (11) and (10) shows that the critical state model and the basic model differ only be the simple expression 1:7 À e in the nonlinear term. Yet, this simple expression makes big difference in the model performance.
The four parameters in the basic model (10) are determined at critical state with E i =r c ¼ 170, m i ¼ 0:2, / crt ¼ 30 and w crt ¼ 0 . Equation (11) is used to simulate triaxial compression tests with four different initial densities (see Fig. 9). The densities with e 0 ¼ 0:5 and e 0 ¼ 0:6 are denser than critical, while the density e 0 ¼ 0:8 is looser than critical.
Some salient behavior of sand is well reproduced by the critical state model, e.g., strain softening and dilatancy for dense sand, strain hardening and contractancy for loose sand. At large strain, the shear strength at critical state (characterized by the friction angle / crt ) and the void ratio at critical state (characterized by e crt ) are approached.
It should be pointed out that the critical state model (11) is only for demonstration purpose. Some improvements can be readily made. For example, both e crt and e min in (5) are known to depend on stress.

Stress response envelope
Until now, only numerical simulations of triaxial tests along certain stress paths are considered. The general properties of constitutive model be best appreciated by the so-called response envelope. The response envelope at a Void ratio (e) e=0.5 e=0.6 e=0.7 e=0.8 Fig. 9 Numerical simulation of triaxial compression test with different initial densities (w ¼ 10 ), stress-strain curves, volumetric strain curves and void ratio curves given stress state r 0 , called initial stress, is the surface spanned by all stresses r ¼ r 0 þ Dr. The stress increments Dr are obtained from the probes of strain rates of all possible directions with the same magnitude jj _ jj. The response envelopes can be depicted in the three-dimensional stress space [27]. For visual inspection, however, the Rendulic plane with r 2 ¼ r 3 is often the better choice.
Some response envelops of the basic model (10) are shown in Fig. 10. The following material parameters are used: E i =r c ¼ 170, m i ¼ 0:2, / ¼ 30 and w ¼ 10 . The initial stress states of the five response envelopes are on the same deviatoric plane. The dotted line represents hydrostatic stress states. The upper and lower straight line represents the failure criterion for compression and extension. Any point on the response envelop is associated with a certain strain increment with the ensuing stress increment. Since the strain rates have the same magnitude, the distance from any point on the response envelope to the initial stress state signifies the directional stiffness as shown in (3).
The following observations are shown in Fig. 10. The response envelope in the middle is symmetric with reference to the hydrostatic axis, which is dictated by the requirement of isotropy. The two outmost response envelopes are for the stress states on the failure surface. Apparently, the initial stress on failure surface lies on the response envelop. A closer look at these two response envelops shows that a small part of the response envelopes are outside the failure surface, which brings us to the derivation of the failure surface and bound surface.

Failure surface
Unlike plasticity theory with a priori defined yield and failure surface, the failure surface in hypoplasticity can be derived as outcome of the constitutive equation. The two failure lines in Fig. 10 are the traces of the failure surface on the Rendulic plane. We proceed to derive the failure surface for the basic equation (10) by requiring that the stress rate at failure vanishes, i.e., r ¼ 0. Moreover, the failure surface will be derived at critical state as a major concern is to build critical state model based on the basic model.
Equation (10) can be separated into a spherical part and a deviatoric part. Let us first consider the spherical part, which can be obtained by taking the trace of both sides of Eq. (10) Note that tr r ¼ 0 and tr _ ¼ 0 at critical state. By making use of the relationship trr _ ¼ trr Ã _ Ã ¼ jjr Ã jjjj _ Ã jj cos h, with h being the angle between r Ã and _ Ã , the following equation can be obtained by letting tr r ¼ 0 Let r c denote the stress ratio jjr Ã jj=trr at critical state. The following expression can be obtained from the above equation The term cos h in the above equation represents the flow direction with reference to stress at critical state. As will be shown thereafter, we have cos h ¼ 1. In this case, we are left with the following relationship between c 3 and c 4 The above expression shows that the failure surface is that of Drucker-Prager. The stress ration r c is independent of the first and second terms in Eq. (10) and depends only on the ratio c 4 =c 3 Now, let us turn our attention to the deviatoric part of the constitutive equation (10) The above equation is quadratic in r c . Combining (15) and Eq. (18), we can solve for c 3 to obtain The parameter c 3 can be set into Eq. (15) to get c 4 To sum up, the failure surface from constitutive equation (10) can be expressed in the following form with r c ¼ Àc 4 =c 3 . Obviously, the failure surface (21) has a conical shape, which likes that of Drucker-Prager. If the same friction angle u crt is used for Drucker-Prager model and basic hypoplastic model, the failure surfaces of these two models will coincide. The trace of the failure surface on a p-plane with trr ¼ 1 is a circle with the radius r c . The flow rule is associated at critical state. Figure 10 shows that for initial stresses on the failure surface, some strain rates may lead to stresses outside the failure surface. We proceed to derive the bound surface, which enclose all accessible stresses following the procedure proposed by Wu and Niemunis [27]. Assume that there exists a bound surface bðrÞ ¼ 0 with bðrÞ being an isotropic function of stress. Consider the stress r b which happens to be on the bound surface so that bðr b Þ ¼ 0. The outward normal to the bound surface at r b is per definition

Bound surface
The bound surface can be found by solving the following equation system where L is the fourth-order tensor and N the second-order tensor from Eq. (3). For constitutive equation (10), solving the above equation system yields the following equation 9c 2 1 r 2 1 þ r 2 2 þ r 2 3 À Á À 2c 2 4 r 2 1 þ r 2 2 þ r 2 3 À r 2 r 3 À r 1 r 2 r 1 r 3 where r 1 , r 2 and r 3 are the principle stresses. The above equation represents a right circular cone with its apex at the origin. The failure surface and the bound surface of constitutive equation (10) are depicted on a p-plane in Fig. 11. The following material parameters are used: E i =r c ¼ 170, m i ¼ 0:2, / crt ¼ 30 and w crt ¼ 0 .
The bound surface ensures that the stresses from the constitutive equation (10) remain bounded. In numerical calculations, however, large step size may lead to stresses beyond the failure surface. In elastoplasticity, some return mapping algorithms are needed to bring the stress back to the yield surface [21]. The existence of bound surface can be regarded as a self-correcting mechanism for the hypoplastic constitutive model.

Conclusions
The hypoplastic model with critical state by Wolffersdorff is based on a basic model with three terms. The number of parameters is reduced from four in previous models to three. However, this basic model with three terms is severely flawed in that the model provides decent predictions only when the parameters are calibrated for vanishing initial Poisson ratio, which is at odds with experimental observations. The model becomes corrupted when realistic Poisson ratios between 0.2 and 0.3 are used.
A basic model is developed by introducing a new term, which vanishes at critical state. Its simplicity and robustness make it an ideal basic model to build more sophisticated models, e.g., critical state model and high-order  [5,8,9,13,30,31]. The failure surface is derived from the constitutive equation and likes that of Drucker-Prager. An extension to the failure surface of Matsuoka and Nakai is straightforward.