Deformations and stability of granular soils: Classical triaxial tests and numerical results from an incremental model

The results presented in this article are based on experiments carried out in a classical triaxial apparatus manufactured by GDS Instruments [see 25, 47]. Vertical and lateral local displacements are measured by special gauges, which use the Hall effect. The gauges are installed directly on the sand sample. Vertical stress is controlled by moving the table on which the sample is placed, and lateral stresses are applied through the water pressure in the chamber. In the case of saturated or fully drained samples, the pore pressure and volume of water in the sample are read from a pneumatic controller. Figure 1 presents the preparation of a sample for testing.

Figure 1

The whole system for triaxial testing consists of a tri-axial chamber, three pressure/volume controllers, an electronic control system and a computer. The apparatus has a Bishop and Wesley chamber [see 3] with a built-in force sensor and a lower load cell, which imposes a vertical load on the sample. The table is moved by increasing the pressure in the lower load cell. With a rigid piston in the upper part of the chamber, a vertical compression force is generated.

The controllers are used for controlling the pressure in the main chamber of the apparatus, for controlling the load in the lower cell and for measuring the pore pressure in the sand sample. They can operate in a pressure change mode or a water volume change mode, so it is possible to perform stress- or strain-controlled tests.

This article presents the results of experiments performed on samples with a diameter of 38 mm and a height of 80 mm. It is important to ensure the homogeneity of the soil sample when it is being formed. The soil samples were prepared in a membrane-lined split moulder. Loose samples were formed by the moist tamping method and dense samples by water pluviation.

All experiments presented in this article were performed on coarse quartz Skarpa sand characterised by the following parameters: mean particle size D₅₀ = 0.42 mm, uniformity coefficient U = 2.5, minimum void ratio e_min = 0.432, maximum void ratio e_max = 0.677, angle of internal friction for loose/dense sand = 34°/41°. In all analyses, the soil mechanics sign convention is applied, in which compression is positive.

The original incremental model proposed by Sawicki and Świdziński [40, 41] has the form of semi-empirical incremental equations and describes the deformation of soils before the limit state is reached. The equations are given as follows: (1)dεv=Mdp′+Ndqd\varepsilon _v = Mdp^\prime + Ndq(2)dεq=Pdp′+Qdpd\varepsilon _q = Pdp^\prime + Qdp The invariants used in Eqs. (1) and (2) are the following: p′=1/3(σ1'+2σ3')p^\prime = 1/3\left( {\sigma _1^\prime + 2 \sigma _3^\prime } \right) , mean effective stress; q=σ1'−σ3'q = \sigma _1^\prime - \sigma _3^\prime , deviatoric stress; ε_v = ε₁ + 2ε₃, volumetric strain; and ε_q = 2/3 (ε₁ − ε₃), deviatoric strain, where σ1'\sigma _1^\prime , ε₁ are the normal effective vertical stress and vertical strain, respectively, and σ3'\sigma _3^\prime , ε₃ are the normal horizontal effective stress and horizontal strain, respectively.

The functions M, N, P, Q are determined from analytical approximations of selected experimental results [see 40]. The definitions of loading and unloading assumed in this model are the following: (3)dp′>0 is the spherical loading,dp^\prime > 0\;{\rm{is}}\;{\rm{the}}\;{\rm{spherical}}\;{\rm{loading}},(4)dp′<0 is the spherical unloading,dp^\prime < 0\;{\rm{is}}\;{\rm{the}}\;{\rm{spherical}}\;{\rm{unloading}},(5)dq>0 or dη>0 is the deviatoric loading,dq > 0\;{\rm{or}}\;d\eta > 0\,{\rm{is}}\;{\rm{the}}\;{\rm{deviatoric}}\;{\rm{loading}},(6)dq<0 or dη<0 is the deviatoric unloadingdq < 0\;{\rm{or}}\;d\eta < 0\,{\rm{is}}\;{\rm{the}}\;{\rm{deviatoric}}\;{\rm{unloading}} where η is a dimensionless stress ratio defined as (7)η=qp′\eta = {q \over {p^\prime}}

In Sawicki’s original model, the signs of dq or dη determine whether loading or unloading is applied (Eqs. 5 and 6), but this generally leads to contradictions. One has to consider the sign of dη, because the information on the sign of dq is not sufficient to clearly determine whether loading or unloading is applied, which is shown by predictions for stress paths for q = const. Therefore, to definite deviatoric loading and unloading, it was necessary to include a condition in which the sign of dη is specified. When this fact is taken into account, the model gives good predictions of soil deformations corresponding to the stress path for q = const or stress paths characteristic of the liquefaction phenomenon [see 41]. However, such definitions lead to some contradictions, and there is a part of the stress plane in which whether the case considered represents deviatoric loading or unloading is unknown. To prove it, let us consider a stress path originating at a point K in the p′ − q space located within the grey zone in Figure 2 – for example, the path KL. It is evident that the path KL can be described by conditions dq > 0 and dη < 0. Therefore, in this case, conditions (5) and (6) are unclear, and the kind of loading cannot be clearly defined. Consequently, it is unclear how precisely to build equations (1) and (2). In Sawicki’ s model, this problem is solved experimentally. It is simply checked, which sign (dq or dη) results in a better prediction. However, with this approach, the model still cannot be algorithmised, and, therefore, the present authors propose new definitions of deviatoric loading and unloading and, hence, a modification of the main equations.

Figure 2

The revised definitions have the following forms: (8)dp′>0 is the spherical loading,dp^\prime > 0\;{\rm{is}}\;{\rm{the}}\;{\rm{spherical}}\;{\rm{loading}},(9)dp′<0 is the spherical unloading,dp^\prime < 0\;{\rm{is}}\;{\rm{the}}\;{\rm{spherical}}\;{\rm{unloading}},(10)dη>0 is the deviatoric loading,d\eta > 0\;{\rm{is}}\;{\rm{the}}\;{\rm{deviatoric}}\;{\rm{loading}},(11)dη<0 is the deviatoric unloading,d\eta < 0\;{\rm{is}}\;{\rm{the}}\;{\rm{deviatoric}}\;{\rm{unloading}}, The definitions of spherical loading and unloading are the same as in the previous case (Eqs. 3 and 4), whereas the definitions of deviatoric loading and unloading are modified and only the sign of dη is considered. Thus, the incremental equations have the following forms: (12)dεv=Mdp′+Ndηd\varepsilon _v = Mdp^\prime + Nd\eta (13)dεq=Pdp′+Qdηd\varepsilon _q = Pdp^\prime + Qd\eta

The functions M, N, P, Q are determined from analytical approximations of selected experimental results. However, their forms are different from those of M, N, P, Q in Eqs. (1) and (2). The basic stress paths that are needed to determine these functions are different too.

To illustrate the physical meaning of the variable η, the possible directions and signs of its increments are shown in the p′ − η plane in Figure 3. The proposed incremental equations (12) and (13) are built based on the p′ and η invariants, and, therefore, this space is used in the description. When the stress path is defined by the condition dη > 0, the soil is approaching the limit state. The condition dη < 0 means that the soil is ‘moving away’ from collapse, whereas at a constant η, the soil is stable.

Figure 3

The following units are adopted in the model: 10⁵ N/m², corresponding to 100 kPa, for stress and 10⁻³, corresponding to 0.1%, for strain. The values of all stresses and strains used in the proposed model follow this convention, and the experimental data are presented in these units.

In the modified equations, the functions M, P correspond to deformations for the stress paths η = const, whereas the functions N, Q correspond to deformations for p′ = const. In each of the functions M, N, P, Q, a distinction has to be made regarding the sign of dp′ and dη, which results in two types of each function. The upper indices ^l and ^u correspond, respectively, to loading and unloading. The definitions of these functions are the following: (14)Ml=∂εvOA∂p′, Mu=∂εvAO∂p′M^l = {{\partial \varepsilon _v^{OA} } \over {\partial p^\prime}},\;\;\;\;\;M^u = {{\partial \varepsilon _v^{AO} } \over {\partial p^\prime}}(15)Nl=∂εvAB∂η, Nu=∂εvAB∂ηN^l = {{\partial \varepsilon _v^{AB} } \over {\partial \eta }},\;\;\;\;\;N^u = {{\partial \varepsilon _v^{AB} } \over {\partial \eta }}(16)Pl=∂εqOA∂p′, Pu=∂εqAO∂p′P^l = {{\partial \varepsilon _q^{OA} } \over {\partial p^\prime}},\;\;\;\;\;P^u = {{\partial \varepsilon _q^{AO} } \over {\partial p^\prime}}(17)Ql=∂εqAB∂η, Qu=∂εqBA∂ηQ^l = {{\partial \varepsilon _q^{AB} } \over {\partial \eta }},\;\;\;\;\;Q^u = {{\partial \varepsilon _q^{BA} } \over {\partial \eta }}

To effectively describe the results of experiments and to determine the functions M, N, P, Q, the idea of a ‘common’ curve is adopted, see [40]. To obtain the functions M, P, a series of experiments for different slopes of the stress paths OA and AO are conducted by the authors. The experiments are performed for two cases, namely, for contractive and dilative ‘Skarpa’ sands. Both these states of soil are defined by the position with respect to the steady-state line, see [40].

3.3.1

Functions for contractive sand

To investigate the behaviour of contractive soil and to determine volumetric and deviatoric strains, experiments in a triaxial apparatus are performed for a set of different stress paths for η = const. Each test consists of two phases. In the first stage, the soil sample is subjected to anisotropic compression through the confining pressure, so, in this case, the mean effective stress p′ increases, and the specimen is subjected to pure spherical loading (OA path). In the second phase, a return path is applied, with the same value η = const and a decreasing p′ (AO path). Therefore, the sample is subjected to pure spherical unloading. In this way, the soil behaviour during single loading–unloading loop is tested. Thus a series of experiments are performed to analytically determine the functions M^l, P^l, M^u, P^u describing the behaviour of volumetric and deviatoric strains.

Figure 4 shows the stress paths tested. Figure 5 shows the development of volumetric strain corresponding to contractive soil behaviour for each stress path in Figure 4 for spherical loading paths. Figure 6 shows the development of volumetric strain for spherical unloading paths. Table 1 presents a list of experiments performed on contractive sand samples and their results for a single spherical loading–unloading loop. ID0I_D^0 denotes the initial index of density of each phase. A_v, AvuA_v^u are coefficients needed to determine the functions M^l and M^u. Each test is numbered and corresponds to a selected stress path η = const. The least squares method is used for all approximations of baseline tests in this article.

Figure 4

Figure 5

Figure 6

Table 1

List of experiments performed on contractive sand samples and the corresponding results for a single spherical loading–unloading loop.

Test	η = const	dp′ > 0		dp′ < 0
		ID0I_D^0	AvuA_v^u	ID0I_D^0	AvuA_v^u
LH75d	η = 0.64	0.077	10.75	0.239	5.10
LH77d	η = 0.97	0.146	8.82	0.280	4.15
LH79d	η = 0.44	0.23	8.69	0.356	5.07
LH80d	η = 1.12	0.113	9.02	0.248	4.17
LH81d	η = 1.27	0.082	10.26	0.204	4.49
LH81dd	η = 1.22	0.057	10.42	0.201	4.45
LH82d	η = 0.23	0.094	8.3	0.181	4.54
LH84d	η = 0	0.056	8.4	0.178	4.72

Figure 5 shows the experimental results and an analytical fit. An analytical approximation of these curves can be described using the following formula: (18)εv=Avp′=9.33p′\varepsilon _v = A_v \sqrt {p^\prime} = 9.33\sqrt {p^\prime} According to definition (14), (19)Ml=∂εvOA∂p′=Av2p′=9.332p′M^l = {{\partial \varepsilon _v^{OA} } \over {\partial p^\prime}} = {{A_v } \over {2\sqrt {p^\prime} }} = {{9.33} \over {2\sqrt {p^\prime} }}

Numbers 1 and 3 in Figure 5 denote a limit and an analytical approximation of experimental results, and number 2 denotes an averaged approximation of all approximations, whose equations for each test are given in Table 1. The limit values of the coefficient A_v are 8.3 and 10.75, whereas the average value is 9.33. In order to formulate the function M^l, it was assumed that its character does not depend on the value of η.

Similarly, analytical approximations corresponding to spherical unloading are presented in Figure 6 and described in detail in Table 1. The boundary values of the coefficient AvuA_v^u are 5.1 and 4.15, whereas the average value is 4.59. Therefore, when creating the function M^u, it was assumed that its character does not depend on the value of η.

The averaged approximation is (20)εv=Avup′=4.59p′\varepsilon _v = A_v^u \sqrt {p^\prime} = 4.59\sqrt {p^\prime} According to definition (14), (21)Mu=∂εvAO∂p′=Avu2p′=4.592p′M^u = {{\partial \varepsilon _v^{AO} } \over {\partial p^\prime}} = {{A_v^u } \over {2\sqrt {p^\prime} }} = {{4.59} \over {2\sqrt {p^\prime} }}

To determine the functions P^l and P^u, it is necessary to introduce a ‘common’ curve to approximate experimental results describing the development of deviatoric strains. To determine the function P^l for the set of experimental results presented in Figure 7, the following steps need to be performed. An analytical approximation of the experimental data can be defined using the following formula: (22)εq=Aqp′\varepsilon _q = A_q \sqrt {p^\prime}

Figure 7

Figure 7 shows that each stress–strain relation can be approximated with very good accuracy. The data presented correspond to stress paths given in Figure 6. In Figure 7, on the right side of each set of experimental data, the corresponding η value is shown. Table 2 presents the values of the parameter A_q corresponding to the selected stress path and the initial relative density of soil samples.

Table 2

Series of experiments performed to calibrate the model for contractive sand.

Test	η = const	First phase: dp′ > 0
		ID0I_D^0	Aq
LH75d	η = 0.64	0.077	2.94
LH77d	η = 0.97	0.146	6
LH79d	η = 0.44	0.23	1.3
LH80d	η = 1.12	0.113	9.76
LH81d	η = 1.27	0.082	28.3
LH81dd	η = 1.22	0.057	20.2
LH82d	η = 0.23	0.094	0.21
LH84d	η = 0	0.056	−2.63

The coefficient A_q (similar to the coefficient A_v for volumetric strain) represents the magnitude of deviatoric strain for a given stress path η = const. The higher the value of A_q, the greater the value of deviatoric strain for the corresponding anisotropic stress path and for a given value of the mean effective stress p′. The greatest value of A_q = 28.3 corresponds to the stress path η = 1.27, whereas the smallest A_q = −2.63 corresponds to the stress path η = 0. Figure 7 and Table 2 present that the value of the coefficient A_q increases monotonically with an increasing value of η.

Figure 8 presents the function A_q (η) obtained from experimental data and its analytical approximation. A_q (η) is given using the following formula: (23)Aq=a+bη4=−0.68+9.7η4A_q = a + b\eta ^4 = - 0.68 + 9.74\eta ^4 On the basis of Eqs. (22) and (23), the final averaged approximation has the following form: (24)εq=Aqp′=(a+bη4)p′=(−0.68+9.74η4)p′\matrix{ {\varepsilon _q } \hfill & { = A_q \sqrt {p^\prime} = \left( {a + b\eta ^4 } \right)\sqrt {p^\prime} } \hfill \cr {} \hfill & { = \left( { - 0.68 + 9.74\eta ^4 } \right)\sqrt {p^\prime} } \hfill \cr }

Figure 8

The curve described using Eq. (24) is a ‘common’ curve containing information about the behaviour of deviatoric strain for every anisotropic consolidation test. Despite the complex behaviour of soil in this case, the approximation contains only two coefficients: a and b.

According to definition (16), the function P^l is defined as (25)pl=(a+bη4)2p′=(−0.68+9.74η4)2p′p^l = {{\left( {a + b\eta ^4 } \right)} \over {2\sqrt {p^\prime} }} = {{\left( { - 0.68 + 9.74\eta ^4 } \right)} \over {2\sqrt {p^\prime} }}

To determine the function P^u, it is necessary to analyse the experimental results given in Figure 9, where η = const. The case considered is spherical unloading. The function P^u describes deviatoric strains, so according to (17), for such a case, (26)Pu=0P^u = 0 It indicates that there is no change in deviatoric strain during spherical unloading for stress η = const.

Figure 9

The results of experiments needed to determine the functions N^l, N^u, Q^l, Q^u are presented in [40]. However, the form of these functions is different from that in the original model. According to Eqs. (15) and (17), the functions for contractive Skarpa sand are given as follows: (27)Nl=4c1η4p′,c1=2.97;N^l = 4c_1 \eta ^4 \sqrt {p^\prime} ,c_1 = 2.97;(28)Nu=avp′,av=−0.87;N^u = a_v \sqrt {p^\prime} ,a_v = - 0.87;(29)Ql=b1b2exp⁡ (b2η)p′,b1=0.023, b2=6.245\matrix{ {Q^l = b_1 b_2 \exp (b_2 \eta )\sqrt {p^\prime} }, \hfill \cr {b_1 = 0.023,\;\;\;\;b_2 = 6.245} \hfill \cr } (30)Qu=bqp′, bq=0.76.Q^u = b_q \sqrt {p^\prime} ,\;\;\;b_q = 0.76.

3.3.2

Functions for dilative sand

Similar to contractive sand, a series of experiments are performed to determine the ‘basic functions’ for dilative sand. The qualitative character of the functions M^l, M^u, P^l, P^u is the same as in the previous case, but their quantitative nature is different. Table 3 lists and describes all experiments performed on dilative sand samples. ID0I_D^0 denotes the initial relative density of each phase. A_v and AvuA_v^u are coefficients needed to determine the functions and M^u. Each test is numbered and corresponds to a different stress path η = const.

Table 3

List of experiments performed on contractive sand samples and the corresponding results for a single spherical loading–unloading loop.

Test	η = const	First phase:			Second phase: dp′ < 0
		ID0I_D^0	Av	Aq	ID0I_D^0	AvuA_v^u
SH73d	η = 0.96	0.702	3.3	1	0.749	2.59
SH75d	η = 1.23	0.701	3.28	2.22	0.749	3.04
SH76d	η = 0.64	0.669	3.72	0.52	0.722	3.04
SH77d	η = 0.34	0.703	3.52	0.62	0.753	2.9
SH78d	η = 0	0.701	3.45	0	0.750	2.96
SH80d	η = 1.35	0.632	3.27	6.5	0.678	3.35

For the case considered, the average values of the parameters A_v and AvuA_v^u are 3.42 and 2.98, respectively. The following functions M^l, M^u can be presented (similarly as for the contractive sand): (31)Ml=Av2p′=3.422p′M^l = {{A_v } \over {2\sqrt {p^\prime} }} = {{3.42} \over {2\sqrt {p^\prime} }}(32)Mu=Avu2p′=2.982p′M^u = {{A_v^u } \over {2\sqrt {p^\prime} }} = {{2.98} \over {2\sqrt {p^\prime} }} In this case as well, the coefficient A_q is a monotonic function of η, see Figure 10.

Figure 10

Similar to that for the contractive sand, it can be seen that the higher the value of A_q, the greater the value of deviatoric strain for the corresponding anisotropic stress path and a given value of the mean effective stress p′. The greatest value of corresponds to the stress path and the smallest value (A_q = 0) corresponds to the stress path η = 0. Figure 10 and Table 3 show that the value of the coefficient increases monotonically with an increase in the value of η. Figure 10 presents the function obtained from experimental data. An analytical approximation of this function is given by (33)Aq=bη4=1.62η4.A_q = b\eta ^4 = 1.62\eta ^4 . On the basis of Eqs. (22) and (23), the final averaged approximation is written as follows: (34)εq=Aqp′=bη4p′=1.62η4p′.\varepsilon _q = A_q \sqrt {p^\prime} = b\eta ^4 \sqrt {p^\prime} = 1.62\eta ^4 \sqrt {p^\prime} .

The curve defined using Eq. (34) is a ‘common’ curve containing information about the development of deviatoric strain for every anisotropic consolidation test. Despite the complex response of soil in this case, the approximation function has only one coefficient b.

According to definition (16), the function P^l is (35)Pl=bη42p′=1.62η42p′.P^l = {{b\eta ^4 } \over {2\sqrt {p^\prime} }} = {{1.62\eta ^4 } \over {2\sqrt {p^\prime} }}. The experimental results show that deviatoric strain does not develop in dilative sand for stress paths η = const when dp′ < 0. Therefore, (36)Pu=0.P^u = 0. The experiments needed to determine the functions N^u, Q^l, Q^u for dilative sand are presented (similar to those for the contractive sand) in [40]. However, the form of these functions is different from that original model. According to Eqs. (15) and (17), the functions for dilative Skarpa sand are given as follow: (37)Nl=4c1η4P′, c1=2.97;N^l = 4c_1 \eta ^4 \sqrt {P^\prime} ,\;\;\;c_1 = 2.97;(38)Nu=avp′, av=−0.39;N^u = a_v \sqrt {p^\prime} ,\;\;\;a_v = - 0.39;(39)Ql=b1b2exp⁡ (b2η)p′,b1=0.00035, b2=6.648;\matrix{ {Q^l = b_1 b_2 {\it exp} (b_2 \eta )\sqrt {p^\prime} }, \hfill \cr {b_1 = 0.00035,\;\;\;b_2 = 6.648;} \hfill \cr } (40)Qu=bqp′, bq=0.4.Q^u = b_q \sqrt {p^\prime} ,\;\;\;b_q = 0.4. An approximation of the experimental results for volumetric strain needed to determine the function N^l in [40] is given as follows: (41)εvp′={a1η2+a2η for 0≤η≤η′(a3η2+a4η+a5) (exp⁡ (a6η)−1)for η′≤η≤η″.{{\varepsilon _v } \over {\sqrt {p^\prime} }} = \left\{ {\matrix{ {a_1 \eta ^2 + a_2 \eta \;{\rm{for}}\;0 \le \eta \le \eta ^\prime} \hfill \cr {\left( {a_3 \eta ^2 + a_4 \eta + a_5 } \right)\;(\exp \;(a_6 \eta ) - 1)} \hfill \cr {{\rm for}\;\eta ^\prime \le \eta \le \eta ''.} \hfill \cr } } \right. However, approximation (41) is too complex and contains six coefficients. Therefore, a simplified function proposed by the present authors has the following form: (42)εvp′=a1exp⁡ (a2η2)+a3η3.{{\varepsilon _v } \over {\sqrt {p^\prime} }} = a_1 {\it exp} \;\left( {a_2 \eta ^2 } \right) + a_3 \eta ^3 .

Figure 11 presents the experimental results for pure shearing of dilative sand in the form of a single ‘common’ curve. The function given by Eq. (42) has a simpler form than Eq. (41), as it is characterised by only three coefficients.

Figure 11

It should be noted that the behaviour of dilative sand during pure shearing is different from that of contractive sand. For contractive sand, volumetric strain increases in all stages of pure shearing until the sample is destroyed. For dilative sand, volumetric strain increases slightly and then, when the dimensionless stress η reaches η′, it begins to decrease rapidly. Hence, the function N^l for dilative sand is qualitatively different from that for contractive sand (see Eq. 15) and has the following form: (43)Nl=(2a1a2exp⁡(a2η2)+3a3η2)p′,a1=−0.000074,a2=6.39,a3=0.848.\matrix{ {N^l = \left( {2a_1 a_2 {\it exp} \left( {a_2 \eta ^2 } \right) + 3a_3 \eta ^2 } \right)\sqrt {p^\prime} ,} \hfill \cr {a_1 = - 0.000074,a_2 = 6.39,a_3 = 0.848.} \hfill \cr }

Strains corresponding to any stress path are defined using Eqs. (12) and (13). The information regarding the state of soil (contractive or dilative) is required to use the appropriate set of the functions M, N, P, Q. Table 4 presents all modified functions necessary to build Eqs. (12) and (13). Predictions for complex tests (other than basic stress paths) are always a combination of approximations of the designated M, N, P, Q in accordance with Eqs. (12) and (13).

Let us consider any point V₀ in the p′ − η plane and draw two lines p′ = const and η = const running through this point. The plane p′ − η will thus be divided into four regions, see Figure 12. Any stress path starting at the point V₀ can be located in one of the regions delimited by the lines p′ = const or η = const. Depending on the region in which it is located, appropriate functions will be used to build Eqs. (12) and (13). In this simple way, the complex behaviour of soil can be described for every linear loading path, focusing especially on deformations, as shown in Figure 12, on stability. The forms of base functions and appropriate coefficients are presented in Table 4.

Figure 12

To examine the numerical results of the modified model and to show the unstable behaviour of contractive sand under certain conditions, let us consider three experiments performed in a classical triaxial apparatus. Each experiment is conducted in three stages. First, the soil sample is subjected to isotropic consolidation to 200 kPa. Then, each sample is subjected to pure shearing at different levels of deviatoric stress. Finally, deviatoric stress was kept constant, and the mean effective stress is being reduced. Spherical unloading was performed along three different stress paths described using Eq. (52), where α = 0.5, 1, 1.5 and β = 0, see Figure 13.

Figure 13

According to Eqs. (12) and (46), the second-order work for this kind of stress paths has the following form: (58)d2W=dεvdp′=(dp′)22p′(Avu−8c1η4).d^2 W = d\varepsilon _v dp^\prime = {{\left( {dp^\prime} \right)^2 } \over {2\sqrt {p^\prime} }}\left( {A_v^u - 8c_1 \eta ^4 } \right).

Generally, the sign of the function describing the second-order work is the same as the sign of the function describing an increase in volumetric strains. The volumetric strains are given by (59)εv=εv0+Avup′+87c1q4p′−3.5.\varepsilon _v = \varepsilon _v^0 + A_v^u \sqrt {p^\prime} + {8 \over 7}c_1 q^4 {p^\prime}^{ - 3.5} . The experimental results and numerical calculations for average values of the parameters (see Table 4) are presented in Figure 14.

Figure 14

The qualitative agreement with the experimental data is achieved, and the quantitative agreement is sufficient. During the first stage of spherical unloading, dilation takes place and the volume strain decreases. After it reaches the minimum, compaction occurs. The soil loses its stability at the moment when the sign of the volumetric strain increment changes.

The value of η at the moment whenchanges its sign is given as follows: (60)η=Avu8c14.\eta = \root 4 \of {{{A_v^u } \over {8c_1 }}} . For “Skarpa” sand, η = 0.66 and, according to Hill’s criterion, when this value is reached, the soil loses stability.

Saturated soil exhibits traits of a macroscopically solid body. When pore water drainage is prevented, soil can be liquefied, and then it exhibits liquid characteristics. Soil liquefaction is an undesirable phenomenon, often causing catastrophic damage (e.g. buildings sinking in the ground or pipelines being dislodged from the seabed). The phenomenon of liquefaction has been investigated in many monographs, see for example [12, 13, 15, 16, 23, 26, 29, 47].

Undrained conditions can also be investigated in the laboratory, for example, in a triaxial apparatus. First, the sample is isotropically loaded, and then it is sheared under undrained conditions. Only a contractive soil can liquefy. The proposed incremental model can also be used to describe the behaviour of saturated sand under undrained conditions, bearing in mind that liquefied soil behaves like an incompressible body, so it is characterised by a null volumetric strain: (61)dεv=Mudp′+Nldη=0.d\varepsilon _v = M^u dp^\prime + N^l d\eta = 0. After substituting the appropriate functions from Table 4, one obtains (62)dεv=Avu2p′dp′+4c1η3p′dη=0.d\varepsilon _v = {{A_v^u } \over {2\sqrt {p^\prime} }}dp^\prime + 4c_1 \eta ^3 \sqrt {p^\prime} d\eta = 0. Solving Eq. (62) with the initial condition for η = 0, one obtains the following solution: (63)q=p′Avu2c1ln⁡p0′p′4.q = p^\prime\root 4 \of {{{A_v^u } \over {2c_1 }}\ln {{p^\prime_0} \over {p^\prime}}} .Equation (63) describes the effective stress path leading to liquefaction.

The second-order work defined in terms of stress invariants [see 32] is the following: (64)d2W=dεvdp′+dεqdq,d^2 W = d\varepsilon _v dp^\prime + d\varepsilon _q dq, and under undrained conditions (dε_v = 0): (65)d2W=dεqdq=b1b2exp⁡ (b2η)1p′(dq−ηdp′)dq.d^2 W = d\varepsilon _q dq = b_1 b_2 {\it exp} \;(b_2 \eta ){1 \over {p^\prime}}\left( {dq - \eta dp^\prime} \right)dq. In the first stage of the liquefaction process, the stress deviator increases to the maximum and then begins to decrease. The soil loses its stability when the sign of the deviator changes. Figure 15 shows the liquefaction process: the experimental, see [48], and numerical results for an initial mean effective stress of 200 kPa. The continuous line represents the numerical results (see Eq. 62) for the average values of the coefficients (Table 4), and the dotted line represents the numerical results for slightly modified values of the parameters: Avu=6.59A_v^u = 6.59 , c₁ = 1.75. The empty dots indicate experimental data. The qualitative agreement with the experimental data is achieved. The quantitative agreement for the slightly altered values of the parameters is good. Changes in all parameter values in the model proposed correspond to the changes described in [33] and their consequences.

Figure 15

The signs of the coefficients b₁, b₂ and the increment dη are positive, so expression (65) changes sign when dq = 0. Figure 15 shows that q = const only when function (63) reaches an extreme: (66)qmax⁡=p′Avu8c14.q_{\max } = p^\prime\root 4 \of {{{A_v^u } \over {8c_1 }}} . Therefore, the instability line η_IL is given by the equation (67)ηIL=Avu8c14,\eta _{IL} = \root 4 \of {{{A_v^u } \over {8c_1 }}} , and for the average values of the coefficients, η_IL = 0.66.

The behaviour of sheared dilative sand under undrained conditions is qualitatively different from that of contractive sand. The mean effective stress p′ decreases to the minimum and then starts to increase. The main general description of stresses is given by Eq. (61), and in the case of dilative sand, (68)dεv=Avu2p′dp′+(2a1a2ηexp⁡(a2η2)+3a3η2)p′dη=0.\matrix{ {d\varepsilon _v } \hfill & { = {{A_v^u } \over {2\sqrt {p^\prime} }}dp^\prime} \hfill \cr {} \hfill & { +\; \left( {2a_1 a_2 \eta {\it exp} \left( {a_2 \eta ^2 } \right) + 3a_3 \eta ^2 } \right)\sqrt {p^\prime} d\eta = 0.} \hfill \cr } The solution of the above equation at the initial condition η(p′₀) = 0 is the following: (69)a1ηexp⁡(a2η2)+a3η3=Avu2ln⁡(p′0p′).a_1 \eta {\it exp} \left( {a_2 \eta ^2 } \right) + a_3 \eta ^3 = {{A_v^u } \over 2} {\it ln} \left( {{{p^\prime_0 } \over {p^\prime}}} \right).

Figure 16 shows the experimental and numerical results for an initial mean effective stress of 200 kPa. The continuous line represents the numerical results (see Eq. 69) for the average values of the coefficients (Table 4), whereas the dotted line represents the numerical results for slightly different values of the parameters a₂ = 5.8, a₃ = 0.6. The empty dots indicate experimental data. A satisfactory qualitative agreement with the experimental data is achieved, and the quantitative agreement for the slightly altered values of the parameters is very good.

Figure 16

The equation describing the phase transformation line η_PTL is given as follows: (70)ηPTL=Avu3a33.\eta _{PTL} = \root 3 \of {{{A_v^u } \over {3a_3 }}} . The sign of the second-order work d²W, see Eq. (64), is always positive, so, according to Hill’s criterion, the soil does not lose stability in this case.