A simple soil model for low frequency cyclic loading

Computers and Geotechnics 84 (2017) 225–237 Contents lists available at ScienceDirect Computers and Geotechnics journal homepage: www.elsevier.com/locate/compgeo Research Paper A simple soil model for low frequency cyclic loading Yazen Khasawneh a,⇑, Antonio Bobet b, Robert Frosch b a b Geosyntec Consultants, Ann Arbor, MI 48105, USA Lyles School of Civil Engineering, Purdue University, West Lafayette, IN 47907, USA a r t i c l e i n f o Article history: Received 23 July 2016 Received in revised form 31 October 2016 Accepted 5 December 2016 Keywords: Integral abutments Soil-structure interaction Constitutive model Numerical simulation Physical model Large-scale test a b s t r a c t A three-dimensional elastoplastic soil constitutive model capable of capturing the response of granular soils under low-frequency cyclic loading is introduced and verified. The model is piecewise linear with a hyperbolic stress-strain relationship. The size of the hysteresis loop is controlled using different scaling factors with a shift in the backbone curve at load reversal. The model introduces a new algorithm to better capture the soil’s response upon reloading for plane strain. Model verification with experimental results at different scales shows that the model has good capabilities in capturing the response of granular soils under low frequency cyclic loading. Ó 2016 Elsevier Ltd. All rights reserved. 1. Introduction In classical formulations of geotechnical problems (such as bearing capacity), the displacements to reach plastic equilibrium are not considered, and the changes in stresses and stiffness with deformation are not generally included. As such, when displacements of a geotechnical structure are of concern or when the relative stiffness of the soil with respect to surrounding structures controls design, it becomes necessary to mathematically formulate the soil response to changes in stress or strain through a soil constitutive model. The mathematical formulation of the soil response requires understanding of the behavior and coupled response of stress changes and volumetric/pore pressure changes. Due to the degradation of the soil modulus with increasing stress/strain, the nonlinear response of the soil usually exhibits a hyperbolic relationship in the stress-strain space. The first formulation of a hyperbolic stressstrain curve for soils was based on Kondner’s argument that soil response to drained triaxial loading in compression can be approximated by a hyperbola [1]. The first implementations of a hyperbolic response of soils were based on either stress dependent stiffness [2] or strain dependent shear modulus [5,6]. The formulations of a hyperbolic relationship utilized the small strain modu- ⇑ Corresponding author. E-mail addresses: [email protected] (Y. Khasawneh), [email protected] (A. Bobet), [email protected] (R. Frosch). http://dx.doi.org/10.1016/j.compgeo.2016.12.003 0266-352X/Ó 2016 Elsevier Ltd. All rights reserved. lus; subsequent modifications recognized the dependency of the small strain modulus on the mean effective confining stress [7]. Formulation of a soil constitutive model that is capable of capturing the soil’s response under any loading combination, strain rate with various drainage conditions, and boundaries is a challenging task. The challenge arises from the need to model the multiphase nature of the soil structure (solids, voids with water and/or air) in a continuum formulation [10,1]. In general, soil constitutive models fall under three categories: The first category is simple elastic-perfectly plastic models. In elastic-perfectly plastic models, the initial loading curve follows Hooke’s law up to the yield stress. At yield, the soil obeys a yield criterion, e.g. Mohr-Coulomb. The advantage of simple elasticperfectly plastic models is the limited number of required parameters (E, t, /, c). In these models, however, the stress path dependency of the soil response cannot be captured; the associated volumetric deformation with changes of stress and the dependence of the stiffness on stress level are not considered [10]. At the other end of the spectrum of soil constitutive models (the third category) are advanced models that accurately capture the soil behavior regardless of the stress path and couple volumetric change with changes of stress level and with changes in soil stiffness. The advanced constitutive models either implement bounding surface plasticity such as MIT-S1 [12] or Multi Yield Surface Plasticity models [18]. The disadvantages of such advanced models are the large number of parameters that require specialized testing and the difficulty in implementing such models in a numerical framework by practitioners, limiting their use to research applications. Pestana 226 Y. Khasawneh et al. / Computers and Geotechnics 84 (2017) 225–237 et al. [13] suggested that seven (7) tests ranging from hydrostatic compression, to triaxial, to resonant column tests were required to obtain the thirteen (13) parameters needed to capture Toyoura sand behavior with the MIT-S1 model. As such, the need arises for relatively simpler soil constitutive models that are capable of capturing the soil response to specific loading combinations and for specified drainage and boundary conditions (the second category): Such models, however, may be only applicable to the particular conditions for which they are developed. The extension of constitutive models to capture cyclic behavior adds to model complexity as the response of the soil to cyclic loading exhibits hysteretic behavior. The well-known Massing rules developed in 1926 and subsequent modifications have been implemented in many constitutive models to capture the response of soils to irregular cyclic loading such as those generated by earthquakes [15,17,8]. In this study, a simple elastoplastic constitutive model is proposed to capture the soil response to cyclic loading from the expansion and contraction of Integral Abutment Bridges (IAB); that is, under small to moderate strains. The intent of the model is to reasonably capture the response of granular soils while maintaining the number of parameters to a minimum under these specific loading conditions. The proposed constitutive model is a modified version of what was proposed by Jung [8] to study the response of retaining walls to seismic loading. The modifications of the constitutive model were implemented based on the observed behavior of backfill and foundation soils of integral abutment bridges [9,4]. The modifications made to the Jung [8] model are: extension of the two dimensional model to three dimensions; rotation of the backbone curve to account for the increase in soil pressure from abutment movement with number of cycles; and capture of the soil response upon reloading. The model has been implemented in AbaqusÒ standard (UMAT) and explicit (VUMAT). Verification of the proposed constitutive model is performed using test results from element and physical model tests and a large-scale test of a laterally loaded pile. Verification and calibration of the model from largescale tests of an integral abutment bridge and from a full scale instrumented bridge will be presented in a future paper. 2. Development of a model for IAB structures The proposed model is developed as a piecewise linear rate independent model in the general elastoplasticity framework. The proposed model is isotropic with dependency of the shear modulus and small strain shear modulus on the octahedral shear strain and the mean effective stress, respectively. The strain increments are decomposed into elastic and plastic increments. The DruckerPrager (D-P) yield criteria with unassociated flow rule is adopted to identify the state of stresses at which plastic strains evolve. The relationship between stresses and strains is given based on the following equation (Eq. (1)): drij ¼ C eijkl deekl ð1Þ where drij = incremental stress tensor, C eijkl = elastic modulus tensor, and d = incremental elastic strain tensor. K¼ E 3ð1 2lÞ ; G¼ E 2ð1 þ lÞ _ total ¼ _ eij þ _ pij ij ð4Þ where d_ total = incremental total strain tensor, ij d_ eij = elastic strain tensor increment, and d_ pij = the plastic strain tensor increment. Hardin and Drnevich [5] identified several factors that control the degradation of the modulus of a granular soil during loading, namely strain amplitude, mean stress, void ratio, and others such as cementation. In addition, Hardin and Drnevich [6] proposed a modified hyperbolic stress-strain relation to better capture soil response which was achieved by distorting the strain scale by defining a hyperbolic relation for the strain as shown in Eqs. (5) and (6). G¼ ds 1 ¼ Go dc ð1 þ ch Þ2 ð5Þ Defining ch as: ch ¼ c ½1 þ a exp ð bðc=cr Þފ cr ð6Þ Defining cr as: cr ¼ sf ð7Þ G0 where G0 = small strain shear modulus, c = shear strain, cr = reference shear strain, a and b = fitting parameters, and sf = shear stress at failure. The small strain shear modulus dependency on the mean effective stress is given as follows by Eq. (8): Go ¼ Gref o r0m r0m;ref !a ð8Þ where r0m = mean effective stress, Gref o = reference small strain shear modulus, r0m;ref = reference effective mean stress, and ð2Þ The small strain shear modulus can be measured either directly by means of geophysical methods (e.g. cross hole), laboratory testing (resonant column), or estimated from empirical correlations. For example, Hardin [7] proposed the following relation for granular soils (PI = 0) with 0.4 < e < 1.2: ð3Þ Gref o ¼ 625 The elastic constants are written as shown in Eq. (2):
 E El dil djk þ dik djl þ 2ð1 þ lÞ ð1 þ lÞð1 During elastoplastic response, the total strain increment is equal to the summation of the elastic strain increment and the plastic strain increment, as shown in Eq. (4): a = constant (material-dependent). e ij C eijkl ¼ where E = Young’s modulus, l = Poisson’s ratio, K = bulk modulus, G = shear modulus, and dij = Kronecker delta and stands for 0 when i – j and for 1 when i = j. 2lÞ dij dkl qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 Pa r0m;ref 0:3 þ 0:7e2 ð9Þ Y. Khasawneh et al. / Computers and Geotechnics 84 (2017) 225–237 where e = void ratio, and Pa = atmospheric pressure. Upon load reversal (changing the direction of loading), the stress-strain relationship follows ‘‘rules” that control the initial stiffness upon loading/unloading and the size of the stress-strain hysteresis loop. The existing rules are mostly modifications of the original Masing rules of 1926. Modifications on the Masing rules were made primarily to overcome some of the shortcomings of the original rules. The Masing rules call for a symmetric backbone curve and for the initial stiffness for unloading/reloading equal to the initial stiffness of the backbone curve as shown in Eqs. (10a) and (10b) [15,17]. f bb ðcÞ ¼  f bb ðcÞ c>0 1  f bb ðcÞ c < 0 1 Gu ¼ Gr ¼ Go  1þ ðc 2 crev Þ ncr ð10aÞ ð10bÞ where fbb(c) = backbone curve function, c = shear strain, Gu and Gr = stiffness as a function of the shear strain of the unloading and reloading curves, respectively, crev = shear strain at the load reversal, cr = reference shear strain, Eq. (7), and n = 2 to satisfy Masing rules. Application of the Masing rules, as described in Eqs. (10a) and (10b), results in a symmetrical backbone curve and a stressstrain curve that is independent of the number of cycles. Furthermore, the Masing rules do not result in following the original backbone curve if the reloading/unloading is beyond the reversal point. In addition, if the unloading/reloading occurs within a cycle, then upon reloading, the model may predict a stress higher than the shear strength. Pyke [15] introduced the constant ‘‘n” in Eq. (10b), based on the ratio of the stress level to the shear strength. The estimate of ‘‘n” allows extension of the Masing rules to irregular cyclic loading (e.g. earthquakes) and this extension is referred to as the extended Masing rules. Using the same ‘‘n” for both loading and unloading results in cycles that will always be tangential to the backbone curve at the first reversal point. Tatsuoka et al. [17] suggested, based on the results of a series of cyclic loading tests on dry sand [11], a set of external and internal rules, and a ‘‘drag” rule for the backbone curves. The set of rules control the development of the hysteresis loop with unsymmetrical backbone curve for loading and unloading. The use of different scaling factors ‘‘n” for the unloading and reloading seems to be reasonable and provides a better prediction of behavior, as suggested by Jung [8], and is shown in Eqs. (11)–(15). G ¼ Go
C¼ 1 1 þ 1n jC Crev j ð11Þ 2  i coct h 1 þ a exp bðcoct =coct;r Þ coct;r Crev ¼  i coct;rev h 1 þ a exp bðcoct rev =coct;r Þ coct;r 2 3 coct ¼ ½ð11 ð12Þ ð13Þ 22 Þ2 þ ð22 33 Þ2 þ ð33 11 Þ2 0:5 þ6ð212 þ 223 þ 231 ފ ð14Þ 8 > <1 n ¼ n1 > : n2 9 Initial Loading > = unloading > ; reloading 227 ð15Þ where n = hysteresis loop scaling factor, coct;rev = octahedral shear strain at the point of load reversal, coct;r = reference octahedral shear strain, and  @u i ij ¼ 12 @u þ @xij , where ui are the displacements along the xi @xj axis. The Jung [8] formulation with different ‘‘n” for unloading/ reloading curves results in unloading/reloading curves that are not tangential to the original backbone curves, offset of the unloading/reloading curves from the reversal points that is symmetrical, and the difference in stress magnitude between the reversal points from cycle (n 1) and cycle (n) that is always constant and independent of the number of cycles. In the current study, modifications to the Jung [8] formulation are implemented to capture the behavior of granular soils under cyclic loading where it is observed that the deviatoric stress increase in the first few cycles in displacement controlled tests (controlled stress escalation with number of cycles) and that the soil response upon reloading, in plane strain, is initially close to a linear relationship in the stress-strain space [11,9,4]. Modifications include the implementation of a modified ‘‘drag” rule (after Tatsuoka et al., [17]) and forcing the initially ‘‘close to linear response” upon reloading. The modified ‘‘drag” rule consists of an offset of the backbone curve as a function of the incremental plastic strain and cycle number. The slope of the initially close to linear response during reloading is a function of the cycle number. The performance of the model with the proposed modifications is verified, as shown later, through a series of experiments conducted at two different scales and under different confining pressures. The concept of backbone shifting is illustrated in Fig. 1, which shows schematically the response to displacement controlled dynamic loading for the first two cycles, along with the offset of the backbone curve for the second, third, and fourth load reversals. For example, at the second load reversal, as shown in Fig. 1, the backbone curve is shifted in the unloading direction along the strain axis (identified by 2nd shift in the figure). Because the shear modulus degradation, Eqs. (11)–(13), is based on the octahedral shear strain, the offset of the backbone curve is also based on the octahedral shear strain. The formulation is given by Eq. (16). coct@load rev ersal 8 9 no shifting > > < 0; = ¼ 0 þ f ðe:p ; NÞ; unloading > : 0 f ðe: ; NÞ; reloading > ; p ð16Þ where e:p = plastic strain increment, and N = number of cycles. Because the function controlling the offset of the backbone curve is dependent of the direction of loading (unloading or reloading), a robust identification of the load reversal is required. Load reversal is recognized by the multiplication of two subsequent octahedral shear strain increments. If the result is negative, (c_ oct;i c_ oct;iþ1 < 0), where i represents the increment number, then there is as a load reversal. Once this is done, it is determined whether there is unloading ðc_ oct;iþ1 < 0Þ or reloading ðc_ oct;iþ1 > 0Þ to use the correct offset equation. 228 Y. Khasawneh et al. / Computers and Geotechnics 84 (2017) 225–237 Fig. 1. Offsetting of the backbone curve (2nd indicates the second shift corresponding to the 2nd load reversal, etc.). The second modification of the original formulation of Jung’s [8] model is inclusion of a ‘‘close to linear” response during the initial portion of the reloading curve (up to zero strain during reloading). It should be noted that this initial ‘‘close to linear response” upon reloading is dependent on the stress path. For example, this quasi-linear response is not observed in triaxial or direct simple shear cyclic tests of dry sand [14,19]. On the other hand, the quasi-linear response is observed in the plane strain cyclic tests on dry granular soils by Masuda et al. [11] (Fig. 2a). A similar observation was made from cyclic tests of dry granular soils using a laboratory apparatus that simulates bridge abutment movement [9,4] (Fig. 2b). The quasi-linear response may be attributed to the formation of the active wedge at the first unloading; upon reloading the energy imposed by the load is converted to work against pushing Fig. 2. Linear and hyperbolic response upon reloading. Y. Khasawneh et al. / Computers and Geotechnics 84 (2017) 225–237 the active wedge on the previously formed active failure surface (no new failure surfaces are formed). After the initial position is reached (displacement is zero), the energy from loading is converted against the attempt to form a new failure surface (passive wedge), thus kinetic work is utilized against displacement (shearing) of the sand matrix, which will result in a close to a hyperbolic response, as shown in Fig. 2a and b. It should be noted that when the strain upon reloading is greater than zero, the stress-strain relation becomes close to hyperbolic as shown in Fig. 2. The quasi-linear portion of the reloading curve can be captured well with a constant slope, as a fraction of the initial small strain shear modulus (Go). Eq. (17) present the relationships used for reloading. 8 1 Go 2 > > < ð1þ1njC Crev jÞ G¼ Go =C 1 > > 1 : Go 2 1 ð1þnjC Crev jÞ 9 Without initial linear response > > = With initial linear response; e < 0 > > With initial linear response; e > 0 ; ð17Þ where C1 = constant, e = strain in the loading direction, and Go ; C; and Crev are calculated from Eqs. (8), (12), and (13), respectively. In implementation of the proposed constitutive model, the use of linear or hyperbolic response upon reloading is determined by the user through an option selection in the model. The Druker-Prager (D-P) yield criterion with unassociated flow rule was selected to identify the state of stresses at which soil yields. The yield surface equation in the principal stress space (r0 1-r0 2-r0 3) is given by Eq. (18) in terms of the first stress invariant I1 (Eq. (19)), and the second deviatoric stress invariant J2 (Eq. (20)). F¼ 1 I1 tan a 3 pffiffiffiffiffiffiffi 3J 2 ð18Þ ð19Þ 1 Sij Sij 2 ð20Þ where a = D-P friction angle (/), j = D-P cohesion, and Sij = deviatoric stress tensor. The plastic potential is a function of the first stress invariant (I1), the second deviatoric stress invariant (J2), and the dilation angle (w), as shown in equation Eq. (21): Q¼ The model has eleven input parameters, five of which require calibration. The parameters are: (1) reference small strain shear modulus, ‘‘Gref o ”, (2) reference effective mean stress, ‘‘r0m;ref ”, (3) Poisson’s ratio, ‘‘t”, (4) friction angle, ‘‘/”, (5) cohesion, ‘‘c”, (6) and (7) the two modulus degradation parameters ‘‘a” and ‘‘b”, (8) and (9) the two scaling factors ‘‘n1”and ‘‘n2”, (10) ‘‘C1” the parameter that defines the quasi-linear reloading for plane strain loading, and (11) dilation angle, ‘‘w”. The a, b, n1 and n2 constants are fitting parameters that control the shape of the stress-strain curve. The constants selection is based on model calibration against laboratory testing. The recommended values for the constant ‘‘a” and ‘‘b” are between 0.5 and 0.1 and between 0.1 and 0.7, respectively. The choice of the scaling factors depends on the size of the hysteresis loop, with values ranging from about 1.7 to 2.3. It should be noted, however, that based on the calibration results, the selected constants could fall outside the given ranges. The model was implemented in a user-defined subroutine for AbaqusÒ Standard (UMAT) and for AbaqusÒ Explicit (VUMAT). For applications with low frequency cyclic loading (e.g. those resulting from thermal expansion and contraction of structures in contact with soil), the use of AbaqusÒ Standard is appropriate because inertia forces are relatively small. For moderate to high frequency (e.g. seismic loading) AbaqusÒ Explicit should be used. 3. Model verification j¼0 I1 ¼ rkk J2 ¼ 229 pffiffiffiffi 1 J 2 þ I1 tan w 3 The model is verified by comparing predictions from finite element simulations with results from different tests across different scales (10 1 ft through 102 ft). In this paper, the model predictions are compared with test results from element tests (10 1 ft), tests conducted on a physical model (100 ft), and with a pile load test (101 ft). Model verifications at other, larger, scales can be found in Khasawneh [9]. The following provides details of the comparisons. In general, the model shows acceptable performance in predicting stresses and displacements for all the tests, which brings confidence in the model and in its implementation in AbaqusÒ. ð21Þ The plastic strain increment is calculated from the equivalent plastic strain increment (depl ) and the derivative of the plastic potential function (Q) with respect to the deviatoric stress tensor as follows: pl Depl ij ¼ de dQ ¼ depl dSij ! pffiffiffi 3 1 1 pffiffiffiffi Sij þ tan wdij 2 3 J2 ð22Þ The equivalent plastic strain increment (depl ) is given as: 8 9 0 No yield > > > > < pl = depl ¼ De1 ðAll loading paths except direct shearÞ other wise > > > dcpl12 > : ; pffiffi ðShear loading pathÞ 3 ð23Þ Fig. 3. Model dimensions, applied confinement stress, and applied displacement of the cyclic drained triaxial test. 230 Y. Khasawneh et al. / Computers and Geotechnics 84 (2017) 225–237 Table 1 Soil model parameters-cyclic triaxial test simulation. Go,ref psf (N/m2) 5 rm,ref psf (N/m2) 7 9.4 * 10 (4.5 * 10 ) 3 5 6.3 * 10 (3.0 * 10 ) t a (deg) j psf (N/m2) 0.3 35 0.02 (1) a 0.06 b n1 n2 1.3 2.05 2.4 Two cyclic drained tests on granular soils were selected from the literature: a triaxial cyclic test by [19] and a direct simple shear cyclic test by Pradhan et al. [14]. Drained triaxial cyclic tests by [19] were conducted on Ottawa sand. One test with a relative density of approximately 40% and a confining stress of 3.0 kg/cm2 (6300 psf) was selected for comparison. Simulation was performed using the AbaqusÒ explicit software as a displacement controlled test with a VUMAT subroutine for the soil constitutive model. The dimensions of the model are 1 in. in diameter and 2 in. in length. The confining stress was applied in the first step, and then the deviatoric stress was applied incrementally by applying vertical displacement. Model dimensions, applied confinement, and displacement time history are presented in Fig. 3 which conform to the test conducted by [19]. The model was discretized using 8-node linear brick, general stress elements. The model parameters are either directly taken from the reference [19] or estimated based on accepted correlations (from Eq. (9)). Table 1 presents the model parameters used in the simulations. A comparison between the stress ratio (r1 =r3 ) versus strain data from the triaxial test conducted by [19] and the numerical simulation is presented in Fig. 4. It can be seen that the model captures well the stress-strain response of Ottawa Sand subjected to cyclic drained triaxial loading. Fig. 4. Comparison between the model predictions and the results from laboratory triaxial cyclic test on Ottawa Sand by [19]. Fig. 6. Comparison between the model and results from laboratory cyclic direct simple shear on Toyoura sand by Pradhan et al. [14]. 3.1. Element test Fig. 5. Model dimensions, applied confinement stress, and applied displacement for the cyclic direct simple shear test. Table 2 Soil model parameters-cyclic direct simple shear test simulation. Go,ref psf (N/m2) 5 rm,ref psf (N/m2) 7 8.8 * 10 (4.2 * 10 ) 3 5 2.1 * 10 (1.0 * 10 ) t a (deg) j psf (N/m2) 0.3 40 0.02 (1) a 0.5 b n1 n2 0.65 1.66 1.82 Y. Khasawneh et al. / Computers and Geotechnics 84 (2017) 225–237 and a confining stress of 1.0 kg/cm2 (2100 psf). Standard AbaqusÒ software was used with the UMAT subroutine. Dimensions of the model are 2 in. in length, width, and height. Confining stress was applied in the first step, and then the shear stress was applied incrementally by subjecting the sample to horizontal displacements. The DSS model dimensions, applied confining stress, and applied displacement-time history duplicate those of Pradhan et al. [14] and are presented in Fig. 5. Similar to the triaxial test simulation, the model was discretized using 8-node linear brick, general stress elements. The physical parameters were estimated from published correlations based on the initial void ratio and effective confining stress as reported by Pradhan et al. [14]. The selected parameters for the analysis are listed in Table 2. The stress-strain curve obtained from the simulation is plotted in Fig. 6 along with the test results as reported by Pradhan et al. [14]. As shown, the model captures well the initial backbone curve, the initial slopes of unloading and reloading, the increase of shear stress with number of cycles, and the modulus degradation with strain. Fig. 7. Testing apparatus [9]. Table 3 Laboratory test results for 430 Wedron Sand. a 231 3.2. Physical model test Property 430 Wedron Sand Ottawa Sanda Cu emax emin 1.7 0.77 0.49 1.7 0.78 0.48 El Mohtar et al. [3]. A Cyclic Direct Simple Shear (DSS) test by Pradhan et al. [14] was conducted on Toyoura sand. Simulation, as in the original test, was performed on a sample with a relative density of about 40% To investigate the response of backfill soils behind a bridge abutment to displacement cycles induced by the thermal expansion and contraction of an integral abutment bridge, a physical model was designed and built [9]. The specimen dimensions were selected based on numerical simulations to ensure that boundary effects on the specimen response were minimized (further details can be found in [9]). The soil specimens were 3.0 ft long and 1 ft high and wide and were prepared by placing three layers of sand that were compacted by tamping such that a relative density of Fig. 8. Test results from the physical model with smooth walls, with 0° and 45° skew. 232 Y. Khasawneh et al. / Computers and Geotechnics 84 (2017) 225–237 Fig. 9. Boundary conditions and contacts used in the numerical simulations. Table 4 Model parameters for the plates. Material E psf (N/m2) t Stainless steel 4.2 * 109 (2.0 * 1011) 0.3 30–35% was achieved. The testing apparatus is shown in Fig. 7. The physical model consisted of moving elements to simulate the abutment wall and fixed elements to simulate the ‘‘far field” boundary. Contacts between the moving and fixed elements were coated with TeflonÒ to minimize friction so that the applied load was transferred to the soil. As shown in Fig. 7, displacement was imposed on the plate that simulates an abutment wall using a hydraulic jack. Displacements were recorded by two Linear Position Transducers (LPTs) and the resulting force was recorded by a load cell. Deformation of the backfill soil surface was captured using two cameras and then analyzed using three dimensional Digital Image Correlation (DIC) software. The testing apparatus allows for various angles of the front plate (simulating different bridge skew angles), different roughness of the front plate (smooth or rough), and different loading-unloading sequences and amplitudes. To minimize friction between the soil and the walls of the testing apparatus, TeflonÒ sheets were installed on the inside faces of the plates. On the front plate, either TeflonÒ or sand paper was attached to achieve either a smooth or a rough surface, respectively. The backfill soil consisted of 430 Wedron (from Wedron, IL) uniform silica sand. This sand has physical and mechanical properties Table 5 Soil model parameters. Wall angle (0°) Go,ref psf (N/m2) rm,ref psf (N/m2) t a (deg) j psf (N/m2) a 0 45 1.0 * 105 (4.8 * 106) 1.1 * 102 (5.3 * 103) 0.3 41 1 (47.9) 0.25 0.3 b 0.1 0.3 n1 n2 C1 3.0 3.5 2.3 1.85 3.1 3.3 Y. Khasawneh et al. / Computers and Geotechnics 84 (2017) 225–237 233 Fig. 10. Comparison between simulations and laboratory measurements. Force and displacement results for 0° skew and smooth wall. similar to those of standard Ottawa sand (ASTM C778). This was verified through a series of laboratory tests which included Particle Size Gradation and maximum and minimum densities, as presented in Table 3. Results from the physical models tests with smooth walls at 0° and 45° skew in terms of force versus horizontal displacement are presented in Fig. 8 for the first, second, fifth, and tenth cycle. For a 300 ft long bridge with an abutment height of 8 ft, the estimated abutment movement is 0.8 in. for an 80 °F change in temperature, which corresponds to 0.8% movement of the abutment height. In the physical model test, the 0.8% movement of the plate height (1 ft) corresponds to 0.1 in. horizontal displacement, which is the displacement magnitude selected for most of the tests. The displacement field is uniform at the plate (the design of the physical model does not allow for plate rotation), while the stress tensor varies across the plate depending on the skew. It should be noted that the recorded force reflects the average stress on the plate. It can be seen from Fig. 8 that the force evolution versus horizontal displacement relation shows a nonlinear response with hysteresis for both wall configurations. The observed relationship between force and displacement follows a hyperbolic relation during loading and unloading, while during reloading a close to linear relation is observed initially then followed by a hyperbolic relationship. In addition, escalation of the force is observed with increasing cycle number up to the fifth cycle, followed by softening. A similar observation was made by Hassiotis and Xiong [20] who analyzed the data collected from an instrumented bridge in Massachusetts and found that the lateral earth pressure coefficient increased with displacement until it peaked, and then it went through a small decrease. Hassiotis and Xiong suggested that after mobilization of the peak internal friction angle and with more displacement, the sand internal friction angle started approaching the critical internal friction angle. Similarly, Tatsuoka and Ishihara [16] showed, for a drained strain controlled cyclic triaxial test on sand, that stress evolved with increasing number of cycles up to the fifth cycle and then a steady state stress was approached. The two tests discussed and presented in Fig. 8 (smooth wall with 0° and 45° skew) were simulated for verification of the proposed constitutive model. Comparisons between other tests (rough surface and other skew angles) and their simulations are not included as they yield similar conclusions. The results from other comparisons can be found in Khasawneh [9]. The simulations were conducted with AbaqusÒ standard and the UMAT subroutine. The model was discretized using 8-node linear brick, general stress elements. The numerical model dimensions are identical to those of the physical model. In the numerical simulations, the front and back plates of the apparatus (shown in Fig. 7) were modeled as steel plates, while the side and bottom plates were replaced by rollers to simulate the frictionless contact between the soil and the walls. The fixed element shown in Fig. 7 (the ‘‘far field” conditions at the back plate) were simulated as a fixed boundary. The AbaqusÒ frictionless surface-to-surface contact algorithm was used for the interface between the back plate and the soil as well as the soil and the front plate and the soil. The displacement-time history, which replicates the displacement-time history of the test, was applied to the front 234 Y. Khasawneh et al. / Computers and Geotechnics 84 (2017) 225–237 Fig. 11. Comparison between simulations and laboratory measurements. Force and displacement results for 45° skew and smooth wall. Table 6 Displacement history. Cycle # Amplitude (in.) 1 2 3 4 5 6 7 8 9 ±0.25 ±0.5 ±0.75 ±1.0 ±1.5 ±2.0 +2.5 +3.0 2.0/+4.0 Fig. 12. Schematic of the test pile configuration. Table 7 Model parameters for the pile. plate. The boundary conditions, contacts, and applied displacement-time history for the zero and 45-degree skew walls are presented in Fig. 9. The input parameters used for the models are listed in Tables 4 and 5. Note that the front and back plates are taken as linearelastic. The reference small strain shear modulus for the sand was estimated based on Eq. (9). The void ratio of the sand was obtained from the relative density (Dr  32%) and the maximum and minimum void ratios of the tested sand (emax = 0.77 and emin = 0.49). The internal friction angle (/) at Dr = 32% is approximately 31°. Relative density determination and sand characterization (emax, emin, and /) were based on minimum and maximum density tests and direct shear tests (details of the tests can be found in Khasawneh [9] and Frosch et al. [4]). For the simulation, the constant C1 was selected to best capture the observed ‘‘close to linear” response upon reloading as shown in Fig. 2. Part Material E psf (N/m2) t Pile Concrete filled steel tube transformed section EI 4.2 * 109 (2.0 * 1011) 0.2 It should be noted that the physical parameters are the same in the two simulations because the samples for the two tests were prepared in an identical manner (similar relative density). The nonphysical ‘‘fitting” parameters, however, are not the same. The differences between the two simulations (no skew and 45° skew) are explained given the test results from the physical model presented in Fig. 8. The first observation is that the force measured during initial loading is lower for the test with a skewed plate than the test with no skew (Fig. 8a). This observation indicates more stiffness degradation with displacement for the test with the 235 Y. Khasawneh et al. / Computers and Geotechnics 84 (2017) 225–237 Fig. 13. Numerical model dimensions. Fig. 14. Boundary conditions and displacements at the top of the pile. skewed wall than for the test with no skew resulting in different values for the fitting parameters a and b. During the unloading/ reloading hysteresis cycles, the test with a skewed plate demonstrated a larger hysteresis loop during unloading and a smaller loop during reloading than for the test with no skew (Fig. 8c and d). For the skewed test, therefore, the scaling factors are larger during unloading (n1) and smaller during reloading (n2), as compared to the scale factors for the test with no skew. During the close to linear portion of the reloading phase, the force versus displacement curves for both tests are almost parallel resulting in similar values for the constant C1. The simulations consist of two steps: The first step generates the initial stresses in the sand through a geostatic solution, while the second step imposes the displacement time history (10 cycles of approximate amplitude 0.1 in.) to the front plate. The results from the AbaqusÒ simulation along with the results from the test with a smooth 0° skew plate and a smooth 45° skew plate wall are presented in Figs. 10 and 11, respectively. As shown in Fig. 10, the model approximates reasonably well the response of the sand due to horizontal cyclic loading imposed by the front plate. The initial loading curve (initial backbone curve) is captured well (Fig. 10a). In addition, rotation of the backbone curve during successive cycles is reproduced well by the model (Eq. (16)) for all cycles, which is supported by the fact that the model is able to duplicate the maximum and minimum forces on the wall during each cycle. Capturing rotation of the backbone curve enables the model to predict the observed force escalation with increasing number of cycles. As discussed previously, force escalation is only observed in the first five cycles, which is followed by softening of the response. The proposed model, however, does not include any rules to account for the softening experienced after the fifth cycle. The softening observed between the fifth and tenth cycles for the test with a smooth wall and 0° skew was small (<10% of the maximum force recorded in the fifth cycle) and the model is able to capture the response up to the tenth cycle (Fig. 10d). The use of different scaling factors for unloading/reloading (n1 and n2 in Eq. (15)) enables the model to capture well the observed sizes of the hysteresis loops with at most a 5% difference in the area between the observed and predicted hysteresis loops. As shown in Fig. 11, the model predicts well the response of sand to horizontal cyclic loading with a 45° skewed smooth wall up to the fifth cycle (Fig. 8a–c); however, it overpredicts the force at the tenth cycle by about 25%, as shown in Fig. 11d. The performance of the model in capturing the response of the initial loading curve (initial backbone curve), in capturing the force escalation up to the fifth cycle, and in predicting the size of hysteresis loop is similar to that of the simulation of the test with a 0° skew and smooth wall. However, the observed softening in the response between the fifth and tenth cycles is more significant for the case with a 45° skew. Because no rules are included in the model to account for softening, the model continues to escalate the forces. Due to this discrepancy between the actual behavior and the model, the accuracy of the model for a rough wall beyond five cycles reduces. For integral abutment bridges, the conditions in the field are not as severe as those tested in the laboratory since actual measurements of instrumented IABs [9,4] indicate that the backfill response reaches steady state (no increase in the earth pressure behind the IAB) after the first few years of being subjected to thermal loading. As such, the model is capable of capturing the backfill response behind a skewed rough abutment wall during the period of anticipated increase in the earth pressure. 3.3. Large-scale test on a laterally loaded pile A 6 in. diameter, 12.5 ft long concrete filled tube (CFT) pile was driven in a granular soil deposit at the Bowen Laboratory for LargeScale Civil Engineering Research of the Lyles School of Civil Engineering at Purdue University. The steel pipe with a 0.25 in. wall thickness was driven first. A 1 in. diameter steel pipe was also installed at the center of the pipe to allow for installation of instrumentation to record lateral deformations of the pile with depth. Table 8 Soil model parameters for the foundation soils. Go,ref psf (N/m2) 5 rm,ref psf (N/m2) 7 8.5 * 10 (4.1 * 10 ) a 1250 6 * 10 4 t a (deg) j psf (N/m2) 0.3 41 1 (47.9) Not applicable because the linear option upon reloading was not used for the foundation. a 0.35 b N1 N2 C1 0.1 2.0 1.8 NAa 236 Y. Khasawneh et al. / Computers and Geotechnics 84 (2017) 225–237 The 6 in. diameter steel pipe was subsequently filled with 5000 psi concrete. Subsurface conditions at the site were determined based on borings that were drilled for construction of the Laboratory. The site is comprised of an upper layer of medium dense to dense sand with a blow count of 25 blows/ft down to an approximate depth of 20 ft. The sand layer is underlain by medium compact silt (N  30 blows/ft) down to an approximate depth of 42 ft. A very compact silt layer (blow counts greater than 50) is encountered underneath extending down to the explored depth of 51 ft. The water table was encountered at a depth of about 9 ft. A pile cap and a thrust block were cast to allow for lateral loading of the pile. Lateral loads were applied through the use of a hydraulic actuator. Various sensors were installed to monitor displacement of the pile cap, applied pressure to the hydraulic ram and pile deformation with depth. A schematic showing the test configuration, and instrumentation is presented in Fig. 12. Displacements were applied at the centroid of the pile cap by the hydraulic actuator, as shown in Fig. 12. The displacement history is presented in Table 6 with positive values for actuator extension and negative values for actuator retraction. Displacements and pressures were recorded continuously while the deformed shape of the pile was measured at the maximum, minimum, and zero amplitude during each cycle. Details of the test, instrumentation, test procedure, and results can be found in Khasawneh [9] and Frosch et al. [4]. The pile test was simulated to achieve two objectives: (1) verify model performance at a scale larger than the element or laboratory (over a range of about one order of magnitude of effective stresses in the soil), and (2) calibrate the soil constitutive model for a fullscale bridge simulation (to be discussed in a future paper). The pile was simulated in three-dimensional space with AbaqusÒ standard. The soil was assumed to obey the elastoplastic model developed here, which is implemented as a user-defined material model (UMAT) in AbaqusÒ, while the pile was assumed as elastic. The pile and the surrounding soil are discretized using 4-node linear tetrahedral, general stress elements. The section properties of the modeled pile are based on the calculated composite effective area (7.8 in.2) and moment of inertia (24.9 in.4) (Frosch and Lovell, [21]). The modeled pile length was identical to the tested pile (12.5 ft in depth below ground and 2.5 ft of stick up above ground). The numerical domain dimensions were selected such that the boundaries are at least 5.25 pile diameters from the center of the pile as shown in Fig. 13. The applied displacement history was identical to that used in the field test (Table 6). Rollers were assigned to the sides of the numerical domain, and a fixed boundary was applied at the bottom of the model. The boundary conditions and displacement-time history are shown in Fig. 14. The contacts and constraints are considered to replicate what is believed to exist in the actual test. There are few contact algorithms available in AbaqusÒ standard, mostly based on classical Coulomb friction. The selected frictional contact algorithm between the pile shaft and soil allows for detachment at the contact. A tie constraint was applied between the tip of the pile and the foundation soil. Model parameters for both the pile and soil were determined. The selected elastic parameters for the pile are provided in Table 7 while the elastoplastic model parameters for the foundation soils are presented in Table 8. The selected parameters are based on the subsurface conditions at the location of the tested pile. Based on the blow counts of the foundation soils and the effective vertical stress, the relative density was estimated to be approximately 50%. The initial void ratio and the internal friction angle were also estimated as 0.25 and 31°, respectively. The reference small strain shear modulus was estimated based on Eq. (9). The simulation was conducted in two steps: The first step generated the initial stresses in the soil while the second step applied the displacement-time history. The deformed shape of the pile at 0.25 in., 0.5 in., 1.0 in., and 2.0 in. amplitude are plotted along with the test measurements in Fig. 15. As shown in Fig. 15, the model is able to capture the pile response under lateral loading. The point of fixity is captured very Fig. 15. Comparison of the pile’s deformed shape between simulations and experimental results. Y. Khasawneh et al. / Computers and Geotechnics 84 (2017) 225–237 well for all simulated displacement amplitudes. In addition, the deformed shape is reproduced well which leads to a reasonable estimate of both the shear forces and bending moments in the pile. 4. Summary and conclusions A three-dimensional, hyperbolic, elastoplastic soil constitutive model with the Druker-Prager (D-P) yield criterion and modified Masing rules for unloading and reloading was developed to capture the backfill response of integral abutment bridges to thermal loading/unloading cycles. Modifications of the Masing rules include different initial slopes and scaling factors for the unloading and reloading of the soil, with ‘‘shifting” of the backbone curve upon load reversals. The different initial slopes for the reloading include ‘‘a close to linear response” initially for specific reloading conditions. The soil model was implemented as a subroutine in the explicit as well as standard AbaqusÒ. The model’s ability to predict soil behavior at different scales was verified by performing a series of numerical simulations at various scales. Model verification was conducted by comparing the results of the analysis with the results from triaxial and direct simple shear cyclic tests, with the results of cyclic loading tests on backfill soils tested in the laboratory to simulate the movement of the abutment wall of an integral abutment bridge, and with the results of a large-scale test on a pile subjected to lateral loading. Based on the results from the simulations and their comparison with the actual observations, it can be concluded that the proposed constitutive model performs well under the range of stress levels and length scales evaluated here. Acknowledgment This research was supported by the Indiana Department of Transportation (INDOT) through the Joint Transportation Research Program (JTRP), under Grant No. SPR-3318. This support is gratefully acknowledged. References [1] Brinkgreve RB. Selection of soil models and parameters for geotechnical engineering application. In: Yamamuro Jerry A, Kaliakin Victor N, editors. GeoFrontiers 2005. ASCE; 2005. p. 69–98. 237 [2] Duncan JM, Chang CY. Nonlinear analysis of stress and strain in soils. J Soil Mech Found Div ASCE 1970;96(SM5):1629–53. [3] El Mohtar CS, Bobet A, Santagata MC, Drnevich VP, Johnston CT. Liquefaction mitigation using bentonite suspensions. J Geotech Geoenviron Eng 2013;139 (8):1369–80. [4] Frosch RJ, Bobet A, Khasawneh Y. Reduction of bridge construction and maintenance costs through coupled geotechnical and structural design of integral abutment bridges. Joint Transportation Research Program Publication No. FHWA/IN/JTRP-2014/06. West Lafayette, IN: Purdue University; 2014. http://dx.doi.org/10.5703/1288284315500. [5] Hardin BO, Drnevich VP. Shear modulus and damping in soils: measurement and parameter effects. J Soil Mech Found Div ASCE 1972;98(SM6):603–24. [6] Hardin BO, Drnevich VP. Shear modulus and damping in soils: design equations and curves. J Soil Mech Found Div ASCE 1972;98(SM7):667–92. [7] Hardin BO. The nature of stress-strain behavior of soils. In: Proceedings of earthquake engineering and soil dynamics of foundations division, ASCE Pasadena, California, vol. 1; 1978. p. 3–89. [8] Jung C. Seismic loading on earth retaining structures. Thesis, presented to Purdue University at West Lafayette, IN, in partial fulfillment of the requirements for the degree of Doctor of Philosophy; 2009. [9] Khasawneh Y. Soil Structure interaction of integral abutment bridges. Thesis, presented to Purdue University at West Lafayette, IN, in partial fulfillment of the requirements for the degree of Doctor of Philosophy; 2014. [10] Lade P. Overview of constitutive models for soils. In: Yamamuro Jerry A, Kaliakin Victor N, editors. GeoFrontiers 2005. ASCE; 2005. p. 1–34. [11] Masuda T, Tatsuoka F, Yamada S, Sato T. Stress-strain behavior of sand in plane strain compression, extension and cyclic loading tests. Soils Found 1999;39:31–45. [12] Pestana JM, Whittle AJ. Formulation of unified constitutive model for clays and sands. Int J Numer Anal Meth Geomech 1999;23:1215–43. [13] Pestana JM, Whittle AJ, Salvati LA. Evaluation of constitutive model for clays and sands: part I – sand behavior. Int J Numer Anal Meth Geomech 2002;26:1097–121. [14] Pradhan TB, Tatsusoka F, Sato Y. Experimental stress-dilatancy relations of sand subjected to cyclic loading. Soils Found 1989;29(1):45–64. [15] Pyke R. Nonlinear soil models for irregular cyclic loadings. J Geotech Eng Div ASCE 1979;105(GT6):715–25. [16] Tatsuoka F, Ishihara K. Drained deformation of sand under cyclic stresses reversing direction. Soils Found 1974;14(3):51–65. [17] Tatsuoka F, Masuda T, Siddiquee MS, Koseki J. Modeling the stress-strain relations of sand in cyclic plane strain loading. J Geotech Geoenviron Eng ASCE 2003;129:450–67. [18] Yang Z, Elgamal A, Parra E. Computational model for cyclic mobility and associated shear deformation. J Geotech Geoenviron Eng, ASCE 2003;129:1119–27. [19] Yu HS, Khong C, Wang J. A unified plasticity model for cyclic behavior of clay and sand. Mechanics Research Communication 2007;34:97–114. [20] Hassiotis S, Xiong K. Deformation of cohesionless fill due to cyclic loading. The 2005 UTRC Research Initiative Region II Publication No. SPR ID# C-05-03. Hoboken, NJ: Stevens Institute of Technology; 2007. http://www.utrc2.org/ publications/deformation-cohesionless-fill-due-cyclic-loading-0. [21] Frosch RJ, Lovell M. Long term behavior of integral abutment bridges. Joint Transportation Research Program Publication No. FHWA/IN/JTRP-2011/16. West Lafayette, IN: Purdue University; 2011. http://dx.doi.org/10.5703/ 1288284314640.