Development of a constitutive model for rockfills and similar granular materials based on the disturbed state concept

Abstract

Behavior of rockfills was investigated experimentally and theoretically. A series of standard triaxial compression tests were carried out on a quarried rockfill material at different stress levels. It was found that both the stress level and the shear stress ratio, like most of granular materials, controls the behavior of rockfill materials. At lower shear stress ratios the behavior is much more similar to a nonlinear elastic solid. When the shear stress goes further, the stress-strain curve shows an elasto-plastic behavior which suggests using the disturbed state concept to develop a constitutive model to predict the stress-strain behavior. The presented constitutive model complies reasonably with the experimental data.

Keywords

constitutive model / granular material / rockfill / plasticity / disturbed state concept / stress level

Introduction

Rockfills are of important materials with various applications in geotechnical engineering, among them construction of earth and rockfill dams, basement of highways, foundations and other applications can be addressed. Like other granular materials, rockfills behavior can be characterized by making use of the theory of elasticity, or in a more complex manner, application of the theory of plasticity. For granular materials, several constitutive models have been developed and applied in a variety of problems (nonlinear elastic model by Duncan and Chang [1]; elasto-plastic hardening models for non-cohesive soils by Lade and Duncan [2]; Nova and Wood [3]; elasto-plastic hardening model for soil, rock and concrete by Lade and Kim [4]; elasto-plastic model for strain softening materials by Guo and Lee [5]; development and application of a stress level based nonlinear model for sand by Veiskarami et al. [6]). In rockfill materials, the stress-strain behavior was modeled by linear or nonlinear elastic models in the past [7,8]; however, implementation of more advanced elasto-plastic models based on the theory of plasticity and the disturbed state concept [9,10] has been widely developed in the recent years [11,12]. Rockfill materials are often very large in their particles and hence, it is common to reduce the size of materials for testing purpose. There are several methods in scaling down the size of the materials, e.g., scalping (cut-off), replacement of materials, generation of quadratic size distribution curve and parallel gradation techniques [11]. The last method was found more appropriate for rockfill materials [11,13].

Very recently, Bauer et al. [14] studied the influence of pressure and density on the rheological properties of coarse-grained rockfills based on a hypoelastic constitutive model. They concluded that the behavior of rockfills depend highly on the crack propagation within the particles, abrasions and grain breakage. This fact is emphasized in this paper, i.e. the role of particle breakage in the behavior of granular (rockfill) materials. In the current study, the behavior of a quarried rockfill material is studied and the influence of the stress level and deviatoric stress ratio on the material behavior is investigated. Then, an elasto-plastic constitutive model, based on the disturbed state concept, has been developed.

Behind the theory

There are two distinct problems often confronted in geotechnical engineering: stability problems, requiring only the yield or the strength of materials and deformation problems which requires a constitutive model to describe material behavior. Establishment of a suitable constitutive model for geomaterials has been always the major concern of most of researchers, while applications are also many. Classical constitutive models for geomaterials often include an elasticity relationship along with a yield surface which can expand (or contract) by a hardening law. Elastic and nonlinear elastic models can also be regarded as particular cases of the most general elasto-plastic models which do not include irrecoverable strains. Such models are useful for very strong materials, like sound rocks or dense and coarse granular materials, in the range of small strains. When the stresses are very high and the strains are beyond the elastic limits, even strong materials can exhibit elasto-plastic behavior. In such cases, they are no longer “sound” or “strong,” rather, they have their states “changed” from their reference or original to a partly or entirely adjusted one. The following section presents a particular concept in which, a combination of both types of behavior constitutes the foundation of more versatile models for geomaterials.

Modeling granular materials and the disturbed state concept

Granular materials, unlike cohesive ones, exhibit more strength when they are subjected to higher confining pressure. The fact that the behavior of granular materials depends on the level of stress has been a basis in development of pressure dependent yield criteria and constitutive models like those presented by different researchers [2,15-18]. In the recent decades, besides nonlinear elastic and general elasto-plastic models in which, the material behavior is assumed to be uniform throughout the sample undergoing deformation, the concept of the disturbed state has been rather widely employed in modeling engineering materials behavior.

As stated by Desai [10], as the founder of the disturbed state concept (DSC), a deforming body does not exist in a unique composition, i.e. any material element is treated as a mixture of an initial part, also known as the relative intact or RI state, and the self-adjusted part also known as the fully adjusted or FA state. Since the deformation may start from a state which is not an ideally intact state, the term relative intact is used. In the other hand, the fully adjusted state is an idealized behavior which can be reached asymptotically after reasonably large deformation. The intermediate behaviors correspond to a mixture of the reference and the fully adjusted behaviors. In fact, the deformation process can be considered as a transition from a relative intact to a fully adjusted state during which, more and more portion of the material undergoes the self-adjusted state. The fully adjusted state is sometimes considered as the critical state. Once the deformation starts to occur, the transition between the reference and the fully adjusted state is governed and characterized by a so called disturbance parameter, D. This parameter should be determined to describe any intermediate behavior at different states of deformation [10]. It is although very hard to be well characterized, very recently, attempts have been made to simulate and investigate the particle breakage under static and dynamic loading [19,20]. Figure 1(a) shows the transition from the relative intact to fully adjusted states in a deforming body. In the relative intact (initial condition) there is no relative disturbance and hence, the parameter D equals 0. Formation of microcracks or further development of them along with particles relative movement alters the material behavior up to the fully adjusted or the critical state in which, no further disturbance will happen. This state corresponds to D = 1. Figure 1(b) shows the composition of the materials when they undergo deformation. This composition is a mixture of the relative intact and fully adjusted states controlled by the disturbance parameter, D.

According to Refs. [9, 10], the relative intact state can be described by either of an elastic or elasto-plastic behavior. Also, the fully adjusted state can be modeled by an elasto-plastic behavior, or, as an idealized behavior, by a viscous liquid which cannot sustain any shear stress. The disturbed state concept appears to be a versatile tool in capturing the stress-strain behavior of granular materials which can inherently capture the influence of the stress level on the material parameters. It is remarkable that the role of stress level on the geotechnical properties of geomaterials, in particular on their yield strength, has been studied by many researchers [21-24; Gan et al., 1988; Clark, 1998].

Applications of the disturbed state concept are many and cover a variety of problems. Known to the authors, application to failure analysis of rock masses subjected to unloading [25], characterization of the mesoscopic softening behavior of clays [26], calculation of the active and passive earth pressure behind retaining structures [27], applications to the liquefaction potential assessment [28] and even implementation in nonlinear finite element analysis of reinforced concrete frames [29], can be addressed as the most recent advances in applications of the DSC theory.

In the current study, a certain type of rockfill materials has been tested in a large scale triaxial apparatus under different confining pressures to develop a constitutive model based on the disturbed state concept (DSC) and to check the most consistent model for any of these two reference and fully-adjusted behaviors. The main concern is to develop a rather simple constitutive model with the minimum number of required model parameters.

Experimental data

An experimental program was conducted on a quarried angular rockfill material in a large scale triaxial apparatus. Coarse-grained rockfill materials were collected from the Shahriar pit located in Alborz province (in central region of Alborz Mountain range), Iran. It is among the largest barrow pits and a source for quarried materials, ranging from sand-size to cobble-size aggregates for many projects in the central and northern Iran. Consolidated-drained tests were carried out to measure both stress-strain and volumetric-strain behavior of samples at a wide range of confining pressures, as far as the capacity of the apparatus allowed. Supplementary tests were also performed to find out general rockfill geotechnical properties of the materials which are outlined in Table 1. These tests included the particle size distribution (aggregation) test (ASTM D422-63), Los Angeles test (ASTM C-131) and mineralogy (chemical) tests to describe different mineral constituent of material. Figure 2 shows an image of materials under study. It is remarkable that the aggregation curve was developed based on the Fuller’s curve (after Fuller and Thompson [30]) which accounts for the densest state of granular materials.

Prior to the start of the tests, a portion of materials were colored in dark blue to allow further investigation of their mechanical behavior (possible breakage) during the tests. Painting was performed on a small portion (less than say 2% by weight) and on randomly chosen particles. Particle breakage was then possible to be checked and the translation from the relative intact to fully adjusted state, partly because of the particle breakage was justifiable. A schematic view of the large scale triaxial apparatus is shown in Fig. 3 which allows the tests to be conducted under confining pressures up to 1000 kPa.

The triaxial compression tests were conducted at different confining pressures to cover a wide range of stress levels required for further study ranging between 50 and 800 kPa. Standard (conventional) triaxial tests were carried out in accordance with ASTM D7181-11 on disturbed (remolded) specimens. First, the oversize particles were removed. It was done by making use of the parallel gradation technique which was found reasonable by the researchers [11,13]. The maximum size of the sample was therefore 25.1 mm. Samples were prepared using a split mold and each layer was compacted using a 50 Hz vibrator. The sample density achieved by this method was 16.8 kN/m³. The sample dimensions were 153 mm in 302 mm which found to be reasonably large for the rockfill materials in this study. It was then subjected to saturation over a course of 6 h. It is remarkable that the sample mostly comprises quite coarse grained materials and hence, saturation was fulfilled very quickly. The deviatoric stress was applied by a rate of 0.76 mm/min after the first phase, i.e., application of the confining pressure.

To find the best fit to the locus of the yield surface in the principal stress space, a triaxial extension test was also performed. There are different techniques for triaxial extension tests. In the conventional method, the confining pressure increases as the vertical stress is held constant. However, this technique requires a large confining pressure to initiate the failure. In this research, the triaxial extension test was performed by reducing the cell pressure under a constant vertical stress. This test is denoted by RTE (reduced triaxial test) in the literature. Tests results are shown in Fig. 4 in terms of the stress-strain-volume change data whereas in Fig. 5, the normalized results are also plotted. This latter plot reveals two important matters: i) the ultimate state of all samples are almost the same when normalized to the confining pressure and ii) there is no peak strength in samples tested at higher confining pressures which can be a result of particle breakage and, further transition from a relatively intact to fully adjusted state. To check this conclusion, Fig. 6, shows the marked grains before and after the test under 800 kPa confining pressure. It is evident that although the strength of the particles is relatively very high, breakage occurred in at least a quarter of marked grains. A summary of general observations are provided in Table 2 revealing the influence of the stress level on the strength parameters of samples tested under different confining pressures.

Interpretation of the experimental results

The results revealed some major points which are important in investigation of the behavior of rockfills and also required for the purpose of this work. As stated earlier, disturbed state concept suggests a transition from a relatively intact state to a fully adjusted state when the loading advances. In this regard, an insight into the results indicated that the onset of inelastic behavior commences after 1.0 to 2.0% axial strain, below which, strains are almost fully recoverable. The fairly linear trend in the stress-strain behavior of all samples below certain strain levels, i.e., 1.0 to 2.0% can substantiate this conclusion. Therefore, the yield will never occur until sufficient amount of elastic strain in the observed stress-strain curve has been reached. This level, also called the threshold (after Lade and Duncan [2]) may differ from material to material. Since the elastic limit is an inherent material property, therefore the threshold seems to be highly dependent on the type, origin, shape and even the density of materials. A general conclusion on the threshold limit and its dependency on physical properties of a material require a deep study and a wide range of database which is out of the scope of this work. For the purpose of this work and as a practical conclusion, the threshold can be approximated from experimental data; here, it is set to be equal to 2.0% axial strain. This threshold value will be later used to define the initial yield function.

The Lade-Duncan yield surface can be a suitable yield criterion for the material under study regarding the type of the material, different behavior in compression and extension and dependency of the strength to the confining pressure. Table 3, represents the rendered results as required by the Lade-Duncan yield criterion. To check this conclusion, the projections (normalized to √3 to be transformed to the π-plane) of the ultimate stress state for both triaxial compression and extension tests, corresponding to the critical state or, the failure state, has been depicted in Fig. 7 on the π-plane, normalized to the confining pressure of relevant tests. It is evident that the Lade-Duncan yield surface can well capture the location of the ultimate stress state on both compression and extension axes.

It is also remarkable that tested samples were separated into two groups; one for model parameters calibration and the other for verification. These two sets of samples are distinguished in the last column of Table 3.

While the critical state friction angle, although is very scattered, shows a relatively constant value for all test results, variation of the peak friction angle with the confining pressure is important. The stress level dependency of soil friction angle has been relatively well investigated in the literature [21-24,31; Clark, 1998]. Figure 8, shows a plot of the peak friction angle versus the confining pressure. A power-law equation [Clark, 1998; 31] was fitted to the results complying relatively good with the observed results (R2= 0.87). The results indicate that the power-law equation still works reasonably when employed to capture the influence of the confining pressure on the variation of the peak friction angle. This conclusion is one of the simplest results observed in the experimental studies which can be directly used in stability problems, like the determination of the bearing capacity of foundations, when only the peak strength of materials is required and there is no need to the stress-strain relationship. This decrease in the friction angle can be related to the particle breakage, as a result of widening of microcracks occurred during the production process of quarried rockfills. In intact materials, such decrease can be related to the weak particles breakage and other factors preventing the dilatancy of the sample. Generally, an increase in the confining pressure leads to the reduction in tendency to dilative behavior.

Constitutive soil model

The behavior of samples during the tests at different confining pressures suggest that the samples be regarded as a mixture of both intact and fully adjusted parts, based on the disturbed state concept. Generally, a constitutive model consists of the following quadruple ingredients:

• An elastic relationship

• A yield criterion

• A plastic potential function

• A hardening law

Relative intact state

The first step to develop a constitutive model capable of capturing the stress-strain-volume change behavior, when the disturbed state concept is employed, is to find the stress-strain relationship in the relative intact state. According to the range of confining pressures and the definition of the relative intact state, it is convenient to assume a nonlinear elastic behavior. It can be found easily by stress-strain data obtained at the lowest confining pressure, i.e., at 50kPa. This confining pressure is the lowermost value below which, no test is available and hence, can be regarded as the initial state from which, the behavior of all samples commences. Following Varadarajan et al. [11], among others, this nonlinear stress-strain behavior can be expressed by the following equation:

(1){\mathit{\sigma }}_{\mathrm{1}}-{\mathit{\sigma }}_{\mathrm{3}}={\displaystyle \frac{{\mathit{ϵ}}_{\mathrm{1}}}{a+b{\mathit{ϵ}}_{\mathrm{1}}}},

In this equation, σ₁ and σ₃ are the major and minor principal stresses, ϵ₁ is the major principal strain and a and b are model parameters. These model parameters can be determined from data obtained experimentally at the early stages of loading under the lowermost confining pressure. In fact, a equals to 1/E_i where E_i is the initial (tangent) modulus of elasticity, i.e., the initial slope of the stress-strain curve under the lowermost confining pressure. The parameter b can be found by curve-fitting to the data from the beginning of the loading up to the peak strength. It should be remarked that it is assumed that the behavior of sample in those exhibiting peak strengths, is mostly comprises of a nonlinear elastic one, i.e., a relative intact behavior. This assumption is true for those samples tested at relatively low confining pressure since there is neither a particle breakage nor a reasonably large particles rearrangement/reorientation before the peak strength is reached; moreover, strains are small enough for an elastic behavior (practically smaller than 1% to 2%). Attention should be taken to the fact that the elastic behavior, although is associated with the behavior of all samples, changes from a confining pressure to another. It can be included in the development of the relative intact behavior by assuming the change of the initial modulus of elasticity as a function of the confining pressure via the following power-law equation:

(2)E_i=E_{\mathrm{ref}}{\left({\displaystyle \frac{{\mathit{\sigma }}_{\mathrm{3}}}{P_{\mathrm{a}}}}\right)}^M,

where, the exponent M is a parameter, determined experimentally by plotting the initial tangent modulus of all tests on a semi-logarithmic plot and E_ref is the initial tangent modulus corresponding to the unit confining pressure which can also be determined by extrapolation of the data to the unit confining pressure. Also, P_a is the atmosphere pressure to make the equation dimensionless. It should be noted that the Poisson’s ratio can also be easily determined from the axial strain-volume change diagram in the course of loading at the start of the application of the deviatoric stress change.

Fully adjusted state

The stress-strain behavior of the fully adjusted state was assumed to be a hardening elasto-plastic behavior which can be described by a linear elastic-hardening model with appropriate assumptions. As a central assumption in development of the DSC based constitutive model, it was assumed that the stress-strain behavior of samples is governed by the FA state when the sample is tested at the highest confining pressure and at reasonably large strains which guarantee an inelastic behavior to commence (e.g., beyond 2% to 3% axial strain). It is remarkable that the model will depend on the test results corresponding to the highest confining pressure. However, it should be noted that the fully adjusted state as well as the relative intact state are ideal states and in practice, this highest pressure should be chosen in a way to cover the entire range of the stress levels the material may experience in the actual problem. Therefore, such assumption may not be too far from the actual behavior of materials and seems to be fairly suitable. Model parameters corresponding to this state can be calibrated by making use of the data obtained for the test at the uppermost value of the confining pressure. In this stage, the yield criterion is of the highest importance since it will be needed to find the plastic strain increments. Several attempts have been made to find the most general form of a yield surface for coarse-grained geomeaterials. In particular for rockfill materials, Varadarajan et al. [11,12] following Desai [10] suggested using a yield function which can be degenerated to all renowned yield criteria suggested earlier. In fact, a proper selection of the parameters of this yield function gives yield criteria of Mohr-Coulomb (1776), Lade and Duncan [2] or some others found in the literature [10]. Although this yield criterion covers almost all possible yield conditions for geomaterials, but it involves a relatively large number of parameters which must be determined experimentally. To reduce the complications within this yield function and based on the experimental observations, also consistency between the results for coarse-grained materials, the yield criterion of Lade and Duncan has been chosen as the yield function for the further development of the constitutive model. This model is expressed in terms of the first and third stress invariants, I₁ and I₃ and a model parameter, k₁, expressed through the following equation:

(3)f(I_{\mathrm{1}},I_{\mathrm{3}},k_{\mathrm{1}})=I_{\mathrm{1}}^{\mathrm{3}}-k_{\mathrm{1}}{I}_{\mathrm{3}}=\mathrm{0},

(4a)I_{\mathrm{1}}=tr[{\mathit{\sigma }}_{ij}]={\mathit{\sigma }}_{ii},

(4b)I_{\mathrm{3}}=Det[{\mathit{\sigma }}_{ij}]=e_{ijk}{\mathit{\sigma }}_{i\mathrm{1}}{\mathit{\sigma }}_{j\mathrm{2}}{\mathit{\sigma }}_{k\mathrm{3}},

where, e_ijk is the permutation symbol and [σ_ij] is the matrix representation of the stress tensor. Note that Einstein summation convention holds and k₁ is the yield function parameter which is always greater than 27. A threshold value for k₁ should be found experimentally beyond which, plastic strains start to develop. The plastic potential function, required to derive the incremental plastic strains, was assumed to take the form analogous to the yield surface as suggested by Lade and Duncan [2]:

(5)g(I_{\mathrm{1}},I_{\mathrm{3}},k_{\mathrm{2}})=I_{\mathrm{1}}^{\mathrm{3}}-k_{\mathrm{2}}{I}_{\mathrm{3}},

Again, k₂ is a model parameters which is yet undetermined.

Moreover, the elastic part of the FA state is assumed to obey a linear elastic stress-strain relationship with the model parameters (elastic modulus) determined for the relative intact state. Now, comes the instant to define a relationship accounting for the model parameters k₁ and k₂. As stated earlier, the development of the DSC based constitutive model was started by assuming a nonlinear elastic stress-strain behavior for the relative intact state. Hence, if the model parameters are found, the FA state can be decomposed of the observed stress-strain behavior. It should be noted that yields only occurs in the fully adjusted state and hence, once the elasto-plastic stress-strain was extracted from the observed stress-strain data, the model parameters for the elasto-plastic state (i.e., the FA state) can be determined. It was done the results were used for further study of the relationship between the model parameters and the observed data. Since, k₁ varies during the test, it must be governed by some hardening parameter. There are several assumptions describing the hardening parameter. For example, the accumulated plastic strain [32], accumulated plastic work [2], the volumetric plastic strain rate [33,34].

To find the relationship between the model parameter, k₁ at different stages of loading, i.e., to establish the hardening rule, the plastic strains have been first decomposed from the observed stress-strain behavior of samples tested at the highest confining pressure, i.e., at 400 kPa and 800 kPa. Then, the trajectory of the plastic strains have been computed and plotted versus the computed values of k₁ at different stages of loading. This trajectory, ξ_D, is defined as:

(6){\xi }_D={\displaystyle \int {(\mathrm{d}{ϵ}_{ij}^{\mathrm{p}}\mathrm{d}{ϵ}_{ij}^{\mathrm{p}})}^{\mathrm{1}/\mathrm{2}},}

Prior to establish the governing equation, there are two important points which should be noticed. First, the yield does not immediately start with the deviatoric loading; rather, it starts from some threshold which should be determined experimentally. Therefore, plastic strains will start to occur once this threshold is reached. Second, the relationship between the total plastic strain trajectory and the k₁ value is valid only for the fully adjusted state since there are only elastic strains in the RI state, as mentioned earlier.

Therefore, the relationship between ξ_D and k₁ is focused particularly for those samples which have the most tendencies to show the FA state, i.e., those tested at the highest confining pressure. Although it is possible to find such a relationship for samples with lower tendency toward a FA state, i.e., at lower confining pressures, it was found that it would result in higher errors in capturing the stress-strain behavior of samples. It should be noted that decomposition at this stage was done by assuming that the material behaves, after some amount of strains, completely in a fully adjusted state. It corresponds to setting the value of the disturbance parameter, D, equal to unity when decomposing the plastic strains from the total strains. However, the value of D is not actually equal to 1.0 and hence, the model parameters found this way, should be refined in an iterative manner, when the model is calibrated, to obtain better results. For the time being, the initial values of model parameters relating k₁ to ξ_D have sufficient accuracy for the purpose of the work. Figure 9 shows the variation of the total plastic strain trajectory with k₁ for samples tested at 200 kPa, 400 kPa and 800 kPa confining pressures.

A hyperbolic equation of the following form:

(7)k_{\mathrm{1}}={\displaystyle \frac{{\xi }_D}{{\mathit{\chi }}_{}+{\mathit{\eta }}_{}{\xi }_D}},

was found to be reasonably fitted to the test data with parameters χ₁ and η₁.

By the same analogy, a similar form was adopted for the variations of k₂, the model parameter accounting for the plastic potential function:

(8)k_{\mathrm{2}}={\displaystyle \frac{{\xi }_D}{{\mathit{\chi }}_{}+{\mathit{\eta }}_{}{\xi }_D}},

Having known the plastic strains from the tests on samples at the highest confining pressure, it is possible to find the unknown parameters χ₁ and η₁ simply by a curve-fitting scheme. To determine the relationship between the value of k₂ and the total plastic strain trajectory however, it is required to compute k₂ at different stages of loading. To do so, following Lade and Duncan [2] the parameter, ν_p, is defined as the ratio of the lateral to the vertical plastic strain increment:

(9){\nu }_p=-{\displaystyle \frac{\mathrm{d}{ϵ}_{\mathrm{3}}^{\mathrm{p}}}{\mathrm{d}{ϵ}_{\mathrm{1}}^{\mathrm{p}}}},

This parameter is computed over the entire range of the observed stress-strain data (when decamped to give only the plastic strains, having known the nonlinear elastic RI state) after some threshold, and then applied to find the parameter k₂ by the flow rule as follows [32]:

(10)\mathrm{d}{ϵ}_{ij}^{\mathrm{p}}=\mathrm{d}\lambda {\displaystyle \frac{\partial g}{\partial {\sigma }_{ij}}},

In this equation, dλ is the plastic multiplier and \mathrm{d}{ϵ}_{ij}^{\mathrm{p}} is the plastic strain increment. Having known the equations for the plastic potential function one can make the following substitutions:

(11a)\mathrm{d}{ϵ}_{\mathrm{1}}^{\mathrm{p}}=\mathrm{d}\lambda {\displaystyle \frac{\partial g}{\partial {\sigma }_{\mathrm{1}}}}=\mathrm{d}\lambda \left({\displaystyle \frac{\partial g}{\partial I_{\mathrm{1}}}}{\displaystyle \frac{\partial I_{\mathrm{1}}}{\partial {\sigma }_{\mathrm{1}}}}+{\displaystyle \frac{\partial g}{\partial I_{\mathrm{3}}}}{\displaystyle \frac{\partial I_{\mathrm{3}}}{\partial {\sigma }_{\mathrm{1}}}}\right)=\mathrm{d}\lambda (\mathrm{3}{I}_{\mathrm{1}}^{\mathrm{2}}-k_{\mathrm{2}}{\sigma }_{\mathrm{3}}^{\mathrm{2}}),

(11b)\mathrm{d}{ϵ}_{\mathrm{3}}^{\mathrm{p}}=\mathrm{d}\lambda {\displaystyle \frac{\partial g}{\partial {\sigma }_{\mathrm{3}}}}=\mathrm{d}\lambda \left({\displaystyle \frac{\partial g}{\partial I_{\mathrm{1}}}}{\displaystyle \frac{\partial I_{\mathrm{1}}}{\partial {\sigma }_{\mathrm{3}}}}+{\displaystyle \frac{\partial g}{\partial I_{\mathrm{3}}}}{\displaystyle \frac{\partial I_{\mathrm{3}}}{\partial {\sigma }_{\mathrm{3}}}}\right)=\mathrm{d}\lambda (\mathrm{6}{I}_{\mathrm{1}}^{\mathrm{2}}-\mathrm{2}{k}_{\mathrm{2}}\sigma {\mathrm{1}}_{}{\sigma }_{\mathrm{3}}),

Hence, the parameter, ν_p, can be found as the ratio of these two expressions for \mathrm{d}{ϵ}_{\mathrm{1}}^{\mathrm{p}}and \mathrm{d}{ϵ}_{\mathrm{3}}^{\mathrm{p}}, eliminating dλ:

(12){\nu }_p=-{\displaystyle \frac{\mathrm{6}{I}_{\mathrm{1}}^{\mathrm{2}}-\mathrm{2}{k}_{\mathrm{2}}{\sigma }_{\mathrm{1}}{\sigma }_{\mathrm{3}}}{\mathrm{3}{I}_{\mathrm{1}}^{\mathrm{2}}-k_{\mathrm{2}}{\sigma }_{\mathrm{3}}^{\mathrm{2}}}},

This latter expression is different from that of Lade and Duncan [2] because of the difference in definition of the volumetric strain in their work and in the current study. Here, the volumetric strain was taken to be equal to ϵ_v = ϵ₁ + 2ϵ₃ as the case in most of continuum mechanics and plasticity texts. This value can be computed based on experimental data to find the value of k₂ during the test. Eventually, a plot of the variations of k₂ with ξ_D can be used to find the model parameters, χ₂ and η₂ for the plastic potential function.

The next step is to derive the required equations for the elasto-plastic stress-strain relationship in the FA state. To do it, the theory of plasticity is employed to find the necessary equations. Beginning with the first assumption, the strain increment consists of elastic and plastic parts, i.e.,

(13)\mathrm{d}{ϵ}_{ij}=\mathrm{d}{ϵ}_{ij}^{\mathrm{e}}+\mathrm{d}{ϵ}_{ij}^{\mathrm{p}},

The elastic strain increment is related to the stress increment by the elastic constitutive tensor, E_ijkl:

(14)\mathrm{d}{\sigma }_{ij}=E_{ijkl}\mathrm{d}{ϵ}_{kl}^{\mathrm{e}},

The plastic strain increment is related to the plastic potential function through the well known flow rule [32]:

(15)\mathrm{d}{ϵ}_{ij}^{\mathrm{p}}=\mathrm{d}\lambda {\displaystyle \frac{\partial g}{\partial {\sigma }_{ij}}},

Since during any elasto-plastic deformation, the yield criterion must hold, the consistency condition [35] requires:

(16)\mathrm{d}f(I_{\mathrm{1}},I_{\mathrm{3}},k)={\displaystyle \frac{\partial f}{\partial {\sigma }_{ij}}}\mathrm{d}{\sigma }_{ij}+{\displaystyle \frac{\partial f}{\partial k_{\mathrm{1}}}}\mathrm{d}{k}_{\mathrm{1}}=\mathrm{0},

Taking the dependency of the parameter k₁ to the plastic strain increment through the total plastic strain trajectory requires:

(17)\mathrm{d}f(I_{\mathrm{1}},I_{\mathrm{3}},k)={\displaystyle \frac{\partial f}{\partial {\sigma }_{ij}}}\mathrm{d}{\sigma }_{ij}+{\displaystyle \frac{\partial f}{\partial k_{\mathrm{1}}}}{\displaystyle \frac{\partial k_{\mathrm{1}}}{\partial {ϵ}_{ij}^{\mathrm{p}}}}\mathrm{d}{ϵ}_{ij}^{\mathrm{p}}=\mathrm{0},

However, the plastic strain increment is related to the stress increment by the flow rule and hence:

(18)\mathrm{d}f(I_{\mathrm{1}},I_{\mathrm{3}},k)={\displaystyle \frac{\partial f}{\partial {\sigma }_{ij}}}\mathrm{d}{\sigma }_{ij}+{\displaystyle \frac{\partial f}{\partial k_{\mathrm{1}}}}{\displaystyle \frac{\partial k_{\mathrm{1}}}{\partial {ϵ}_{ij}^{\mathrm{p}}}}\left(\mathrm{d}\lambda {\displaystyle \frac{\partial g}{\partial {\sigma }_{ij}}}\right)=\mathrm{0},

On the other hand, the elasticity equations can be used to derive the stress increment-elastic strain increment as follows:

(19)\mathrm{d}{\sigma }_{ij}=E_{ijkl}\mathrm{d}{ϵ}_{kl}^{\mathrm{e}}=E_{ijkl}(\mathrm{d}{ϵ}_{kl}-\mathrm{d}{ϵ}_{kl}^{\mathrm{p}})=E_{ijkl}\left(\mathrm{d}{ϵ}_{kl}-\mathrm{d}\lambda {\displaystyle \frac{\partial g}{\partial {\sigma }_{kl}}}\right),

This latter equation can be substituted in the consistency equation and further rearrangement will result in the equation for the plastic multiplier as follows:

(20)\mathrm{d}\lambda ={\displaystyle \frac{{\displaystyle \frac{\partial f}{\partial {\sigma }_{ij}}}{E}_{ijkl}\mathrm{d}{ϵ}_{kl}}{{\displaystyle \frac{\partial f}{\partial {\sigma }_{ij}}}{E}_{ijkl}{\displaystyle \frac{\partial g}{\partial {\sigma }_{kl}}}-{\displaystyle \frac{\partial f}{\partial k_{\mathrm{1}}}}{\displaystyle \frac{\partial k_{\mathrm{1}}}{\partial {ϵ}_{ij}^p}}{\displaystyle \frac{\partial g}{\partial {\sigma }_{ij}}}}},

Components of this equation have been derived based on the shape of the yield and the plastic potential surfaces and presented in the Appendix of this paper.

Disturbance parameter

The final parameter which has been left undetermined is the disturbance parameter, D, which is required to establish the constitutive model. The equation for the disturbance parameter suggested by Varadarajan et al. [11] requires relatively complicated elaboration in determination of supplementary and auxiliary parameters. In the current work, the original form of the disturbance parameter represented by Desai and Toth [9], following Weibull (1951) and Lewis (1982), has been employed with three parameters collected into the following equation:

(21)D=D_u\left[\mathrm{1}-{\left(\mathrm{1}+{\left[{\displaystyle \frac{{\xi }_D}{h}}\right]}^w\right)}^{-s}\right],

where w, s and h are parameters which control the width, shape and height of the disturbance function when plotted against ξ_D. D_u is a parameters which governs the actual value of D when the fully adjusted state is reached. This means, the fully adjusted state cannot be reached unless in an idealized condition; in practice, at relatively large strain, the portion of the fully adjusted state in the deforming matter is close to but not equal to unity. This portion is denoted by D_u. To find those parameters contributing in the disturbance function, further assumption was made: in accordance with Varadarajan et al. [11] and formerly, Desai [10], it is possible to make the following assumption:

(22)\mathrm{d}{ϵ}_{ij}^a=\mathrm{d}{ϵ}_{ij}^i=\mathrm{d}{ϵ}_{ij}^c,

It means that there is no relative motion between the RI and FA states, that is, the observed strain increments are the same as the strains in the RI and FA parts. Such assumption, although leads to errors in volume change prediction, simplifies the determination of the model parameters. Once the RI state and its parameters were found experimentally, it is possible to decompose the stress-strain relationship into the RI and FA states and compute the value of D at the observed points. Then a curve-fitting scheme could be used to determine the parameters w, s and h. An increase in the number of tests at different confining pressures improves the results. It is also possible to employ other procedures as outlined by Desai and Toth [9] to determine D experimentally.

It is remarkable that once the disturbance function has been found, the stress state at each stage of the loading can be computed as follows [10]:

(23)\mathrm{d}{\sigma }_{ij}^a=(\mathrm{1}-D)\mathrm{d}{\sigma }_{ij}^i+D\mathrm{d}{\sigma }_{ij}^c+\mathrm{d}D(\mathrm{d}{\sigma }_{ij}^c-\mathrm{d}{\sigma }_{ij}^c),

In this equation, \mathrm{d}{\sigma }_{ij}^a, \mathrm{d}{\sigma }_{ij}^i and \mathrm{d}{\sigma }_{ij}^c are stress increments in the observed, RI and FA states respectively.

Predictions

The model parameters have been determined for the experimental data and are presented in Tables 4 and 5 for relative intact and fully adjusted states respectively. It is remarkable that a threshold value for k₁ was also found to be roughly equal to 38. This value was determined in accordance with the threshold value of the axial strains equal to 2%. Scattering in the stress-strain data points prevents better approximation in determination of model parameters, which is although an inherent property of all constitutive models developed for rocks and rockfill materials. It is mainly due to lack of certainty on the observed properties in experimental tests.

These model parameters were implemented into a developed computer code in MATLAB to predict the stress-strain behavior of the rockfill materials. Figure 10 shows the predicted versus measured data points in terms of stress-strain and axial strain-volumetric strain plots for all samples comprising those separated for model calibration and those put aside for further verification. It is obvious that the developed constitutive model can predict the stress-strain-volume change behavior with reasonable accuracy. It should be noted that the model parameters, after they had been initialized by procedure outlined earlier, were refined by further calibration of the predicted to measured data. Better results can be obtained by further refinements to get more accurate model parameters.

The sensitivity of the predictions to variations in the model parameters was also checked. Among all parameters, it was found by inspection that the model showed an insignificant sensitivity to change in the model parameters η₁ and η₂. It is evident from plotted data on Fig. 11 that±33 change in η₁ (i.e., 33% change in the basic value used in predictions) has no significant effect on the predictions for samples tested at 400 kPa and 800 kPa confining pressures respectively. On the other hand, both model parameters η₁ and η₂ are fairly close to each other. It could be expected from the behavior of k₁ and k₂; since they obey the same trend. Therefore, model parameters η₁ and η₂ can be assumed to be equal or used interchangeably for practical purposes. Moreover, selection of the Poisson’s ratio as it varies in a limited range for rockfill materials (say 0.2 to 0.3) has also a little effect on the predicted behavior of the materials by the developed model. However, the response to the model parameters χ₁ and χ₂ was found to be considerable and they must be chosen with care. This is the same for the parameters involved in determination of the disturbance parameters. Therefore, by making use of the abovementioned simplifying assumptions, the number of model parameters can be reduced to some extents.

Conclusions

An experimental study was conducted to investigate the behavior of rockfill materials in triaxial compression test. A series of large scale triaxial compression tests were performed on a sample of coarse-grained quarried rockfill material obtained from Shahriar pit located in Iran. Experimental data showed that there are primarily some important points in the stress-strain behavior of rockfill materials. First of all, under relatively high confining pressures, particles start to break. This breakage can be emerged from the widening of former microcracks occurred during the production process or existence of weak materials. This particle breakage decreases the tendency of materials to dilate which in turn results in the peak friction angle reduction. A power-law equation, formerly suggested by researchers, was found to be still useful in prediction of the peak friction angles in the course of the confining pressures applied in this study.

The disturbed state concept (DSC) of Desai [10] was then employed to develop a relatively simple constitutive model for rockfill materials; this time, with making use of the Lade-Duncan yield criterion and a nonlinear elastic model for the relative intact state. Observations indicated that both assumptions worked properly and the developed constitutive model can capture the stress-strain behavior of the rockfill materials with relatively good accuracy. Model parameters within the developed constitutive model can be found experimentally, but further refinements are required to achieve the most reasonable results. Predictions were made after the model parameters were found, on the observed data in the experimental program. A computer code in MATLAB was developed to compute the stress-strain-volume change behavior of samples based on the developed constitutive model. The developed constitutive model was found to be capable of capturing both dilative and contractive behaviors. The advantage of this model over former models is its fewer model parameters which allow the model parameters to be determined in a more convenient way. Sensitivity of the model to the model parameters was checked and it was found that it is possible to reduce the number of model parameters when some simplifying assumptions are made.

Appendix: Derivation of the required equations for the constitutive model

Required terms in computing dλ:

Recalling that:

(A-1)\mathrm{d}\mathit{\lambda }={\displaystyle \frac{{\displaystyle \frac{\partial f}{\partial {\mathit{\sigma }}_{ij}}}{E}_{ijkl}\mathrm{d}{\mathit{ϵ}}_{kl}}{{\displaystyle \frac{\partial f}{\partial {\mathit{\sigma }}_{ij}}}{E}_{ijkl}{\displaystyle \frac{\partial g}{\partial {\mathit{\sigma }}_{kl}}}-{\displaystyle \frac{\partial f}{\partial k_{\mathrm{1}}}}{\displaystyle \frac{\partial k_{\mathrm{1}}}{\partial {\mathit{ϵ}}_{ij}^{\mathrm{p}}}}{\displaystyle \frac{\partial g}{\partial {\mathit{\sigma }}_{ij}}}}},

First one should note that:

(A-2){\displaystyle \frac{\partial f}{\partial {\sigma }_{ij}}}\mathrm{d}{\mathit{\sigma }}_{ij}={\displaystyle \frac{\partial f}{\partial I_{\mathrm{1}}}}{\displaystyle \frac{\partial I_{\mathrm{1}}}{\partial {\mathit{\sigma }}_{ij}}}\mathrm{d}{\mathit{\sigma }}_{ij}+{\displaystyle \frac{\partial f}{\partial I_{\mathrm{3}}}}{\displaystyle \frac{\partial I_{\mathrm{3}}}{\partial {\mathit{\sigma }}_{ij}}}\mathrm{d}{\mathit{\sigma }}_{ij},

(A-3){\displaystyle \frac{\partial g}{\partial {\sigma }_{ki}}}\mathrm{d}{\mathit{\sigma }}_{ki}={\displaystyle \frac{\partial g}{\partial I_{\mathrm{1}}}}{\displaystyle \frac{\partial I_{\mathrm{1}}}{\partial {\mathit{\sigma }}_{ki}}}\mathrm{d}{\mathit{\sigma }}_{ki}+{\displaystyle \frac{\partial g}{\partial I_{\mathrm{3}}}}{\displaystyle \frac{\partial I_{\mathrm{3}}}{\partial {\mathit{\sigma }}_{ki}}}\mathrm{d}{\mathit{\sigma }}_{ki},

Moreover:

(A-4a){\displaystyle \frac{\partial f}{\partial I_{\mathrm{1}}}}{\displaystyle \frac{\partial I_{\mathrm{1}}}{\partial {\mathit{\sigma }}_{ij}}}\mathrm{d}{\mathit{\sigma }}_{ij}={\displaystyle \frac{\partial f}{\partial I_{\mathrm{1}}}}{\displaystyle \frac{\partial I_{\mathrm{1}}}{\partial {\mathit{\sigma }}_{\mathrm{1}}}}\mathrm{d}{\mathit{\sigma }}_{\mathrm{1}}+{\displaystyle \frac{\partial f}{\partial I_{\mathrm{3}}}}{\displaystyle \frac{\partial I_{\mathrm{3}}}{\partial {\mathit{\sigma }}_{\mathrm{1}}}}\mathrm{d}{\mathit{\sigma }}_{\mathrm{1}},

(A-4b){\displaystyle \frac{\partial f}{\partial I_{\mathrm{3}}}}{\displaystyle \frac{\partial I_{\mathrm{3}}}{\partial {\mathit{\sigma }}_{ij}}}\mathrm{d}{\mathit{\sigma }}_{ij}={\displaystyle \frac{\partial f}{\partial I_{\mathrm{1}}}}{\displaystyle \frac{\partial I_{\mathrm{1}}}{\partial {\mathit{\sigma }}_{\mathrm{3}}}}\mathrm{d}{\mathit{\sigma }}_{\mathrm{3}}+{\displaystyle \frac{\partial f}{\partial I_{\mathrm{3}}}}{\displaystyle \frac{\partial I_{\mathrm{3}}}{\partial {\mathit{\sigma }}_{\mathrm{3}}}}\mathrm{d}{\mathit{\sigma }}_{\mathrm{3}},

(A-5a){\displaystyle \frac{\partial g}{\partial I_{\mathrm{1}}}}{\displaystyle \frac{\partial I_{\mathrm{1}}}{\partial {\mathit{\sigma }}_{kl}}}\mathrm{d}{\mathit{\sigma }}_{kl}={\displaystyle \frac{\partial g}{\partial I_{\mathrm{1}}}}{\displaystyle \frac{\partial I_{\mathrm{1}}}{\partial {\mathit{\sigma }}_{\mathrm{1}}}}\mathrm{d}{\mathit{\sigma }}_{\mathrm{1}}+{\displaystyle \frac{\partial g}{\partial I_{\mathrm{3}}}}{\displaystyle \frac{\partial I_{\mathrm{3}}}{\partial {\sigma }_{\mathrm{1}}}}\mathrm{d}{\mathit{\sigma }}_{\mathrm{1}},

(A-5b){\displaystyle \frac{\partial g}{\partial I_{\mathrm{3}}}}{\displaystyle \frac{\partial I_{\mathrm{3}}}{\partial {\mathit{\sigma }}_{kl}}}\mathrm{d}{\mathit{\sigma }}_{kl}={\displaystyle \frac{\partial g}{\partial I_{\mathrm{1}}}}{\displaystyle \frac{\partial I_{\mathrm{1}}}{\partial {\mathit{\sigma }}_{\mathrm{3}}}}\mathrm{d}{\mathit{\sigma }}_{\mathrm{3}}+{\displaystyle \frac{\partial g}{\partial I_{\mathrm{3}}}}{\displaystyle \frac{\partial I_{\mathrm{3}}}{\partial {\mathit{\sigma }}_{\mathrm{3}}}}\mathrm{d}{\mathit{\sigma }}_{\mathrm{3}},

The partial derivatives are:

(A-6)\begin{array}{lll}{\displaystyle \frac{\partial f}{\partial I_{\mathrm{1}}}}=\mathrm{3}{I}_{\mathrm{1}}^{\mathrm{2}}& \mathrm{a}\mathrm{n}\mathrm{d}& {\displaystyle \frac{\partial f}{\partial I_{\mathrm{3}}}}=k_{\mathrm{1}}\end{array},

(A-7)\begin{array}{lll}{\displaystyle \frac{\partial g}{\partial I_{\mathrm{1}}}}=\mathrm{3}{I}_{\mathrm{1}}^{\mathrm{2}}& \mathrm{a}\mathrm{n}\mathrm{d}& {\displaystyle \frac{\partial g}{\partial I_{\mathrm{3}}}}=k_{\mathrm{2}}\end{array},

(A-8){\displaystyle \frac{\partial I_{\mathrm{1}}}{\partial {\mathit{\sigma }}_{\mathrm{1}}}}=\mathrm{1},{\displaystyle \frac{\partial I_{\mathrm{3}}}{\partial {\mathit{\sigma }}_{\mathrm{1}}}}={\mathit{\sigma }}_{\mathrm{3}}^{\mathrm{2}},{\displaystyle \frac{\partial I_{\mathrm{1}}}{\partial {\mathit{\sigma }}_{\mathrm{3}}}}=\mathrm{2},{\displaystyle \frac{\partial I_{\mathrm{3}}}{\partial {\mathit{\sigma }}_{\mathrm{3}}}}=\mathrm{2}{\mathit{\sigma }}_{\mathrm{1}}{\mathit{\sigma }}_{\mathrm{3}},

Also:

(A-9){\displaystyle \frac{\partial k_{\mathrm{1}}}{\partial {\mathit{ϵ}}_{ij}^{\mathrm{p}}}}={\displaystyle \frac{\partial k_{\mathrm{1}}}{\partial {\mathit{ϵ}}_{\mathrm{1}}^{\mathrm{p}}}}+{\displaystyle \frac{\partial k_{\mathrm{1}}}{\partial {\mathit{ϵ}}_{\mathrm{3}}^{\mathrm{p}}}}={\displaystyle \frac{\partial k_{\mathrm{1}}}{\partial {\xi }_D}}{\displaystyle \frac{\partial {\xi }_D}{\partial {\mathit{ϵ}}_{\mathrm{1}}^{\mathrm{p}}}}+{\displaystyle \frac{\partial k_{\mathrm{1}}}{\partial {\xi }_D}}{\displaystyle \frac{\partial {\xi }_D}{\partial {\mathit{ϵ}}_{\mathrm{3}}^{\mathrm{p}}}},

(A-10){\displaystyle \frac{\partial k_{\mathrm{1}}}{\partial {\xi }_D}}={\displaystyle \frac{\partial }{\partial {\xi }_D}}\left({\displaystyle \frac{{\xi }_D}{{\mathit{\chi }}_{\mathrm{1}}+{\mathit{\eta }}_{\mathrm{1}}{\xi }_D}}\right)={\displaystyle \frac{{\mathit{\chi }}_{\mathrm{1}}}{{({\mathit{\chi }}_{\mathrm{1}}+{\mathit{\eta }}_{\mathrm{1}}{\xi }_D)}^{\mathrm{2}}}},

(A-11){\displaystyle \frac{\partial {\xi }_D}{\partial {\mathit{ϵ}}_{\mathrm{1}}^{\mathrm{p}}}}={\displaystyle \frac{\partial }{\partial {\mathit{ϵ}}_{\mathrm{1}}^{\mathrm{p}}}}{\displaystyle \int \sqrt{{(\mathrm{d}{\mathit{ϵ}}_{\mathrm{1}}^{\mathrm{p}})}^{\mathrm{2}}+\mathrm{2}{(\mathrm{d}{\mathit{ϵ}}_{\mathrm{3}}^{\mathrm{p}})}^{\mathrm{2}}}}={\displaystyle \frac{\partial }{\partial {\mathit{ϵ}}_{\mathrm{1}}^{\mathrm{p}}}}{\displaystyle \int \sqrt{\mathrm{1}+\mathrm{2}{(-{\mathit{\nu }}_p)}^{\mathrm{2}}}}\mathrm{d}{\mathit{ϵ}}_{\mathrm{1}}^{\mathrm{p}}=\sqrt{\mathrm{1}+\mathrm{2}{(-{\mathit{\nu }}_p)}^{\mathrm{2}}},

(A-12){\displaystyle \frac{\partial {\xi }_D}{\partial {\mathit{ϵ}}_{\mathrm{3}}^{\mathrm{p}}}}={\displaystyle \frac{\partial }{\partial {\mathit{ϵ}}_{\mathrm{3}}^{\mathrm{p}}}}{\displaystyle \int \sqrt{{(\mathrm{d}{\mathit{ϵ}}_{\mathrm{1}}^{\mathrm{p}})}^{\mathrm{2}}+\mathrm{2}{(\mathrm{d}{\mathit{ϵ}}_{\mathrm{3}}^{\mathrm{p}})}^{\mathrm{2}}}}={\displaystyle \frac{\partial }{\partial {\mathit{ϵ}}_{\mathrm{3}}^{\mathrm{p}}}}{\displaystyle \int \sqrt{\mathrm{2}+{(-{\mathit{\nu }}_p)}^{-\mathrm{2}}}}\mathrm{d}{\mathit{ϵ}}_{\mathrm{3}}^{\mathrm{p}}=\sqrt{\mathrm{2}+{(-{\mathit{\nu }}_p)}^{-\mathrm{2}}}.

Nomenclatures

References

[1] Duncan J M, Chang C Y. Nonlinear analysis of stress and strain in soils. Journal of Soil Mechanics and Foundations Division, 1970, 96(SM5): 1629-1653

[2] Lade P V, Duncan J M. Elastoplastic stress-strain theory for cohesionless soil. Journal of the Geotechnical Engineering Division, 1975, 101(GT10): 1037-1053

[3] Nova R, Wood D M. A constitutive model for sand in triaxial compression. International Journal for Numerical and Analytical Methods in Geomechanics, 1979, 3(3): 255-278

[4] Lade P V, Kim M K. Single hardening constitutive model for soil, rock and concrete. International Journal of Solids and Structures, 1995, 32(14): 1963-1978

[5] Guo R, Li G. Elasto-plastic constitutive model for geotechnical materials with strain-softening behaviour. Computers and Geotechnics, 2008, 34: 14-23

[6] Veiskarami M, Jahanandish M, Ghahramani A. Prediction of foundations behavior by a stress level based hyperbolic soil model and the ZEL method, Computational Methods in Civil Engineering (CMCE). University of Guilan Press, 2010, 1(1): 37-54

[7] Kulhawy F H, Duncan J M. Stresses and movements in oroville dam. Journal of Soil Mechanics and Foundation Engineering Division, 1972, 98(7): 653-665

[8] Escuder I, Andreu J, Rechea M. An Analysis of Stress-Strain Behaviour and Wetting Effects on Quarried Rock Shells. Canadian Geotechnical Journal, 2005, 42(1): 51-60

[9] Desai C S, Toth J. Disturbed state constitutive modeling based on stress-strain and nondestructive behavior, International Journal of Solids and Structures, 1996, 33(11): 1619-1650

[10] Desai C S. Mechanics of Materials and Interfaces-the Disturbed State Concept, CRC Press, 2001

[11] Varadarajan A, Sharma K G, Venkatachalam K, Gupta A K. Testing and modeling two rockfill materials. Journal of Geotechnical and Geoenvironmental Engineering, 2003, 129(3): 206-218

[12] Varadarajan A, Sharma K G, Abbas S M, Dhawan A K. Constitutive model for rockfill materials and determination of material constants. International Journal of Geomechanics, 2006, 6(4): 226-237

[13] Ramamurthy T, Gupta K K. Response paper to how ought one to determine soil parameters to be used in the design of earth and rockfill dams. In: Proceedings of Indian Geotechnical Conference, New Delhi, India, 1986, 2: 15-19

[14] Bauer E, Fu Z, Liu S. Influence of pressure and density on the rheological properties of rockfills. Frontiers of Structural and Civil Engineering, 2012, 6(1): 25-34

[15] Coulomb C A. Essai Sur une Application des règles des Maximis et Minimis à Quelques Problemes de Statique Relatifs à l’architecture. Paris: De l'Imprimerie Royale, 1776, 5, 7

[16] Von Mises R. Mechanik der Festen Koerper im Plastisch-Deformablen Zustand. Nachrichten von der Gesellschaft der Wissenschaften zu Göttingen, Mathematisch-Physikalische Klasse, 1913: 582-592

[17] Drucker D C, Prager W. Soil mechanics and plastic analysis or limit design. Quarterly of Applied Mathematics, 1952, 10: 157

[18] Matsuoka T, Nakai K. Stress-deformation and strength characteristics of soil under three different principal stresses. Proceedings of Japan Society of Civil Engineers, 1974, 232: 59-70

[19] Ezaoui A, Lecompte T, Di Benedetto H, Garcia E. Effects of various loading stress paths on the stress-strain properties and on crushability of an industrial soft granular material. Granular Matter, 2011, 13(4): 283-301

[20] Bouteiller C, Naaim M. Aggregate breakage under dynamic loading. Granular Matter, 2011, 13(4): 385-393

[21] Meyerhof G G. The bearing capacity of sand. Dissertation for the Doctoral Degree, London: University of London, 1950

[22] Bolton M D. The strength and dilatancy of sands. Geotechnique, 1986, 36(1): 65-78

[23] Maeda K, Miura K. 1999, Confining stress dependency of mechanical properties of sands. Soils and Foundations. Japanese Geotechnical Society, 1999, 39(1): 53-67

[24] Kumar J, Raju K V S B, Kumar A. Relationships between rate of dilation, peak and critical state friction angles. Indian Geotechnical Journal, 2007, 37 (1): 53-63

[25] GanJ K M, FredlundD G, RahardjoH. Determination of the shear strength parameters of an unsaturated soil using the direct shear test. Canadian Geotechnical Journal, 1988, 25: 500–510

[26] ClarkJ I. The settlement and bearing capacity of very large foundations on strong soils: 1996 R.M. Hardy keynote address. Canadian Geotechnical Journal, 1998, 35: 131–145

[27] Zheng J Y, Wu A L. Mesoscopic analysis of the utilization of hardening model for a description of softening behavior based on disturbed state concept theory. Journal of Zhejiang University, Science A, 2008, 9(9): 1167-1175

[28] Wu G, Zhang L. Studying unloading failure characteristics of a rock mass using the disturbed state concept. International Journal of Rock Mechanics and Mining Sciences, 2004, 41(3): 437

[29] Zhu J F, Xu R Q, Li X R, Chen Y K. Calculation of earth pressure based on disturbed state concept theory. Journal of Central South University of Technology, 2011, 18(4): 1240-1247

[30] Park K, Choi J, Park I. Assessment of liquefaction potential based on modified disturbed state concept. KSCE Journal of Civil Engineering, 2012, 16(1): 55-67

[31] Akhaveissy A H, Desai C S. Application of the DSC model for nonlinear analysis of reinforced concrete frames. Finite Elements in Analysis and Design, 2012, 50: 98-107

[32] Fuller W, Thompson S E, The laws of proportioning concrete. Transactions of the American Society of Civil Engineers, 1907, 1053: 67-143

[33] Veiskarami M, Jahanandish M, Ghahramani A. Prediction of the Bearing Capacity and Load-Displacement Behavior of Shallow Foundations by the Sstress-Level-Based ZEL Method. Scientia Iranica, 2011, 18(1): 16-27

[34] Hill R. The Mathematical Theory of Plasticity, Oxford: Clarendon Press, 1950

[35] Schofield A N,Wroth C P. The Critical State Soil Mechanics. London: McGraw-Hill, 1958

[36] Yu H S. CASM: A unified state parameter model for clay and sand. International Journal for Numerical and Analytical Methods in Geomechanics, 1998, 22(8): 621-653

[37] Prager W. Recent developments of mathematical theory of plasticity. Journal of Applied Physics, 1949, 20(3): 235

[38] WeibullW A. A statistical distribution function of wide applicability. Applied Mechanics, 1951, 18: 293–297

[39] LewisC D. Industrial and Business Forecasting Methods. London: Butterworth Scientific, 1982

RIGHTS & PERMISSIONS

Higher Education Press and Springer-Verlag Berlin Heidelberg