Complexity
Volume 2018 (2018), Article ID 2829873, 10 pages
https://doi.org/10.1155/2018/2829873

Research Article

Rigorous Solution of Slopes’ Stability considering Hydrostatic Pressure

Chengchao Li, Pengming Jiang, and Aizhao Zhou

Department of Civil and Architecture Engineering, Jiangsu University of Science and Technology, Zhenjiang, Jiangsu 212003, China

Correspondence should be addressed to Chengchao Li; justlcc2015@163.com

Received 10 December 2017; Revised 23 April 2018; Accepted 8 May 2018; Published 12 June 2018

Academic Editor: Rafał Burdzik

Copyright © 2018 Chengchao Li et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract.

According to characteristics of soils in failure, a sliding mechanism of slopes in limit state is divided into five parts, for building a slip line field satisfying all possible boundary conditions. An algorithm is built to obtain the rigorous solution approaching upper and lower bound values simultaneously, which satisfies the static boundary and the kinematical boundary based on the slip line field, while stress discontinuity line and velocity discontinuity line are key points. This algorithm is copared with the Spencer method to prove its feasibility with a special example. The variation of rigorous solution, including an ultimate load and a sliding belt the rigid body sliding along rather than a single slip surface for friction-type soils, is achieved considering hydrostatic pressure with soil parameters changing.

1. Introduction

The stability of slopes has been regarded as a classic and difficult problem for engineers because of less boundary constraints, compared with the earth pressure of retaining wall and the bearing capacity of foundation. Over the past few years, many investigators have evaluated slope stability, thereby developing many methods for meeting engineering requirements, such as the limit equilibrium method (LEM) [1, 2], the finite element method (FEM) [3, 4], and the limit analysis method (LAM) [5–7]. The LEM captures the static equilibrium of rigid blocks on a particular slip surface, while not considering the plastic deformation of soils. The equilibrium equation accounts for the whole slice but not guarantees each point in soils. The strength reduction method (SRM) [8–10] is the main finite element slope stability method currently employed, by which stress field and displacement of soils in slopes can be calculated with an elastic-plastic constitutive model to get the safety factor. Although the displacement mutation, the numerical calculation of nonconvergence, and other criteria can determine the slope instability, the slope displacement calculated by the SRM is not the actual displacement of soils and the safety factor required is an approximation. For the slope stability, it is sometimes not necessary to get the variation of stress field and strain field, only to get the ultimate load. Based on the extremum principal [11], the lower bound (LB) [12, 13] solution can be got by static analysis for limit equilibrium problems and the upper bound (UB) [14–17] solution can be got by dynamic analysis. If the lower bound solution satisfies all the kinematic conditions or the upper bound solution satisfies all the static conditions, the solution will be the rigorous one. Compared with obtaining the safety factor, the solution calculated by the upper and lower bound theorems is closer to the real condition because the ultimate load approaches the upper bound solution and the lower bound solution at the same time, satisfying all possible boundary conditions.

When a slope is on the verge of collapse, for type, a sliding belt is emerged within the slope due to the friction between soils. And the sliding belt is not a single slip surface, but a thin shear zone. This paper builds the slip line field satisfying the static boundary condition and the velocity boundary condition according to the characteristics of the stress discontinuity line and the velocity discontinuity line and compiles an algorithm to gain the distribution of the sliding belt and the ultimate capacity of the slope considering the hydrostatic pressure. It also indicates the variation of the ultimate bearing capacity and the sliding belt with different parameters. The application of this algorithm is proved by comparing with the Spencer’s method.

2. Slip Line Field

To analyze the slope stability based on the LAM, the slip line field is constructed according to the stress boundary conditions; one of which is a noncharacteristic line stress boundary with normal stress and shear normal stress, and the other one is the interface of a rigid region and a plastic region or a stress discontinuity surface [18]. Under the limit state, only the noncharacteristic line boundary can be used to construct the slip line field, while the other is unknown. Basic coordinate and noncharacteristic line stress boundary is shown in Figure 1. is the angle between the direction of the maximum principal stress , and the -axis and is the noncharacteristic line stress boundary. is the angle between the tangential direction of the boundary and the -axis, so the angle between and is . According to and , two Mohr’s circles tangent to the Coulomb failure line can be drawn (Figure 2). When soils are in the extreme state, , , and on the boundary are expressed, respectively, as follows:

where , , and are normal stress, shear stress, and tangent normal stress on the boundary , respectively. is average stress; is radius of Mohr’s circle; is soil cohesion; and is internal friction angle. is cohesive internal stress, which is given by: .

Figure 1: Basic coordinate system and noncharacteristic line stress boundary.

Figure 2: Mohr’s circles on noncharacteristic line stress boundary.

Because is a point on the Mohr’s circle, can be rewritten as follows:

According to (1.4) and (2):

Solving the unary quadratic equation of :

When , according to (1.3):

When , (1.3) should be rewritten as follows:

where and are the maximum principal direction angles.

According to the above analysis, the analytic expressions of noncharacteristic line stress boundary conditions can be uniformly written as follows:

With on the noncharacteristic line stress boundary, (6) can be used to draw the stress line field near the stress boundary:

where and are the angles between stress characteristic lines through the point and the -axis; is the angle between the stress characteristic line and the direction of major principal stress and given as follows:

In this case, stress boundary conditions are achieved. In order to facilitate the following calculations, new variables and are used instead of , , and (Figure 3). is the effective stress, and is the direction angle of the maximum principal stress, written as the (9). In addition, and of each point in slip line field can be calculated by solving the basic boundary value problems. In the theory of hyperbolic partial differential equations, there are three basic boundary value problems [19, 20], of which Cauchy problem is the key issue the others can be solved by. Here, only the Cauchy problem is introduced and the rest of the problems (Riemann problem and mixed boundary value problem) can be referred to the literature [19]. Taking into account of the hydrostatic pressure, the soils above the ground water take unit natural weight and the soils below ground water take unit floating weight .

Figure 3: Effective stress of Mohr’s circle.

As shown in Figure 4(a), is a smooth and continuous stress boundary line, a nonslip line, in which the coordinates of each point, effective stress , and the direction angle are known. By solving the Cauchy problem, the distribution of slip line field in the triangular surrounded by the lines and through the points of and , respectively, could be obtained.

(a)

(b)

Figure 4: Cauchy boundary problem considering the impact of the ground water table.

Such as , , , and on the point 22 can be obtained by the points of 21 and 32, along line :

along line :

when soils are above ground water, . When soils are below ground water, .

If the ground water table is between two points, it should be simplified to calculate (Figure 4(b)). Setting an allowable value , if , the of the points 21 and 22 all take . If , the of the points 21 and 22 all take . If the and , of the points 21 and 22 take and , respectively. In fact, the ground water table is not a straight line. The smaller the grid in the slip line field is, the more similar to a straight line the water line is.

The above two sets of nonlinear equations ((10.1) and (10.2)) need to be solved by an iterative method. It can be assumed that the initial values of point 22 is

Substituting and in (11) into (10.2), the approximations of and can be obtained. Similarly, substituting and in (11) into (10.1), the approximations of and can be obtained. Then, substituting the first approximations of and into (10.1), the second approximations of and can be gained, so repeatedly until the accuracy of the calculation to meet the precision. And , , , and of each point can be achieved in the entire . Finally, according to the (9), , , and in the slip line field could be gained.

3. Stress Discontinuity Line and Velocity Discontinuity Line

When soils reach the limit state, the stress discontinuity field and the velocity discontinuity field will be generated under complex boundaries [21]. With different boundary conditions, the stress discontinuity line and the velocity discontinuity line will also change. Therefore, the distribution and the characteristic of the two discontinuity lines are significant for the slope stability.

3.1. Characteristics of Stress Discontinuity Line

3.1.1. Stress Equation

The stress discontinuity line is actually a thin transition zone, where the stress changes rapidly. It divides the element of the discontinuity line crossing through two plastic regions ① and ②, while only the tangential stress can be interrupted and the normal stress and the shear stress stay the same on both sides (Figure 4). In the plastic regions ① and ②, stresses on the discontinuity line should obey the Mohr-Coulomb yield criterion.

where and are the maximum principal direction angles of the plastic regions ① and ②; is the tangential direction angle of the discontinuity line; and and are the angles between the maximum principal stress direction and tangential direction of stress discontinuity line and are given as follows (Figure 5):

Figure 5: Stress discontinuity line.

The relationship of equivalent stresses on both sides of the discontinuity line is

where and are equivalent stresses on both sides of the discontinuity line; is given as follows:

3.1.2. Geometric Condition

As shown in Figure 6, the angle of the principal stress element is . The angle between the stress discontinuity line and the principal stress plane in the zone ① is , and the angle in the zone ② is . From the geometric relationship shown in Figure 6, it can be obtained:

Figure 6: Stress elements on both sides of stress discontinuity line.

According to the geometric relationship shown in Figure 2, (12) can be rewritten as follows:

The relationship of on both sides of a stress discontinuity line can be obtained eliminating and from (19).

Substituting (18.1) and (18.2) into (20):

The position of the stress discontinuity line can be determined according to (15) and (21) in the slip line field.

3.1.3. Kinematic Equation

Because the stress discontinuity line is emerged by the rigid region reducing to the limit state, the tangential velocity on the discontinuity line remains unchanged [21]:

where and are velocity components in the and directions.

Assuming that and are the velocity components in the direction of the characteristic line, a rigid kinematic equation of the stress discontinuity line will be obtained:

Substituting (22.2) into (22.1):

3.2. Characteristics of Velocity Discontinuity Line

Obeying the associated flow rule, the velocity discontinuity line is a slip line or an envelope of slip line. The boundary between a rigid region and a plastic region is a velocity discontinuity line. So the velocity discontinuity line is the slip line at both ends of the stress discontinuity line, dividing the slope into rigid zone and plastic zone. For the Mohr-Coulomb material, the angle between the velocity direction and the slip line is [22, 23].

4. Algorithm

For slope stability problems, the equilibrium equation and the plastic flow equation should be settled simultaneously to obtain the rigorous solution, which is difficult to achieve in mathematics and only depends on numerical methods. If a statically admissible stress field has been got, the strain rate field and the kinematically admissible velocity field will be obtained by the stress field based on the associated flow rule. If the strain rate field and the velocity field are nothing less than the kinematical field, in this case, the plastic region corresponding to such a statically admissible stress field and a kinematically admissible velocity field must be the sliding belt in the extreme state and the external load must be the ultimate bearing capacity. Based on the upper and lower bound theorems, the numerical algorithm is established to get the distribution of the sliding belt and the ultimate bearing capacity considering the hydrostatic pressure.

When a slope is on the verge of collapse (shear failure), not all the points in sliding mass reach to the yield state, but the mixture of the plastic region and the rigid region is emerged. As shown in Figure 7, the entire slope is divided into five areas. The stress discontinuity line divides the slope foot into the strong plastic region and the weak plastic region, passive, and active, respectively. The boundary between two regions is the velocity discontinuity line, which is the slip line at both ends of the stress discontinuity line. So the stress discontinuity line and the velocity discontinuity line determine the distribution of the plastic region and the rigid region. The stress and deformation in each region are different as well as kinematic features. Only to meet all the possible kinematically admissible conditions, the solution will be exact.

Figure 7: Calculation model.

The calculation process (shown in Figure 8) is consisted of two parts: a statically admissible field and a kinematically admissible field. The specific calculation process is as follows.

Figure 8: Flow chart of rigorous solution for slope stability.

4.1. Static Calculation

(1)The maximum principal stress on the noncharacteristic line stress boundary is 0, and its direction angle is , which is perpendicular to the slope surface. From the boundary , the slip line of the active region is calculated in a very dense grid by solving the Cauchy problem.(2)To get the stress discontinuity line In the active region, the specific approach is to firstly select the point in the grid. The stress discontinuity line is drawn by (15) and (21). A group of the maixmum principal stress direction angles of the stress discontinuity line in the active region is assumed as . Then, the direction angle , the tangential direction angle of the discontinuity line , and the coordinate of each point in the line in the passive region can be obtained by the (13) and (15). And the slip line field of the region can be obtained by solving the Cauchy problem.(3)According to the direction angle of the maximum principal stress artificially assumed and the coordinate of each point in the line , the effective stress is obtained. From the lines and , the slip line field of region can be got by solving the Riemann problem. On the basis of maximum principal stress direction angle , the slip line field of region can be got by solving the mixed boundary value problem.(4)The numerical integration is carried out along the slip lines and to get the forces in the and directions and the moments of each force to the point . Then the vertical stress on the boundary is calculated.(5)At this point, the static calculation is over. Then, check whether the force and the moment of the rigid block satisfy the equilibrium condition and whether the points within the rigid block satisfy the yield condition. If not satisfied, the positon of point is modified to recalculate. The whole calculation process is shown in Figure 9.

Figure 9: Flow chart of static calculation.

4.2. Velocity Field Analysis

According to the associated flow rule, the dynamic analysis is carried out on the basis of the slip line field got by the static calculation.

(1)The absolute velocity of point on slope top is assumed as 1 m/s. The velocity component and in the slip line ABCD is obtained according to the velocity characteristic equation.

(2)The velocity field of region II can be obtained by solving the mixed boundary value problem for the condition that the tangential velocity on the discontinuity line remains unchanged (23). From the slip lines and , the velocity field in regions III and V and the velocity of each point in the line can be got by solving the Riemann problem. Then, the velocity and of each point in the line in the passive region II is translated into the velocity and in the active region I. Finally, the velocity field in the region I and the velocity in line are obtained by solving the Cauchy problem from the line .(3)At this point, the velocity analysis is over. Then, check whether the velocity in the lines and satisfy the kinematically rigid condition of the sliding block IV. If not satisfied, the input values of on the slip line and the position of the point are modified to recalculate all.

5. Example Verification

Geometric parameters and material properties are shown in Figure 10. The Fortran 95 compiler is utilized to compute the calculation program for the rigorous solution satisfying the static equilibrium and the kinematical requirement.

Figure 10: Parameters of a calculation model.

When a slope reaches the limit state, the ultimate bearing capacity is composed of two parts (Figure 11). While the load on the rigid region reaches 77.54 kPa, the other on the plastic region changes from 102.05 kPa to 112.45 kPa, considering the hydrostatic pressure. The rigid body rotates around point , slipping along the sliding belt. The velocity of each point in the velocity discontinuity line increases gradually from top to foot (Figure 12). The angle between a velocity discontinuity line and a slip line is . The position of points and calculated determines the general distribution of the sliding belt. In order to facilitate the appearance, the figures (Figures 11 and 12) show only part of the slip lines and the actual existence is hundreds of groups.

Figure 11: Distribution of the sliding belt and ultimate loads under the limit state.

Figure 12: Velocity field of the slope under the limit state.

The limit state can be seen as a status with safety factor equal to 1.0, and the external load at this time is the ultimate bearing capacity. When considering the hydrostatic pressure, the soils below ground water take floating weight resulting in the slide force falling and the slide resistance remains unchanged. From the Table 1, the ultimate load has become lager compared with the load without considering the ground water. As the depth increases causing the slide force to rise, the ultimate loads in both regions decrease, of which the one on the rigid region is influenced greater than the other by the ground water. While changing the depth, the distribution of sliding belt hardly alters. As the water level rises, the width of the sliding belt increases just a little.

Table 1: Variation of the ultimate loads and the distribution of the sliding belt with the depth .

Ground water depth (m)
Ultimate load (kPa)
(m)
(m)

Rigid region
Plastic region (average)

0
77.80
107.69
2.32
0.87

0.5
77.54
107.25
2.32
0.87

1.0
76.69
106.74
2.31
0.87

1.5
73.34
105.39
2.31
0.86

2.0
70.78
104.38
2.30
0.86

2.5
64.09
101.67
2.29
0.86

3.0
59.83
99.93
2.28
0.85

4.0
43.89
93.51
2.27
0.85

No ground water
33.85
86.98
2.25
0.85

From Table 2, when weight of soils increases gradually, the ultimate load on the rigid region decreases corresponding to the other on plastic region increasing. The distribution of sliding belt undergoes great changes due to the reducing width of . For kN/m³, the rigid region in the slope disappears under the limit state and the whole slip field becomes plastic. Changes of cohesive only bring about the increase of ultimate loads on both regions, and the sliding belt has a little change. Internal friction angle has a great impact on the load on the plastic region and the sliding belt. With friction angle rising, the load on plastic region increases sharply and the width of sling belt fills out gradually. When the angle reduces to 0 resulting the plastic region disappeared, the sliding belt transforms into a signal slip surface and the bearing capacity of a slope depends on the cohesive. At this time, the ultimate load calculated by the algorithm is considered as the external load on the slope top and this status can be calcualted by the Spencer method to get the slip surface and the safety factor. As shown in Figure 13, the critical slip surface obtained by this paper’s method is basically consistent with the Spencer method. The safety factor calculated by the Spencer method is 0.99, when the ultimate load is 65.53 kPa. Generally speaking, the comparison proves the feasibility of the algorithm proposed by this paper.

Table 2: Variation of the ultimate loads and the distribution of the sliding belt with different soil parameters.

(kPa), (°), (kN/m³), (m)
Ultimate load (kPa)
(m)
(m)

Rigid region
Plastic region

83.23
84.86
2.25
1.14

79.99
93.97
2.29
0.99

77.54
107.25
2.32
0.87

71.91
125.62
2.34
0.77

53.52
78.61
2.33
0.85

77.54
107.25
2.32
0.87

110.57
191.98
2.31
0.89

198.88
331.25
2.30
0.92

54.22
75.89
2.39
0.79

77.54
107.25
2.32
0.87

108.52
198.55
2.28
1.03

163.35
387.97
2.26
1.12

Figure 13: Comparison between this paper’s method and the Spencer method.

With the hydrostatic pressure, the difference between the two loads on both regions becomes smaller compared with the state of no ground water. The vairable amplitude of the ultimate loads and the sliding belt with the different parameters is also diminished.

6. Discussion

During the process of searching the rigorous solution of slopes stability, it is necessary to find a slip line field satisfying all possible static and kinematic conditions. Although the calculation is so difficult that the iterative computation takes thousands of times, the author used the Fortran 95 compiler to compile the calculation program, which only takes less than one minute to obtain the rigorous solution and greatly simplified the calculation. Since the grid does not meet the requirements of infinite subdivision, it is not each point on the boundary but the intersection of the grid and the boundary that satisfies the stress boundary condition and the velocity boundary condition.

With the calculation model, the program is compiled based on the algorithm developed to obtain the ultimate bearing capacity and the distribution of the sliding belt. The influence of the ground water table on the ultimate bearing capacity and the sliding belt of a slope is discussed by an example. The results show that the ultimate load increases considering the hydrostatic pressure, while decreasing with the descending water level. However, the sliding belt remains constant. The soil weight and the internal friction angle have a tremendous impact on the belt, especially for the later. The ultimate load is influenced by the three parameters, of which the friction angle mainly controlled the load on the plastic region. When the internal friction angle reduces to zero, the sliding belt will translate into a traditional slip surface. Considering the ultimate load in this status as the external load on the slope top, the critical slip surface obtained by the Spencer method is basically consistent with the one got by this paper and the safety factor obtained by the Spencer method is 0.99 very close to 1.0. So the feasibility of this algorithm is verified by the specific example. Due to the hydrostatic pressure, the difference of loads on the two regions and the variation of sliding belt become smaller compared with not considering the ground water.

On the basis of this study, the next research focus is to add the seepage field into the slip field, which are not completely coincide. The authors hope that others can be motivated to consider adding the penetration force into the equilibrium equation of slip line to gain the numerical solution and the influence of seepage on the sliding belt and the ultimate load.

7. Conclusion

Based on the upper and lower bound theorems, the rigorous solution satisfying all static boundaries and kinematic boundaries is obtained, which approaches lower and upper bound solutions at the same time and is considered as the actual solution. An algorithm for calculating the ultimate bearing capacity and the distribution of sliding belt has been developed, and the computer program has been accomplished. The application of this algorithm is verified by a specific example (Figure 13).

Considering the hydrostatic pressure, the variation of the ultimate load and the sliding belt under different parameters has been gained. The ultimate load decreases with the rising water level and become the minimum with no ground water. The sliding belt is controlled by the internal friction angle and the soil weight, while cohesive has a obvious impact on the ultimate loads. The difference of ultimate loads on the two regions and the variation of sliding belt become smaller because of the soils taking the floating weight.

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Acknowledgments

This study was financially supported by the National Natural Science Foundation of China (Grant no. 51579119).

References

X.-L. Chen, C.-G. Liu, Z.-F. Chang, and Q. Zhou, “The relationship between the slope angle and the landslide size derived from limit equilibrium simulations,” Geomorphology, vol. 253, pp. 547–550, 2016.

D. Dong-ping, L. Liang, W. Jian-feng, and Z. Lian-heng, “Limit equilibrium method for rock slope stability analysis by using the Generalized Hoek–Brown criterion,” International Journal of Rock Mechanics and Mining Sciences, vol. 89, pp. 176–184, 2016.

H. Zheng, D. F. Liu, and C. G. Li, “Slope stability analysis based on elasto-plastic finite element method,” International Journal for Numerical Methods in Engineering, vol. 64, no. 14, pp. 1871–1888, 2005.

S.-H. Jiang, D.-Q. Li, L.-M. Zhang, and C.-B. Zhou, “Slope reliability analysis considering spatially variable shear strength parameters using a non-intrusive stochastic finite element method,” Engineering Geology, vol. 168, pp. 120–128, 2014.

D.-Q. Li, T. Xiao, Z.-J. Cao, K.-K. Phoon, and C.-B. Zhou, “Efficient and consistent reliability analysis of soil slope stability using both limit equilibrium analysis and finite element analysis,” Applied Mathematical Modelling, vol. 40, no. 9-10, pp. 5216–5229, 2016.

W. F. Chen and X. L. Liu, Limit Analysis in Soil Mechanics, Elsevier, Amsterdam, Netherlands, 1990.

F.-T. Liu, Y.-H. Fan, and J.-H. Yin, “The use of QP-free algorithm in the limit analysis of slope stability,” Journal of Computational and Applied Mathematics, vol. 235, no. 13, pp. 3889–3897, 2011.

W. Fu and Y. Liao, “Non-linear shear strength reduction technique in slope stability calculation,” Computers and Geotechnics, vol. 37, no. 3, pp. 288–298, 2010.

Y. Tu, X. Liu, Z. Zhong, and Y. Li, “New criteria for defining slope failure using the strength reduction method,” Engineering Geology, vol. 212, pp. 63–71, 2016.

Q. Xu, H. Yin, X. Cao, and Z. Li, “A temperature-driven strength reduction method for slope stability analysis,” Mechanics Research Communications, vol. 36, no. 2, pp. 224–231, 2009.

D. C. Drucker and W. Prager, “Soil mechanics and plastic analysis or limit design,” Quarterly of Applied Mathematics, vol. 10, no. 2, pp. 157–165, 1952.

D. Li and Y. Cheng, “Lower bound limit analysis using nonlinear failure criteria,” Procedia Earth and Planetary Science, vol. 5, pp. 170–174, 2012.

J. Zhou, Q. Chen, and J. Wang, “Rigid block based lower bound limit analysis method for stability analysis of fractured rock mass considering rock bridge effects,” Computers and Geotechnics, vol. 86, pp. 173–180, 2017.

S. W. So and W. A. Strauss, “Upper bound on the slope of steady water waves with small adverse vorticity,” Journal of Differential Equations, vol. 264, no. 6, pp. 4136–4151, 2018.

N. Deusdado, A. N. Antão, M. V. d. Silva, and N. Guerra, “Application of the upper and lower-bound theorems to three-dimensional stability of slopes,” Procedia Engineering, vol. 143, pp. 674–681, 2016.

I. B. Donald and Z. Chen, “Slope stability analysis by the upper bound approach: fundamentals and methods,” Canadian Geotechnical Journal, vol. 34, no. 6, pp. 853–862, 1997.

K. Krabbenhoft, A. V. Lyamin, M. Hjiaj, and S. W. Sloan, “A new discontinuous upper bound limit analysis formulation,” International Journal for Numerical Methods in Engineering, vol. 63, no. 7, pp. 1069–1088, 2005.

M.-c. Yang and Y.-r. Zhen, “The analysis expression of non-characteristic line stress boundary condition in the slip line field theory for geotechnical materials,” Rock and Soil Mechanics, vol. 22, no. 4, pp. 395–398, 2011.

N. Louat, “Solution of boundary value problems in plane strain,” Nature, vol. 196, no. 4859, pp. 1081-1082, 1962.

S. Momani, S. Abuasad, and Z. Odibat, “Variational iteration method for solving nonlinear boundary value problems,” Applied Mathematics and Computation, vol. 183, no. 2, pp. 1351–1358, 2006.

R. T. Shield, “Stress and velocity fields in soil mechanics,” Journal of Mathematics and Physics, vol. 33, no. 1–4, pp. 144–156, 1954.

O. C. Zienkiewicz, C. Humpheson, and R. W. Lewis, “Associated and non-associated visco-plasticity and plasticity in soil mechanics,” Geotechnique, vol. 25, no. 4, pp. 671–689, 1975.

H. Liu, “Unified sand modeling using associated or non-associated flow rule,” Mechanics Research Communications, vol. 50, pp. 63–70, 2013.