Seismic bearing capacity of a strip footing on rock media

Obaidur Rahaman, Jyant Kumar

Department of Civil Engineering, Indian Institute of Science, Bangalore, Karnataka, 560012, India

Keywords:Earthquake loads Limit analysis Lower bound Power cone programming (PCP)Rock mass Strip footing

ABSTRACT The bearing capacity factors for a rough strip footing placed on rock media,which is subjected to pseudostatic horizontal earthquake body forces, have been determined using the lower bound finite element limit analysis in conjunction with the power cone programming (PCP).The rock mass is assumed to follow the generalized Hoek-Brown(GHB)yield criterion.No assumption needs to be made to smoothen the GHB yield criterion and the convergence is found to achieve quite rapidly while performing the optimization with the usage of the PCP.While incorporating the variation in horizontal earthquake acceleration coefficient (kh), the effect of changes in unit weight of rock mass (γ), ground surcharge pressure (q0) and the associated GHB material shear strength parameters (geological strength index(GSI), yield parameter (mi), uniaxial compressive strength(σci)) on the bearing capacity factors has been thoroughly assessed.Non-dimensional charts have been developed for design purpose.The accuracy of the present analysis has been duly checked by comparing the obtained results with the different solutions reported in the literature.The failure patterns have also been examined in detail.

1.Introduction

A number of studies have been performed in the literature to determine the bearing capacity for different types of foundations resting over rock media (Serrano and Olalla, 1994, 1996; Serrano et al., 2000; Yang et al., 2003; Yang and Yin, 2005; Saada et al.,2008, 2011; Merifield et al., 2006; Yang, 2009; Keshavarz et al.,2016; Kumar and Mohapatra, 2017; Keshavarz and Kumar, 2021).Using the stress characteristics method, Serrano and Olalla (1994)computed the bearing capacity of a strip footing placed over weightless rock media,with horizontal ground surface,on the basis of the original Hoek-Brown failure criterion (Hoek and Brown,1980).The same methodology was later employed by Serrano and Olalla (1996) to compute the bearing capacity of a strip footing placed over weightless rock media but with sloping ground surface.Serrano et al.(2000) obtained the bearing capacity of a strip footing over rock media to consider the modified version of the Hoek-Brown failure criterion (Hoek et al., 1992).Using the lower bound (LB) limit analysis but with an assumption of the failure mechanism and using the modified Hoek-Brown criterion,Yang et al.(2003)computed the bearing capacity of a strip footing resting over weightless rock media.The upper bound(UB)solution for computing the bearing capacity based on the generalized Hoek-Brown (GHB) yield criterion (Hoek et al., 2002) was obtained by Yang and Yin(2005),using a generalized tangential technique and with an assumption of the multi-wedge translation failure mechanism.In this method, instead of using the actual GHB failure criterion,a linearized Mohr-Coulomb failure criterion,which remains always tangential to the GHB failure envelope, was used.This procedure was,however,found to overestimate the actual ultimate bearing capacity as demonstrated by Saada et al.(2008).Merifield et al.(2006) employed the lower and UB finite element limit analysis (FELA) technique and bracketed the ultimate bearing capacity of a strip footing on rock masses approximately within 2%.The effects of the inclination of a rock slope,ranging from 0?to 30?,and the horizontal seismic coefficient kh,ranging from 0 to 0.2,on the bearing capacity factor of a strip footing placed on the edge of a rock slope, have been examined by Yang (2009), using the UB theorem of the limit analysis,based on the methodology proposed by Yang and Yin (2005).Saada et al.(2011) implemented the kinematic approach of the limit analysis to determine the reduction in the bearing capacity of a spread footing placed on rock slopes with an inclusion of the seismic forces.The ultimate bearing capacity of a strip footing placed on weightless rock media in the presence of seismic forces has also been assessed by Keshavarz et al.(2016)using the stress characteristic method and with the usage of the earlier version of the Hoek-Brown yield criterion (Hoek andBrown, 1980).Keshavarz and Kumar (2021) also computed the bearing capacity of a ring foundation using the stress characteristics method and the obtained solution was also compared with that computed on the basis of the FELA.

Amongst the numerous methods for solving different geomechanics stability problems,the FELA(Lysmer,1970;Sloan,1988)is the most rigorous one in terms of accuracy and computational efficiency.It inherits a number of advantages over the other conventional stability methods (Sloan, 2013).Unlike the limit equilibrium method, this method does not require any kind of assumption associated with the failure mechanism.Rather, this method generates the actual collapse mechanism from the analysis.Like the displacements based elasto-plastic finite element method,it can handle irregular geometry,complicated boundary condition,complex loading condition, material anisotropy and inhomogeneity.In the elasto-plastic finite element method, the definition of the complete material constitutive relationship is necessary.On the other hand,the FELA considers the problem only on the verge of failure/collapse.Therefore, in the FELA, it requires only the shear strength parameters for performing the stability analysis.Unlike the conventional elasto-plastic finite element method, the need of any step-by-step incremental analysis is completely avoided in the FELA which eventually results in saving of the computational cost.

Over the years, the efficiency of this methodology has been significantly improved especially with an introduction to the conic programming techniques such as the second order cone programming (SOCP) (Makrodimopoulos and Martin, 2006) and the semidefinite programming(SDP)(Krabbenhoft et al.,2008).For the rock mass obeying the Hoek-Brown yield criterion, Kumar and Mohapatra (2017) applied the SDP technique which is computationally more efficient than the nonlinear programming (NLP), to determine the bearing capacity factors for strip and circular footing based on the LB FELA.The same technique was later implemented by Ukritchon and Keawsawasvong (2018) for dealing with various three-dimensional (3D) stability problems.The proposed SDP formulation,however,overestimates the bearing capacity factors as it requires a slight modification in the GHB criterion where the value of the exponential factor needs to be kept equal to 0.5.This limitation was later overcome by Kumar and Rahaman (2020) by implementing the power cone programming (PCP) while employing the LB FELA.This technique was also proved to be computationally more efficient.

Although a number of studies with different methodologies are available to estimate the ultimate bearing capacity of a strip footing under static condition, only a few investigations have been performed to account for the consideration of the seismic forces.Most of the existing studies related to the seismic bearing capacity factors are based on the UB analysis(Yang,2009;Saada et al.,2011)or the slip line method(Keshavarz et al.,2016);in both the methods a predefined failure mechanism has been assumed.Moreover, these studies do not consider the most recent GHB failure criterion(Hoek and Brown,2019)in its true form.The present study examines the effect of the pseudo-static earthquake inertial body forces while computing the bearing capacity factors for a strip footing resting over rock media and subjected to surcharge pressure.The current analysis is based on the application of the LB FELA and the PCP.The analysis involves the implementation of the GHB yield criterion in its true form without requiring any kind of assumption associated with the smoothing of the GHB yield criterion.

2.The yield criterion for rock mass

Among the various failure criteria proposed in the literature,the GHB failure criterion as proposed by Hoek et al.(2002) was considered to be the most acceptable basis for modelling the shear strength of rock mass.This yield criterion was developed through a series of laboratory tests which covered numerous types of undisturbed and disturbed rock specimens over a wide range of confining stresses.This criterion is used extensively by practitioners for a variety of rock engineering projects(Ulusay,2014).By Using this criterion, a number of stability problems have been solved:for instance,(i)the determination of the bearing capacity of footings(Merifield et al.,2006;Chakraborty and Kumar,2015);(ii)the stability of slopes (Li et al., 2008); and (iii) the stability of tunnels(Suchowerska et al.,2012;Ukritchon and Keawsawasvong,2019; Rahaman and Kumar, 2020).The latest GHB criterion (Hoek et al., 2002; Hoek and Brown, 2019), which is followed in this study, has been written in the following form:

Note that in the above expression,there is a change in the sign convention as compared to failure criterion reported by Hoek and Brown (2019).In this expression, the tensile normal stress is considered to be positive.The variables σ1and σ3imply the major and minor principal stresses at failure,respectively;and σcidefines the uniaxial compressive strength of the rock specimen.The strength parameters mb, s and α are described by the following equations:

Note that all these parameters are a function of the geological strength index(GSI),disturbance factor(D)and the yield parameter(mi).The parameter GSI signifies a measure of the rock mass strength under different geological conditions.The range of the GSI typically varies from 10 for an extremely poor rock to 100 for an intact rock mass.The factor D is a disturbance coefficient that incorporates the disturbance of the rock mass which may occur for instance due to blast damage, impact loading and sudden stress relaxation.It varies from 0 for an undisturbed rock mass to 1 for very disturbed rock mass.The parameter miwhich can be obtained from triaxial compression test data is analogous to the frictional strength of the intact rock mass.The approximate values of mifor different types of rocks are described in detail by Hoek(1990).It has an approximate range of 5-35.

3.Description of the problem

A strip footing of width B is placed on rock media with horizontal ground surface,as shown in Fig.1a.The footing is assumed to be perfectly rough.The rock mass has a unit weight γ and the ground surface is subjected to vertical uniform surcharge pressure of q0.The rock mass is assumed to follow the GHB yield criterion.This study employs the pseudo-static approach to consider the influence of seismic inertial forces in the system.The rock media and the overlying superstructure are subjected to horizontal seismic acceleration of khg,where khis the acceleration coefficient,and g is the acceleration due to gravity.The load on the footing comes from the superstructure weight and the horizontal pseudostatic earthquake inertial forces.If the vertical load on the footing per unit length at collapse is Qu, the footing will be subjected to ahorizontal shear force of khQualong the interface of the footing and underlying rock mass.Similarly,if the ground surface is loaded with a surcharge pressure of q0, it will also be subjected to a uniform shear stress of khq0.The rock media will be subjected to body forces per unit volume of γ and khγ in vertical and horizontal directions,respectively.It is to determine the magnitude of Qufor various specified values of material parameters defining the GHB yield criterion, and different values of γ, q0and kh.

Fig.1.(a)Problem definition and stress boundary conditions,and(b)A typical chosen mesh.

It should be mentioned that the best way to analyze the stability of any structure due to occurrence of an earthquake is to perform a dynamic analysis with the consideration of the actual/design acceleration-time history.However, this kind of approach often requires a complex numerical analysis which in turn needs quite high computational effort.In most of the geotechnical stability analyses in the presence of seismic forces, the conventional pseudo-static method is often employed to determine the stability of any structure in the event of occurrence of an earthquake.In this approach, the seismic inertial force is replaced by an equivalent static force which is calculated by multiplying the unit weight of the medium with the seismic acceleration coefficient.The magnitude of the seismic acceleration coefficient is determined from the knowledge of the peak ground acceleration (PGA) for the design earthquake.The analysis following the pseudo-static method is generally considered to be conservative(Li et al.,2009;Yang,2009;Saada et al.,2011).Despite its limitations,the conventional pseudostatic method is widely employed due to its simplicity and ease in implementation.In the presence of seismic forces,this method has been applied by various researchers for finding the solutions of numerous geotechnical stability problems, such as finding (i) the bearing capacity of foundations resting on horizontal ground surface (Kumar and Rao, 2002; Keshavarz et al., 2016), (ii) the bearing capacity of foundation on slopes (Kumar and Rao, 2003;Saada et al., 2011; Kumar and Chakraborty, 2013), and (iii) the stability of slopes(Li et al.,2009),tunnels(Sahoo and Kumar,2012;Saada et al., 2013), and retaining walls (Kumar 2001; Conte et al.,2017; Chehade et al., 2021).The pseudo-static approach is also used to countercheck the results of a more sophisticated dynamic analysis.Therefore,this approach has been employed in the current research.

4.Analysis

In order to solve the stability problem which involves layers,complicated geometry and complex boundary conditions,the FELA is proved to be a powerful tool.This method involves(i)the LB and UB limit theorems of the plasticity,(ii)the discretization technique of the finite elements,and(iii)an optimization technique.Although both the LB and UB solutions help to bracket the true limit load in a bound form,in practice,however,the LB analysis predicts a collapse load which is always less than or equal to the true answer.In the present study,the LB FELA is employed for the computation of the safe design load.

4.1.Finite elements

The very first step in the formulation of the LB FELA is to discretize the selected domain with a mesh of chosen finite elements.The present study employs the three noded linear stress triangular elements.The chosen problem domain and a typical triangular element with the unknown nodal stresses (σx,σyand τxy) are presented in Fig.1a.The stresses(σx,σyand τxy)at any point within the element are assumed to vary linearly according to the following expressions:

where σx,i, σy,iand τxy,iare the stresses at the node i, and the parameter Nirepresents the linear shape function.If the coordinate of the node i is denoted by (xi,yi) and the three nodes of a triangular element are numbered in a counterclockwise direction with 1,2 and 3,then the expression for the shape function associated with the node 1 is given as

where Δ represents the elemental area,i.e.Δ = 0?5|(x1-x3)(y2-y3) - (x3- x2)(y3- y1)|.For rest of the two nodes, the shape functions can further be defined by following the same pattern.

The statically admissible stress discontinuities in between adjoining elements are permitted in the LB FELA.The discontinuous stress field generates additional degrees of freedom which help to improve the LB solution significantly(Sloan,1988; Pastor, 2003).

Fig.1b presents a typical mesh used in the analysis.Note that the sizes of the elements gradually decrease as it approaches towards the edges of the strip footing.To simulate the singularity of stress at the footing edge,a fan type of the mesh is chosen along the edges of the footing.The horizontal and vertical extents of the semi-circular domain (UVR), denoted by Lhand Lv, respectively, are kept sufficiently large so that the proximity of the stresses near the edge of the boundaries remains far away from the yield.It is also ensured that the magnitude of the collapse load remains unchanged even if the chosen boundaries of the domain are extended further.

4.2.LB formulation

As per the theory of the plasticity, the LB solution is achieved from a statically admissible stress field which satisfies equilibrium conditions everywhere in the domain and stress boundary conditions and it nowhere violates the yield criterion (Sloan,1988).

4.2.1.Equilibrium equation

As the rock media are subjected to horizontal pseudo-static seismic acceleration kh, the equilibrium of each element must be satisfied using the following equations:

Note that the pseudo-static horizontal seismic body forces are incorporated in the analysis with the application of Eq.(5a).For kh= 0, the term on the right hand side of equal sign in Eq.(5a)becomes simply equal to zero.Substituting Eq.(3) into Eq.(5)will lead to the following equality constraint(s):

where Aequiand bequirepresent the matrix and vector associated with the equilibrium conditions,respectively;andcomprises the global vector of unknown stresses at the nodes.

4.2.2.Stress discontinuities

The continuities of normal and shear stresses are maintained along the interfaces of all the elements.Along the interface of any two adjacent elements i and j,with the nodal pairs(1,2)and(3,4),the following conditions must hold good:

where subscripts n and tn represent the normal and tangential directions, respectively, along the discontinuity line.The nodes 1 and 3 are associated with the element i,and the nodes 2 and 4 form a part of the element j.This will result in the generation of the following equality constraints:

where Asdand bsdrepresent the matrix and vector associated with the continuity of the stresses along the interface (discontinuity line) of the two adjoining elements.

It should be mentioned that the state of stress for any plane strain problem can be defined with the usage of three stress variables.The continuity of the normal and shear stresses only ensures the equal values of the two of the three stress variables on either side of the stress discontinuity.The difference in the value of the third stress variable makes the stress state different on either side of the stress discontinuity.

4.2.3.Boundary conditions

The boundary conditions applicable for this problem are shown in Fig.1a,in which ST denotes the base of the footing.Note that at failure, the vertical load on the footing is Quand the horizontal shear force is khQu.The extent of the semi-circular periphery(UVR)is kept sufficiently large, and it is ensured, like everywhere within the problem domain, the GHB yield condition is nowhere violated along this boundary.

(1) Boundary conditions along the ground surface

A vertical surcharge of q0is applied over ground surface either side of the footing, i.e.along the boundaries RS and TU.Therefore,the vertical normal stress (σn) along these boundaries will be simply equal to the applied surcharge pressure (q0).Due to horizontal inertial seismic force, a shear stress of khq0will develop along these boundaries.Hence, the following boundary stresses need to be enforced along the boundaries RS and TU:

(2) Boundary conditions along the footing-rock interface

Since the footing is assumed to be rough, the yield strength of the footing-rock interface is assumed to be the same as that of the yield strength of the rock medium.Therefore, no exclusive constraint on the shear stress needs to be imposed for any node along the footing-rock interface.However, the effect of the horizontal inertial force due to the overlying mass of the superstructure is incorporated by introducing the following equality constraint along the footing-rock interface:

The imposition of the stress boundary conditions can be specified by the following condition:

where Abcand bbcrepresent the matrix and vector in connection with the applied boundary conditions, respectively.

4.2.4.Imposition of the yield condition

The GHB yield criterion,as defined by Eq.(1),can be expressed as

where p = -mb(-σci)(1-α)/α and r = s(-σci)1/α.

Introducing a new variable t such that

Since σ1>σ3and t ≥0, hence, Eq.(12) can be re-written as

For a plane strain problem, we have

Therefore,Eq.(14),after substituting Eqs.(15a)and(15b),can be written as

Since t ≥0,Eq.(16)can easily be expressed in the standard form of a second order conic constraint:

where x1= t, x2= (σx-σy) and x3= 2τxy.

For a plane strain case problem,using Eq.(15a),the value of t in Eq.(13) can be written in terms of σx,σyand τxy:

Since the value of the exponent α is always positive, using Eqs.(16) and (18), it can be written as

If x4, x5and x6are defined by the following expressions:

Then, Eq.(19) will take the form of a power conic constraint as defined by

Note that the variable x4needs to be always positive.

Accordingly, the optimization problem in the LB FELA can be solved by imposing the GHB failure criterion as the summation of one quadratic conic constraint (Eq.(17)) and one power conic constraint (Eq.(21)).

It should be mentioned that an n-dimensional power cone can be defined as

In the present formulation,the following form of the 3D power cone has been considered:

This form of the GHB yield constraint has already been indicated in Eq.(21) with u3= x6, u1= x4and u2= x5.

4.3.The optimization problem

Finally, all the constraints were assembled to maximize the objective function, i.e.the magnitude of the collapse load, by forming the following canonical form of the optimization problem:

where c is the vector containing the coefficients of the objective function; Aeqis the matrix containing the coefficients of all the equality constraints as expressed in Eqs.(6), (8) and (11); beqis corresponding right known vector of the equality constraints;Aconeis the matrix containing the coefficients of all the conic constraints,and bconeis the corresponding right known vector of the conic constraints.If the total number of nodes is denoted by NN, the vectors σ and xconeare defined as

The computer code to perform the LB FELA with the usage of the GHB was written in MATLAB.For this large-scale conic optimization problem, the solver MOSEK, which is based on the primal-dual interior-point method (Andersion et al., 2003) and has already been recommended by a number of researchers (Martin and Makrodimopoulos,2008;Tang et al.,2014;Kumar and Mohapatra,2017; Ukritchon and Keawsawasvong, 2018) because of its robustness and computational efficiency,was used in the present analysis.The computations were carried out on a desktop computer (Intel Core i7-7700 K CPU @ 4.20 GHz,16 GB RAM) in the Windows 10 operating-based system.

5.Results and discussions

5.1.Bearing capacity equation

Similar to Terzaghi (1943)’s bearing capacity equation for a footing resting on a soil medium,the following equation is used to determine the ultimate bearing capacity (qu= Qu/A) of a strip footing placed on rock mass:

where Nσ and Nqrefer to the non-dimensional bearing capacity factors associated with the self-weight of rock mass and surcharge pressure,respectively.Deriving the bearing capacity factor(Nσ)due to unit weight of the rock mass requires surcharge pressure to be zero (q0= 0), thus it is calculated as

If the rock mass is considered to be weightless,then the bearing capacity factor Nσ is termed as Nσ0.

The bearing capacity factor due to surcharge (Nq) can be expressed as

The bearing capacity factors for different seismic coefficients(kh),ranging from 0.0 to 0.5(Terzaghi,1950),have been computed and the final results are presented in both graphical and tabular forms for all the practical range of the Hoek-Brown parameters as mentioned in the earlier section.The effects of γ and q0on bearing capacity factors have been examined in the form of dimensionless parameters σci/(γB)and q0/σci.The value of σci/(γB)is varied from 100 to infinity (inf); the infinite value of σci/(γB) implies a zero value of the unit weight.The value of q0/σciwas varied from 0 to 1.It should be noted that the rock mass in this study is considered to be an undisturbed (D = 0) one.

5.2.Validation of the present analysis

The present analysis has been validated by comparing the obtained solution with that reported in the literature.The comparison of the results has been carried out for the static case as well with the consideration of the pseudo-static inertial forces.On account of non-availability of the results, at present,no comparison could be,however,made for the true dynamic case;this will form the scope of the future study.

5.2.1.Comparison for static case

The accuracy of the present solutions is verified by comparing the bearing capacity factor Nσ0for a weightless rock medium with the available results of Kulhawy and Carter (1992) using the LB analysis but with an assumption of the collapse mechanism,Serrano et al.(2000)on the basis of the slip line method,Merifield et al.(2006)using the average of the LB and UB solutions using the nonlinear optimization, and the LB solution of Chakraborty and Kumar (2015).The comparisons of all these results are shown in Table 1 for kh= 0.The bearing capacity factors proposed by Kulhawy and Carter (1992) from the simplified LB solutions are always found to be remarkably lower than all the other reported solutions.The slip-line solutions estimated by Serrano et al.(2000)are found to be quite close to the present FELA solutions and the difference between the two solutions is found to be very nominal.The present LB solutions are also found to match well with the average of the LB and UB solutions of Merifield et al.(2006),and the LB solutions of Chakraborty and Kumar (2015).The studies ofMerifield et al.(2006)and Chakraborty and Kumar(2015)are based on the FELA technique using the nonlinear programming.

5.2.2.Comparison for pseudo-static case

For γ = 0, q0= 0 and kh= 0-0.2, Table 2 provides a comparison of the present results with the existing solution of Saada et al.(2011) on the basis of the kinematic limit analysis approach with the incorporation of the pseudo-static inertial forces.It can be noted that the values of Nσ0 reported by Saada et al.(2011) are found to be marginally greater than the present LB solution; it is quite an expected outcome since the present analysis is based on the application of the LB limit theorem.

Table 2A comparison of the bearing capacity factor for γ = 0 and q0 = 0 with kh ≥0.

5.3.The variation of the bearing capacity factors with kh

The variation of the bearing capacity factor Nσ0for a weightless rock medium with respect to changes in kh,ranging from 0 to 0.5,is illustrated in Fig.2.This figure displays the variation of Nσ0for seven different values of miranging from 5 to 35,with an interval of 5,and six different values of GSI,i.e.10,20,40,60,80 and 100.It can be clearly noted that the magnitude of Nσ0decreases continuously with an increase in the value of kh.For GSI=10 and mi=5,with an increase in khfrom 0 to 0.2,the factor Nσ0decreases by 30.14%,and from 0.2 to 0.4, the factor Nσ0decreases by 35.96%.Similarly, for GSI=10 and mi=30,with an increase in khfrom 0 to 0.2,the factor Nσ0decreases by 33.77%, and from 0.2 to 0.4, the factor Nσ0decreases by 38.77%.Furthermore, for GSI = 80 and mi= 5, with an increase in khfrom 0 to 0.2,the factor Nσ0decreases by 27.54%,and from 0.2 to 0.4, the factor Nσ0 decreases by 33.89%.

Figs.3 and 4 display the variation of the bearing capacity factor Nσ with the changes in khfor σci/(γB)=100 and 1000,respectively.Four different sub-plots are meant for four different values of GSI,i.e.20,40,60 and 80.It can be seen that in all the cases,the value of the factor Nσreduces invariably with an increase in the value of kh.

Fig.2.Variation of Nσ0 with kh for different mi values and (a) GSI =10, (b) GSI =20, (c) GSI =40, (d) GSI =60, (e) GSI = 80, and (f) GSI = 100.

To examine the effect of surcharge (q0)on the bearing capacity factor,several computations have been performed and these results have been presented in Figs.5-13 in terms of the variation of the factor Nqwith an increase in the value of kh.Similar to Nσ, a significant reduction in Nqwith an increase in khis observed.For example,in case of σci/(γB)= inf, q0/σci=0.25 and mi= 20,if the value of khchanges from 0 to 0.3,the corresponding magnitude of Nqdecreases by 49.22% and 45.06% corresponding to the value of GSI equal to 20 and 60,respectively.Fig.5 has been drawn to bring out the effect of σci/(γB)on the variation of Nqwith changes in kh.Three different values of σci/(γB),i.e.100,1000 and infinity(γ=0)have been used.Note that the magnitude of Nqincreases marginally with an increase in σci/(γB).The difference in the values of Nqfor σci/(γB) = 1000 and inf has been found to be much smaller than the corresponding difference between the values of Nqfor σci/(γB) = 100 and 1000.The variation of Nqfor σci/(γB)=100,q0/σci=0.25,with respect to khfor different values of mibetween 5 and 35, and GSI between 20 and 80 is shown in Fig.6.Likewise,Figs.7-9 present the results for σci/(γB)=100 corresponding to the value of q0/σciequal to 0.5, 0.75 and 1, respectively.Note that the magnitude of Nqdecreases marginally with an increase q0/σci, it should be mentioned that the product q0Nqalways increases with an increase in q0/σci.The magnitude of Nqincreases continuously with increases in miand GSI.Figs.10-13 provide the variation of Nqwith an increase in khfor σci/(γB)= inf corresponding to q0/σciequal to 0.25, 0.5,0.75 and 1,respectively.As indicated earlier,themagnitude of Nqfor σci/(γB)=inf becomes marginally greater than that with σci/(γB) = 100.

Fig.3.Variation of Nσ with kh for different mi values with (a) GSI = 20, (b) GSI = 40, (c) GSI = 60, and (d) GSI = 80.

Fig.4.Variation of Nσ with kh for σci/(γB) = 1000 for different mi values with (a) GSI = 20, (b) GSI = 40, (c) GSI = 60, and (d) GSI = 80.

Fig.5.Effect of σci/(γB)on Nq with changes in kh for(a)mi=15,GSI=20 and q0/σci=0.25;(b) mi=15,GSI=80 and q0/σci=0.25;(c)mi=30,GSI=20 and q0/σci=0.25;(d) mi=30,GSI=80 and q0/σci= 0.25; (e) mi=15, GSI=20 and q0/σci= 0.75; and (f) mi=15, GSI=80 and q0/σci= 0.75.

5.4.Failure patterns

In order to investigate how the failure patterns are influenced by the variation in the seismic coefficient with and without the presence of surcharge pressure,the proximity of the stress state to failure was defined in terms of a non-dimensional ratio,a/b,where a = (σ1-σ3)and b = (pσ1+r)α.This ratio is defined in a way such that it will attain a maximum value of 1 when the point is in the state of yield,otherwise,it will have a value smaller than unity for a non-yielded stress state.The failure patterns are shown in Figs.14 and 15.With σci/(γB) = 100, Fig.14 illustrates the failure patterns,without any surcharge pressure,for different values of kh,i.e.0,0.02,0.05 and 0.1 with GSI=30,mi=5,and GSI=80,mi=15.For kh= 0,the failure patterns become symmetrical about the center line passing through the footing base.In the presence of kh, the pattern becomes unsymmetrical.With an increase in kh,the extent of the plastic zone on one side of the footing,both horizontally and vertically, reduces continuously.The extent of the plastic zone for larger values of GSI and mibecomes continuously greater.Fig.15 illustrates the effect of surcharge on the shape of the yielded zone for different values of GSI, miand khand q0/σci=0.5 and 1.Once again, an introduction to the seismic force produces a non-symmetrical yielded zone, and the size of the yielded zone decreases continuously with an increase in surcharge pressure.Since the footing-rock interface is considered to be perfectly rough,in all the cases, a small non-plastic wedge is found to exist invariably below the footing base.

Fig.6.Variation of Nq with kh for σci/(γB)=100, q0/σci= 0.25 for different mi values with (a) GSI = 20, (b) GSI = 40, (c) GSI = 60, and (d) GSI = 80.

Fig.7.Variation of Nq with kh for σci/(γB)=100 and q0/σci= 0.5 for different mi with (a) GSI = 20, (b) GSI = 40, (c) GSI = 60, and (d) GSI = 80.

Fig.8.Variation of Nq with kh for σci/(γB)=100 and q0/σci= 0.75 for different mi with (a) GSI = 20, (b) GSI = 40, (c) GSI = 60, and (d) GSI = 80.

Fig.9.Variation of Nq with kh for σci/(γB)=100 and q0/σci= 1 for different mi with (a) GSI = 20, (b) GSI = 40, (c) GSI = 60, and (d) GSI = 80.

Fig.10.Variation of Nq with kh for σci/(γB)= inf and q0/σci= 0.25 for different mi values with (a) GSI = 20, (b) GSI = 40, (c) GSI = 60, and (d) GSI = 80.

Fig.11.Variation of Nq with kh for σci/(γB)= inf and q0/σci= 0.5 for different mi values with (a) GSI = 20, (b) GSI = 40, (c) GSI = 60, and (d) GSI = 80.

Fig.12.Variation of Nq with kh for σci/(γB)= inf and q0/σci= 0.75 for different mi values with (a) GSI = 20, (b) GSI = 40, (c) GSI = 60, and (d) GSI = 80.

Fig.13.Variation of Nq with kh for σci/(γB)= inf and q0/σci= 1 for different mi values with (a) GSI = 20, (b) GSI = 40, (c) GSI = 60, and (d) GSI = 80.

Fig.14.The collapse mechanisms for σci/(γB) = 100 and q0/σci=0 with(a)mi=5,GSI=30 and kh=0.0;(b)mi=5,GSI=30 and kh=0.02;(c)mi=5,GSI=30 and kh=0.05;(d)mi=5,GSI =30 and kh=0.1; (e) mi=15, GSI =80 and kh=0.0; and (f) mi=15, GSI =80 and kh=0.1.

6.Conclusions

The LB FELA in conjunction with the PCP has been employed to compute in an accurate fashion the seismic bearing capacity factors for a strip footing lying on rock media.The yielding of rock media has been modelled by the GHB criterion.No approximation is needed to make any assumption either about the geometry of the failure mechanism or smoothing of the yield criterion.The effects of the equivalent inertial forces on the bearing capacity factors, induced due to an occurrence of earthquake, are accounted for by means of equivalent pseudo-static horizontal seismic forces with the usage of the earthquake acceleration coefficient (kh).It has been observed that the bearing capacity factors Nσand Nqdecrease continuously with an increase in the magnitude of the earthquake acceleration coefficient (kh).The magnitudes of Nσ and Nqincrease continuously with an increase in the values of GSI and mi.The failure patterns in the presence of the earthquake acceleration become always un-symmetrical.The size of the plastic zone increases with an increase in the values of GSI and mi.The results provided in this study will be useful for designing strip footings on rock media in a seismically active zone.

Fig.15.The collapse mechanisms for σci/(γB) = 100 with(a)mi=5,GSI=30,q0/σci=0.5 and kh=0.0;(b) mi=5,GSI=30,q0/σci=0.5 and kh=0.02;(c) mi=15,GSI=80,q0/σci=0.5 and kh=0.0; (d) mi=15, GSI =80, q0/σci= 0.5 and kh=0.02; (e) mi=15, GSI =80, q0/σci=1 and kh=0.0; and (f) mi=15, GSI =80, q0/σci=1 and kh=0.1.

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

List of symbols

B Width of foundation

q0Surcharge pressure

γ Unit weight of rock mass

khHorizontal acceleration coefficient

g Acceleration due to gravity

QuVertical collapse load

σ1Major principal stress

σ3Minor principal stress

σciUniaxial compressive strength

mb, mi, s,α Hoek-Brown material constants

GSI Geological strength index

D Disturbance factor

LhHorizontal extents of the circular domain

LvVertical extents of the circular domain

σxNormal stress on x-plane

σyNormal stress on y-plane

τxyShear stress with respect in x-y plane

NiShape function

σnNormal stress

τ Shear stress

AequiMatrix containing the left hand side of all the equilibrium equations

AsdMatrix containing the left hand side of all the discontinuity equations

AbcMatrix containing the left hand side of all the boundary conditions

AeqMatrix containing the coefficients of all the equality constraints

AconeMatrix containing the coefficients of conic constraints

bequiVector containing the right hand side of all the equilibrium equations

bdisVector containing the right hand side of all the discontinuity equations

bbcVector containing the right hand side of all the boundary conditions

beqVector containing the right hand side of the equality constraints

bconeRight known vector of the conic constraints

Global unknown vector of stress variables

σ Global unknown vector of stress and slack variables

t An auxiliary variable

c Vector containing the objective function

NN Number of nodes

quUltimate bearing pressure

NσBearing capacity factor due to unit weight

Nσ0Bearing capacity factor of weightless media

NqBearing capacity factor due to surcharge

a, b Parameters which defines the state of stress

?nn-dimensional power cone

Rnn-dimensional real number