The bearing capacity factor Nγ of strip footings on c–ϕ–γ soil using the method of characteristics

Background The method of characteristics (also called as the slip-line method) is used to calculate the bearing capacity of strip footings on ponderable soil. The soil is assumed to be a rigid plastic that conforms to the Mohr–Coulomb criterion. The solution procedures proposed in this paper is implemented using a finite difference method and suitable for both smooth and rough footings. By accounting for the influence of the cohesion c, the friction angle ϕ and the unit weight γ of the soil in one failure mechanism, the solution can strictly satisfy the required boundary conditions. Results The numerical solution of Nγ are consistent with published complete solutions based on cohesionless soil with no surcharge load. The relationship of Nγ between smooth and rough foundations is discussed which indicates that the value of Nγ for a smooth footing is only half or more of that for a rough footing. The influence of λ (λ = (q + ccot ϕ)/γB) on Nγ is studied. Finally, a curve-fitting formula that simultaneously considers both ϕ and λ is proposed and is used to produce a series of Nγ versus λ curves. Conclusions The surcharge ratio λ and roughness of the footing base both have significant impacts on Nγ. The formula for the bearing capacity on c–ϕ–γ soil can be still expressed by Terzaghi’s equation except that the bearing capacity factor Nγ depends on the surcharge ratio λ in addition to the friction angle ϕ. Comparisons with the exact solutions obtained from numerical results indicate that the proposed formula is able to provide an accurate approximation with an error of no more than ±2 %.

q and γ, respectively. Equation (1) was proposed by Terzaghi (1943) and assumes that the factors N c , N q and N γ can be obtained by superposition. The soil is treated as weightless when computing N c and N q (i.e., q ≠ 0, c ≠ 0, γ = 0) and as ponderable but having no cohesion or surcharge when calculating N γ (i.e., q = 0, c = 0, γ ≠ 0). The errors caused by this superposition have been discussed by many researchers (Bolton and Lau 1993;Davis and Booker 1971;Griffiths 1982). The superposition approach has been concluded to lead to over-conservative results, which are on the safe side of a design. However, the bearing capacity determined by superposition method superimposes two failure mechanisms, which is different from the real situation. Therefore, it is necessary to calculate the bearing capacity with one failure mechanism to obtain the exact results.
When the bearing capacity is computed on general c-φ-γ soil without superposition and the result is still written in the form of Eq. (1), some researchers have found that the value of N γ relates to not only the soil friction angle φ but also to other parameters, such as q, c, γ and B. Cox (1962) revealed that the parameters associated with stress characteristic equations are φ and a dimensionless parameter G(G = γB/2c) for a smooth footing without surcharge. Chen (1975) introduced a foundation depth and width ratio, D/B, and computed the changes in N γ with the D/B for different internal friction angles. Xiao et al. (1998) calculated the bearing capacity using the method of characteristics (MOC) and revealed that q, c, γ and B all affect N γ , and N γ is only affected by φ and γB/(c + qtan φ) when the load is vertical. Michalowski (1997) and Silvestri (2003) studied the influence of c/γB and q/γB on N γ using the limit analysis method and the limit equilibrium method, respectively. Their research demonstrated that for a given φ, the value of N γ significantly changes with c/γB or q/γB. Zhu et al. (2003) showed that N γ is not only related to the friction angle φ but also to the surcharge ratio λ (λ = (q + ccot φ)/γB). Sun et al. (2013) noted that there are two types of failure mechanisms for rough footings, and whether the trapped non-plastic wedge traverses the footing edge depends on the surcharge ratio λ. Sun et al. (2013) also studied the variation of N γ with φ when λ equals the critical surcharge ratio λ c . The researchers above studied different factors influencing N γ , but none of these studies proposed a formula to calculate N γ .
The MOC is one of main methods applied in the bearing capacity issue which has been discussed by many researchers (Bolton and Lau 1993;Lundgren and Mortensen 1953;Martin 2003;Sokolovskii 1965). The classical bearing capacity factor N γ when q = 0, c = 0, γ ≠ 0 by MOC is calculated to a high degree of precision and is proven to be exact by checking for coincident lower and upper bounds, and by extending the lower bound stress field throughout the semi-infinite soil domain (Martin 2005;Smith 2005). The bearing capacity on general c-φ-γ soil is calculated in one failure mechanism and can be therefore treated as exact solution. In this paper, the MOC is employed to calculate the bearing capacity of strip footings on general c-φ-γ soil and is implemented with a self-coded finite difference method program. The computation of the bearing capacity is carried out with one failure mechanism instead of using superposition approximation, which avoids assuming the shape of the slip lines. The present procedures for computing both smooth and rough footings are unified, which satisfies all of the requirements of the boundary conditions and the symmetric conditions of the surface footings. The numerical results of N γ are compared with other published results, and the sources of errors in the other results are discussed. The bearing capacity factor N γ on general c-φ-γ soil is found to be a function of not only the friction angle φ but also the surcharge ratio λ. Then, the effects of the footing base roughness and the surcharge ratio λ on N γ are investigated. Finally, a curve-fitting formula for N γ that considers φ and the surcharge ratio λ is proposed based on the numerical results.

Characteristic equations at any point in the limit equilibrium state
With the coordinate system depicted in Fig. 1, if a notation is adopted similar to that of Sokolovskii (1965), the normal stresses σ x and σ y and the shear stress τ xy at any point M satisfy the following equilibrium equations: If the soil satisfies the Mohr-Coulomb criterion, the point M in the limit equilibrium state also satisfies the equation where σ 1 and σ 3 are the major and minor principal stresses, respectively. As shown in Fig. 2, the relations of the principal, normal and shear stresses in the limit equilibrium state can be expressed as the following: where σ = (σ 1 + σ 3 )/2 + ccot φ is defined as a characteristic stress, and η is the angle between the major principal stress direction and the horizontal axis ox.
According to Eqs.
(2) and (4), the differential equations along the α and β characteristic lines can be obtained as follows: (2)

Derivation of the finite difference equation
A theoretical formula for the bearing capacity of strip foundations that considers the weight of the soil is not available. However, a finite difference method is usually used to solve Eqs. (5) and (6). Figure 3 displays the α and β characteristic lines through the point M. If A and B are two points on the α and β characteristic lines near M and the state of these two points is known, the dx, dy, dη, and dσ between M and A along the α characteristic line can be approximately expressed as dx = x − x A , dy = y − y A , dη = η − η A and dσ = σ − σ A . By substituting the above formulas into Eq. (5) and taking η = η A and σ = σ A , the finite difference equations are written as follows: Likewise, the finite difference equations along the β characteristic line at point M are (6) along the β characteristic line dy = tan(η + µ)dx dσ + 2σ tan φdη = γ (dy + tan φdx)  The expressions for x, y, η and σ at M can be determined from Eqs. (7) and (8).
The accuracy of the solutions derived from Eqs. (9) to (12) depends on the spacing of the characteristic lines and the approximation error. In this paper, the soil beneath half of the footing base is divided into more than 1000 elements to attain results with good accuracy. To reduce the error due to approximating η = η A and σ = σ A along the α characteristic line and η = η B and σ = σ B along the β characteristic line, an iterative algorithm is used. Set x ′ = x, y ′ = y, η ′ = η and σ ′ = σ after the first calculation, and replace η and σ with η = (η ′ + η A )/2 and σ = (σ ′ + σ A )/2 in Eq. (5) and η = (η ′ + η B )/2 and σ = (σ ′ + σ B )/2 in Eq. (6). Recalculate the solutions using the updated formulas, and repeat the above procedure until all of the four components converge. Convergence is achieved when where Er is the allowable error in the calculation and is assumed to be 10 −15 .
The state of the field is determined point by point, and the final point exists at the bottom of the footing, for which y and η are known. As a result, x and σ at the final point can be directly solved using Eq. (7) without iteration.

Computation procedure
The solution procedure for computing strip footing using the MOC is depicted in Fig. 4. In each figure, OE is the half width of the base, and EG is the center line of the foundation and the axis of symmetry. Note that the edge of footing O is a singularity that can be taken as an α characteristic with zero length. There are a considerable number of x, y, η and σ at point O. The number depends on the divisions of the angle of the transition area DOC. The abscissa x at any point on the free surface OA is known, and the corresponding y and η equal 0. According to Eq. (4), σ = (q + ccot)/(1 − sinφ). Therefore, all of the four components (i.e., x, y, η and σ) at the free surface are known. The computation starts from the free surface OA. The solutions for the internal points can be calculated sequentially along the α characteristic. For a smooth footing, the characteristics end at , the base. The point at the smooth base has η = π/2 and y = 0. The calculation is terminated when x = −B/2. However, for a rough footing, η equals 3π/4 + φ/2 if the characteristic exists at the footing base. It is impossible for all of the soil under the footing base to be in the yield state because of the roughness. Therefore, the critical problem for rough foundations is determining the boundary of the non-plastic wedge and the yield zone. The α characteristics are assumed to progress to the footing at the beginning, and all the α characteristics thus start from the free boundary and end at the footing base. According to the symmetry requirements, there is no shear stress at the center of the foundation. If the last α characteristic and the centerline of the foundation intersect at point I, the point I will have the properties of x = −B/2 and η = π/2 at the same time. The β characteristic noted as FI in Fig. 4b is the boundary of the non-plastic wedge and the yield zone. The area between FI and the footing is the non-plastic wedge in which the dotted lines in the figure represent nonexistent characteristics. The region FIAO is the yield zone, and the characteristics in this region are real.
From the construction of the characteristic field above, the computational procedure for both smooth and rough footings can be unified. The computation initiates with α characteristics progressing from the free surface to the footing and terminates at the point in which x = −B/2 and η = π/2 simultaneously. The area enclosed by the β b a Fig. 4 Solution procedure for the MOC. a Characteristics of a smooth footing, b characteristics of a rough footing characteristic passing through the terminal point and the footing base is the non-plastic zone. The smooth footing is simply the special case presented in Fig. 4b in which the terminal point I coincides with point F and the middle point E of the footing base, which indicates that there is no non-plastic wedge, as shown in Fig. 4a.
The bearing capacity of the smooth and rough footings is given by where σ y and τ xy are the normal stress and shear stress in the y direction, respectively, and are obtained from Eq. (4); dW is the differential weight of the soil wedge EFI along the curve FI. For a smooth foundation, the two points I and F coincide at the point E, and Eq. (14) is simplified to Martin (2003) and Sun et al. (2013) found that there are two types of failure mechanisms for rough footings when constructing the stress field. In one type, no α characteristics progress to the footing base; therefore, a complete non-plastic wedge exists under the footing base. Whereas in the other type, the α characteristics enter the region beneath the footing and result in a partial non-plastic wedge. Sun et al. (2013) stated that the type of failure mechanism depends on φ and λ. The proposed construction method of the characteristics in this paper satisfies all the boundary requirements, and the type of failure mechanism is automatically determined using the computation. Moreover, the computation of the bearing capacity of smooth and rough footings is unified with the same termination condition.

Equivalent solution of the bearing capacity problem
It is widely accepted that the bearing capacity factors N c and N q are correlated with each other through the following formula: When the bearing capacity is computed on general c-φ-γ soil without superposition and the result is still written in the form of Eq. (1), the bearing capacity factor N γ is not the value that computed by superposition method.
Combining Eqs. (16) and (1) gives By dividing both sides of Eq. (17) by γB and setting p u = (q u + ccot φ)/γB and λ = (q + ccot φ)/γB, the bearing capacity formula is further transformed to (18) is a general solution of the bearing capacity of strip footings that is equivalent to Eq. (1). The N γ values deduced by exact bearing capacity equal to those by superposition method only when λ = 0. To obtain an exact solution of N γ on general c-φ-γ soil, it is necessary to calculate the bearing capacity in the real failure mechanism using a method other than the superposition approximation proposed by Terzaghi. Zhu et al. (2003) computed the bearing capacity factor N γ of rough strip foundations using the critical slip field method. p u and λ are defined as the normalized bearing capacity and the surcharge ratio, respectively. The value of N γ was found to be influenced not only by φ but also by the surcharge ratio λ. However, Zhu et al. (2003) assumed that the inclined angle of the active wedge underneath the footing was π/4 + φ/2 with respect to the horizontal line, leading to discrepancies between the calculations and the exact solutions. The proposed method in this paper avoids this assumption and results in better numerical results. Moreover, the improved computation extends the application to both a smooth footing and a rough footing. This approach is helpful in attaining a better fitting formula for N γ based on the exact numerical results.
From Eq. (18), N γ can be written as The numerical results of Zhu et al. (2003) implied that p u is constant with fixed values for φ and λ. The calculations using the MOC in this paper also confirm this conclusion. Consequently, N γ is influenced by φ and λ as long as N q is a function of φ or φ and λ. The bearing capacity factor N q is typically regarded as not being influenced by the soil weight with a theoretical formula given by Reissner (1924) as Shield (1954) studied the bearing capacity of strip footings using plastic theory and reached the conclusion that the well-known, closed-form expressions of N q and N c given by Reissner (1924) and Prandtl (1921) are exact solutions for weightless soil regardless of the footing roughness. It is most straightforward to use the closed form solutions for N q and N c derived for weightless soil and use N γ to account for all the effects of self weight and its interaction with q and c, by using the surcharge ratio λ.

Comparisons of N γ with other known results
The results of N γ have been computed by many investigators based on cohesionless soil with no surcharge load (that is, λ = 0 in this paper). The present results of N γ for smooth footings when λ = 0 are listed in Table 1. It is noted that the results of N γ corresponding to λ = 0 are computed when λ = 10 −10 . The results calculated by other researchers are also presented in Table 1.
The present values when λ is 0 in Table 1 are equivalent to complete solutions given by Martin (2005) and Smith (2005), which indicates that the results by present method can be treated as exact solutions. The computations provided by Bolton and Lau (1993) or Kumar (2009) using the MOC have little difference compared to present results. Moreover, the results determined by the proposed method are between the upper and lower   bounds given by Hjiaj et al. (2005). The N γ values determined by Woodward and Griffiths (1998) using the finite element method are consistent with those obtained using the MOC. Compared to the calculations given by Hjiaj et al. (2005) and Smith (2005), the results from Frydman and Burd (1997) using fast lagrangian analysis of continua(FLAC) exceed the upper bound when φ equals 35° and are below the lower bound when φ equals 40° or 45°, although the errors are small compared with the calculations in this paper. As inferred in previous section, the value of N γ depends only on λ at a determined φ if the bearing capacity is calculated without superposition on general c-φ-γ soil. The values of N γ are given in Table 1 as well when λ equals to 0.1, 1, 10, 100 and 10 4 . The computations in Table 1 show that the value of N γ increases with the growth of λ if the friction angle is determined. When λ equals to 10 4 or even larger, the value of N γ is found to approach the theoretical upper bound given by Chen (1975) in the Hill mechanism: where u = π/4 + φ/2, f = tan φ.
The calculations of N γ in Table 1 when equals 10 4 have errors of no more than 0.1 % compared to the solutions by Eq. (21). When the weight of soil decreases to 0, the surcharge ratio λ will approach ∞. In this case, the failure surface computed by MOC is consistent with the Hill mechanism. So the upper bound of N γ in Eq. (21) can be treated as the exact theoretical solution when λ = ∞. Table 2 contains the N γ results for rough footings. Similar to the results for a smooth footing, λ is treated as λ = 10 −10 when λ = 0. Some results associated with λ = 0 that were published by other researchers are also listed in Table 2.
The present values of N γ when λ = 0 are basically treated as exact solutions because they are equal to complete solutions computed by Martin (2005) and Smith (2005). The N γ results obtained by Bolton and Lau (1993) are much greater than the present values and even exceed the maximum values corresponding to λ = 10 4 when φ equals 5° or 10°. The large errors are mainly ascribed to the assumption that the trapped wedge beneath the foundation has a base angle of π/4 + φ/2. Kumar (2009) abandoned this assumption and determined the partly trapped wedge by computation and, consequently, obtained better results. Following Terzaghi's assumptions, Kumbhojkar (1993) achieved a numerical solution for N γ , and the results are in agreement with Terzaghi's calculations. Zhu et al. (2001) determined the base angle of the active wedge when N γ is a minimum using the method of triangular slices, and the corresponding N γ results are better than those determined by Kumbhojkar (1993). Hjiaj et al. (2005) meshed fine finite elements to determine the yield zones instead using an arbitrary assumption. The errors do not exceed 3.42 % between the rigorous lower and upper bound solutions, and the results are in good agreement with the present calculations.
When λ equals to 0.1, 1, 10, 100 and 10 4 , the values of N γ are also given in Table 2. Similar to smooth footings, the value of N γ approaches the upper bound in the Prandtl mechanism when λ equals to 10 4 . The exact theoretical solution of the upper bound is also given by Chen (1975), which is twice of the value calculated with Eq. (21). The values of N γ in Table 2 when λ equals 10 4 are basically equal to the theoretical solutions with the errors less than 0.1 %.

Ratio of N γ for smooth and rough footings
The numerical results in Tables 1 and 2 indicate that the N γ values have large differences for smooth and rough foundations for given φ and λ values. This implies that the roughness of the footing base has a large impact on N γ . Accordingly, the ratio of the N γ values for smooth and the rough foundations is defined as R N , i.e., where N s γ and N r γ are the bearing capacity factors N γ for the smooth and rough foundations, respectively, for given φ and λ values. The curves of R N versus φ with different λ are plotted in Fig. 5, and the results given by Hjiaj et al. (2005) are also marked in the figure. The computations by Hjiaj et al. (2005) in the case of q = 0 and c = 0 are in good agreement with the curve for λ = 0.
The numerical calculations of N γ for smooth footings and rough footings reveal that the N γ value for a smooth footing is only half or more than half of that for a rough footing. Figure 5 also demonstrates that R N becomes less sensitive to φ as λ increases. Equivalent to the solution for a granular soil with zero surcharge, the numerical result of N γ when λ equals 0 is a minimum solution with a determined φ. For a rough foundation, the collapsed surface when λ = ∞ is the same as that in the Prandtl mechanism, and the computational result of N γ equals the closed-form solution deduced in the Prandtl mechanism. Similarly, the N γ for a smooth footing is identical to the theoretical expression in Hill's failure mechanism. As stated by Chen (1975), the N γ in the Prandtl mechanism is exactly twice the value in the Hill mechanism, i.e., R N = 0.5. The relationship of R N and φ when λ = ∞ in Fig. 5 verifies Chen's judgment.

Influence of the surcharge ratio on N γ
The computations of N γ also exhibit large discrepancies when λ = 0 and λ = ∞ at the same φ regardless of having a smooth or rough footing base. To distinguish the N γ for different values of λ, N γ is noted as N γ,min when λ = 0 and as N γ,max when λ = ∞. It can be easily inferred that the N γ,max of a smooth foundation is exactly calculated by Eq. Numerical results of the R N ratio versus φ and the N γ,max of a rough footing is twice of that of a smooth footing. If K N is defined as the ratio of N γ at λ = 0 and λ = ∞, then Figure 6 reveals the relationship of K N to φ for smooth and rough footings. A smaller φ results in a larger difference in K N between smooth and rough footings. The values of K N for smooth and rough foundations are very close when φ is greater than 40°. Furthermore, the numerical results of K N can be approximated using a polynomial expression in the form below Taking n as 4 in Eq. (24), the error between the approximated and numerical values appears to be no more than ±1 % when φ is greater than 2°. Therefore, the suggested expression for smooth footings can be written as follows: Furthermore, the fitting formula for K N for a rough footing is given as follows: The curves of K N versus φ according to Eqs. (25) and (26) are plotted in Fig. 6, and the curves are in good agreement with the numerical computations.

Proposed formula of N γ
The bearing capacity factor N γ is influenced by both λ and φ regardless of Eq. (19) or the numerical calculations using the MOC. The calculations of N γ related to λ with a series of φ are plotted in Figs. 7 and 8 for smooth footings and rough footings, respectively. The where p and A 0 are fitting parameters. As mentioned previously, the N γ,max of a foundation can be exactly calculated by theoretical solutions. However, the exact closed form solution of N γ,min is not available when λ = 0, although plenty of empirical formulas are given by different researchers. The N γ,min value is proposed to be calculated by the solution of N γ,max times K N based on Eqs. (25) and (26).
The parameter p ranges from 0.75 to 0.8 for both smooth and rough footings. The value of p was selected to be 0.75 because the variation of p from 0.75 to 0.8 has little effect on N γ . The factor A 0 is a fitting coefficient that is related to φ and is defined as The value of A 0 has satisfactory accuracy when n is 3. The proposed formulas for A 0 are given below: The curves of N γ versus λ for smooth and rough footings are plotted in Figs. 7 and 8, respectively, based on the fitting formula (27). For both smooth and rough footings, the approximate results agree well with the numerical results within errors of ±2 %. Therefore, Eq. (27) is able to estimate N γ with adequate accuracy.
The N γ data computed by Zhu et al. (2003) are given in Fig. 8 as well and are much greater than the results of the proposed method when λ is less than 10. As mentioned above, the discrepancies in the results are mainly attributed to the assumption that the base angle of the active wedge underneath the footing base equals 45° + φ/2. The error resulting from that assumption rapidly decreases with increasing λ. As seen in Fig. 8, there is little difference between the present values and the results provided by Zhu et al. (2003) when λ is greater than 10. Moreover, the theoretical solution of N γ given by Zhu et al. (2003) is the same as the present value in the case of λ = ∞. Similar to the pattern of the present results, the calculations by Zhu et al. (2003) also have an "S" shape for a fixed φ, which implies that the results can be estimated by expression (27) as well, except N γ,min , A 0 and p differ from the values in this paper.
It should be noted the proposed approximate formula of N γ in Eq. (27) is limited to the classic issue on the bearing capacity of strip footings that the soil is treated as a rigid plastic and obeys Mohr-Coulomb criterion. If the soil beneath the strip footing does not flow Mohr-Coulomb criterion, the proposed method may be no longer applicable. Because the suggested approximate formula is based on the conclusion that the bearing capacity factor N γ depends on the surcharge ratio λ in addition to the friction angle φ. This conclusion is only valid for Mohr-Coulomb soil. Further research is required for smooth footings A 0 = 0.222 tan 3 φ + 0.101 tan 2 φ + 0.102 tan φ + 0.188 (30) for rough footings A 0 = 0.22 tan 3 φ + 0.684 tan 2 φ − 0.042 tan φ + 0.354 whether the conclusion is suitable when the soil meets other yield criteria other than Mohr-Coulomb criterion. A comprehensive research is also needed whether the conclusions and proposed formula in this paper can be extended to circular or rectangular footings.

Conclusions
The MOC is employed to calculate the exact bearing capacity of strip footings. By considering the influence of c, φ and γ with one failure mechanism, the computational procedures for smooth and rough foundations are unified without assuming the failure mechanism. The computations were implemented using a self-coded finite difference program. If the bearing capacity of the footings is calculated using the formula proposed by Terzaghi (1943) and N q and N c are obtained using the theoretical solutions given by Prandtl (1921) and Reissner (1924), the value of N γ is influenced by not only the friction angle φ but also by the surcharge ratio λ. The computations were compared with other published results. The comparisons and analysis indicate the following conclusions: 1. In the case of no overload, the computed N γ values in this paper for a granular soil are treated as exact solutions because the values are consistent with complete solutions given by Martin (2005) and Smith (2005). Some researchers assume failure surfaces or mechanisms that are not the same as the real state; therefore, their results have considerable errors compared with the exact solutions. 2. The roughness of the footing base has a significant impact on N γ . The ratio of the bearing capacity factor N γ for smooth foundations and rough foundations, which is R N , indicates that the value of N γ for a smooth footing is only half or more of that for a rough footing. The curve of R N versus φ with λ = 0 has good agreement with the results given by Hjiaj et al. (2005). A value of R N equal to 0.5 when λ = ∞ supports Chen's statement that the N γ in the Prandtl mechanism is exactly twice the value in the Hill mechanism. 3. The surcharge ratio λ also significantly affects N γ , and the ratio K N defined by N γ,min /N γ,max can be approximately evaluated using a polynomial expression when λ = 0 and λ = ∞. When λ is sufficiently large, the solution of N γ is demonstrated to approach the upper bound that deduced by Chen (1975) in a closed-form solution. Therefore, N γ,max is obtained by exact theoretical formula, and thus, N γ,min can be accurately estimated using K N and N γ,max . 4. The present N γ value in the case of λ = ∞ is exactly the same as the theoretical solution of N γ given by Zhu et al. (2003). However, the calculations of Zhu et al. (2003) have obvious errors compared with the present results when λ is less than 10 primarily due to the assumption that the base angle of the active wedge underneath the footing base equals 45° + φ/2. The values of N γ can be calculated by the approximate formula (27) containing two factors: φ and λ. The discrepancies between the approximate results and the numerical solutions are less than ±2 % for both smooth and rough foundations. Formula (27) is demonstrated to be suitable for evaluating N γ when considering the factor λ.