# GEOTECHNICAL INNOVATION FOR TRANSPORT INFRASTRUCTURES

EDITED BY : Sanjay Shrawan Nimbalkar, Prabir K. Kolay and Yifei Sun PUBLISHED IN : Frontiers in Built Environment

#### Frontiers eBook Copyright Statement

The copyright in the text of individual articles in this eBook is the property of their respective authors or their respective institutions or funders. The copyright in graphics and images within each article may be subject to copyright of other parties. In both cases this is subject to a license granted to Frontiers. The compilation of articles constituting this eBook is the property of Frontiers.

Each article within this eBook, and the eBook itself, are published under the most recent version of the Creative Commons CC-BY licence. The version current at the date of publication of this eBook is CC-BY 4.0. If the CC-BY licence is updated, the licence granted by Frontiers is automatically updated to the new version.

When exercising any right under the CC-BY licence, Frontiers must be attributed as the original publisher of the article or eBook, as applicable.

Authors have the responsibility of ensuring that any graphics or other materials which are the property of others may be included in the CC-BY licence, but this should be checked before relying on the CC-BY licence to reproduce those materials. Any copyright notices relating to those materials must be complied with.

Copyright and source acknowledgement notices may not be removed and must be displayed in any copy, derivative work or partial copy which includes the elements in question.

All copyright, and all rights therein, are protected by national and international copyright laws. The above represents a summary only. For further information please read Frontiers' Conditions for Website Use and Copyright Statement, and the applicable CC-BY licence.

ISSN 1664-8714 ISBN 978-2-88963-841-3 DOI 10.3389/978-2-88963-841-3

#### About Frontiers

Frontiers is more than just an open-access publisher of scholarly articles: it is a pioneering approach to the world of academia, radically improving the way scholarly research is managed. The grand vision of Frontiers is a world where all people have an equal opportunity to seek, share and generate knowledge. Frontiers provides immediate and permanent online open access to all its publications, but this alone is not enough to realize our grand goals.

#### Frontiers Journal Series

The Frontiers Journal Series is a multi-tier and interdisciplinary set of open-access, online journals, promising a paradigm shift from the current review, selection and dissemination processes in academic publishing. All Frontiers journals are driven by researchers for researchers; therefore, they constitute a service to the scholarly community. At the same time, the Frontiers Journal Series operates on a revolutionary invention, the tiered publishing system, initially addressing specific communities of scholars, and gradually climbing up to broader public understanding, thus serving the interests of the lay society, too.

#### Dedication to Quality

Each Frontiers article is a landmark of the highest quality, thanks to genuinely collaborative interactions between authors and review editors, who include some of the world's best academicians. Research must be certified by peers before entering a stream of knowledge that may eventually reach the public - and shape society; therefore, Frontiers only applies the most rigorous and unbiased reviews.

Frontiers revolutionizes research publishing by freely delivering the most outstanding research, evaluated with no bias from both the academic and social point of view. By applying the most advanced information technologies, Frontiers is catapulting scholarly publishing into a new generation.

#### What are Frontiers Research Topics?

Frontiers Research Topics are very popular trademarks of the Frontiers Journals Series: they are collections of at least ten articles, all centered on a particular subject. With their unique mix of varied contributions from Original Research to Review Articles, Frontiers Research Topics unify the most influential researchers, the latest key findings and historical advances in a hot research area! Find out more on how to host your own Frontiers Research Topic or contribute to one as an author by contacting the Frontiers Editorial Office: researchtopics@frontiersin.org

# GEOTECHNICAL INNOVATION FOR TRANSPORT INFRASTRUCTURES

Topic Editors:

Sanjay Shrawan Nimbalkar, University of Technology Sydney, Australia Prabir K. Kolay, Southern Illinois University Carbondale, United States Yifei Sun, Ruhr University Bochum, Germany

Citation: Nimbalkar, S. S., Kolay, P. K., Sun, Y., eds. (2020). Geotechnical Innovation for Transport Infrastructures. Lausanne: Frontiers Media SA. doi: 10.3389/978-2-88963-841-3

# Table of Contents

	- K. Krishnan, Koushik Halder and Debarghya Chakraborty

Sudhir K. Tripathi, Prabir K. Kolay, Vijay K. Puri and Sanjeev Kumar

# Editorial: Geotechnical Innovation for Transport Infrastructures

Sanjay Nimbalkar <sup>1</sup> \*, Prabir K. Kolay <sup>2</sup> and Yifei Sun3,4

*<sup>1</sup> School of Civil and Environmental Engineering, University of Technology Sydney, Sydney, NSW, Australia, <sup>2</sup> Department of Civil and Environmental Engineering, Southern Illinois University, Carbondale, IL, United States, <sup>3</sup> Key Laboratory of Ministry of Education for Geomechanics and Embankment Engineering, Hohai University, Nanjing, China, <sup>4</sup> Faculty of Civil and Environmental Engineering, Ruhr University Bochum, Bochum, Germany*

Keywords: geotechnical engineering, ground improvement, transportation infrastructures, road, railway

**Editorial on the Research Topic**

**Geotechnical Innovation for Transport Infrastructures**

# BACKGROUND

The transportation and transit systems are crucial parts of a nation's economy playing major roles in efficiently managing the conveyance for passengers, mail, or freight. However, the service life performance of these systems is impeded by deterioration of key components including road/rail, base, and subgrade due to ever-growing traffic. Providing safe, efficient, and sustainable transport infrastructure is often challenging owing to complex geotechnical aspects of the ground. This Research Topic is designed to accommodate the interests of engineers and researchers in modeling and designing both geomaterials and geostructures that could improve performance, sustainability, and life cycle of transportation infrastructures. The topic involves a wide coverage of timely issues on technologies and innovations focusing on broad aspects of geotechnical innovations in order to address global grand challenges and UN's sustainable development goals with great social and economic importance. We are proud to present this topic containing 10 peer-reviewed contributions providing insights into geotechnical innovations crucial for sustainable transport in many continents, including Asia, Australia, Europe, and North America.

#### Edited and reviewed by:

*Sakdirat Kaewunruen, University of Birmingham, United Kingdom*

#### \*Correspondence: *Sanjay Nimbalkar Sanjay.Nimbalkar@uts.edu.au*

#### Specialty section:

*This article was submitted to Transportation and Transit Systems, a section of the journal Frontiers in Built Environment*

> Received: *21 April 2020* Accepted: *04 May 2020* Published: *29 May 2020*

#### Citation:

*Nimbalkar S, Kolay PK and Sun Y (2020) Editorial: Geotechnical Innovation for Transport Infrastructures. Front. Built Environ. 6:78. doi: 10.3389/fbuil.2020.00078*

# MAJOR HIGHLIGHTS OF CONTRIBUTIONS

Pile foundations supporting high-rise buildings and transport infrastructures are often subjected to eccentric lateral cyclic load arising from the action of wind, waves, high speed traffic, ship impacts etc. Such lateral load can induce torsion in pile and lead to progressive degradation of the soil strength and the axial pile capacity as reported through boundary element modeling (BEM) approach by Nimbalkar et al.. The results from this study can be utilized for formulating the design criteria for pile subjected to axial and torsional cyclic loads relevant to transport environment.

The granular soils including ballast usually experience dynamic loads in the field. The existing models are complex and can mainly be used to predict the monotonic stress-strain response of granular soils. However, the cyclic fractional constitutive model developed by Li et al. is simple and is able to capture the cyclic (repetitive) loading representative of traffic. By employing the discrete element modeling (DEM) approach in PFC 3D, Dahal and Mishra highlighted the significance of particle breakage in permanent deformation accumulation of ballast layer under cyclic loading. The gravel loss is a major limitation for unsealed roads and effects of three significant factors, material properties, traffic, and weather condition are assessed through extensive field investigations by Pardeshi et al.. Such gravel loss can be reduced by promoting the use of modified gravel than conventionally used material.

Hou et al. determined the material properties of aggregate quarry by-products through laboratory testing and materials characterization, and evaluated the feasibility of their use in various pavement applications. The study promotes the potential use of these materials with fly ash stabilization as sustainable pavement construction strategies. The laboratory study by Tripathi et al. evaluates the liquefaction characteristics of 20– 30 Ottawa sand with various percentages of monofilament polypropylene fibers through a series of stress-controlled cyclic triaxial tests. Test findings revealed a significant improvement in liquefaction resistance when the polypropylene fiber content exceeded beyond 0.075% at 34.47 kPa effective confining stress. In the case of earthquake, bearing capacity assessment of a strip footing placed over an embankment of anisotropic clay is crucial as demonstrated by Krishnan et al.. The adoption of finite element limit analysis (FELA) approach combined with second-order conic optimization technique can enable correct assessment of such a strip footing used for various infrastructure projects including bridge abutments, roadways, and railways on flood plains or riverbanks.

The selection of a wide variety of subgrades, geogrid characteristics and placement locations in the study by Jiang and Nimbalkar have enabled the comprehensive assessment of geogrid-reinforced railroad under monotonic and cyclic loads. The results of finite element modeling (FEM) study in PLAXIS 2D reveals the more prominent role of geogrid for railroad located over the poor quality subgrade. The correct placement of geogrid at the ballast-sub-ballast interface as well as increasing the number of geogrids are found to bear favorable implications on track stability. The inclusion of 3D geoinclusion such as geocell in comparison with planar geogrid is found attractive in mitigation of traffic induced vibrations through finite difference modeling approach in FLAC 3D (Hegde and Venkateswarlu). The efficacy of geocell is more when it is placed at the shallow depth. Sun has investigated the potential use of recycled rubber tires for railroad using a 3D FEM implemented in ABAQUS. The rubber tire is found to radially confine the infilled materials effectively, and thus reduce excessive lateral spreading and vertical displacement.

# SUMMARY

The salient contributions amply demonstrate the useful innovations in the area of geotechnical engineering, constitutive, and numerical modeling as well as ground improvement relevant to transportation infrastructures around the globe. We, editors, hope that you enjoy reading these articles and the Research Topic, as a whole.

# AUTHOR CONTRIBUTIONS

SN prepared the draft. PK and YS edited the draft.

**Conflict of Interest:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2020 Nimbalkar, Kolay and Sun. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Piles Subjected to Torsional Cyclic Load: Numerical Analysis

Sanjay Shrawan Nimbalkar <sup>1</sup> , Piyush Punetha<sup>2</sup> , Sudip Basack <sup>3</sup> \* and Mehdi Mirzababaei <sup>4</sup>

*<sup>1</sup> Faculty of Engineering and Information Technology, School of Civil and Environmental Engineering, University of Technology Sydney, Ultimo, NSW, Australia, <sup>2</sup> Geotechnical Engineering Division, CSIR-Central Building Research Institute, Roorkee, India, <sup>3</sup> Department of Civil Engineering, Kaziranga University Jorhat, Jorhat, India, <sup>4</sup> School of Engineering and Technology, Central Queensland University, Melbourne, VIC, Australia*

Pile foundations supporting large structures (such as high-rise buildings, oil drilling platforms, bridges etc). are often subjected to eccentric lateral load (in addition to the vertical loads) due to the action of wind, waves, high speed traffic, and ship impacts etc. The eccentric lateral load, which is usually cyclic (repetitive) in nature, induces torsion in the pile foundation. This paper presents a numerical model based on boundary element approach to study the performance of a single pile subjected to the torsional cyclic load. The model is initially validated by comparing it with the experimental data available from the literature. Thereafter, the model has been utilized to conduct a parametric study to understand the influence of the torsional cyclic loading parameters on the axial pile capacity. The results indicated that the model is able to capture the degradation in the axial pile capacity due to the torsional cyclic loading with a reasonable accuracy. Moreover, the parametric study showed that the frequency, amplitude and number of cycles play a significant role in the torsional cyclic response of the pile. The present study is essential for the development of design guidelines for pile foundations subjected to torsional cyclic load.

#### Edited by:

*Sakdirat Kaewunruen, University of Birmingham, United Kingdom*

#### Reviewed by:

*Rims Janeliukstis, Riga Technical University, Latvia Zhang Qian, Qingdao University, China*

> \*Correspondence: *Sudip Basack basackdrs@hotmail.com*

#### Specialty section:

*This article was submitted to Transportation and Transit Systems, a section of the journal Frontiers in Built Environment*

> Received: *20 December 2018* Accepted: *13 February 2019* Published: *08 March 2019*

#### Citation:

*Nimbalkar SS, Punetha P, Basack S and Mirzababaei M (2019) Piles Subjected to Torsional Cyclic Load: Numerical Analysis. Front. Built Environ. 5:24. doi: 10.3389/fbuil.2019.00024* Keywords: piles and piling, lateral load, torsion, mathematical modeling, foundations

# INTRODUCTION

The construction of large structures such as high-rise buildings, oil drilling platforms, electrical transmission towers, wind turbines, bridges, and railway embankments over a soft compressible clay poses a serious challenge to the design engineers. Therefore, a cost-effective foundation system with an acceptable degree of safety is required and consequently, these structures are usually supported by pile foundations. These piles are often subjected to cyclic (repetitive) lateral loads in addition to the vertical loads, during their service life due to wind, sea waves, high speed train traffic and ship impacts etc. (Arshad and O'Kelly, 2016; Haiderali and Madabhushi, 2016). The cyclic lateral load often acts eccentrically and generates cyclic torsion on the pile foundation (Barker and Puckett, 1997). The inadequate design of these piles against the cyclic torsional loading may affect the safety and serviceability of the foundation, leading to disastrous consequences (Vickery, 1979; Barker and Puckett, 1997). Moreover, in case of foundations supporting offshore platforms, the amplitude (magnitude) of the cyclic torsional load is high and is usually accompanied with low frequency vibrations, while the reverse is true for onshore structures including railways (Nimbalkar and Indraratna, 2016).

The cyclic load initiates the reversal of soil-pile interface shear stress, thus producing a progressive deterioration in strength and stiffness of the surrounding soil and consequently, pilesoil interactive performance undergoes significant degradation (Basack, 2015). Such degradation may reduce the load carrying capacity of the pile and increase the pile head displacement (settlement). The primary reason for such degradation is the gradual reorientation and rearrangement of the soil particles adjacent to the pile-soil interface. The other reasons include the generation of excess pore water pressure and the development of irrecoverable plastic deformation of soil adjacent to the pile soil interface (Basack and Dey, 2012).

Several theoretical, laboratory and field-based investigations on the pile foundations subjected to vertical and lateral cyclic loads have been conducted in the past (El Naggar et al., 1998; Cavey et al., 2000; Taciroglu et al., 2006; Basack, 2010a; Jardine and Standing, 2012). Moreover, the behavior of pile subjected to static torsion has been investigated by several researchers (e.g., Stoll, 1972; Poulos, 1975; Kong and Zhang, 2008; Misra et al., 2014; Chen et al., 2016). However, the studies pertaining to the influence of torsional cyclic load on the behavior of pile foundation are rather limited. It has been found that the torsional load can significantly influence the load carrying capacity of the pile. The skin friction mobilized due to the axial load interacts with the circumferential pile-soil interface shear stress (induced due to the torsional load) and consequently, the axial pile capacity decreases and the axial displacement increases (Basack and Sen, 2014a). Therefore, the investigation of the behavior of pile subjected to cyclic torsional loading becomes imperative for the safe design and satisfactory long-term performance of the pile foundations supporting the large offshore or onshore structures.

The behavior of pile foundation subjected to monotonic and cyclic axial loading can be investigated by using several numerical and analytical techniques available in the literature. These include the dynamic response analysis, cyclic stability analysis, finite element and boundary element analysis etc. to name a few (Poulos, 1982, 1988; Bea, 1992; Basack and Dey, 2012; Fatahi et al., 2014). In the present paper, a novel numerical methodology based on boundary element modeling (BEM) is proposed to capture the response of a single, vertical, floating pile subjected to combined axial, and torsional cyclic loads. The non-linear stress-strain response of the soil is incorporated in the model by using a hyperbolic function. Moreover, the constitutive behavior of the pile material is assumed to be elastic-perfectly plastic. The effect of progressive degradation of soil strength and stiffness under interface shear stress reversal (or simply cyclic loading) has been incorporated by using an exponential correlation and a semi-logarithmic rate function. The analysis could also have been conducted using finite element method which might have required the use of three-dimensional meshes to represent the pile, the interface and the surrounding soil mass (Lebeau, 2008; Kim and Jeong, 2011). However, an enormous computational effort is required in finite element analysis to incorporate the non-linear stress-strain response of soil and progressive slippage at the interface. Moreover, the studies conducted by Poulos (1989), Basile (2010), and Fattah et al. (2012) underlines several advantages of the use of boundary element method as an alternative of FEM to solve the pile-soil interaction problems.

The paper is presented in the following sequence: first, a model is formulated using BEM approach followed by its validation by comparing the BEM computed results with the available field data. Such comparison indicates reasonable accuracy of the proposed numerical solutions. Thereafter, a prototype parametric study is conducted using the developed model to analyse the influence of cyclic loading parameters on the soil-pile interactive performance. Finally, the normalized soil-pile interface shear stress profiles (predict using the proposed solution) are presented to show the distribution of stresses along the length of pile for static and post-cyclic condition.

# NUMERICAL MODELING: MATHEMATICAL FORMULATIONS

# Problem Definition

**Figure 1A** shows a single, vertical, floating pile of diameter D with embedded depth of L, subjected to an axial static load of V<sup>t</sup> and a two way symmetrical cyclic torsional load (given by Equation 1):

$$
\pi\_{\text{cyc}} = \pi\_{\text{cyc}}^{\text{max}} \sin 2\pi ft \tag{1}
$$

Where, τ max cyc and f are the amplitude and frequency of the cyclic torsional load, respectively. The imposed axial and torsional loads induce the stresses τ<sup>b</sup> and σ<sup>b</sup> at the base of the pile, and shear stresses τ t z and τ v z at the soil-pile interface along the length of the pile in the horizontal and vertical directions, respectively (see **Figure 1A**). The stresses τ v z and σbare primarily induced due to the axial load and are static, whereas, the stresses τ t z and τ<sup>b</sup> are induced due to the torsional load and are cyclic. The primary aim of this study is to evaluate these unknown interface shear stress components and subsequently, determine the axial pile capacity after completion of a certain number of load cycles (N).

The pile material in the present study has been idealized as elastic-perfectly plastic. The stress-strain behavior of the soil in shear is assumed to be non-linear up to the peak shear stress (τu), followed by a perfectly plastic post-peak response (Basack and Sen, 2014a). The non-linear pre-peak behavior has been represented using a hyperbolic equation (Equation 2) with an initial tangent modulus of G<sup>i</sup> and a reduction factor of R<sup>f</sup> (Duncan and Chang, 1970). The value of the reduction factor, R<sup>f</sup> usually varies in the range of 0.8–1.0 (Randolph, 2003).

$$
\pi = \frac{\mathcal{V}}{\frac{1}{G\_i} + \frac{R\_f}{\mathfrak{r}\_\mu}\mathcal{V}} \tag{2}
$$

Where, τ andγ are the shear stress and shear strain, respectively. The post peak response can be mathematically represented as:

$$
\mathfrak{r} = \mathfrak{r}\_u \tag{3}
$$

## Boundary Element Modeling

The pile-soil interaction in the present study has been analyzed using the Boundary Element Modeling (BEM) following the methodology of (Basack and Sen, 2014a,b). The pile is longitudinally discretized into n cylindrical elements of equal height (thickness), δ (see **Figure 1B**). The unknown pile-soil interface shear stress components in the horizontal and vertical directions in the i th element have been denoted as τ t i and τ v i ,respectively. Moreover, the displacement components at the central nodal plane (i.e., the central plane of each pile element) have been denoted as ρ<sup>i</sup> and θ<sup>i</sup> corresponding to the vertical and torsional modes, respectively.

Initially, the analysis has been performed for static loading with subsequent extension for the cyclic loading by using appropriate parameters for simulating the degradation of soil strength and stiffness (under cyclic loading), and the influence of the cyclic loading parameters. First, the static axial and torsional loadings are analyzed separately, followed by a coupled analysis to arrive at specific solutions. The governing differential equations for the static torsion (Equation 4) and static axial load (Equation 5) are given as (Basack and Sen, 2014a,b):

$$\frac{d^2\theta}{dz^2} = \frac{\pi D^2}{2J\_\mathcal{P} G\_\mathcal{P}} \mathbf{r}\_z^t \text{ (For static torsion)}\tag{4}$$

$$\frac{d^2\rho}{dz^2} = \frac{1}{E\_\rho} \left(\frac{4\pi\_z^\nu}{D} - \nu\_\ ^\nu\right) \quad \text{(For static axial load)}\tag{5}$$

Where, θ is the angle of twist; ρ is the vertical displacement; z is the depth; J<sup>p</sup> is the polar moment of inertia of pile cross section; G<sup>p</sup> is the modulus of rigidity of the pile; E<sup>p</sup> is the Young's modulus of the pile; γ ′ <sup>p</sup>is the unit weight of the pile. For static torsion, the governing differential equation is solved using the finite difference technique by developing correlations between the twist angle (θ) and some functions of D, δ, τ t i , Jp, G<sup>p</sup> and pile head torque (Tt). The correlations are then compiled together (in a matrix form) and the resultant matrix is given by:

$$[\![\![\!M]\!]\!] (\theta) = \{a\} \tag{6}$$

Where, [CM] is a coefficient matrix of order (n+1) x (n+1), {θ} is a column vector of order (n+1) x 1 and {a} is augment vector of order (n+1) x 1. The elements of the coefficient matrix, column vector, and augment vector can be found in Basack and Sen (2014a). There are two sets of unknown quantities in Equation 6, namely θ<sup>i</sup> andτ t i . Therefore, an initial no-slip condition is assumed and a correlation between θ<sup>i</sup> and τ t i (given by Randolph, 1981) is used (Equation 7) to reduce the number of unknown quantities in Equation (6).

$$\theta = \frac{\mathbf{r}\_i^t}{2G\_s} \tag{7}$$

Where, G<sup>s</sup> is the soil secant modulus. Using the provided correlation (Equation 7), the Equation (6) is modified to:

$$[DM]\{\pi\_i^t\} = \{b\} \tag{8}$$

Where, [DM] and {b} are coefficient matrix and augment vector, respectively and the elements of these matrices/vectors are functions of D, δ, Jp, G<sup>s</sup> , Gp, and T<sup>t</sup> . The initial values of the unknown horizontal torsional interface shear stress (τ t i ) are computed using Equation (8). The values of τ t i for each element are then compared with the limiting values of interface shear stress (τu) for each element. The elements are assumed to have slipped if the value of τ t i exceeds theτu. The initial value of τ t i for slipped elements is replaced by τ<sup>u</sup> and appropriate adjustments are made to the initial values of G<sup>s</sup> to incorporate the soil nonlinearity. The entire computation procedure is then repeated for the rest of the non-slipped elements until the desired convergence is achieved.

The governing differential equation for the static axial load is solved by using a similar procedure as described above i.e., by utilizing the finite difference technique and a correlation between ρ and τ v i given by Randolph and Wroth (1978) (for no-slip condition). After evaluating the unknown parametersτ t i andτ v i from the separate analyses for axial and torsion loads, a coupled analysis is conducted by evaluating the resultant interface shear stress for each element, through vector addition of τ t i andτ v i (see **Figure 2**), given by:

$$\mathbf{r}\_{i} = \begin{bmatrix} \mathbf{r}\_{i}^{t} + \mathbf{r}\_{i}^{\nu} \end{bmatrix}^{0.5} \tag{9}$$

The resultant interface shear stress is then compared with the ultimate shear stress and the values of τ t i are recomputed (or adjusted). The procedure is repeated until the desired convergence is achieved. Finally, the values of twist angles are computed for each element using Equation 6. The detailed formulations for the analysis of static load are published elsewhere (Basack and Sen, 2014a,b; Basack and Nimbalkar, 2017).

The analysis for the torsional cyclic loading has been performed by using a quasi-static method with a peak torsion of τ max cyc , wherein the elemental stress and displacement components have been adjusted after the completion of a desired number of load cycles. The application of cyclic loading influences the strength and stiffness of the soil to a large extent. On one hand, the cyclic loading leads to the degradation of shear strength and stiffness of soil due to generation of excess pore-water pressure, generation of irrecoverable plastic strain in soil around the pile and rearrangement of soil particles in the vicinity of the pile

(Poulos, 1981; Basack, 2015). However, on the contrary, the soil strength and stiffness increase with an increase in the loading rate (or frequency) (Poulos, 1989; Rodriguez and Alvarez, 2008). The following mathematical expression (Basack and Nimbalkar, 2017) addresses the combined effects of these two phenomena on the soil-pile interaction:

$$\mathbf{D}\_i^s = \left[ 1 + F \log\_{10} \left( \frac{2f \mathcal{D} \theta\_i}{\lambda\_r} \right) \right] \mathbf{N}^{-\frac{\mathcal{V} \cdot \mathbf{C}}{A + B \cdot \mathbf{y}}} \tag{10}$$

Where, D s i is the nodal soil degradation factor, F is a nondimensional rate factor, λ<sup>r</sup> is a datum loading rate, γ<sup>c</sup> is the peak nodal shear strain and A and B are the nondimensional cyclic soil parameters. The soil degradation factor (D s i ) is defined as the ratio of the post cyclic to pre-cyclic values of soil strength and stiffness. The derivation and the details of the parameters for Equation 9 can be found elsewhere (Basack and Nimbalkar, 2017).

The post cyclic axial pile capacity is then evaluated, based on the degraded values of the shaft and end bearing resistance, using the following expression:

$$\mathbf{Q}\_{\mathfrak{u}}^{\varepsilon} = \pi D \delta \sum\_{i=1}^{n} D\_{i}^{\varepsilon} \mathfrak{r}\_{\mathfrak{u}i} + \frac{\pi D^{2}}{4} \sigma\_{b\mathfrak{u}} - W\_{\mathfrak{P}} \tag{11}$$

Where, Q c u is the post-cyclic axial pile capacity, δ is the height of pile elements, τui is the elemental ultimate soil strength, σbu is the ultimate base restraint and W<sup>p</sup> is the self-weight of pile.

Finally, the pile degradation factor (Dp) is evaluated. It is defined as the ratio of the post-cyclic to static axial pile capacities (Equation 12):

$$\mathcal{D}\_{\mathfrak{p}} = \frac{Q\_{\mathfrak{u}}^{\mathfrak{c}}}{Q\_{\mathfrak{u}}^{\mathfrak{s}}} \tag{12}$$

where, Q s u is the ultimate static axial pile capacity. A user-friendly computer program PTCYC is written in FORTAN 90 language to conduct the desired computations. **Figure 3** shows the flowchart of the computer program.

# MODEL VALIDATION

The proposed solutions have been validated using the field test results available in literature. The field investigation on piles under composite torsional cyclic and axial static loading is extremely difficult and expensive. In absence of such research work, the authors have used the available laboratory and field test data to validate their numerical model. Stoll (1972) conducted a full-scale torsional load test on two concrete filled steel pipe piles embedded in a soil with a linearly increasing soil modulus. Guo and Randolph (1996) developed the analytical and numerical solutions for the torsional response of piles embedded in heterogeneous soil and validated the model with the field test results of Stoll (1972). The soil-pile interface shear stress profiles computed using the present model have been compared with the theoretical results of Guo and Randolph (1996) to validate the model (refer to **Figure 4**). From **Figure 4**, it can be observed that the results obtained using the present model are in good agreement with the solutions developed by Guo and Randolph (1996), with an average variation of about 21%.

Stuedlein et al. (2016) conducted a field study to evaluate the transfer of torsional load to soil along the soil-pile interface for two drilled concrete shafts in a site consisting of silty clay overlying a silty sand deposit. In addition to static torsion, cyclic torsion tests were also conducted on these shafts under one-way displacement-controlled mode with 20 displacement cycles at a frequency of 0.57 cycles per minute (cpm). The values of the resulting peak torsion evaluated from the present BEM computation have been compared with the field observations [for both the torsion test drilled shaft with production base (TDS) and torsion test drilled shaft with frictionless base (TDSFB)] of Stuedlein et al. (2016) (**Figure 5**). It is clear from this figure that the computed results are in acceptable agreement with the field test data with an average deviation of about 6%. Thus, the present numerical solutions are able to capture the soil-pile interaction under both static and cyclic torsional loads with an acceptable accuracy.

# PARAMETRIC STUDIES: ANALYSIS AND INTERPRETATION

The present boundary element model has been used to predict the response of a prototype vertical concrete floating pile embedded in soft clay subjected to a combined axial and cyclic torsional loading. **Table 1** shows the properties of the soil and pile used (adopted from Basack, 2010a). The soil unit cohesion and the initial tangent shear modulus at the ground surface are 30 kPa and 300 MPa, respectively, which are assumed to increase linearly with depth at a rate of 3 kPa/m and 30 GPa/m, respectively. Moreover, the key input parameters to account for the strength and stiffness degradation (or improvement) under cyclic loading are: A = 4.5, B = 2.5, F = 0.1, λr = 0.11 mm/s. The axial load on the pile is assumed to be 0.4 times the axial pile capacity (for pure axial load condition) i.e., the load ratio (Vt/Vu0) is 0.4. The number of pile elements are fixed at 100 after conducting a sensitivity check (Basack and Nimbalkar, 2017).

In the present study, the variation of the pile degradation factor (Dp) with cyclic loading parameters namely, number of cycles (N), frequency (f) and cyclic loading level (Lc) has been investigated. The cyclic loading level L<sup>c</sup> is defined as the ratio of peak cyclic torsion to the static ultimate torsional pile capacity. The number of cycles, frequency and cyclic loading level have been varied in the range of 10–1,000 cycles, 5–30 c.p.m., and 15– 30%, respectively. Moreover, the analysis has been conducted for TABLE 1 | Input parameters of soil and pile for the parametric study.


two values of reduction factor (0.85 and 0.95) to investigate its influence on the soil-pile interactive performance.

**Figures 6A–C** show the variation of the D<sup>p</sup> with N, f and L<sup>c</sup> , respectively. It can be observed that the parameter D<sup>p</sup> decreases with an increase in the number of loading cycles. However, the trend shows an asymptotically stabilizing tendency, i.e., the parameter D<sup>p</sup> becomes almost constant after a certain number of loading cycles. This may be attributed to the exponential degradation of the soil strength and stiffness with N (Idriss et al., 1978). Moreover, D<sup>p</sup> increases with an increase in the loading frequency. This is because the strength and stiffness of the soil increases logarithmically with f (Poulos, 1989). Furthermore, the pile degradation factor (Dp) decreases with an increase in L<sup>c</sup> following a curvilinear pattern with an increasing slope. This observation is reasonable because an increase in the value of the torsional cyclic amplitude or cyclic loading level is likely to cause the failure of more number of elements (which is initiated by the rapid yielding of soil adjacent to the soil-pile interface) (Basack, 2010b). It must be noted that the reduction factor R<sup>f</sup> shows insignificant effect on the soil-pile interactive performance.

**Figure 7** shows the profiles for the soil-pile interface shear stress (normalized by the product of unit cohesion and adhesion factor) for both static and post cyclic condition. It can be observed that the normalized shear stress attains a maximum value at the pile head and reduces to a minimum value at the base, in a curvilinear manner. It is interesting to note that the values of the post-cyclic shear stress have decreased as compared to the relevant value under static loading. Moreover, the percentage reduction is higher in the pile head (15%) as compared to the pile base (11%). The normalized shear stress values at the pile head and base are 0.74 and 0.45, respectively for the static condition while reduced to 0.63 and 0.4 after the application of the cyclic loading. This reduction is due to the degradation of the strength and stiffness of soil. Since, there is an overall reduction in the pile capacity, therefore, the influence of soil stiffness and strength degradation due to cyclic loading on the pile capacity is much higher as compared to the influence of loading frequency (which tends to enhance the soil strength and stiffness).

# PRACTICAL APPLICATIONS

The present study attempts to investigate the effect of the cyclic loading parameters namely, frequency, number of cycles and the amplitude, on the axial load carrying capacity of the pile foundation subjected to axial and torsional cyclic load. The results from the present study show that the axial loading capacity of the pile decreases with an increase in the number of loading cycles up to a particular number of cycles, beyond which the capacity becomes constant. Moreover, the pile capacity decreases with an increase in the amplitude of the cyclic torsional loading or the cyclic loading level. These findings can be used to predict the in-situ load carrying capacity of the pile installed below the offshore structures or transportation embankments in the real field (for the pile and soil conditions used in the present analysis) for the known amplitude and number of load cycles (which can be evaluated through physical observation). Similarly, the long-term performance of the pile foundations subjected to axial and cyclic torsional loadings can be predicted for the known values of the cyclic loading parameters. Nevertheless, the pile foundations must be designed with adequate factor of safety against the ultimate failure and acceptable displacement at the pile head (serviceability criterion) (Basack, 2010a). For the known cyclic loading parameters and ground conditions, the factor of safety against ultimate failure can be evaluated by computing the degraded pile capacity. Similarly, the variation of

## REFERENCES


twist angle with the loading parameters can be evaluated using the present model and a permissible limit could be established.

Another practical aspect of the present study is the proper assessment of the in-situ load carrying capacity of the pile foundation subjected to cyclic torsional loading. This assessment might help in the design of a suitable ground improvement measure such as electro-osmosis and high-voltage electro-kinetics, which could significantly improve the load carrying capacity of the pile foundations subjected to cyclic torsional loading.

# CONCLUSIONS

In the present study, a numerical solution based on the boundary element modeling has been developed for a single floating pile subjected to combined axial and torsional cyclic loads. The numerical model is successfully calibrated using appropriate values of the key input parameters and is validated against the field data published in the literature. The validation of the computed results with available field data exhibits the accuracy of the proposed solution. The results of numerical analysis indicate that the cyclic loading parameters, viz. number of load cycles, frequency, and cyclic loading level, significantly influence the degradation of the axial pile capacity due to torsional cyclic loading. Moreover, the interface shear stress has been found to decrease in a curvilinear pattern from a maximum value at the ground surface to a minimum value at the pile base. Furthermore, the proposed numerical solution can be used to evaluate the post-cyclic factor of safety relevant to the ultimate pile capacity. Thus, the results of the present parametric studies (conducted to investigate influence of key design parameters) can be utilized for formulating the design criteria for pile subjected to axial and torsional cyclic loads.

# AUTHOR CONTRIBUTIONS

All authors listed have made a substantial, direct and intellectual contribution to the work, and approved it for publication.

# ACKNOWLEDGMENTS

The Authors acknowledge the in-kind support from the School of Civil and Environmental Engineering, University of Technology, Sydney, Australia. The assistance received from Mr. Sankhasubhra Sen, former postgraduate student of Bengal Engineering and Science University, India in developing the computer program is acknowledged as well.


**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2019 Nimbalkar, Punetha, Basack and Mirzababaei. This is an openaccess article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# LIST OF NOTATIONS (BASIC SI UNITS ARE SHOWN IN PARENTHESES)

A & B Non-dimensional cyclic soil parameters (Dimensionless) D Diameter of pile (m) D s <sup>i</sup> Nodal soil degradation factor (Dimensionless) D<sup>p</sup> Pile degradation factor (Dimensionless) E<sup>p</sup> Young's modulus of the pile (N/m<sup>2</sup> ) F Non-dimensional rate factor (Dimensionless) f Loading frequency (Dimensionless) G<sup>i</sup> Initial tangent modulus (Dimensionless) G<sup>p</sup> Modulus of rigidity of the pile (N/m<sup>2</sup> ) G<sup>s</sup> Secant modulus (N/m<sup>2</sup> ) J<sup>p</sup> Polar moment of inertia of the pile (m<sup>4</sup> ) L Embedded depth of the pile (m) L<sup>c</sup> Cyclic loading level (N) N Number of load cycles (Dimensionless) n Number of pile elements (Dimensionless) Q c <sup>u</sup> Post-cyclic axial pile capacity (N) Q s <sup>u</sup> Ultimate static axial pile capacity (N) R<sup>f</sup> Reduction factor (Dimensionless) t Time (s) T<sup>t</sup> Pile head torque (N.m) V<sup>t</sup> Axial static load (N) Vu0 Pile capacity for pure axial load (N) W<sup>p</sup> Self-weight of pile (N) z Depth (m) δ Height of pile elements (m) θ<sup>i</sup> Twist on the ith element (V) ´ λ<sup>r</sup> Datum loading rate (Dimensionless) ρ<sup>i</sup> Vertical displacement at the ith element (m) τ Shear stress (N/m<sup>2</sup> ) τ<sup>b</sup> Shear stress at the base of the pile (N/m<sup>2</sup> ) τ max cyc Amplitude of the cyclic torsional load (N) τui Elemental ultimate soil strength (N/m<sup>2</sup> ) τ t z Interface shear stress component in horizontal direction at depth z (N/m<sup>2</sup> ) τ t i Interface shear stress component in horizontal direction for ith element (N/m<sup>2</sup> ) τ v z Interface shear stress component in vertical direction at depth z (N/m<sup>2</sup> ) τ v i Interface shear stress component in vertical direction for ith element (N/m<sup>2</sup> ) τ<sup>u</sup> Peak shear stress (N/m<sup>2</sup> ) σ<sup>b</sup> Normal stress at the base of the pile (N/m<sup>2</sup> ) σbu Ultimate base restraint (N/m<sup>2</sup> ) γ Shear strain (Dimensionless) γ<sup>c</sup> Peak nodal shear strain (Dimensionless) γ ′ <sup>p</sup> Unit weight of pile (N/m<sup>3</sup> ) γ ′ <sup>s</sup> Unit weight of soil (N/m<sup>3</sup> )

# Cyclic Fractional Plastic Model for Granular Soils

Ye Li 1,2, Yifei Sun<sup>2</sup> \* and Wen Ju<sup>3</sup>

*<sup>1</sup> China Guodian Corporation, Guodian Science and Technology Research Institute, Nanjing, China, <sup>2</sup> Key Laboratory of Ministry of Education for Geomechanics and Embankment Engineering, Hohai University, Nanjing, China, <sup>3</sup> China Nuclear Industry Huaxing Construction Co. Ltd., Nanjing, China*

Granular soils, e.g., sand, ballast, and rockfill, usually experience dynamic loads in the field. Traditional constitutive models for monotonic loading conditions cannot be used for advanced characterization of the complex loading behavior of granular soils. In this study, a simple fractional plastic model is developed, based on the generalized fractional plastic flow rule which considers the loading and unloading differences under triaxial compression and extension conditions. The model is further validated against a series of cyclic loading behavior of different granular soils, where a good predicting performance is observed.

Keywords: left-sided fractional derivative, right-sided fractional derivative, fractional plasticity, fractional flow rule, cyclic loads, granular soils

#### Edited by:

*Sakdirat Kaewunruen, University of Birmingham, United Kingdom*

#### Reviewed by:

*Trung Ngoc Ngo, University of Wollongong, Australia Ado Farsi, Imperial College London, United Kingdom*

> \*Correspondence: *Yifei Sun sunny@hhu.edu.cn*

#### Specialty section:

*This article was submitted to Transportation and Transit Systems, a section of the journal Frontiers in Built Environment*

> Received: *19 December 2018* Accepted: *08 March 2019* Published: *29 March 2019*

#### Citation:

*Li Y, Sun Y and Ju W (2019) Cyclic Fractional Plastic Model for Granular Soils. Front. Built Environ. 5:40. doi: 10.3389/fbuil.2019.00040* INTRODUCTION

According to many experimental (Aursudkij et al., 2009; Nimbalkar et al., 2012; Nimbalkar and Indraratna, 2016; Sun et al., 2017a, 2018e) and numerical (de Bono and McDowell, 2014; McDowell and Li, 2016; Li et al., 2018) studies, it is widely recognized that granular soils, including sand, ballast and rockfill, usually exhibit complex strength and deformation behavior, such as contraction accompanied by strain hardening and dilation accompanied by strain softening, when subjected to static and dynamic loads. Correct representation of such complex behavior of granular soils is the key factor for the design and safe operation of engineering facilities, for example, railroad and retaining wall (Nimbalkar and Choudhury, 2008; Nimbalkar et al., 2018). Hence, a number of different approaches, e.g., elastoplasticity, viscoplasticity, damage-plasticity, and bounding surface plasticity has been proposed. The models can be divided into three categories: (i) models that did not consider state dependence; (ii) models that cannot reflect non-associated flow; (iii) models that were not suitable for cyclic loading. For category (i), a number of constitutive methods can be found, for example, the disturbed state concept by Desai and Toth (1996) which was adopted for constitutive modeling of rockfills (Varadarajan et al., 2006) under different confining pressures; the damage-plasticity concept (Einav et al., 2007; Sun et al., 2015a), the diffuse failure models by Daouadji et al. (2011), the micromechanical models by Yin et al. (2010, 2017), the fractional cumulative models (Sun et al., 2015b, 2016a,b), the bounding surace models (Dafalias, 1986; Sun et al., 2014, 2017b; Sun and Shen, 2017). Some other models captured the state dependence but fall within category (ii), for example, the elastoplastic models by Sun et al. (2013, 2018c,f). In addition, most monotonic models cannot appropriately simulate the cyclic behavior of soil without modification on material flow and hardening, and thus fall within category (iii), for example, the state-dependent plastic models (Gajo and Muir Wood, 1999b; Li and Dafalias, 2000; Einav et al., 2007; Sun and Xiao, 2017; Sun et al., 2018d, 2019a).

To capture the state-dependent non-associated plastic behavior of sand, Been and Jefferies (1985) firstly suggested the concept of state dependence for sands, which was then promptly developed in various constitutive models (Li and Dafalias, 2000; Yang and Li, 2004; Sun et al., 2018a,b,c, 2019a). However, most of these models usually assumed a plastic flow rule that was only suitable for monotonic loads. For more complex loading conditions, for example, the cyclic loading case, more work needs to be carried out. In addition, classical constitutive models (Gajo and Muir Wood, 1999b; Imam et al., 2005; Collins et al., 2010) were usually dependent on the assumption of an additional plastic potential function, in order to correctly capture the non-associated stress-strain behavior of granular soils under either monotonic or cyclic loads. Is there a possible way to capture the state-dependent non-associated behavior of granular soils without using additional plastic potential? Sumelka (2014a,b) suggested the incorporation of fractionalorder derivatives into classical associated plasticity, from which a non-associated plastic flow rule can be achieved, without using additional plastic potential. However, the approach did not consider the physical practice in geomaterials, thus cannot reflect the underlying strength and deformation mechanisms of granular soil, unless proper modifications have been made. By connecting the state-dependent stress-dilatancy phenomenon with fractional plastic flow rule, a family of fractional plasticity models has been developed and applied in characterizing the non-associated stress-strain behavior of granular soil (Sun et al., 2017c, 2018b,c,d,f, 2019b,c). As demonstrated, the directions of plastic flow are no longer necessarily normal to the yielding or potential surfaces. However, these models did not consider the flow and hardening differences between loading, unloading and reloading, thus cannot properly capture the cyclic behavior of granular soils.

In fact, unlike the integer order derivative usually used in classical plasticity, the fractional order derivatives are nonlocal which induces a strong memory of the loading stress σ. Therefore, by using fractional derivative, the plastic flow of a material point was not only determined by the current stress state but also by its loading history. The size of the history and collection of yielding information are dependent on the employed definition of the derivative, for example, the Riesz-Caputo definition used by Sumelka and Nowak (2016) and the Riemann-Liouville definition used by Sun et al. (2016b). The choice of whatever fractional operator should rely on the specific physical issues is to be taken into consideration.

However, the present fractional plasticity models were mainly based on monotonic fractional flow rule, and cannot be used for capturing the cyclic behavior of granular soils, without possible modifications. Unlike monotonic models (Gajo and Muir Wood, 1999a; Li, 2002; Yang and Li, 2004), the plastic flow behavior and hardening/softening behavior loading, unloading and reloading should be different. Thus, the previous fractional plasticity model cannot properly capture the cyclic behavior of granular soils. To extend the fractional plasticity approach, a possible extension of the previous fractional plastic flow rule and hardening modulus for monotonic loading should be carried out. In contrast to previous works (Sun et al., 2018d, 2019c), this paper presents a simplified fractional-order elastoplastic model for granular soils subjected to cyclic loads, based on a general fractional plastic flow rule suitable for both monotonic and cyclic loads, which is the main innovation of this study. The study is structured as follows: section Definition of Fractional Derivative presents the basic definition of the fractional derivatives, while section Constitutive Model presents the basic constitutive relations; section Calibration of Model Parameters discusses how to identify each model parameters, while section Model Performance validates the proposed model; section Conclusions concludes the study by summarizing several main findings.

## DEFINITION OF FRACTIONAL DERIVATIVE

To begin, it is necessary to present a clear definition of the fractional derivative being used. According to the theory of fractional calculus, the fractional derivatives of the loading and unloading dynamic processes are different. In previous studies (Sun and Shen, 2017; Sun and Xiao, 2017), only the triaxial behavior of granular soils subjected to monotonic loading was considered. Therefore, only the left-sided Caputo fractional derivative with moving upper terminal (loading stress, σ ′ ) and fixed lower terminal (initial stress state, σ ′ <sup>0</sup>) was introduced in constitutive modeling. However, in this study, both the loading and unloading behaviors in triaxial tests will be addressed. The left-sided fractional derivative is suitable for describing the loading and unloading states of a material under triaxial compression where the deviator stress q > 0, while the rightsided fractional derivative is an operator performed on the loading and unloading states of a material subjected to triaxial extension where the deviator stress q < 0. Therefore, in this study, both the left-sided and right-sided Caputo's fractional derivatives (Agrawal, 2007) of a function f are used, such that:

$$\,\_{a}D\_{x}^{\alpha}f(\boldsymbol{x}) = \frac{1}{\Gamma(n-\alpha)} \int\_{a}^{\boldsymbol{x}} \frac{f^{(n)}(\boldsymbol{\tau})d\boldsymbol{\tau}}{(\boldsymbol{x}-\boldsymbol{\tau})^{\alpha+1-n}}, \qquad \boldsymbol{x} > a \tag{1}$$

$$\_xD\_b^\alpha f(\mathbf{x}) = \frac{\left(-1\right)^n}{\Gamma(n-\alpha)} \int\_x^b \frac{f^{(n)}(\mathbf{r})d\mathbf{r}}{\left(\pi-\mathbf{x}\right)^{\alpha+1-n}}, \qquad b > \infty \tag{2}$$

where D means derivation; Ŵ denotes the gamma function, defined as Ŵ(x) = R <sup>∞</sup> 0 e −τ τ <sup>x</sup>−1dτ . a and b are the lower and upper terminals used for integration. x is an independent variable and is designated as loading stress, σ, in this study. The fractional order, α , ranges from n−1 to n, where n = 1 or 2. Clearly, the fractional derivative is defined on an interval that is contrary to the integer order differential operators defined at a single point. Therefore, due to this non-locality, fractional approach in this study would intrinsically memorize the cyclic loading history.

## CONSTITUTIVE MODEL

In this study, only homogeneous and isotropic materials loaded under triaxial stress conditions are under consideration, where compressive stress and strain are regarded as positive. All the stresses used for derivation are effective stresses unless otherwise specified. Therefore, the following triaxial stress notations can be given:

$$\mathbf{s}^{\epsilon} = \begin{bmatrix} \varepsilon\_{\nu}^{\epsilon}, \varepsilon\_{s}^{\epsilon} \end{bmatrix}^{T} \tag{3}$$

$$\mathbf{e}^{\mathcal{P}} = \begin{bmatrix} \boldsymbol{\varepsilon}\_{\boldsymbol{\nu}}^{\mathcal{P}}, \boldsymbol{\varepsilon}\_{\boldsymbol{s}}^{\mathcal{P}} \end{bmatrix}^{T} \tag{4}$$

where ε indicates strain tensor for triaxial loading while the superscripts, e and p, indicate elastic and plastic components, respectively. The total volumetric (εv) and shear (εs) strains can be, respectively, defined as:

$$
\varepsilon\_{\nu} = \varepsilon\_1 + 2\varepsilon\_3 \tag{5}
$$

$$
\varepsilon\_s = 2(\varepsilon\_1 - \varepsilon\_3)/3 \tag{6}
$$

where ε<sup>1</sup> and ε<sup>3</sup> are the first and third principal strains, respectively. The total strain tensor ε = ε <sup>e</sup> + ε p , while the effective stress tensor, σ ′ , can be expressed as:

$$\boldsymbol{\sigma} = \begin{bmatrix} p', q \end{bmatrix}^T \tag{7}$$

where the mean effective principal and deviator stresses can be, respectively, defined as:

$$p' = \left(\sigma'\_1 + 2\sigma'\_3\right)/3\tag{8}$$

$$q = \sigma'\_1 - \sigma'\_3 \tag{9}$$

in which σ ′ <sup>1</sup> and σ ′ <sup>3</sup> are the first and third effective principal stresses, respectively. In this study, calculation of the elastic stress and strain response is based on the traditional theory of elasticity. However, the plastic strain is determined by using fractional plasticity (Sun and Xiao, 2017). A fundamental difference between classical and fractional plasticity is the calculation of the incremental plastic strain. In classical plasticity, the incremental plastic strain tensor, ε˙ p , is obtained by:

$$
\dot{\mathfrak{e}}^{\mathcal{P}} = \Lambda D^1 f(\sigma) \tag{10}
$$

where a superimposed dot indicates increment. 3 is the plastic multiplier and f often denotes the plastic potential or yielding function, depending on the flow rule that is chosen. Due to the equivalence between the first-order derivative and the classical differentiation of a point, the derivative interval in Equation (10) is omitted for clarity. As can be expected from Equation (10), the direction of plastic flow is fixed at each stress point once the function f is given. However, the plastic flow of granular soil varies with soil types. A plastic flow rule without considering the varied plastic flow was usually not able to unified modelling of the constitutive behavior of granular soils. Therefore, a flow rule for granular soils modified by using fractional-order differential operator was proposed (Sun and Shen, 2017; Sun and Xiao, 2017) where the incremental plastic strain tensor, ε˙ p , can be determined by:

$$
\dot{\mathfrak{s}}^{\rho} = \Lambda D^{\alpha} f(\sigma) \tag{11}
$$

where the modified Cambridge (Schofield and Wroth, 1968) relation (f) is used in this study for simplification:

$$f = \left(2p' - p'\_{\,0}\right)^2 + \left(\frac{2q}{M}\right)^2 - p'\_{\,0}^2 = 0\tag{12}$$

where p ′ 0 is the intercept of f with the abscissa. M = 6 sin ϕc/(3t−sin ϕc), is the critical-state stress ratio under triaxial compression. t = +1 for compressive loading whereas t = −1 for extensive loading. According to Sun and Xiao (2017), the plastic strain can be expressed as:

$$
\dot{\mathbf{e}}^{\mathcal{P}} = \frac{1}{\Pi} \mathbf{m}^T \mathbf{n} \dot{\boldsymbol{\sigma}} \tag{13}
$$

where the generalized flow direction (**n**) can be expressed as:

$$\mathbf{n} = \begin{bmatrix} n\_{\nu}, n\_{s} \end{bmatrix}^{T} \tag{14}$$

in which n<sup>v</sup> and n<sup>s</sup> are the flow directions induced by compression and shearing, respectively. Due to the distinct formulae of the fractional derivatives in describing loading and unloading, a generalized plastic flow rule suitable for both cases can be derived by substituting Equation (12) into Equations (1) and (2) (Sun et al., 2018a), such that:

$$n\_{\nu} = t \frac{8 \left[ p' - (2 - \alpha) p'\_0 / 2 \right] p'^{1 - \alpha}}{\Gamma(3 - \alpha) \left\| D^{\alpha} f(\sigma') \right\|} \tag{15}$$

$$m\_3 = \operatorname\*{tr} \frac{8M^{-2} \left| q \right|^{1-\alpha}}{\Gamma(3-\alpha) \left\| D^\alpha f(\sigma') \right\|} \tag{16}$$

where the gradient, D α f(σ ′ ) , is defined as:

$$\left\| D^{\alpha} f(\sigma') \right\| = \frac{8}{\Gamma(3-\alpha)} \sqrt{\left( \frac{q^{1-\alpha}}{M^2} \right)^2 + \left[ p' - \left( 1 - \frac{\alpha}{2} \right) p'\_0 \right]^2 p'^{2-2\alpha}} \tag{17}$$

The stress-dilatancy relationship for loading and unloading can be therefore obtained by using Equations (16) and (17), where the classical modified Cam-clay stress-dilatancy relationship can be also achieved by using α = 1. Detailed derivations of Equations (15) and (16) can be found in Sun et al. (2018a) and thus not repeated here for simplicity. In addition, the loading tensor **m** is assumed to be the same as the flowing tensor, such that:

$$m\_{\nu} = t \frac{d}{\sqrt{1 + d^2}}\tag{18}$$

$$m\_s = t \frac{1}{\sqrt{1 + d^2}} \tag{19}$$

where the dilatancy ratio d can be defined as:

$$d = \frac{M^2 - (1 - \alpha/2)[\eta^2 + M^2]}{t|\eta|^{2 - \alpha}}\tag{20}$$

The hardening modulus 5 should satisfy the following conditions according to Li (2002): (i) 5 = +∞ when the stress ratio η = q/p ′ = 0, (ii) 5 = 0 when specimens are at the critical and drained peak stress states. Therefore, the loading plastic modulus 5 can be expressed as:

$$
\Pi = h\_L G \frac{t\eta\_\rho - \eta}{\eta} e^{m\psi} \tag{21}
$$

where m is a model constant; the shear modulus (G) is defined as:

$$G = G\_0 \frac{(2.97 - e)^2}{1 + e} p\_a \sqrt{\frac{p'}{p\_a}} \tag{22}$$

in which G<sup>0</sup> is the elastic modulus of the material and

$$h\_L = (h\_1 - h\_2 \psi) \tag{23}$$

where h<sup>1</sup> and h<sup>2</sup> are hardening parameters. The state parameter ψ is defined as

$$
\psi = e - e\_{\mathfrak{c}} \tag{24}
$$

where e is the current void ratio; e<sup>c</sup> is the critical-state void ratio that can be expressed as (Li and Wang, 1998):

$$e\_{\mathfrak{c}} = e\mathfrak{r} - \lambda \left(\mathfrak{p}'/p\_a\right)^{\sharp} \tag{25}$$

where λ , eŴ, and ξ are critical state parameters, describing the critical state line in the e − p ′ξ plane. It should be noted that the current model does not explicitly consider the particle breakage behavior of granular soil. But, as evidenced by laboratory tests (Bandini and Coop, 2011; Ghafghazi et al., 2014; Yu, 2017a,b), the influence of particle breakage on the mechanical behavior of granular soil was reflected by shifting the critical state line in the e − ln p ′ plane. To capture this behavior, numerous critical state lines incorporating particle breakage in the e − ln p ′ plane were proposed. However, as suggested in Li and Wang (1998), a linear representation of the critical state line in the e − p ′ξ plane can be simply used to implicitly consider the particle breakage, which had been used by a number of researchers, including Gajo and Muir Wood (1999a,b) and Dafalias and Taiebat (2016), etc.

Moreover, a slight difference of the hardening parameter is found between the current and previous study (Sun and Xiao, 2017) where h<sup>L</sup> was correlated to the initial ψ<sup>0</sup> and e<sup>0</sup> rather than the evolution of the current ψ and e. However, one may not be able to know the initial ψ<sup>0</sup> and e<sup>0</sup> of soils that have been already sheared. In such cases, the modification shown in Equation (23) works. η<sup>p</sup> is the virtual peak stress ratio, which can be correlated to the state parameter as:

$$
\eta\_{\mathcal{P}} = M e^{-m\psi} \tag{26}
$$

For reloading, the plastic modulus can be further expressed as:

$$
\Pi = h\_L h\_C G \frac{t \eta\_p - \eta}{\eta} e^{m\psi} \left( 1 + \frac{tM}{\eta} \right)^n \tag{27}
$$

where h<sup>C</sup> = e −rε<sup>v</sup> , is the densification factor, considering accumulation of strain caused by loading and reloading; and n and r are material constants, which are only used for capturing cyclic loading. For undrained loading condition, ε<sup>v</sup> = 0, the plastic modulus is not affected by h<sup>C</sup> , as suggested by Ling and Yang (2006). Hence, the dependence of material hardening on the material density occurs via its dependence on ψ .

To characterize the unloading behavior of granular soils, the unloading plastic modulus by Ling and Yang (2006) is used:

$$
\Pi = h\_U h\_C G \left(\frac{M}{\eta\_U}\right)^\vartheta \tag{28}
$$

where h<sup>U</sup> and ϑ are model parameters. η<sup>U</sup> is the stress ratio at which the unloading occurs. For tests with |M/ηU| ≤ 1, ϑ = 0

To consider the elastic deformation during loading and unloading, the hyperelasticity is used to determine the elastic incremental strain, such that:

$$
\dot{\mathbf{e}}^{\varepsilon} = \mathbf{C}^{\varepsilon} : \dot{\mathbf{o}} \tag{29}
$$

where **C** e is the elastic compliance tensor:

$$\mathbf{C}^{\epsilon} = \begin{bmatrix} 1/K \\ & 1/(3G) \end{bmatrix} \tag{30}$$

where the bulk modulus (K) can be defined as:

$$K = \frac{\langle 2 + 2\nu \rangle}{\Im(1 - 2\nu)} G \tag{31}$$

where ν is the Poisson's ratio, describing the lateral deformation capability of the material. Thus, a complete description of the elastoplastic stress-strain behavior of granular soils under cyclic loads can be achieved, by using Equations (13) and (29).

## CALIBRATION OF MODEL PARAMETERS

There are in total 14 model parameters, i.e., four critical state parameters (ϕ<sup>c</sup> , λ , eŴ, ξ ), one fractional order (α ), seven hardening parameters (h1, h2, m, hU, r, n, ϑ ) and two elastic constants (G0, ν ). Details of how to determine the model parameters are described as follows.

The critical state parameters (ϕ<sup>c</sup> , λ , eŴ, ξ ) define the critical state of the material, which can be determined by fitting the critical state points in the p ′ − q and e − p ′ planes. For most granular soils, the critical state friction parameters are independent of the loading state.

The fractional order, α , determines the plastic flow directions of the material. Therefore, it can be determined by using the leastsquares method to fit the stress-dilation relationship as shown in **Figures 1**, **2**. To be compatible with the critical state soil mechanics, α is equal to unit when the critical state void ratio is reached.

The hardening parameters, h<sup>1</sup> and h2, determine the hardening and softening behavior of the material, which can be determined by fitting the ε<sup>s</sup> − η relationship of specimens

under different initial monotonic test conditions, as discussed in Sun and Shen (2017). The peak failure constant, m, can be calibrated from the stress points at peak failure state by using:

$$m = \frac{1}{\psi\_{\mathcal{P}}} \ln \frac{M}{\eta\_{\mathcal{P}}} \tag{32}$$

where ψ<sup>p</sup> and η<sup>p</sup> are two values of ψ and η at the peak stress state. The peak stress decreases as k increases. There are four hardening parameters (hU, θ , r, n) for describing the unloading/reloading behavior of the material. h<sup>U</sup> can be determined by fitting the slope of the first unloading stress-strain curve while ϑ can be determined from the rate of change of the slope of the first unloading curve. r is determined by fitting the hysteretic loops in the stress-strain curve. n can be obtained by fitting the first reloading stress-strain curve of the material. Detailed discussions on determining the hardening parameters for cyclic loads can be found in Ling and Yang (2006) and thus not repeated here.

The elastic constant, G0, mainly determines the elastic characteristics of the material, which can be obtained by rearranging Equation (19):

$$G\_0 = \frac{(1+e)G}{\left(2.97-e\right)^2 \sqrt{p'p\_a}}\tag{33}$$

The Poisson ratio, ν , usually ranges between 0.05 and 0.35 for most granular soils. It defines the lateral deformation ability that can be determined by:

$$\upsilon \approx \frac{9\varepsilon\_s - 2\varepsilon\_\nu}{18\varepsilon\_s + 2\varepsilon\_\nu} \tag{34}$$

Detailed values of the model parameters of each material simulated in this study can be found in **Table 1**.

## MODEL PERFORMANCE

In this section, the proposed fractional order elastoplastic model is validated by simulating the drained and undrained triaxial behaviors of different granular soils, including sand and rockfill. Specifically, **Figures 1**–**3** present the model simulations of the monotonic and cyclic behavior of Fuji River sand (Ishihara et al., 1975). The model simulations of the drained triaxial behavior of rockfill (Li, 1988; Fu et al., 2014) are shown in **Figures 4**–**6**.

**Figures 1**–**3** present the model predictions of the drained and undrained triaxial behavior of Fuji River sand (Ishihara et al., 1975). The material primarily consisted of sub-angular aggregates with d<sup>50</sup> of 0.22 mm and C<sup>u</sup> of 2.21. Samples of 50 mm in diameter and 100 mm in height were prepared by pluviating fresh sand into the molds which were filled with deaired water. The e<sup>0</sup> used for simulating undrained tests are 0.740, 0.731, and 0.718, with the corresponding σ ′ <sup>3</sup> equal to 98, 196, and 294 kPa, respectively. The e<sup>0</sup> for simulating drained monotonic tests are 0.750, 0.747, and 0.751, with the corresponding σ ′ 3 equal to 98, 196, and 294 kPa, respectively. It can be observed in **Figures 1**, **2** that the undrained and drained monotonic stress-strain relationship of Fuji River sand with various initial conditions can be well-captured by using α = 0.95. The initial

TABLE 1 | Model parameters related to monotonic loading.


deviator stress vs. predicted axial strain.

contraction and the subsequent dilation of the samples (**Figure 2**) are well-simulated, highlighting the rationality of the adopted fractional flow rule. Moreover, the simulated deviator stress increases rapidly until reaching a critical state, which agrees well with the experimental results shown in **Figure 2**. The undrained cyclic performance of Fuji River sand with e<sup>0</sup> = 0.737 and σ ′ 3 = 206.5 kPa is simulated in **Figure 3**, by using the additional model parameters: h<sup>L</sup> = 0.5, h<sup>U</sup> = 0.1, r = 130, n = 1, ϑ = 4. It is found that the model simulation of the stress path is in reasonable agreement with the corresponding test results prior to liquefaction. However, the simulation result of the variation of deviator stress vs. axial strain is less favorable as strain increases.

Monotonic and cyclic test results of Xiaolangdi rockfill reported by Fu et al. (2014) and Li (1988) are simulated in **Figures 4**–**6**. The material primarily consisted of slightly weathered sandstones. The e<sup>0</sup> equal to 0.199 was used for all the tests. **Figure 4** shows the model prediction of the monotonic stress-strain behavior of Xiaolangdi rockfill (Fu et al., 2014). Concordance between the model simulations and corresponding test results can be observed by using α = 0.93. More specifically, the material hardening and softening accompanied by volumetric contraction at high confining pressure and dilation at relatively low confining pressure are well-characterized. **Figures 5**, **6** show the model simulations of the cyclic behavior of rockfill. The initial σ ′ <sup>3</sup> = 0.5 and 1 MPa while the additional model parameters used for model simulation are: h<sup>L</sup> = 0.2, h<sup>U</sup> = 8, r = 10, n = 2, ϑ = 3.3. It can be observed that the loading and unloading

FIGURE 4 | Model prediction of the drained monotonic behavior of Xiaolangdi rockfill (Fu et al., 2014): (A) deviator stress vs. axial strain, (B) volumetric strain vs. axial strain.

stress-strain responses can be reasonably simulated by the proposed fractional-order model.

# CONCLUSIONS

A fractional-order plastic flow rule was suggested in previous studies. However, due to the limitations induced by the mathematical definitions of the (left-sided) fractional derivative,

the suggested fractional-order flow rule can only be applied to model the monotonic stress-strain response of granular soils. To solve this problem, a new cyclic fractional plasticity model was developed in this study. The main conclusions can be drawn as:


# AUTHOR CONTRIBUTIONS

All authors listed have made a substantial, direct and intellectual contribution to the work, and approved it for publication.

## REFERENCES


# ACKNOWLEDGMENTS

We would like to express our sincere gratitude to Prof. Wen Chen in Hohai University, for his lifelong inspiration on fractional mechanics.


**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2019 Li, Sun and Ju. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Finite Element Modeling of Ballasted Rail Track Capturing Effects of Geosynthetic Inclusions

#### Yajun Jiang<sup>1</sup> and Sanjay Nimbalkar <sup>2</sup> \*

*<sup>1</sup> Key Laboratory of Transportation Tunnel Engineering, Ministry of Education, Southwest Jiaotong University, Chengdu, China, <sup>2</sup> School of Civil and Environmental Engineering, University of Technology Sydney, Sydney, NSW, Australia*

This paper presents a two dimensional finite element (FE) approach to investigating beneficial aspects of geogrids in the railway track. The influences of different factors including the subgrade strength, the geogrid stiffness, the placement depth of geogrid, the effective width of geogrid, the strength of ballast-geogrid interface and the combination of double geogrid layers were investigated under the monotonic loading. The results indicated the role of geogrid reinforcement is more pronounced over the weak compressible subgrade. A stiffer geogrid reduces ballast settlement and produces a more uniform stress distribution along a track. The placement location of a geogrid is suggested at the ballast-sub-ballast interface to achieve better reinforcement results. Although the width of a geogrid layer should be sufficient to cover an entire loaded area, excessive width does not guarantee additional benefits. Higher interface strength between a ballast and a geogrid is beneficial for effective reinforcement. Increasing the number of geogrid layers is an effective way to reinforce the ballast over weak subgrades. The results of the limited cyclic FE simulations revealed the consistency of the reinforcement effect of the geogrids under monotonic and cyclic loads.

#### Edited by:

*Sakdirat Kaewunruen, University of Birmingham, United Kingdom*

#### Reviewed by:

*Javad Sadeghi, Iran University of Science and Technology, Iran Trung Ngoc Ngo, University of Wollongong, Australia*

> \*Correspondence: *Sanjay Nimbalkar Sanjay.Nimbalkar@uts.edu.au*

#### Specialty section:

*This article was submitted to Transportation and Transit Systems, a section of the journal Frontiers in Built Environment*

> Received: *28 February 2019* Accepted: *07 May 2019* Published: *22 May 2019*

#### Citation:

*Jiang Y and Nimbalkar S (2019) Finite Element Modeling of Ballasted Rail Track Capturing Effects of Geosynthetic Inclusions. Front. Built Environ. 5:69. doi: 10.3389/fbuil.2019.00069* Keywords: finite element modeling, geogrid, railway ballast, reinforcement mechanism, settlement

# INTRODUCTION

Ballasted rail tracks progressively deform in vertical and lateral directions, which cause deviations from the desired geometry under repeated traffic loading (Selig and Waters, 1994). A ballast can deteriorate due to the breakage of angular sharp corners under cyclic loading, which accelerates the settlement of a track (Indraratna et al., 2005). The differential settlement of a track reduces the level of safety of the track and causes significant risk to trains. Therefore, track maintenance, including tamping and replacement of ballasts, is periodically conducted to correct geometry faults (Esveld, 2001).

The ballast layer is often contaminated by fouling agents such as downward coal spillage from moving wagons, upward clay migration and internal particle breakage (Nimbalkar et al., 2012; TolouKian et al., 2018). This leads to rapid deterioration along with settlements often combined by the lateral spread owing to less lateral restraints (Nimbalkar and Indraratna, 2016; Sun et al., 2017). The regular track maintenance operations are needed which attracts huge funds to maintain track in operations (Sadeghi et al., 2018, 2019). In recent years, studies have focused on the improving confining pressure on a ballast to reduce the settlement of railway tracks during operation (Indraratna et al., 2005, 2013; Lackenby et al., 2007).

Various types of geosynthetics have been extensively employed in practice to improve or modify the behaviors of soil, aggregate and other construction materials in geotechnical engineering (Ghosh and Madhav, 1994; Moghaddas-Nejad and Small, 1996; Haeri et al., 2000; Raymond and Ismail, 2003; Palmeira, 2009; Indraratna et al., 2015). In current practical applications, the insertion of a geosynthetics layer in track substructures is an effective and economic method for increasing the confining pressures on a ballast and conducting the renewal of a track (Amsler, 1986; Raymond, 1999; Lieberenz and Weisemann, 2002; Indraratna and Salim, 2003; Fernandes et al., 2008; Indraratna et al., 2010, 2014a).

Geogrids have been utilized to reinforce railway track substructures since 1980s (Coleman, 1990). Previous laboratory studies (Bathurst and Raymond, 1987; Gobel and Weisemann, 1994; Raymond, 2002; Shin et al., 2002; Brown et al., 2007; Indraratna et al., 2011a) and field studies (Indraratna et al., 2010, 2014a,b) have demonstrated the successful application of geogrids in track reinforcement. Limited numerical simulations have been conducted to analyze the reinforcement mechanism of geogrids in railway tracks (Indraratna et al., 2007; Jirousek et al., 2010; Chen et al., 2012; Indraratna and Nimbalkar, 2013). The discrete element models are usually applied to capture the particle-scale mechanism and interface behavior using geosynthetics (Han et al., 2011; Tutumluer et al., 2012; Ngo et al., 2017; Gao and Meguid, 2018). However, these models are not able to capture the overall track response (Li et al., 2018). A generally accepted method or standard for the design or renewal of railway tracks with geogrids is lacking.

In contrast to laboratory and field investigations, the numerical method is a convenient and economical method for investigating the mechanical behaviors of geogrids and the surrounding components of track substructures. Thus, the main objective of this study was to employ a finite element (FE) model to improve the comprehensive understanding of the geogrid reinforcement mechanisms of railway ballasts. Influencing factors such as the subgrade strength, the geogrid stiffness, the placement depth of a geogrid layer inside a ballast, the effective width of a geogrid layer, the strength reduction factor of the geogrid/ballast interface and the effect of double geogrid layers were considered. Subsequently, the reinforcement performance of a geogrid and the influences of these factors were illustrated based on the simulation results. The simulation results of cyclic loading was compared with field data to explore ballast behavior under repeated train loads.

# ESTABLISHMENT OF FE MODEL

The FE method has been proven to be a useful and effective method for analyzing the behaviors of track structures (Selig and Waters, 1994; Esveld, 2001). The authors have performed a series of fundamental numerical simulations, which focused on the behaviors of ballasts with and without geosynthetics, using plane-strain track models (Indraratna et al., 2007; Indraratna and Nimbalkar, 2013). Thus, this study can be considered as an extension of the authors' previous studies.

## Modeling

PLAXIS has demonstrated its success in the limit analysis of geotechnical problems (PLAXIS, 2007). In this study, a twodimensional plane-strain finite element model of a composite multi-layer track system, including the rail, sleeper, ballast, subballast and subgrade (as shown in **Figure 1**), was numerically simulated using PLAXIS2D.

Only half of the rail track was simulated in the model due to the symmetry of the track. The width and height of the subgrade is 6 and 3 m, respectively. The heights of the subballast, ballast and sleeper are 150, 300, and 200 mm, respectively. The gauge length of the track is 1.68 m, the length of the half sleeper is 1.25 m, and the side slope of the rail track embankment is 1:2. The rail is simplified as a rectangle with a width of 160 mm and a height of 50 mm to ensure that the total area is equivalent to the Australian 60 kg/m rail (Doyle, 1980).

The nodes along the bottom boundary of the section were considered to be standard and fixed. The left and right boundaries were restrained in horizontal directions to represent smooth contact in the vertical direction. The rail was restrained in the horizontal direction to simulate the connection between the sleeper and the rail.

The geogrid in the model was created with the Geogrid Element provided by PLAXIS2D, and the interface elements along each side of the geogrid were created to simulate the interaction relationship between ballast and geogrid (PLAXIS, 2007).

The geogrid was initially placed at the ballast/sub-ballast interface in the model; however, different placement locations were also considered at different phases during the simulations.

## Parameters

#### Material Parameters of Main Components

The material parameters and constitutive models that were employed for the main components of the model are listed in **Table 1**, and the ballast was simulated using a hardening soil model, which is suitable to describe the behavior of a ballast under wheel loading (Indraratna et al., 2012; Indraratna and Nimbalkar, 2013). The parameters are chosen based on the actual measurement data of ballast behavior as observed in the laboratory tests (Indraratna and Nimbalkar, 2013) and field trials (Nimbalkar and Indraratna, 2016). The numerical model is an extension of previous works (Indraratna and Nimbalkar, 2011; Nimbalkar and Indraratna, 2014) and results are in agreement with these numerical studies. The accuracy and reliability of the model is thus demonstrated.

#### Subgrade Strength

The strength parameters of different typical subgrades, including poor, fair, good and rocky subgrades, are listed in **Table 2**; they were considered in the simulations to explore the influence of subgrade strength.

#### Axial Geogrid Stiffness

The axial geogrid stiffness always differs by the factors of the aperture dimensions, the tensile strength of the manufacture

TABLE 1 | Material parameters and constitutive models in the finite element analysis (Indraratna et al., 2012).


*E* = *elasticity modulus, Eref <sup>50</sup>* = *secant stiffness at 50% strength for loading conditions, E ref ur* = *triaxial unloading/reloading stiffness, Eref oed* = *tangent stiffness for primary oedometer loading,* γ = *unit weight,* ν = *Poisson's ratio for loading conditions,* ν*ur* = *Poisson's ratio for unloading/reloading conditions, c* = *effective cohesion,* φ = *effective friction angle,* ψ = *dilatancy angle, pref* = *reference confining pressure, m* = *stress dependent stiffness factor, knc* <sup>0</sup> = *coefficient of earth pressure at rest for normal consolidation, R<sup>f</sup>* = *failure ratio and Rinter* = *strength reduction factor in the interface.*

material and the rib thickness. Thus, the EA values of the geogrids in the simulations were varied for comparison purposes (i.e., EA = 275, 525, 775, 1,025, and 1,275 kN/m).

TABLE 2 | Material parameters and constitutive models for subgrades (adapted from Khordehbinan, 2009).


### Placement Depth of Geogrid Layer

The placement depth of the geogrid layer inside the ballast was initially set at the ballast/sub-ballast interface in the model, i.e., 300 mm below the sleeper. Other placement depths inside the ballast were also considered in the simulations (i.e., 50, 100, 150, 200, and 250 mm below the sleeper).

#### Effective Width of Geogrid Layer

The effective width of the geogrid layer was initially set to 2.2 m when the geogrid layer was placed at the ballast/sub-ballast interface in the model. To investigate the influence of width, effective widths of the geogrid layer of 0.4, 0.8, 1.2, 1.6, and 2.0 m were also considered.


*a "P, F, G, R" stands for Poor, Fair, Good and Rocky, respectively.*

*b including the values of 275, 525, 775, 1,025, and 1,275.*

*c including the values of 50, 100, 150, 200, 250, and 300.*

*d small differences of the width of geogrid layer were caused by changing the placement of geogrid layer inside the ballast due to the side slope of the rail track embankment.*

*e including the values of 0.4, 0.8, 1.2, 1.6, 2.0, and 2.2.*

*f including the values of 0.4, 0.6, 0.8, and 1.0.*

*<sup>g</sup>double geogrid layers with placement depths of 300 and 200 mm.*

#### Strength Reduction Factor of Geogrid Interfaces

The strength reduction factor Rinter defines the behaviors of the interfaces between the geogrid and the adjacent components. This parameter is influenced by the complex interaction relationship between the aperture size of the geogrids, the mean particle size of the soil materials, the physical properties of the geogrids and the mechanical behaviors of the soils (Goodhue et al., 2001; Brown et al., 2007; Indraratna et al., 2012). Therefore, the values of Rinter, including 0.4, 0.6, 0.8, and 1.0 (rigid), were considered to investigate the influence of the interaction between the geogrid and the ballast on the reinforcement effect.

## Load Determination

The equivalent dynamic wheel load P<sup>d</sup> is empirically expressed as the static wheel load P<sup>s</sup> times the dynamic impact factor ø (Jeffs and Tew, 1991; Sun et al., 2015). The American Railroad Engineering Association (AREA) formula was applied in the calculation of the impact factor ø (Jeffs and Tew, 1991)

$$
\phi = 1 + 5.21\nu/D \tag{1}
$$

where v = vehicle speed (km/h) and D = wheel diameter (mm). A P<sup>d</sup> value of 178.71 kN was considered to simulate a train speed of 80 km/h and a wheel diameter of 970 mm.

## Calculation Cases

The following factors were investigated in this study: the subgrade strength, the geogrid stiffness (EA), the placement depth of the geogrid layer inside the ballast (D), the effective width of the geogrid layer (B), the strength reduction factor of the geogrid/ballast interfaces (Rinter) and the effect of double geogrid layers. To avoid the complex interaction effects among the factors, the simulation investigation was divided into six calculation cases. Details of the parameter settings of these factors for each case are presented in **Table 3**.

TABLE 4 | Results of Case 1 with respect to the factor of subgrade strength.


*a this negative result of R over Rocky subgrade might be caused by the compression of the geogrid layer in the model.*

# ANALYSIS AND DISCUSSION OF MONOTONIC LOADING RESULTS

The reinforcement effect of a geogrid can be evaluated by the reinforcement rate R, which was calculated by

$$R = \frac{\text{S}\_{tu} - \text{S}\_{tr}}{\text{S}\_{tu}} \times 100\% = \frac{\text{S}\_{t}}{\text{S}\_{tu}} \times 100\% \tag{2}$$

where Stu and Str = the settlement on the top of ballast (below the rail) without and with geogrid reinforcement, and 1S<sup>t</sup> = the settlement difference on the top of the ballast with the geogrid reinforcement.

## Subgrade Strength

The results for Stu, Str, 1S<sup>t</sup> , and R of Case 1 are presented in **Table 4**. Improved geogrid reinforcement performance was observed over the weaker subgrades, which can be explained by the geogrid reinforcement mechanisms. The subgrade settlement determined the deformations of other track components and subsequently mobilized the tensile/interlocking capacity of the geogrid (Kwon and Penman, 2009).

To illustrate the influence of subgrade strength, the vertical deformations of the subgrade surface over the different subgrades in the reinforced models in Case 1 are plotted in **Figure 2**. The maximum axial forces on the horizontal direction in the geogrids (Fag )max over different subgrades in the reinforced models are plotted in **Figure 3**. As expected, the larger surface deformation of the subgrade appeared over the weaker subgrade. Consequently, the different axial forces in geogrids were induced by the corresponding subgrade deformations. The different extents of the interlock, confinement and load distribution effects were subsequently observed. Therefore, the subgrade strength is an important factor of geogrid reinforcement performance. Better geogrid reinforcement performance can be expected over the weak subgrades due to the greater compressibility of the weak subgrades.

## Geogrid Stiffness

The results for R and (Fag )max of Case 2 are plotted in **Figure 4**. An approximately linear relationship between R and EA was observed, i.e., increasing the geogrid stiffness increased the enhanced reinforcement rate. (Fag )max almost linearly increased by increasing EA. Subsequently, additional reductions of the deformations and more uniform stress distributions in the surrounding soil layers can be expected with a stiffer geogrid.

This may be attributed to the fact that the geogrid with a higher stiffness value produced greater axial force and less tensile deformations under the same vertical load. Consequently, the reinforcement effects were enhanced due to the higher geogrid stiffness. The influence of subgrade strength on the geogrid reinforcement was observed, i.e., better reinforcement effects were achieved over the weaker subgrades. In addition to the geogrid stiffness, the performance of the geogrid reinforcement is also dependent on the subgrade strength.

# Placement Depth of Geogrid Layer

The results for R and (Fag )max of Case 3 are plotted in **Figure 5**. The reinforcement rates over the fair and poor subgrades increased with an increase in the placement depth of the geogrid layer inside the ballast. The peak values of the reinforcement rates over poor and fair subgrades occurred within the placement depth range of 250–300 mm, which is located near the bottom of the ballast layer. (Fag )max increased with an increase in D, and the maximum value of (Fag )max was obtained over the poor subgrade due to the maximum subgrade deformation. Therefore,

this study recommends that geogrids should be placed near the bottom of a ballast, which is consistent with practical solutions in the field, in which the geogrid layer is usually placed at the ballast/sub-ballast interface to protect geogrids from damage due to tamping maintenance. Some research results also suggested that the optimum placement depth for the geosynthetics inside a ballast should be varied for different conditions (Coleman, 1990). Brown et al. (2007) suggested the bottom of a ballast layer as an optimum placement location for geogrids.

# Effective Width of Geogrid Layer

The results for R of Case 4 are plotted in **Figure 6A**. The reinforcing function of the geogrid over fair and poor subgrades was not activated until the effective width of the geogrid layer exceeded 1.2 m, which was almost equal to the length of the half model sleeper. Subsequently, higher reinforcement rates were observed in the width range of 1.6–2.0 m for the geogrid layer. However, the peak value of the reinforcement rate was not observed at the maximum width of 2.2 m. Therefore, a sufficient width of the geogrid layer is essential for the reinforcement; however, excessive width does not provide significant benefit.

To investigate the loaded area of the location where the geogrid was placed, the normal strains (εn) on the top of

the sub-ballast over different subgrades of unreinforced models are plotted in **Figure 6B**. The loaded areas on the top of the sub-ballast over different subgrades were almost located in the depth range of 1.0–1.5 m under the unreinforced conditions. The excessive width of the geogrid layer cannot provide additional benefits when the width extends beyond the range of the loaded area as the reinforcement benefit by interlock was generated within the loaded area and the geogrid did not have to be anchored beyond the loaded area (Guido et al., 1987).

# Strength Reduction Factor

The results for R of Case 5 are plotted in **Figure 7A**. The distinct trend is an increase in the reinforcement performance with an increase in Rinter. A constant value of the reinforcement rate was attained when Rinter varied from 0.8 to 1.0. As expected, a larger Rinter produced better performance. Minor settlement reduction was observed over the rocky subgrade and the good subgrade. The geogrid reinforcement cannot be employed in the railway ballast over the stiff subgrades, which is consistent with field measurements that indicate an increased effectiveness of the reinforcing geogrids on the soft subgrade (Indraratna et al., 2014a).

The value of Rinter influencesthe interaction between a geogrid and a ballast. Generally, tighter interaction between a geogrid and a ballast may produce a higher axial force in the geogrid. Thus, the (Fag )max values of Case 5 are plotted in **Figure 7B** to show the influence of Rinter. The values of (Fag )max over the rocky and good subgrades did not significantly change, i.e., the reinforcement effect of the geogrid over these stiff subgrades can be disregarded. Thus, the variation in Rinter did not provide any pronounced improvement of the settlement of the ballast. Distinct increments of (Fag )max were not observed when Rinter varied in the range of 0.6–1.0 over the poor and fair subgrades. This trend is consistent with the trends of the reinforcement rates that are presented in **Figure 10**, which reveals that the Rinter may have less influence on the reinforcement performance than the remaining factors and additional laboratory tests and field investigations are needed to comprehensively understand the interaction between a geogrid and a ballast.

# Double Geogrid Layers

The poor subgrade conditions were considered in Case 6. The corresponding reinforcement rates for each single geogrid layer and the double layers are plotted in **Figure 8A**. In this figure, "S" denotes the condition of a single geogrid layer in the model; "D" denotes double layers; and "200" and "300" denote the placement

depth of the geogrid layers inside the ballast. Better reinforcement performance was observed when the geogrid layer was placed at the ballast/sub-ballast interface when the single geogrid layer was considered in the simulation. The reinforcement rate increased with an increase in geogrid stiffness. These observations are consistent with **Figures 4**, **5**. However, better reinforcement performance was achieved by the double geogrid layers. Thus, an increase in the number of geogrid layers is an effective method for reinforcing the ballast over weak subgrades besides increasing the geogrid stiffness.

The (Fag )max values of Case 6 are plotted in **Figure 8B**. The value of (Fag )max in the geogrid of "S-300" was larger than the value of (Fag )max in the geogrid of "S-200," which indicated that better reinforcement results can be expected when the same geogrid layer was placed at the ballast/sub-ballast interface instead of inside the ballast over the poor subgrade. This observation corresponds to **Figure 5** and demonstrates the feasibility of (Fag )max in analyzing the geogrid reinforcement effect of a ballast.

As a result of the double geogrid layers, the reinforcement rate significantly increased. However, the reductions of (Fag )max in each geogrid layer of "D-200" and "D-300" are distinct compared with the reductions of (Fag )max in each geogrid layer of "S-200" and "S-300." As both of the geogrid layers contributed to the

load bearing, the axial forces in each geogrid layer reduced under the same loading, which is also beneficial for the load bearing capacity of the ballast and the geogrids. Thus, double or multiple geogrid layers are more efficient for the reinforcement of a ballast if cost is not considered.

# COMPARISON OF FIELD BEHAVIOR WITH NUMERICAL PREDICTIONS

## Modeling of Cyclic Loading Simulation

The calculation results of the monotonic loading revealed no distinct enhancements of the reinforcement effects over rocky and good subgrades and the possibility of collapsing failure for the poor subgrade under long-term repeated loads due to the relatively weak strength. Therefore, the track model over the fair subgrade was adopted in the cyclic loading to investigate the reinforcement effect of the geogrid under repeated train loads.

The sinusoidal cyclic load (as shown in **Figure 9**) was applied. The loading frequency was set to 15 Hz to simulate a train travel speed of 80 km/h (Indraratna et al., 2007). The minimum value of cyclic load was set to zero because the overburden pressure caused by the rail-sleeper assembly was already considered in the gravity generation by PLAXIS. The peak value of the cyclic load was set to 125 kN, which was equivalent to P<sup>s</sup> , as the impact caused by the different train travel speeds was already considered in the simulation by the loading frequency. The maximum loading number of cyclic load was 100,000 as the ballast usually enters the stable zone and the settlement of the ballast insignificantly changed after 100,000 cycles (Shin et al., 2002; Indraratna and Nimbalkar, 2013).

The geogrid stiffness of 525 kN/m, which was adopted in the cyclic simulation, was obtained from the actual mechanical properties of the biaxial geogrid with an aperture size of 40 mm (Indraratna et al., 2007). Cubical tests by Indraratna and Nimbalkar (2013) also reported an EA of 360 kN/m for a biaxial geogrid with an aperture size of 40 mm. The conditions of the unreinforced track model, the track models reinforced by a single geogrid layer (at the ballast/sub-ballast interface) and reinforced by two geogrid layers (similar to Case 6 in **Table 3**, i.e., D = 200 and 300 mm), were considered in the simulation.

# Results of Cyclic Loading Simulation

The calculation results of Stu, Str and R during the 100,000 load cycles are plotted in **Figure 10** including the result of load cycle of 1, 10, 100, 1,000, 10,000, 50,000, and 100,000. The settlement on the top of the ballast layer kept increasing during loading. The stable zone of the unreinforced track model cannot be observed within the 100,000 load cycles. Conversely, the settlement on the top of the ballast layer of the track model, which was reinforced by the single geogrid layer, attained a constant value after 10,000 load cycles. Better reinforcement performance is observed under the conditions of the track model that was reinforced by double geogrid layers, and a very slight settlement difference occurred during loading.

The logarithmic trend lines for the two series of reinforcement rates are shown in **Figure 10B**, and the regression equations are displayed in **Figure 10B**. Indraratna et al. (2011b) suggested a logarithmic equation for determining the rail track vertical deformation based upon the laboratory findings

$$S\_N = c + d(\ln N) \tag{3}$$

where S<sup>N</sup> is the settlement of the ballast accumulated after N load cycles, c is the settlement after the first load cycle, and d is an empirical constant that is dependent on the initial density, ballast and reinforcement type and moisture content. The simulation results in this study indicate that the reinforcement rate of the geogrid can be described by the similar equation

$$R\_N = a + b(\ln N) \tag{4}$$

where R<sup>N</sup> is the reinforcement rate of the geogrid layers after N load cycles and a and b are the empirical constants that depend on the mechanical properties of the ballast, the subgrade strength, the geogrid stiffness and the number of geogrid layers.

# Comparison With Field Data

The simulation results of the cyclic loading is compared with the measured data from a field trial at a site near Bulli in New South Wales, Australia (Indraratna et al., 2010). The settlements under the rail at the section of the fresh ballast with and without geocomposite reinforcement during the first 120,000 load cycles are plotted in **Figure 11A**. The reinforcement rate of the geocomposite, its logarithmic trend line and the regression equation are displayed in **Figure 11B**.

The comparison of **Figure 10** with **Figure 11** indicates that the simulation results for the ballast settlement of the unreinforced track over the fair subgrade was relatively larger than the field measurement data. The reinforcement rates for the geogrid layers, which were obtained by the cyclic loading simulations, were more pronounced than the field measurement data. This finding can be explained by the different mechanical properties of the track components, the physical properties of the geosynthetics (the EA value of the geocomposite was approximately 390 kN/m according to Indraratna and Nimbalkar, 2013) and the cyclic parameters of the train loads between the simulations and the field conditions. Thus,

additional studies are needed to improve the understanding of the complex relationship between the behavior of a ballast and these factors.

However, similar trends of the reinforcement rates for the geosynthetics are shown in **Figures 10**, **11**, i.e., both the results of the cyclic loading simulations and the field measurements can be successfully described by Equation (4). The parameters a and b in Equation (4) should be varied in every case due to the variation in the mechanical properties of the ballast, the subgrade strength, the geogrid stiffness and the number of geogrid layers. The results indicate the successful application of PLAXIS for simulating and predicting the behaviors of reinforced railway tracks under cyclic loads.

The double geogrid layers were more effective than the single geogrid layer and were beneficial for the load bearing capacity of the ballast and geogrids. The simulation results of the monotonic loading indicate that double or multiple geogrid layers are more efficient for the reinforcement of a ballast if cost is not considered.

# LIMITATIONS OF THIS STUDY

The variations of the parameters of subgrade strength were considered in the simulations, which revealed the important role of subgrade in the reinforcement of geogrids. However, the effect of the strength parameters and the constitutive models of ballast and sub-ballast warrant evaluation. The reduction in the

degradation of the ballast aggregates caused by the confinement of the geogrids has not been considered.

The strength reduction factor Rinter is influenced by the complex relationship between the aperture size of geogrids, the mean particle size of a ballast, the physical properties of geogrids and the mechanical behaviors of a ballast. Therefore, additional experimental research is needed to correctly define the behaviors of the interfaces in the simulations.

A complex relationship exists between the ballast behavior and cyclic loads. The degradation and breakage of a ballast were not considered in the simulation. Therefore, additional studies are needed to investigate the ballast behaviors and the reinforcement effect of geogrids under cyclic loads of varied amplitude and varied frequency. Three-dimensional simulations and Discrete Element Method (DEM) may provide a better understanding of the long-term geogrid reinforcement effects of railway ballast under cyclic train loads.

# CONCLUSIONS

The FE simulations that were conducted in this study investigated the geogrid reinforcement effect on a railway ballast that was subjected to monotonic loads. The reinforcement mechanisms of a geogrid and the influences of different factors, including the subgrade strength, the stiffness, placement depth and effective width of a geogrid; the strength reduction factor of the interfaces between a geogrid and a ballast; and the combination of double geogrid layers, were discussed based on the numerical simulation results. The following conclusions are obtained from this study:


# REFERENCES


beneficial to reinforcement performance. However, the effect of Rinter is relatively insignificant compared with the results of the factor of subgrade strength.


# DATA AVAILABILITY

All datasets generated for this study are included in the manuscript and/or the supplementary files.

# AUTHOR CONTRIBUTIONS

SN elaborated the paper guidelines and corrected the paper. YJ prepared the paper draft.

# ACKNOWLEDGMENTS

This study was made possible by the 2011 Endeavour Research Award (Endeavour Research Fellowship) from the Australian Government Department of Education, Employment and Workplace Relations (DEEWR), and help from Prof. Buddhima Indraratna and Associate Professor Cholachat Rujikiatkamjorn of University of Wollongong. The authors acknowledge their support.


element modelling. Granular Matter 19, 54: 51–16. doi: 10.1007/s100 35-017-0743-4


**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2019 Jiang and Nimbalkar. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Seismic Bearing Capacity of a Strip Footing Over an Embankment of Anisotropic Clay

#### K. Krishnan<sup>1</sup> , Koushik Halder <sup>1</sup> and Debarghya Chakraborty <sup>2</sup> \*

*<sup>1</sup> Department of Civil Engineering, Indian Institute of Technology, Kharagpur, India, <sup>2</sup> Department of Civil Engineering, Indian Institute of Technology, Kharagpur, India*

The present study adopts lower bound finite element limit analysis technique in association with non-linear optimization to compute the seismic bearing capacity of a strip footing placed over an embankment of anisotropic clay. The influence of seismic loading is incorporated in terms of a horizontal seismic acceleration coefficient. The bearing capacity factor *N<sup>c</sup>* of the strip footing is obtained for various combinations of (i) slope angle, (ii) footing setback distance, (iii) horizontal seismic acceleration coefficient, and (iv) anisotropic strength ratio of clay. The bearing capacity factor of the strip footing increases with the increasing setback distances. However, beyond a certain setback distance, hardly any influence is found. For a particular footing position and constant value of horizontal seismic acceleration coefficient, the magnitude of bearing capacity factor increases with the increase in anisotropic strength ratio. However, the percentage increment in the bearing capacity reduces when the value of anisotropic strength ratio becomes greater than one. In addition, the bearing capacity reduces with the increasing slope angle and horizontal seismic acceleration coefficient. Failure patterns are obtained to understand the failure mechanism.

Keywords: seismic bearing capacity, embankment, lower bound limit analysis, non-linear optimization, clay anisotropy

# INTRODUCTION

In recent times, shortage of construction land compels the people to build structures on the embankment; even on its edge or side face. Construction of shallow foundations for short storeyed buildings and bridge abutments, roadways, and railways on flood plains or riverbanks are the perfect examples of such structures. The presence of slope on both the sides of the footing reduces the bearing capacity of the footing, which significantly enhances the vulnerability of the structures (Georgiadis, 2009; Chakraborty and Kumar, 2013; Chakraborty and Mahesh, 2015; Leshchinsky, 2015; Halder and Chakraborty, 2018; Raj et al., 2018; Halder et al., 2019). The risk and vulnerability of the structures increases further under the seismic loading, generated during an earthquake (Kumar and Mohan Rao, 2003; Choudhury and Subba Rao, 2006; Yang, 2009; Kumar and Chakraborty, 2013; Chakraborty and Mahesh, 2015; Halder et al., 2018; Halder and Chakraborty, 2019). The literature review reveals that although there are many studies available for the computation of the bearing capacity of footing on the embankment under seismic loading, no one has considered the effect of anisotropy of clay on the bearing capacity of a strip footing placed on an embankment and subjected to seismic loading. The effect of anisotropy is very much significant for undrained clay (Jakobson, 1955; Lo, 1965; Bishop, 1966; Duncan and Seed, 1966; Davis and Christian, 1971).

#### Edited by:

*Sanjay Shrawan Nimbalkar, University of Technology Sydney, Australia*

#### Reviewed by:

*Subhadeep Banerjee, Indian Institute of Technology Madras, India Anindya Pain, Central Building Research Institute (CSIR), India*

> \*Correspondence: *Debarghya Chakraborty debarghya@civil.iitkgp.ac.in*

#### Specialty section:

*This article was submitted to Transportation and Transit Systems, a section of the journal Frontiers in Built Environment*

Received: *29 June 2019* Accepted: *31 October 2019* Published: *15 November 2019*

#### Citation:

*Krishnan K, Halder K and Chakraborty D (2019) Seismic Bearing Capacity of a Strip Footing Over an Embankment of Anisotropic Clay. Front. Built Environ. 5:134. doi: 10.3389/fbuil.2019.00134* Hence, the introduction of anisotropy in numerical analysis increases the accuracy of the analysis towards field condition. Davis and Christian (1971) formulated an anisotropic failure criterion (Equation 2) for cohesive soil in the form of an ellipse by using the anisotropic strength ratio (b/a) which is a function of undrained shear strengths su0, su45, and su90. In the present analysis, the effect of anisotropy of clay is incorporated using the formulation proposed by Ukritchon and Keawsawasvong (2018) using the yield function of Davis and Christian (1971). From here, throughout the article, DC indicates the yield criterion of Davis and Christian (1971). The advantage of DC yield function is associated with the consideration of three parameters (suo,

su90, and su45); whereas, most of the other anisotropic yield criterions are based on two parameters (suo and su90). The parameters mentioned above are the undrained shear strengths of clay measured in different orientations (δ): suo (plane strain compression, δ = 0 ◦ ), su<sup>90</sup> (plane strain tension, δ = 90◦ ), and su45(simple shear, δ = 45◦ ). It is also indicated by Ukritchon and Keawsawasvong (2018) that the three-parameter model simulates the actual failure envelope in a better way with respect to other models. The influences of other parameters investigated in the study are slope inclination angle (β), footing setback distance (S), and seismic acceleration coefficient (kh). With the consideration of all the parameters mentioned above, the present study computes the bearing capacity of a strip footing situated over an embankment of anisotropic clay.

# STATEMENT OF THE PROBLEM

The objective of the present study is to determine the bearing capacity factor N<sup>c</sup> of a rigid, rough and surface strip footing of width B placed over an embankment of anisotropic clay, in the presence of horizontal earthquake acceleration (khg) using the lower bound finite element limit analysis method. A non-linear optimization technique based on the second-order cone programming (SOCP) method is used. The embankment faces are assumed to incline at an angle β with the horizontal axis. The classical bearing capacity equation of Terzaghi (1943) reduces to Equation (1), as there is no presence of frictional property of soil and overburden pressure in the ground.

$$q\_{\mu} = Q\_{\mu}/B = s\_{\mu\alpha\nu} N\_c \tag{1}$$

where, Q<sup>u</sup> is the ultimate collapse load determined from lower bound limit analysis and suav = (su<sup>0</sup> + su90) is the average undrained shear strength of clay.

# PROBLEM DOMAIN, MESH DETAILS AND STRESS BOUNDARY CONDITIONS

The problem domain used in the present study is represented by **Figure 1A**. The domain is varied up to a length and depth of 7.5B−15B in the horizontal and vertical directions. **Figure 1A** also shows the stress boundary conditions employed on different boundary edges. With the absence of any overburden pressure, the normal and shear stresses are zero along the ground surface (MI and JN) and inclined faces (RM and NO) of the embankment. In order to account the footing roughness, an

equation: τxy <sup>≤</sup> <sup>s</sup>uav is adopted. The semi-infinite and perfectly plastic soil mass is assumed to follow associated flow rule and DC yield criterion. The effect of pseudo-static lateral earthquake force is considered in the study by introducing a lateral shear force (Q<sup>h</sup> = khQu) along the footing-soil interface IJ.

The finite element mesh is first created using ABAQUS software<sup>1</sup> . Then the coordinates of the mesh are transferred to MATLAB<sup>2</sup> . **Figures 1B,C** show the ABAQUS generated finite element mesh, which is later used in MATLAB in order to obtain the collapse load of footing. In **Figures 1B,C**, Ne, E, and Ndisc denote the number of nodes, elements and discontinuity edges, respectively.

# LOWER BOUND FORMULATION

The present study uses the lower bound finite element limit analysis formulation of Sloan (1988). The plane strain lower bound formulation of Sloan (1988) has certain unique qualities such as: (i) it requires only the shear strength parameters for solving any stability problem (ii) it gives safe and conservative solution of the collapse load (iii) the assembly of elements is easier due to the distinct node numbers for each element. As per the lower bound theorem, any statically admissible stress field that satisfies the (i) element equilibrium condition, (ii) stress boundary condition, and (iii) yield criteria, gives a safe solution of the true collapse load. Among the above three, the first two conditions are the equality constraints and the third one is an inequality constraint. In the finite element lower bound limit analysis formulation, the elements are assigned distinct node numbers. Hence, there arises a discontinuity between the edges of shared elements. The discontinuity of the stresses is considered by imposing an additional equality constraint. To solve any stability problem, the above four conditions need to be satisfied.

Equation (2) gives the DC yield condition for anisotropic clay below.

$$F = \sqrt{\left(\frac{\left(\sigma\_\circ - \sigma\_\text{x}\right)/2 - h}{a}\right)^2 + \left(\frac{\tau\_{\text{xy}}}{b}\right)^2} = 1\tag{2}$$

where, σ<sup>x</sup> and σ<sup>y</sup> are the normal stresses on x and y directions, respectively. h = (−su<sup>0</sup> + su90) 2 , a = su<sup>0</sup> + su<sup>90</sup> 2 , and b = √ asu45 su0su90 are

<sup>1</sup>ABAQUS 6.14 [Computer software]. Dassault Systemes, Vélizy - Villacoublay. <sup>2</sup>MATLAB 8.5 [Computer software]. MathWorks, Massachussets, USA.

the center, major semi-axes, and minor semi-axes of the elliptical yield surface formulated by DC yield criteria. The b/a ratio: the ratio between the minor and major axes, respectively, is used as a factor for anisotropic strength.

The objective function is generated by integrating the normal stresses along the footing-soil interface. The canonical form of the optimization scheme is given in Equation (3). A MATLAB code is written to form the canonical structures of the formulation; whereas, an optimization toolbox MOSEK<sup>3</sup> is used to carry out the non-linear optimization.

$$\mathbf{Maximize} : \mathbf{c}^T \boldsymbol{\sigma} \tag{3}$$

$$Subjected to: [A\_{\mathfrak{K}}] \{ \sigma \} = \{ B\_{\mathfrak{K}} \} \tag{3a}$$

$$[A\_{1SOP}]\{\sigma\} + [A\_{2SOP}]\{\xi\} = \{B\_{SOP}\} \tag{3b}$$

Where A<sup>g</sup> and B<sup>g</sup> are the global matrix and vector, respectively constituting equality constraints. A1SOCP, A2SOCP, and BSOCP are the yield constraint matrices and vector, respectively. {σ} is the nodal stress vector and {ξ} is the nodal slag variable vector.

# VARIATION OF BEARING CAPACITY FACTOR N<sup>C</sup> FOR DIFFERENT KH, β, S/B, AND B/A VALUES

The variation of bearing capacity factor N<sup>c</sup> with respect to k<sup>h</sup> is determined for various b/a ratios ranging from 0.8 to 1.2. The **Figures 2**–**4** show the results obtained for different β (=10◦ , 20◦ , 30◦ , and 40◦ ) and S/B (0, 2, and 4) values. The results indicate that as the footing moves away from the slope edge, the bearing capacity increases. However, after a particular setback distance of footing, there is no further increment in the value of bearing capacity factor. The reason behind this is the reduction in the influence of plastic zone over the slope when setback distance is increased. Similar results are also obtained in Chakraborty and Kumar (2015) and Halder and Chakraborty (2018). Results also illustrate that the bearing capacity increases with the increase in anisotropic strength ratio (b/a). It is noted that the increment in N<sup>c</sup> value is higher for b/a < 1. However, the increment becomes constant when b/a > 1. For a slope of β = 30◦ , S/B = 0, and k<sup>h</sup> = 0, the increase in N<sup>c</sup> value is 57.81% with the increase in the value of b/a from 0.8 to 1.0. Whereas, when b/a ratio is increased from 1.0 to 1.2, the percentage increase in the N<sup>c</sup> value is reduced to 9.41%. The decrease in N<sup>c</sup> value with the increase in β also is evident from the results. For example, when b/a ratio is reduced

<sup>3</sup>MOSEK ApS version 8.0 [Computer software]. MOSEK, Copenhagen, Denmark.

FIGURE 6 | Failure contours for *S/B* = 4 and *b/a* = 1 with (A) β = 30◦ , *k<sup>h</sup>* = 0; (B) β = 40◦ , *k<sup>h</sup>* = 0; (C) β = 30◦ , *k<sup>h</sup>* = 0.1; (D) β = 40◦ , *k<sup>h</sup>* = 0.1; and failure contours for S/B = 4 and *b/a* = 0.8 with (E) β = 30◦ , *k<sup>h</sup>* = 0; (F) β= 40◦ , *k<sup>h</sup>* = 0; (G) β = 30◦ , *k<sup>h</sup>* = 0.1; (H) β = 40◦ , *k<sup>h</sup>* = 0.1.

from 1.0 to 0.8 for a footing with S/B = 0, the reduction in N<sup>c</sup> value becomes higher from 13.75 to 45.43% for the increase in β value from 10 to 40◦ , respectively. In the same way when setback distance is varied from S = 0B−4B for a slope angle of 30◦ , the magnitude of N<sup>c</sup> value reduces from 36.63 to 5.85% for the reduction in b/a ratio from 1.0 to 0.8.

TABLE 1 | Comparison of available experimental results with the present study.


# FAILURE PATTERNS

The state of stress at any point with respect to collapse is determined in terms of F as given in Equation (2). As per the DC failure criterion, the soil fails or enters plastic state when F ≥ 1. When F < 1, the soil will remain in a nonplastic state. The failure mechanism of an embankment for S/B = 0 and S/B = 4 are shown in **Figures 5A–H**, **6A–H**, respectively. It is evident from **Figure 5** that when S/B = 0, the failure surface extends to the slope face easily. It is also clear that when b/a ratio is reduced from 1.0 to 0.8, the failure surface extends to a greater depth. The influence of seismic force is well-understood from the **Figures 5C,D,G,H**. With the consideration of earthquake loading, the failure surface shifts more toward the slope face in the direction of seismic force. From **Figures 6A–D**, it is seen that the plastic zone is developed locally within the horizontal surface when b/a = 1 (isotropic condition). However, when the b/a ratio is decreased to 0.8 as in **Figures 6E–H**, additional to the local failure due to the footing on horizontal surface, the plastic zone is developed in the slope as well.

# COMPARISONS

The values of N<sup>c</sup> calculated from the present study for isotropic condition (b/a = 1) are compared with the results of Chakraborty and Mahesh (2015) as shown in **Figure 7A**. The present results are slightly higher than that reported by Chakraborty and Mahesh (2015). It is due to the use of different optimization technique in the present study. Chakraborty and Mahesh (2015) determined the seismic bearing capacity of an embankment of isotropic soil using lower bound limit analysis with linear optimization. The present study is carried out with the utilization of lower bound limit analysis in combination with non-linear optimization. However, to verify the results obtained for the anisotropic condition, the N<sup>c</sup> values of strip footing are determined for b/a ranging from 0.5 to 1.2 over a horizontal ground surface (β = 0 ◦ ) using the DC yield criterion as per Ukritchon and Keawsawasvong (2018). The present results match closely with the results of Ukritchon and Keawsawasvong (2018) as shown in **Figure 7B**. Additionally, the comparative study of Davis and Christian (1971) with the experimental results of Jakobson (1955), Lo (1965), Duncan and Seed (1966), and Bishop (1966) is considered to validate the numerical precision of the present study. Davis and Christian (1971) calculated the N<sup>c</sup> values using the b/a ratios corresponding to the experimental shear strength results of Jakobson (1955), Lo (1965), Bishop (1966), and Duncan and Seed (1966). The N<sup>c</sup> values of the present work are compared with the N<sup>c</sup> values of Davis and Christian (1971) for different b/a ratios in **Table 1**. The comparison shows that the present study match well with the calculated N<sup>c</sup> values of Davis and Christian (1971).

# CONCLUSIONS

By incorporating lower bound finite element limit analysis with second-order conic optimization technique, the seismic bearing capacity of a strip footing over an embankment of anisotropic clay is determined. The variation of bearing capacity factor N<sup>c</sup> is presented as a function of slope angle (β), setback distance of footing (S), seismic coefficient (kh), and anisotropic strength ratio (b/a). Design charts are proposed based on the results of Nc . The design charts can be used by the practicing engineers to analyse the seismic bearing capacity of a strip footing over an embankment of anisotropic clay using the anisotropic strength ratio (b/a).

It is found that the bearing capacity factor N<sup>c</sup> decreases with the increase in the value of β and kh. The trend of the reduction in the value of N<sup>c</sup> with the increase in k<sup>h</sup> for a b/a value greater than or equal to one is almost linear. However, when b/a < 1, the trend of the reduction in the value of N<sup>c</sup> becomes non-linear.


# REFERENCES


than the major axis(a), making the direction of anisotropy to change. For that reason, whenever the anisotropic strength ratio is lesser than 1, the strength reduces greatly.

• When anisotropic condition is coupled with earthquake force, the N<sup>c</sup> values of strip footing situated over an embankment reduces drastically. Hence, more emphasis on anisotropy needs to be given while analyzing the seismic bearing capacity of strip footing over an embankment.

# DATA AVAILABILITY STATEMENT

The datasets generated for this study are available on request to the corresponding author.

# AUTHOR CONTRIBUTIONS

KK developed the codes in MATLAB for carrying out all the numerical simulations required in the present study and prepared the initial draft. KH and DC proposed the project. KH modified the initial draft. DC supervised all the works carried out in this project and given a final shape of the paper.


**Conflict of Interest:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2019 Krishnan, Halder and Chakraborty. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# NOTATIONS


# Mitigation of Traffic Induced Vibration Using Geocell Inclusions

#### Amarnath Hegde\* and Hasthi Venkateswarlu

*Department of Civil and Environmental Engineering, Indian Institute of Technology Patna, Patna, India*

This paper describes the potential use of geocell reinforcement in mitigating the traffic induced vibration. The vibration caused by the vehicular movement was simulated over the unreinforced and geocell reinforced sections using a mechanical oscillator. The displacement amplitude and peak particle velocity were measured to understand the vibration mitigation efficacy of geocell. The effect of depth of placement of geocell on the mitigation of vibration parameters was studied. The inclusion of geocell was found effective in reducing the induced vibration based on the experimental results. The vibration mitigation efficacy of geocell was improved significantly at the shallow depth of placement of geocell mattress. The improvement in elasticity of the subgrade was observed maximum when the geocell was placed at a depth of 0.1*B* from the ground surface. Further, analytical and numerical approaches were used to predict the displacement amplitude vs. frequency response of reinforced soil sections. FLAC3D was used for performing the numerical investigation. The geocell was modeled according to its honeycomb shape to acquire the accurate response of geocell reinforced section. Whereas, mass spring dashpot analogy was followed for the analytical evaluation. In overall, the amplitude response predicted from the numerical and analytical studies were found to be in good agreement with the experimental results.

#### Edited by:

*Sanjay Shrawan Nimbalkar, University of Technology Sydney, Australia*

#### Reviewed by:

*Yajun Jiang, Southwest Jiaotong University, China Anindya Pain, Central Building Research Institute (CSIR), India*

#### \*Correspondence:

*Amarnath Hegde ahegde@iitp.ac.in*

#### Specialty section:

*This article was submitted to Transportation and Transit Systems, a section of the journal Frontiers in Built Environment*

> Received: *28 August 2019* Accepted: *01 November 2019* Published: *19 November 2019*

#### Citation:

*Hegde A and Venkateswarlu H (2019) Mitigation of Traffic Induced Vibration Using Geocell Inclusions. Front. Built Environ. 5:136. doi: 10.3389/fbuil.2019.00136* Keywords: geocell, traffic vibration, amplitude reduction ratio, peak particle velocity, FLAC3D, MSD analogy

# INTRODUCTION

The mitigation of intense levels of ground vibration generated by the rail and road networks is a common concern in urban areas. These vibrations can possibly exhibit adverse effects on the nearby structures, sensitive equipment, technical processes, and inhabitants (Murillo et al., 2009). The intensity of vibration is majorly influenced by both traffic and rail or road performance. The vehicle weight and speed are the major traffic characteristics. Structure/subgrade stability and roughness are the performance parameters. The frequency of vibration emanated from different vehicles is summarized in **Table 1**.

The vibration energy emanated from the traffic is transmitted through the ground in the form of surface and body waves (Woods, 1968; Ujjawal et al., 2019). In an elastic half space, surface waves are exclusively propagate along the surface, and body waves propagate in the form of spherical wave front in all directions. As a result, the body waves attenuate earlier than the surface waves. Miller and Pursey (1955) reported the distribution of total induced energy as 67% Rayleigh waves, 26% shear waves, and 7% compression waves. DIN 4150 (1999) recommended the tolerance limits of vehicle induced vibration for safeguarding the structures and human beings. Similarly, various studies reported the problems associated with the traffic induced vibration. Paolucci et al. (2003) observed the level of discomfort caused by the train induced vibration by comparing the field measurements

TABLE 1 | Frequency of vibration induced by different vehicles (Barneich, 1985).


with the limits suggested in DIN 4150 (1999). The train generated vibrations was found to reach the limits from "troublesome to persons" to the "severe to persons" occasionally.

Various studies reported that the amplitude of vibration is a function of speed of the vehicle and irregularity of the road (Watts and Krylov, 2000; Crispino and D'apuzzo, 2001; D'Apuzzo, 2007). In addition, heavy vehicles was found to produce the most of perceptible levels of vibrations. On the other hand, several methods have been recommended for attenuating the amplitude of induced vibration. The open trench approach was found the first and efficient approach for the mitigation of vibration in active and passive methods (Woods, 1968). However, to prevent the instability of sidewalls of an open trench, it was filled with different geo-materials (Thompson et al., 2016). The bentonite, water, concrete, soil and bentonite mixture, and EPS (expanded polystyrene) geofoam were the examples of such materials.

Alzawi and El Naggar (2011) studied the vibration mitigation efficacy of open and soft-filled trenches through field studies. Çelebi et al. (2009) reported the performance of concrete filled trench based on field measurements. Kim et al. (2000) highlighted the isolation behavior of rubber chips filled open trench. Massarsch (2005) reported the brief review about the screening efficacy of gas cushions. Various researchers evaluated amplitude reduction ratio (ARR) to quantify the screening efficacy of open and infilled trenches (Woods, 1968). The ARR is the ratio between displacement amplitude observed in the presence of barrier system to the displacement amplitude of without barrier system. Similarly, threshold limits were proposed to assess the level of damage experienced by different buildings due to traffic vibration. The threshold limits of ground vibration intensity for different buildings is listed in **Table 2**. Amick and Gendreau (2000) reported the standard criteria for the safe operation of buildings based on field recordings of train induced vibration. It was also stated that the building damage criteria does not essentially indicate that the structure is free from trifling turbulences. Wiss (1981) observed that the threshold range of vibration based on human perception gets "disturbing" at 7 mm/s and "very disturbing" at 25 mm/s. Thus, the ideal vibration mitigation system is aimed to attenuate the vibration levels in order to satisfy two different criteria namely, human discomfort and damage to the buildings.

On the other hand, trench systems are expensive and not a viable approach to protect the structures located in the close proximity to the vibration source. In such cases, subgrade stiffening is an alternative approach of mitigating the ground induced vibration. Few researchers suggested that the subgrade stiffening not only mitigate the traffic vibration but also reduces the track settlement and deflections (Ekanayake et al., 2014; Thompson et al., 2015). The techniques adopted TABLE 2 | Threshold limits of PPV for different structures due to traffic vibrations (Amick and Gendreau, 2000).


for stiffening the subgrade are vibro-compaction, jet grouting, vibro-replacement, excavation and replacement, stabilization, and reinforced earth (Thompson et al., 2015; Venkateswarlu et al., 2018a). Coulier et al. (2015) investigated the efficacy of jet grout columns in mitigating the rail induced vibration through the field and numerical investigation. Various studies reported the benefits of changing the soil profile of the foundation bed in reducing the amplitude of machine vibration (Baidya and Rathi, 2004; Baidya et al., 2006; Mandal et al., 2012). In addition, the significant improvement in natural frequency of the foundation bed was observed in the presence of stiff layer nearer to the surface footing. Few studies evaluated the reduction in ground vibration by the stabilization of soil bed (Mitchell, 1981; Baker, 1982; Welsh, 1986; Saride and Dutta, 2016).

Currently, the reinforced earth technique being popularly used for strengthening the subgrade in order to support the static and cyclic loads (Hegde, 2017). Boominathan et al. (1991) and Haldar and Sivakumar Babu (2009) studied the potential benefits of this method in reducing the vibration induced by industrial machines. In addition, the use of geosynthetics was well-studied for strengthening the track performance (Biabani and Indraratna, 2015; Biabani et al., 2016; Nimbalkar and Indraratna, 2016). However, limited studies have been addressed the vibration isolation ability of geosynthetics reinforced foundation beds. Hegde and Sitharam (2016) reported that the geocell reinforcement enhances the elastic response and natural frequency of the foundation bed. The numerical evaluation of Azzam (2015) found that the presence of confined cell could effectively mitigate the amplitude of vibration through increasing the damping of a subgrade. Venkateswarlu et al. (2018a) reported that increase in elastic response of foundation bed due to the inclusion of geocells. Few studies reported the potential using the geocell in mitigating the lateral spreading of machine induced vibration (Venkateswarlu and Hegde, 2018; Venkateswarlu et al., 2018b).

Based on existing literature, the vibration mitigation efficacy of geocell reinforced subgrade system is not completely understood. The present study is aimed to demonstrate the geocell potential in mitigating the traffic induced vibration. The vibration caused by vehicular movement was generated using oscillator assembly. The influence of depth of placement of geocell on the vibration mitigation efficacy of geocell reinforced subgrade has been investigated. In addition, numerical and analytical approaches have been demonstrated for predicting the amplitude vs. frequency response of different subgrade sections.

# MATERIALS

# Subgrade Soil

The particle size distribution of the soil used for the preparation of different subgrade sections is shown in **Figure 1A**. It consists of 84% coarse fraction and 16% fines content followed by the silt and clay compositions are 11 and 5%, respectively. As per Unified Soil Classification System, soil was classified as silty sand having the group symbol SM. The experimentally determined physical and mechanical properties of SM are listed in **Table 3**.

# Geocell

In this study, the NPA (Novel Polymeric Alloy) geocell mattress with cell pocket having equivalent diameter 250 mm was used. The cell walls were perforated and textured with rhomboidal shape indentation to mobilize friction with infill soil. The maximum tensile load capacity of geocell was determined from the stress vs. strain response as shown in **Figure 1B**. The ultimate load and failure strain were observed as 23.8 kN/m and 12.1%, respectively. The test was performed in accordance with ISO, E. 10319 (2015). The other properties of the geocell reinforcement are listed in **Table 3**.

# EXPERIMENTAL INVESTIGATION

# Experimental Setup

The test setup for inducing the traffic vibration over subgrade sections is shown in **Figure 2**. During the test, sinusoidal dynamic excitation replicating the traffic vibration was generated using mechanical oscillator. It was mounted over the loading plate. The loading plate helps to transfer the induced vibration from the oscillator to the subgrade section. It was made up of concrete and having the dimensions of 600 mm<sup>2</sup> and 200 mm thickness. The combined assembly of plate and oscillator was placed at the center of reinforced section. The 6 HP capacity DC motor was attached with the oscillator to induce the vibration at a required frequency. The frequency was adjusted and measured using speed control device with the use of speed measuring


sensor. The maximum resolution of the sensor was 10,000 RPM. The vibration response was measured in terms of displacement amplitude and peak particle velocity. To record the amplitude of vibration, vibration meter with accelerometer assembly was used. This assembly was designed for the continuous measurement of displacement amplitude. The accelerometer used in this study can measure the acceleration corresponding to vertical mode. The 3D geophone was used to measure the peak particle velocity (PPV) of the applied excitation. The measured PPV is the sum of the particle velocities measured in three orthogonal directions. To record the measurements of geophone, vibration monitoring terminal was used.

# Preparation of Subgrade Sections

Two different subgrade sections, namely, unreinforced and geocell reinforced subjected to traffic induced vibration were studied. Both the sections were prepared with the dimensions of 2 m (length) × 2 m (width) × 0.5 m (depth) at the field. The test sections were prepared in five numbers of layers with thickness of each layer equal to 100 mm. To do so, manual compaction mode was followed using a rammer of 12 kg. The subgrade soil was compacted at optimum moisture content (OMC) pertaining to the standard proctor as per IRC: 37 (2012). Before using the soil for compaction, it was mixed with OMC and allowed to reach the state of maturation. To establish uniform condition of the sections, predetermined amount of blows were applied over the each layer using the height of fall of 500 mm. The number of blows were determined in accordance with the compaction energy of standard proctor. The specifications of IS 2720-29 (1975) were followed to study the dry density variation along the longitudinal and lateral directions of the sections. Total, 9 numbers of soil samples were collected from left (L), center (O), and right (R) locations of the subgrade as shown in **Figure 3a**. **Figures 3b,c** shows the variation of dry density and water content at various locations of the sections. The dry unit weight of the subgrade was found to vary between 17.26 and 17.44 kN/m<sup>3</sup> . Similarly, the variation in optimum moisture content was found in the range of variation was 12.1 ± 4%.

The aforementioned procedure was followed up to the placement of geocell in the preparation of geocell reinforced subgrade section. Over the compacted surface, geocell mattress was positioned and expanded using metallic stacks. Each pocket of geocell was filled using SM material. The total height of geocell was compacted in three layers with the help of tamping rod. In addition, the suitable precautions were followed to protect the cell walls from bending and distortion as reported by Biswas and Krishna (2017) and Hegde (2017). After compacting all the geocell pockets, soil cover with the thickness equal to the depth of placement of geocell was provided. The test section with the partially filled geocell mattress is shown in **Figure 3d**. The dynamic excitation was applied over the test sections corresponding to the frequency of 30 Hz. It replicates the vibration caused by different vehicles as listed in **Table 1**.

# Experimental Results

Primarily, the displacement amplitude and PPV were measured up to the distance of 5 m from the vibration source. Total, 10 monitoring points were selected at an interval of 0.5 m. The efficacy of geocell in controlling the displacement amplitude of vibration was quantified in terms of amplitude reduction ratio (ARR). It is defined as the ratio between amplitude of vibration observed in the reinforced case to the amplitude of unreinforced case. The minimum ARR is generally recommended for better screening of the ground vibration. The ARR vs. normalized distance is shown in **Figure 4**. The normalized distance is obtained by dividing the distance from the vibration source (d) with the width of loading plate (B). The ARR was significantly diminished with the decrease in depth of placement of geocell mattress (U) from the surface of subgrade. The increase in U resulted in the amplification of ARR. Moreover, ARR found to decrease with the increase in normalized distance regardless of the depth of placement of geocell. The increase in ARR was observed beyond the d/B ratio of 5 in all the cases. The variation in soil conditions at the subgrade section and away from the subgrade might be the cause for the amplification of ARR.

The effect of d/B and U on the PPV is shown in **Figure 5**. The efficacy of geocell in attenuating the PPV of the foundation bed was decreased with the increase in depth of placement. The maximum attenuation of PPV was observed more than 48% at the depth of placement of 0.1B. Thus, the presence of geocell nearer to the surface is more beneficial for the effective attenuation of vibration. In general, the subgrade performance under the transit loading conditions is majorly depends on two parameters namely, dynamic shear modulus and elasticity of the bed. Increase in modulus enriches the stiffness of the bed. Whereas, the elasticity controls the deformation of the subgrade under the continuous action of cyclic stresses. Thus, the change in these parameters with the increase in U/B was studied. The elasticity of the bed was quantified using elastic uniform

placement of geocell.

compression (Cu). It was estimated by,

$$C\_u = K \times A \tag{1}$$

$$f\_r = \frac{1}{2\pi} \sqrt{\frac{\text{Kg}}{\text{W}}} \tag{2}$$

where K is the stiffness of the subgrade, A is the contact area between the loading plate and surface of the subgrade, f<sup>r</sup> is the resonant frequency of a system, W is the total weight of vibrating mass, and g is the acceleration due to gravitational force. To measure the f<sup>r</sup> , variation in displacement amplitude was observed for unreinforced and geocell reinforced cases by varying the frequency from 0 to 45 Hz. The displacement amplitude becomes maximum at the f<sup>r</sup> . The observed resonance parameters of different subgrade sections are listed in **Table 4**.

TABLE 4 | Variation of resonance parameters.


Based on f<sup>r</sup> , K and C<sup>u</sup> were evaluated. Similarly, the shear modulus (G) was determined using the following relation as suggested by Timoshenko and Goodier (1970).

$$K = \frac{4Gr\_0}{1 - \vartheta} \tag{3}$$

where r<sup>0</sup> is the equivalent radius of the loading plate. The value of Poisson's ratio (ϑ) was considered as 0.3 while determining the shear modulus (G). The variation of G and C<sup>u</sup> for different cases is shown in **Figure 6**. The maximum improvement in both the parameters was observed when the geocell located at the placement of 0.1B. It was attributed due to the increase in natural frequency of the subgrade system by the additional confinement offered by geocell. Hegde and Sitharam (2016) also reported the improvement in elastic response in the presence of geocell based on the results of cyclic plate load test. Further, the displacement amplitude vs. frequency response of different reinforced sections was studied through the numerical and analytical methods. The description about both the approaches are explained in the subsequent sections.

## NUMERICAL MODELING

Finite difference based three dimensional package FLAC 3D was used to conduct the numerical analysis. The numerical methodology was followed in two major steps. Primarily, the subsurface profile up to a depth of 10 m was simulated using brick element. The subsoil conditions may also have significant influence on the dynamic response of any soil system (Gazetas, 1991). In addition, selection of mesh size, boundary distance, and conditions can significantly affect the accuracy of the model. Thus, the sensitivity analysis was carried out systematically to determine the optimal values of these parameters. Based on the results of sensitivity study, 12B was found an optimum boundary distance from the dynamic actuation for not affecting the dynamic response even at higher frequencies. In addition, mesh density was found to have minimum influence on the results. Hence, the numerical model having the dimensions of 15 × 15 × 10 m was developed. The coarse mesh was selected for discretizing the model. At the ground surface, the subgrade sections of 2 × 2 × 0.5 m was simulated as similar to the experimental study. Venkateswarlu et al. (2018a) reported the details of subsurface profile and their modeling parameters. The loading plate behavior was simulated using linear elastic material. Whereas, the response of subgrade soil was modeled using Mohr-Coulomb yield criteria. The displacement along the bottom plane of the model was restrained in all the three directions (i.e., U<sup>x</sup> = U<sup>y</sup> = U<sup>z</sup> = 0). The vertical faces were restrained only in the horizontal direction. Therefore, the displacement was permitted in the vertical direction (i.e., U<sup>x</sup> = U<sup>y</sup> = 0 and U<sup>z</sup> # 0). In addition, viscous (quiet) boundary conditions were assigned to the vertical faces to prevent wave reflections from the boundary of the model.

**Figure 7A** shows the FLAC3D model for unreinforced case with the boundary conditions. To simulate the geocell-reinforced

section, geocell was created at the required position as shown in **Figure 7B**. The actual honeycomb curvature was considered while developing the geocell. The geocell and infill material were assigned with linear elastic and Mohr-Coulomb constitutive behavior, respectively to simulate the real case conditions. Geocell was modeled using the geogrid structural element. The modeling parameters used for the development of geocell reinforced bed are listed in **Table 5**. In the second stage, the dynamic excitation similar to experimental study was applied over the loading plate using FISH code. To monitor the vibration parameters, zonal points were selected along the longitudinal direction. The variation in dynamic force with the increase in frequency is shown in **Figure 7C**. The figure represents the magnitude of dynamic force acted over the loading plate at different operating frequencies. The dynamic force corresponding to the each operating frequency was applied for the duration of 10 s.

TABLE 5 | Properties of different materials used in modeling for dynamic loading condition.


# ANALYTICAL STUDY

Various studies reported the use of mass spring dashpot (MSD) analogy for predicting the response of pavements and rail track systems. Loizos et al. (2003) carried out the dynamic analysis of pavement subgrade using MSD approach. Choudhury et al. (2008) suggested MSD system with 2-degree of freedom to predict the displacement magnitude of subgrade and ballast layers. The unique advantage of this method is that the consideration of damping and stiffness characteristics during the analysis. In this study, single degree MSD analogy was used to predict the displacement amplitude vs. frequency response of different subgrade sections. Based on MSD model, the subgrade soil was considered as an isotropic and elastic soil medium. It was characterized using linear elastic weight less spring. The dashpot system was used to replicate the damping behavior of the subgrade material. The dynamic excitation applied over the reinforced sections was considered as vertical during the analogy. The loading plate and oscillator assembly was represented with the rigid mass of M. **Figure 8** shows the MSD idealization of subgrade subjected to the vertical mode dynamic excitation.

The generalized equation of motion used to describe the system can be written as,

$$M\ddot{\mathbf{Z}} + \mathbf{C}\dot{\mathbf{Z}} + \mathbf{K}\mathbf{Z} = X\left(t\right) \tag{4}$$

where X(t) is the total vertical dynamic force, Z is the displacement amplitude, Z˙ and Z¨ are the first and second derivatives of Z over time t, C is the damping coefficient, and K is the equivalent stiffness. The variation in displacement amplitude (Z) of various sections at each frequency is determined using

$$Z = \frac{\left(\frac{m\_c e}{M}\right) \left(\frac{\omega}{\omega\_n}\right)^2}{\sqrt{\left(1 - \left(\frac{\omega}{\omega\_n}\right)^2\right)^2 + \left(2D\left(\frac{\omega}{\omega\_n}\right)\right)^2}}\tag{5}$$

where m<sup>e</sup> is the mass of the rotating elements, ω<sup>n</sup> is the natural frequency of the subgrade condition, ω is the operating frequency of dynamic excitation, and D is the damping ratio. The displacement amplitude becomes maximum at the resonance condition (when the ω<sup>n</sup> matches with the ω). The peak displacement amplitude (Zm) of the vibration at the resonance is determined by,

$$\frac{Z\_m M}{m\_e e} = \frac{1}{2D\sqrt{1 - D^2}}\tag{6}$$

Comparison of experimental results with the results of numerical and analytical studies for unreinforced section is shown in **Figure 9A**. The good agreement was noticed among the results of experimental, numerical and analytical studies. However, the numerical results has shown close agreement with the experimental results as compared to the analytical results. The

similar observation was also noticed in the case of geocell reinforced subgrade as shown in **Figure 9B**. The analytical model has exhibited a slight deviation in predicting the post resonance behavior of the geocell reinforced subgrade. It might be due to the lack of consideration of parameters representing the interaction between the soil and geocell reinforcement. In addition, the subgrade is generally considered as a linear elastic weightless spring in the case of analytical study. In spite of these limitations, the analytical model has predicted the resonance response of the geocell reinforced case reasonably well. The reported amplitude vs. frequency is corresponding to the optimum (U = 0.1B) geocell reinforced case. About 7 and 12% variation in the resonant frequency was observed between the experimental and analytical studies for unreinforced and geocell reinforced sections. Whereas, the deviation between numerical and experimental resonant frequency was noticed as 3 and 2% for unreinforced and geocell reinforced subgrade cases.

The influence of geocell properties, namely, elastic modulus and the width of geocell on the ARR and natural frequency (fn) of the subgrade was numerically investigated and presented in **Figures 10A,B**. The geocell modulus was varied as 0.25, 0.5, 1,

and 2 times Young's modulus (E) of geocell used in the present study. The E of the geocell used in the present study was 280 MPa. The width of geocell was considered as 3.3B while studying the effect of geocell modulus. From **Figure 10A**, the reduction in ARR and improvement in f<sup>n</sup> was found with the increase in modulus of the geocell mattress. It was due to the increase in additional confinement offered by the geocell with the increase in its modulus. Further, to study the effect of geocell width, it was increased from 1.6B to 6.7B with an increment of 1.7B. The geocell modulus was considered as 1E while studying the effect of width of geocell mattress. The increase in geocell width resulted in the decrease in ARR and increase in fnof the subgrade as shown in **Figure 10B**. The improvement in f<sup>n</sup> was found marginal beyond the geocell width of 5B.

# CONCLUSIONS

The efficacy of geocell in mitigating the traffic induced vibration was systematically studied in the present study. The vibration caused by traffic was generated over the subgrade sections with and without geocell reinforcement using the oscillator. The reduction in vibration parameters namely, displacement amplitude and peak particle velocity were measured to assess the geocell efficiency. From the experimental results, the reduction in peak particle velocity of subgrade was decreased from 48 to 34% with the change in depth of placement of geocell from 0.1B to 0.5B. Thus, 0.1B was suggested as an optimum depth for the

# REFERENCES


effective mitigation of traffic vibration. At the optimum depth of placement, 1.4 times increase in natural frequency of subgrade was observed. The elasticity of the subgrade was improved by 96%. Similarly, the peak displacement amplitude of the subgrade was reduced by 57% in the presence of geocell reinforcement. The numerical and analytical methods were successfully used to predict the displacement amplitude response of the subgrade sections considered in the present study. The amplitude vs. frequency responses predicted from the analytical and numerical studies have shown good agreement with the experimental results. However, the resonant frequency predicted from the analytical study has exhibited the slightly higher deviation as compared to the numerical results. In overall, the inclusion of geocell not only mitigate the traffic induced vibration but also enriches the dynamic properties of the subgrade section.

# DATA AVAILABILITY STATEMENT

The datasets generated for this study are available on request to the corresponding author.

# AUTHOR CONTRIBUTIONS

HV performed the studies and prepared the manuscript based on the inputs and guidance of the AH. AH and HV reviewed and accepted the final version.


**Conflict of Interest:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2019 Hegde and Venkateswarlu. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# NOMENCLATURE

The following notations are used in the present study.


# Engineering Characteristics and Stabilization Performance of Aggregate Quarry By-Products From Different Sources and Crushing Stages

Wenting Hou<sup>1</sup> , Issam Qamhia<sup>1</sup> , Vincent Mwumvaneza<sup>2</sup> , Erol Tutumluer <sup>1</sup> \* and Hasan Ozer <sup>3</sup>

*<sup>1</sup> Department of Civil and Environmental Engineering, University of Illinois at Urbana-Champaign, Urbana, IL, United States, <sup>2</sup> Dynatest North America, Alpharetta, GA, United States, <sup>3</sup> Arizona State University, Tempe, AZ, United States*

Quarry by-products (QB), usually <1/4 in. (6 mm) in size, are the residual deposits from the production of required grades of aggregate. This paper provides findings of a detailed laboratory study with the objective of characterizing the engineering properties of QB materials produced in the primary, secondary, and tertiary aggregate production stages from four different quarries operating in Illinois. Property tests were conducted for determining aggregate gradation, morphological shape characteristics, compaction properties (moisture-density), chemical composition, and strength properties of QB samples. Since the unconfined compressive strength for QB materials is relatively low, chemical admixture stabilizers such as Portland cement and Class C fly ash were used to improve the strength properties. This study aims at evaluating properties governing the untreated and stabilized strength of QBs such as source variation, compacted density, chemical composition, gradation, particle shape and angularity, as well as the uniformity of distribution and the effectiveness of stabilizer. QB samples treated with 2% cement or 10% Class C fly ash by dry weight were found to be 10–30 times stronger than the virgin QB samples. Such significant increases in the strength of stabilized QB materials observed may indicate suitability of QBs for sustainable pavement applications.

Keywords: aggregate, quarry by-product, engineering properties, admixture stabilization, unconfined compressive strength, sustainable pavement applications, aggregate packing, X-ray fluorescence

# INTRODUCTION

During aggregate production, multiple aggregate quarry processes such as blasting, crushing, and screening of aggregates produce large amounts of by-product mineral fine materials, commonly known as Quarry By-products (QBs). During the crushing stages, QBs are generally produced in three stages: primary crushing, secondary crushing, and tertiary crushing (Petavratzi and Wilson, 2007). Common types of crushers used to produce crushed aggregate materials in the primary, secondary, and tertiary crushing stages include jaw crushers, cone crushers, gyratory crushers, and impact crushers. The quantity of fines (and QB) produced from each crusher type depends on the crushing stage and the operating conditions. Jaw crushers are normally utilized for the primary crushing. Cone crushers are widely used by crushed rock producers due to the relatively

Edited by:

*Yifei Sun, Hohai University, China*

#### Reviewed by:

*Libin Gong, University of Wollongong, Australia Pei Tai, Hong Kong Polytechnic University, Hong Kong*

> \*Correspondence: *Erol Tutumluer tutumlue@illinois.edu*

#### Specialty section:

*This article was submitted to Transportation and Transit Systems, a section of the journal Frontiers in Built Environment*

> Received: *11 September 2019* Accepted: *21 October 2019* Published: *19 November 2019*

#### Citation:

*Hou W, Qamhia I, Mwumvaneza V, Tutumluer E and Ozer H (2019) Engineering Characteristics and Stabilization Performance of Aggregate Quarry By-Products From Different Sources and Crushing Stages. Front. Built Environ. 5:130. doi: 10.3389/fbuil.2019.00130*

**58**

lower operational costs and low fines (QB) generation. Horizontal shaft impact crushers are more commonly used for processing softer materials such as limestone, while vertical shaft impact crushers have more capability of producing cubical aggregates needed for concrete production. With vertical impact crushers, the quantities of fines generated during aggregates production can be limited to acceptable levels if the right operation techniques are implemented (Tarmac Ltd. Partners, 2011).

QBs are typically <¼ in. (6 mm) in size and consist of mostly coarse, medium and fine sand particles, and a small amount of clay/silt fraction. Research conducted in the early 1990s showed that stockpiled fines comprised an average of approximately 12% of the total annual aggregate production of the surveyed companies (Kumar and Hudson, 1992). Hudson et al. (1997) stated that the stockpiling and disposal of aggregate by-products is a major problem facing the aggregate industry. According to the recent NCHRP Synthesis 435 (volume 4), and depending on the type of rock quarried, quarry by-products can be up to 25% of the total aggregate produced (Stroup-Gardiner and Wattenberg-Komas, 2013). A more recent industrial survey distributed to aggregate producers in the State of Illinois in 2013 showed that, among 42 respondents, 90% of the respondents are producing QBs. Among the quarries that produce QBs, 55% of the respondents produce a typical annual amount of QB >100,000 short tons; 26% between 25,000 and 100,000 tons, and 19% <25,000 tons (Mwumvaneza et al., 2015).

Hence, it would be of value to evaluate the properties of QBs and find potential application areas in the construction and rehabilitation of transportation infrastructure. A number of previous researches have focused on laboratory property testing of QBs. Kalcheff and Machemehl (1980) conducted particle size distribution test for different types of QBs. It was reported that screenings generally contain freshly fractured faces, have fairly uniform gradation, and contain fewer plastic fines. The particle distributions of the evaluated QBs followed a similar gradation trend, with particles smaller than sieve No. 200 (0.075 mm) ranging from 6 to 12%. Stokowski (1992) showed that QBs are enriched in CaCO3, SiO2, Al2O3, and Fe2O<sup>3</sup> relative to MgCO3. As a result of this finding, QBs have lower specific gravity and are relatively soft because of calcite (CaCO3) and enrichment of clay minerals (SiO2, Al2O3, and Fe2O3).

According to Dumitru et al. (2001), mineralogical tests such as X-ray diffraction analysis should be conducted to determine the compositions of secondary minerals and to quantify amounts of harmful content that can be detrimental in some applications. Petavratzi and Wilson (2007) concluded that QBs are composed of the same mineral substances as the soil and solid rock from which they are derived, and they are usually inert or nonhazardous. Puppala et al. (2008) reported the plasticity index properties of QBs and also concluded that the compressive strength of untreated (virgin) QBs can be very low. Puppala et al. (2012) conducted one-dimensional vertical free-swelling test of QB. In their study, the QBs were found to be of moderately swelling potential. Mwumvaneza et al. (2015) conducted the unconfined compressive strength test of QBs and found the strength to be relatively low.

Based on laboratory testing results, some researchers utilized chemical stabilization for improving strength of QBs; and based on their findings, different field applications for QBs were recommended. According to Kalcheff and Machemehl (1980), stabilization of QBs with cement developed relatively high rigidity with a small amount of Portland cement compared with granular soil–cement stabilization. The use of low-cement content has the advantage of decreasing the shrinkage cracking. In 1992, Kumar and Hudson examined the unconfined compressive strength, tensile modulus of elasticity, and Poisson's ratio of cement-treated QB. Their study concluded that stabilizing QBs with cement could produce the adequate compressive strength, modulus of elasticity, and tensile strength required for subbase material. They proposed base course material additive, flowable fill, under slab granular fill, and cement-stabilized subbase/base layers as possible pavement applications of QBs. According to the results presented in the study by Wood and Marek (1995), using 3% cement, 8% fly ash, and 89% QBs resulted in a flowable fill with adequate performance.

McClellan et al. (2002) reported engineering backfill as a potential area of QB material use, which was evaluated based on particle size distribution, moisture content, and mineralogy of by-products representing a variety of limestone and dolomitic QB. Naik et al. (2005) examined the use of QBs in Self-Consolidating Concrete (SCC). They found that the addition of QBs minimized the addition of the admixture without reducing the strength of the SCC. Puppala et al. (2008) reported that the addition of 2.3% cement increased the unconfined compressive strength to 174 psi (1,200 kPa). It was concluded that the strength and resilient modulus of the cement-treated QB are similar to those of sandy material with very few fines. Following the laboratory assessment of QBs, field performance tests were also conducted with QBs used as subbase/base material on expansive subgrade treated with lime (Puppala et al., 2012), where no substantial changes in the surface deformation profile were observed during the experimental testing. Lohani et al. (2012) found that the replacement of sand with QBs in concrete improved the properties of the mixture since QBs improved the pozzolanic reaction, micro aggregate filling, and concrete durability. Recently, Qamhia et al. (2018, 2019) utilized aggregate QBs in sustainable pavement applications as a filler material in the voids of large aggregate subgrade rocks and as a stabilized base/subbase material in flexible pavement, and reported satisfactory performance trends from accelerated pavement testing.

# PROBLEM STATEMENT AND OBJECTIVES

Large amount of non-hazardous QB materials are produced each year during various crushing stages and stockpiled over time in quarries. There exists potential for various construction application of QBs given having a good understanding of how these QBs behave alone and with proper stabilizers. Previous studies on QBs have provided many useful background information and insights. However, it is generally agreed that

TABLE 1 | Quarry by-products sources and samples.


TABLE 2 | Experimental plan.


the properties of QBs could not be easily generalized because of the natural variability of parent rock and the different crushing technologies employed (Wood, 1995; Manning, 2004; Stroup-Gardiner and Wattenberg-Komas, 2013).

The main goal of this research is to evaluate physical and engineering properties of wide range of QBs from different sources and different crushing stages in Illinois and identify most feasible paving applications according to laboratory testing results.

To address the research objectives, a thorough laboratory characterization was conducted; which included determining the engineering properties as well as the chemical and compositional characteristics of QBs received from four quarries in Illinois and sampled from multiple crushing stages. Both 2% Portland cement and 10% Class C fly ash were used as stabilizers to investigate trends in strength improvement. Chemical admixture stabilization improvements of representative QB samples were evaluated to provide recommendations for potential sustainable QB applications. Effect of density, gradation, as well as chemical composition on the unconfined compressive strength of stabilized QB was also evaluated.

# MATERIALS AND EXPERIMENTAL PLAN

For simplicity, the QBs used in this study, which were received from four different quarries in the State of Illinois will be referred to as Quarry 1–Quarry 4, or Q1–Q4 for short. Two batches of QB samples were collected from Quarry 1 near the Chicago area and were sampled within a 5-month period. The laboratory testing results for the two batches of materials from Quarry 1 were presented in a recent paper by Mwumvaneza et al. (2015). In this paper, only the average values of the two batches of Quarry 1 are presented. After the first round of the experimental program was completed for the samples collected from Quarry 1, additional QB samples were collected from three other quarries (see **Table 1**). All four quarries generate large quantities of QB annually and are located in different regions of Illinois, in order to evaluate variability in engineering properties of QB from different geological origins and production techniques to better represent QB materials from all across Illinois.

In general, the three crushing and screening stages utilized in quarries (primary, secondary, and tertiary) can produce QBs along with the main aggregate products meeting specifications of various applications. The materials were sampled from three main crushing/screening stages when they were active and contributing to the aggregate production. The materials collected from Quarries 1 and 2 included all three stages however, samples were only collected from the primary and tertiary stages from Quarries 3 and 4 because only two crushing stages were utilized on site at these two locations.

The scope of the experimental program included sieve analysis, moisture density relationship, imaging-based aggregate shape property testing, Modified Methylene Blue (MMB), Atterberg limits, and Unconfined Compressive Strength (UCS) characteristics of QB samples collected from the crushing stages shown previously in **Table 1**. In addition, X-Ray Fluorescence (XRF) technique was also utilized for evaluating chemical composition of QB samples. Both Portland cement and Class C fly ash treated QB samples were also evaluated by standard Proctor test and UCS test. A summary of the experimental framework is illustrated in **Table 2**.

# RESULTS AND DISCUSSION

### Untreated Quarry By-Products Grain Size Distribution, Atterberg Limits, and Soil Classification

Owing to the variety of parent rock types, and different crushers and technologies used in the aggregate production stages, the particle size of aggregate by-products may differ from one another. To quantify these differences, grain size distributions were determined by laboratory dry sieve analysis according to the ASTM C136 method (ASTM International, 2014). Quarry by-products were compared based on the quarry locations and crushing stages. **Figure 1** shows gradation results for all the QB samples obtained from the four quarries also indicating the crushing stages. The aggregate by-products from Quarries 1–3 have a nominal maximum size of 4.75 mm (No. 4 sieve size). About 2% of the aggregate by-products sampled from Quarry 4 tertiary crusher were found to be slightly larger than 4.75 mm.

The percentages of particles finer than 0.075 mm (No. 200 sieve) range from 7 to 15% for all quarry samples. Note that all the aggregate by-products follow a well-graded type of gradation curve except for the QB sample obtained from the tertiary production stage of Quarry 3.

Clay contents of the aggregate by-products were determined from a Modified Methylene Blue (MMB) test, which is a quick assessment method of the amount of harmful clay in the fines portion of the aggregates (Pitre, 2012). For the two batches collected from Quarry 1, almost similar clay contents were observed, and an average was taken to represent the clay content for those samples. Harmful clay contents for samples from Quarries 1 through 4 are summarized in **Table 3**.

In general, harmful clay content comprised < ∼3% of fines in the QB samples. There were some variations between quarries and within quarries. For example, the harmful clay contents

Hou et al. Performance of Aggregate Quarry By-Products



\**Averages of the two batches.*

\*\**Unified soil classification system of the two batches.*

of Quarry 1 samples were lower than those of the samples collected from the other three quarries. Slightly higher clay content was observed in QBs from the primary crushing stage except for Quarry 3, where higher clay contents were observed from samples collected during the tertiary crushing stage. The higher clay content in QBs from the primary crushing stage is possibly the result of weathered rock or overburden material, which are more likely to be found in the quarried rocks in the primary crusher.

Atterberg limit tests were performed in accordance with the ASTM D4318 method (ASTM International, 2017a). All QB samples from the four quarries and from different crushing stages were found to be non-plastic. This indicates that QB samples had a low moisture affinity and a low shrink-swell potential.

Based on the American Association of State Highway and Transportation Officials (AASHTO) Soil Classification System, all the QB samples were classified as A-2-4, which represented silty gravel and sand. Based on the Unified Soil Classification System (USCS), the QB samples were mostly well graded sands (SW) or silty sands as shown in **Table 3**. Classification of silty sand (SM) in general indicates that the percentage of fine materials passing sieve No. 200 was often >12%. Typically, the primary crusher QB samples from Quarries 1–3 produced less fines compared with the secondary and tertiary crusher materials.

#### Morphological Shape Properties

Aggregate particle shape and angularity were studied using an aggregate image analysis system, known as the Enhanced University of Illinois Aggregate Image Analyzer (E-UIAIA). The E-UIAIA is equipped with three high-resolution (1292 × 964 pixels) charge-coupled device progressive-scan color cameras to capture three orthogonal views (front, top, and side) of individual particles to establish the morphological indices of aggregate particle shape, texture, and angularity. More details on the features of the E-UIAIA can be found elsewhere (Moaveni et al., 2013). **Figure 2** shows the E-UIAIA and an example of the three orthogonal views of a particle captured by the camera system.

The Flat and Elongated (F&E) ratio and Angularity Index (AI) were the key image-based indices—measured with the E-UIAIA—for determining the particle shape properties of the different QBs from the different quarries and crushing stages. Approximately 50 particles retained on the No. 8 (2.38 mm) sieve were scanned for QBs from each crushing stage from each quarry to determine trends in particle shape. Average AI values and average F&E ratios for QB materials are presented in **Table 4**. QB samples from the primary crushers generally tended to have lower AI values indicating the least angular particles, except for Quarry 3. The samples from the secondary crusher had slightly higher value than the samples from the tertiary crusher. These findings supported the visual assessment that primary QB samples appeared to have more round particles than the secondary and tertiary QB samples. Another trend observed was that QB from the primary crusher had the lowest F&E ratio and that QB from the secondary and tertiary crushers had very close F&E ratios. These findings indicated that particles from the primary crusher were more cubical and therefore may have better resistance to breakage. Another observation is that QB samples from Quarry 1 tended to be more angular when compared with samples from the other three quarries. Detailed distributions of shape indices for all the particles scanned by the E-UIAIA are presented in **Figure 3**. Like the grain size distribution concept, **Figure 3** presents the cumulative percent of QB vs. certain shape index from different crushing stages and quarries.

#### Standard Proctor Test

In accordance with ASTM D698 (ASTM International, 2012), the moisture-density compaction characteristics of the virgin QB samples from the four quarries sampled in this study were evaluated; the results are summarized in **Table 5**. For all QBs from different quarry locations, the Optimum Moisture Contents (OMC) from standard Proctor tests ranged from 7.9 to 10.4%, whereas the Maximum Dry Densities (MDD) ranged from 129.7 to 142.1 pcf (20.4–22.3 kN/m<sup>3</sup> ). In general, the MDD and OMC pairs varied from one quarry location to another. While there was no unique trend for the variation of OMC, for all quarry locations, higher maximum dry density values were observed in QBs from the primary crushing stages.

#### Unconfined Compressive Strength Test

Unconfined Compressive Strength (UCS) tests were performed on the QB samples from all four quarries. The samples prepared for UCS testing were compacted at the OMC and MDD obtained from the moisture-density characteristics of virgin QB materials. The UCS specimens were compacted in a split mold into cylindrical specimens in three equal lifts using the standard Proctor hammer for compaction. The final specimens were 2.8 in. (71 mm) in diameter by 5.6 in. (142 mm) high. All samples were then cured unsealed in a moisture-controlled room at 100% relative humidity and at room temperature for 7 days. Samples were prepared and tested in accordance with the ASTM D2166 method (ASTM International, 2016). For each QB material obtained from a certain crushing stage, two samples were prepared to conduct UCS tests.

FIGURE 2 | E-UIAIA and the particle front, top, and side views captured (Moaveni et al., 2013).



In **Figure 4**, UCS values are presented with the maximum densities obtained from the standard Proctor tests. Note that all QB samples generally exhibited low UCS values. The average UCS values for the QB materials from Quarry 1 decreased in an orderly fashion in different crushing stages. For samples from Quarry 2, the UCS values decreased from the primary crusher QB to the secondary crusher QB but then increased for the tertiary crusher QB. For samples from Quarry 3, the UCS values also decreased from the primary crusher QB to the tertiary crusher QB. For samples from Quarry 4, the UCS values for both the primary and tertiary crushing stages were the same with an average value of 8.1 psi (55.8 kPa).

In addition, a clear trend can be seen between UCS and the maximum dry densities obtained: for samples from the same quarry source, the higher the maximum dry density, the larger the UCS was achieved. This emphasizes the importance of particle size distribution and packing on performance; and therefore, good compaction in the field is required if QB samples are to be used in future pavement applications. Overall however, the compressive strength values of the QB materials were low, thus

indicating the need for improvement through stabilization depending on the strength requirements for various highway construction applications.


Virgin materials


\**Average values of the two batches.*

# Stabilization of Quarry By-Products

To increase the strength properties of the QB samples, chemical stabilizers were chosen. While chemical stabilization of soil and aggregates improves their physical and engineering properties, this process depends heavily on the chemical reaction between soil/aggregates and the stabilizers. It is therefore important to choose the right stabilizers to effectively improve the strength. The rationale behind the selection of stabilizers is to obtain maximum strength gain with lowest cost and environmental impact.

According to the previous research study on similar materials, stabilization using low-cement contents appeared to be a costeffective alternative successfully used in stabilizing quarry byproducts. In addition, lower cement content has the benefit of reducing shrinkage potential of cemented materials. Therefore, two percent (2%) Type I Portland cement by the oven-dry weight of the aggregate by-products was used as the first stabilizer in the laboratory. The second choice was fly ash, a by-product of a different industry. Like cement, Class C fly ash possesses cementitious and pozzolanic properties that do not depend on the reaction with clay particles to develop strength. Initial trials of 5 and 10% of Class C fly ash were conducted to determine the proper stabilizer amount and effectiveness. The results of initial trials were presented in a previous publication (Tutumluer et al., 2015). Based on initial trials, the percentage of Class C fly ash used in this study was 10% by the oven-dry weight of the aggregate by-products. The Class C fly ash material conformed to ASTM C618 (ASTM International, 2019b) and AASHTO M295 (ASTM International, 2019a) standards.

### Standard Proctor Tests on Stabilized Quarry By-Products

The compaction properties of the cement and Class C fly ash stabilized QB materials were determined as per the ASTM D558 method (ASTM International, 2011). The MDDs and OMCs determined are listed in **Table 6**. Note that for the 2% cement-treated quarry by-products the OMCs ranged from 6.6 to 9.9% and for the 10% Class C fly ash–treated aggregate byproducts, the OMCs ranged from 7 to 8.3%. No clear trend could be established on the effect of quarry location and crushing stage on the MDD and OMC values.

In general, smaller values of OMC and larger values of maximum dry density were observed for the 10% Class C fly ash-treated QB materials when compared to the 2% cement-treated QBs. Such a difference in the characteristics of aggregate by-products stabilized with Portland cement and Class C fly ash was attributed to the differing amounts of the free lime that each stabilizer contributed during the flocculation and agglomeration phase of the QB treatment.

## Unconfined Compressive Tests on Stabilized Quarry By-Products

Unconfined Compressive Strength (UCS) tests were also conducted on the stabilized QB materials to evaluate the benefits of chemical admixture treatment and demonstrate how treated samples of weak soils could be improved to achieve the desired strength.

The maximum dry density and the OMC data obtained from the moisture-density characteristics of treated quarry byproducts were used to prepare samples for the UCS tests. Specimens were prepared as per ASTM D1632 and ASTM D1633 for QB materials treated with each of the two stabilizers (ASTM International, 2017b,c). Significant strength increases were reported for the Quarry 1 QB specimens treated with the two stabilizers in an earlier published study (Mwumvaneza et al., 2015). In this paper, new results for QB materials from Quarries 2–4, together with old results for QB materials from Quarry 1 are presented. For these samples, two specimens from each crushing stage were prepared and tested. All the samples treated with admixtures were cured in a moisture-controlled environment for 7 days at room temperature at 100% humidity. Before UCS testing, all the stabilized samples were soaked for 4 h to evaluate the effects of a harsh, moist environment on strength properties.

**Table 7** lists the average UCS values obtained for all of the QB samples collected from the four quarries at different crushing stages. Significant strength improvements could be achieved with proper admixture treatment. Note that virgin (untreated) aggregate by-products had very low UCS values, with an average of <11 psi (76 kPa). When 10% Class C fly ash was used, strength gains were ∼20–30 times of the original untreated samples. The highest strength values were achieved with 10% fly ash in Quarry 1 samples (more than 300 psi/2,068 kPa), whereas other quarry samples had relatively lower strength values achieved of around 100–300 psi (689– 2,068 kPa). On average, the 2% cement-treated QB specimens from all four quarries and crushing stages had consistently high UCS values above 200 psi (1,379 kPa). Virgin QB samples were strengthened by more than 20 times the initial strength by adding 2% cement. Approximately similar gains were achieved by 2% cement and 10% Class C fly ash treatment of QB samples.

## Factors Affecting Unconfined Compressive Strength of Stabilized Quarry By-Products

#### **Dry density achieved**

Despite the fact that both 2% cement-treated and 10% Class C fly ash–treated QB materials can reach significantly high strength values compared with virgin materials; differences were observed in the stabilized samples. **Figure 5** shows that the fly ash treatment resulted in the highest strength gains for all the QB samples from Quarry 1. However, for all samples from the other quarries, fly ash treatment was less effective than the cement treatment. One possible reason is that materials with 10% Class C fly ash treatment from Quarry 1 achieved significantly higher maximum dry densities when compared to the 2% cement treatment, while the stabilized QB materials from the other quarries did not show such large differences in maximum dry density (see **Figure 5**). Higher density indicates better compaction/packing and higher strengths, which was certainly the case for the rapid nature of tricalcium aluminate (C3A) reaction taking place in this self-cementing 10% Class C fly ash.

Accordingly, one possible reason that would affect the strength properties of stabilized QB samples can be the compaction and density achieved. To take a closer look at the UCS results, strength properties of the QB samples from the same quarry can be related to the densities achieved (six out of the 10 materials). For example, materials with 10% Class C fly ash treatment from Quarry 1 achieved higher maximum dry density than that achieved by the 2% cement treatment, which resulted in higher compressive strength of Class C fly ash treated materials. Similarly, materials from Quarry 2 primary and tertiary crushers treated with 2% cement achieved higher maximum dry densities and higher compressive strengths when compared to those treated with 10% Class C fly ash. Also, Quarry 3 tertiary crusher materials treated with 2% cement and 10% Class C fly ash achieved approximately the same maximum dry density and compressive strength characteristics. However, for the other four materials tested, such as those from secondary crusher from Quarry 2, primary crusher from Quarry 3, and primary and tertiary crushers from Quarry 4, the achieved maximum dry density trends did not correlate well with the compressive strength characteristics, which indicates that unlike virgin QB samples, whose strength mainly depends on density achieved, there are other governing factors affecting the strength properties of stabilized QB.

#### **Gradation and packing**

Another possible factor contributing to the strength of stabilized QB can be the gradation trend and its effect on particle packing. To address the gradation and particle packing concerns, laboratory tests were conducted to determine the effect of QB grain size distribution on the unconfined compressive strength of cement and Class C fly ash stabilized QB. The QBs from Quarry 3 were size separated and engineered to match certain gradation curves that are power exponents of the ratio of the sieve size to the maximum particle size according to the Fuller curve, also commonly known as the Talbot equation: p<sup>i</sup> = Di <sup>D</sup>max <sup>n</sup> , where pi is the percentage passing the ith sieve, D<sup>i</sup> is the ith sieve size in mm, Dmax is the maximum sieve size and is taken to be the No. 4 sieve for sampled QB, and "n" is power exponent of the curve. The "n" values considered in this study were 0.3, 0.4, 0.45, 0.5, and 0.6, to obtain five different gradation curves as shown in **Figure 6A**.

In order to compare the gradation and packing effect, densities and moisture contents for each engineered gradation were set


Tertiary 8.0 21.1 (134.6) 3 Primary 8.2 21.0 (133.9) Tertiary 7.8 20.9 (133.1) 4 Primary 7.2 21.8 (139.0) Tertiary 7.8 21.3 (135.6)

TABLE 6 | Optimum moisture contents and maximum dry densities of stabilized quarry by-products sampled from quarries 1–4.

\**Moisture-density test is conducted on first batch.*

the same as the maximum density and optimum water content for the original gradation obtained from standard Proctor test. Densities achieved were 132.1 pcf (20.7 kN/m<sup>3</sup> ) and 133.5 pcf (21.0 kN/m<sup>3</sup> ) for cement and fly ash stabilized QB samples, respectively. Water contents targeted were 9.2% and 9.4% for cement- and fly ash-stabilized QB samples, respectively. Three replicates were tested for each engineered gradation. More details about this study can be found elsewhere (Qamhia et al., 2016). Accordingly, stabilized QB samples from Quarry 3 had varying strength properties at different engineered gradations. The engineered gradation with the power term "n" equals to 0.45 led to the highest UCS for both the 2% cement and 10% Class C fly ash treatment (see **Figure 6B**).

The virgin QB gradations from quarries 1–4 were then checked with the engineered gradations for the "n" power values of 0.3, 0.4, 0.45, 0.5, and 0.6 to explain some of the behaviors of stabilized QB strength in this study. As stated previously, four out of the 10 materials, i.e., the secondary crusher from Quarry 2, primary crusher from Quarry 3, and the primary and tertiary crushers from Quarry 4, did not exhibit strong correlations between density and strength. For the QB material from the Quarry 3 primary crusher, this can be explained since the QB TABLE 7 | Average UCS for virgin QB materials and stabilized QB materials with 2% cement and 10% class C fly ash.


gradation curve is close to n = 0.6 power curve. Therefore, possible explanations for a lower strength of the fly ash treated materials are insufficient packing and relatively larger particle sizes diminishing the effectiveness of the fly ash reaction with QB particles. In addition, the secondary crusher QB material from Quarry 2 had higher strength compared with the primary crusher and tertiary crusher QB materials even though the secondary crusher QB material had the lowest density achieved. One possible reason for this is due to a better packing of the secondary crusher QB materials from Quarry 2 since its gradation is very close to the n = 0.45 maximum density gradation curve.

#### **Effectiveness of stabilizer**

Effectiveness of the stabilizer can be another factor that affects the strength of stabilized QBs. Note that QB materials from Quarry 1 were tested and characterized for engineering properties first, and then, the rest of the QB materials from the other 3 quarries were tested 10 months later. Class C fly ash has a shelf life that depends on storage conditions, and the effectiveness of fly ash is reduced with time as it hydrates in damp conditions. As a result, the free lime in Class C fly ash decreases over time and the cement hydration process for QB samples from quarries 2– 4 is less effective than for QB samples from quarry 1. This may explain the reduction in the strength gains with three of the QB materials tested afterwards. Thus, it is expected that fly ash was highly reactive when QB materials from Quarry 1 were tested and the strength gains with 10% Class C fly ash were possibly higher than those for the 2% cement treated samples as well as the other QB materials stabilized with 10% Class C fly ash.

#### **Chemical composition of QB**

Another factor that controls the effectiveness of stabilization is the chemical composition of QB. The cementitious and pozzolanic reactions are the major causes of strength increase, which depends highly on the properties of soils and stabilizers. In this study, X-Ray Fluorescence (XRF) technique was utilized to determine the chemical composition of the tested QB materials. XRF is a non-destructive analytical technique to determine

the elemental composition of QB materials. The chemical composition results are shown in **Table 8**. Note that components <1% by composition of QB samples are omitted here.

As seen in **Table 8**, QB samples from the same quarry usually have similar chemical composition since they are composed of the same mineral substances as the parent aggregates from which they are derived. To be more specific, CaO is the major component for QBs from all four quarries. When comparing the percentage of MgO and SiO2, it can be seen that QBs from Quarry 2 and Quarry 3 primary crusher have higher amount of SiO2; QBs from Quarry 1 have higher amount of MgO; while QBs from Quarry 3 have least amount of MgO. However, there still might be substantial differences in composition among different crushing stages due to production processes and weathering condition due to storage of QBs. QBs from primary and tertiary crushers from Quarry 3 are quite different in composition. For example, QB from Quarry 3 primary crusher has significantly lower CaO content, significantly higher SiO<sup>2</sup> content compared with QB from tertiary crusher, and also the highest amount of Al2O<sup>3</sup> among all samples. This might explain why QBs from Quarry 3 primary and tertiary crushers exhibit such different strength behaviors. Unlike all the other samples from Quarries 1, 2, 4, where primary crusher has more strength gain than tertiary crusher, QBs from Quarry 3 primary crushers gains less strength during stabilization. This is probably due to higher CaO amount in QB from tertiary crusher, which accelerates the cement hydration process in the 7-day testing period. Similarly, QBs from Quarry 4 are rich in CaO, which exhibits the largest strength gain during cement stabilization. In terms of long-term strength gain, it is expected that QBs consisting higher amount of SiO<sup>2</sup> and Al2O<sup>3</sup> tend to gain more strength with pozzolanic reaction.

In conclusion, on the basis of the UCS measurements, both 2% cement and 10% Class C fly ash were shown to be effective stabilizers that can be used with the aggregate by-products for strength improvement. For the virgin QB samples, it appears that density is the main factor governing the strength of QB such that high densities correlate to high UCS values for all four quarry materials. On the other hand, for stabilized QB materials, several other factors may control strength, including density, gradation and packing, particle shape and angularity, chemical composition of QB samples as well as the effectiveness of stabilizer.

# SUMMARY AND CONCLUSIONS

This paper presents findings of the Illinois Center for Transportation study on Quarry By-products (QBs), the residual deposits from the production of required grades of aggregate. The objectives of this study were to determine the material properties of aggregate by-products through laboratory testing and materials characterization, and to evaluate the feasibility of using QBs in various pavement applications. A full matrix of laboratory tests i.e., modified methylene blue, dry sieve analysis, moisture density, Atterberg limits, shape properties, XRF, and unconfined compressive strength were performed to explore the feasibility of using QB in pavement applications as an area that can consume higher amounts of QB. The performed tests also helped to evaluate and distinguish QB samples collected from the primary, secondary, and tertiary crushing stages. The overall strength properties of QBs from all the four quarries were found to be low. Accordingly, chemical admixture stabilization alternatives using Portland cement and Class C fly ash were studied as means to improve the strength and durability characteristics of QBs. For that evaluation, moisture density and UCS tests were conducted on treated QB samples. The addition of cement and fly ash considerably increased the strength properties. Several factors that affect the strength

property of stabilized QB samples were discussed and evaluated, including density, gradation and packing, chemical composition, and the effectiveness of stabilizers.

Based on the research findings, the following conclusions can be drawn from the study:


TABLE 8 | Chemical composition of quarry by-products sampled from quarries 1–4.


6. The XRF test on QB showed that chemical composition is another main factor contributing to the strength of the stabilized QB. Further investigation of detailed correlation between chemical composition and strength is needed.

The use of QB materials with fly ash stabilization can provide opportunities to develop sustainable pavement construction strategies. Fly ash is a by-product of coal-burning plants, and using Class C fly ash as an admixture for treating aggregate by-products can be a more cost-effective option for sustainable pavement applications. While the 7-day gain in UCS values is promising, further investigation of the chemical properties of stabilized QB mixes is needed. Recommendations for future work include investigating the 28-day compressive strength properties for stabilized QB mixes, testing more QB sources with different parent materials and geological origins, and testing the dynamic properties of these QB materials with repeated load triaxial testing. Further, field validation is recommended for the full assessment of the performance of the treated QB materials in sustainable pavement applications.

# DATA AVAILABILITY STATEMENT

The datasets generated for this study are available on request to the corresponding author.

# AUTHOR CONTRIBUTIONS

The authors confirm contribution to the paper as follows: WH, IQ, VM, ET, and HO: study conception and design, analysis and interpretation of results. WH, IQ, and VM: data collection. WH and IQ: draft manuscript preparation. All authors reviewed the results and approved the final version of the manuscript.

## ACKNOWLEDGMENTS

This publication is based on the results of ICT-R27-125 project which focused on characterizing quarry by-products in pavement applications. ICT-R27-125 was conducted in cooperation with the Illinois Center for Transportation; the Illinois Department of Transportation, Office of Program Development; and the U.S. Department of Transportation, Federal Highway Administration. The final report of the ICT R27-125 project, entitled, Sustainable Aggregates Production: Green Applications for Aggregate By-Products, can be accessed via https://apps.ict.illinois.edu/projects/getfile.asp?id=3507. The authors would like to acknowledge Pengcheng Wang and Joshua S. Cheung, graduate students at the University of Illinois at Urbana-Champaign, for their help with the laboratory work. Cooperation and support of the Illinois Association of Aggregate Producers is greatly appreciated. The contents of this paper reflect the views of the authors who are responsible for the facts and the accuracy of the data presented herein. This paper does not constitute a standard, specification, or regulation.

# REFERENCES


**Conflict of Interest:** VM is employed by Dynatest Nort America. However, VM was a graduate student at UIUC when he was involved with the research work described in this paper.

The remaining authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2019 Hou, Qamhia, Mwumvaneza, Tutumluer and Ozer. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Discrete Element Modeling of Permanent Deformation Accumulation in Railroad Ballast Considering Particle Breakage

Beema Dahal 1† and Debakanta Mishra<sup>2</sup> \*

*<sup>1</sup> Department of Civil Engineering, Boise State University, Boise, ID, United States, <sup>2</sup> School of Civil and Environmental Engineering, Oklahoma State University, Stillwater, OK, United States*

#### Edited by:

*Sanjay Shrawan Nimbalkar, University of Technology Sydney, Australia*

#### Reviewed by:

*Yang Gui, Hohai University, China Qideng Peter Sun, Coffey International Australia, Australia*

> \*Correspondence: *Debakanta Mishra deb.mishra@okstate.edu*

> > †Present address:

*Beema Dahal, Bailey Engineering Inc., Eagle, ID, United States*

#### Specialty section:

*This article was submitted to Transportation and Transit Systems, a section of the journal Frontiers in Built Environment*

> Received: *16 October 2019* Accepted: *09 December 2019* Published: *14 January 2020*

#### Citation:

*Dahal B and Mishra D (2020) Discrete Element Modeling of Permanent Deformation Accumulation in Railroad Ballast Considering Particle Breakage. Front. Built Environ. 5:145. doi: 10.3389/fbuil.2019.00145* Particle shape and mineralogy have been found to greatly affect railroad ballast response under train loading. In this research effort, the effect of particle breakage on ballast permanent deformation behavior was studied using the Discrete Element Method (DEM). Real ballast particle shapes were digitized using an inexpensive imaging tool, and agglomerates of spheres were used to regenerate those complex shapes. Particle crushing tests were carried out in the laboratory, and discrete element models simulating the crushing test were calibrated based on the laboratory test results. Model parameters established through this calibration process were subsequently used to study the permanent deformation response of ballast layers comprising complex-shaped breakable ballast particles under repeated loading. The effect of particle shape was studied by comparing the permanent deformation behavior for ballast layers comprising simplified ellipsoid particles to those comprising the complex-shaped particles simulated using agglomerates of spheres. The simulation results clearly highlighted the significance of particle breakage in ballast permanent deformation accumulation under cyclic loading. Most particle breakage was observed to occur under the initial load cycles, which also corresponded to the highest rate of permanent strain accumulation. The particle breakage as well as the rate of strain accumulation gradually decreased with increasing load cycles.

Keywords: discrete element method, railroad ballast, ballast shape, ballast breakage, permanent deformation, particle crushing test

# INTRODUCTION

As a part of the railroad track substructure, the ballast layer performs various functions like resisting the vertical, longitudinal, and transverse forces from trains; distributing the high stresses to protect underlying track layers; absorbing the shock from dynamic loads from moving trains, and facilitating the free drainage conditions (Nimbalkar and Indraratna, 2015). As the ballast layer ages, it is progressively fouled with materials that are finer than the initial ballast aggregates (Qian et al., 2014). The main cause of ballast fouling is attributed to ballast degradation and breakdown of the uniformly-graded ballast aggregates under repeated train loading (Selig et al., 1988, 1992). It has been reported that ∼76% of the ballast layer fouling occurs due to ballast breakage (Selig and Waters, 1994). Particle breakage significantly influences the shear strength and deformation behavior of railroad ballast layers (Indraratna et al., 1998; Anderson and Key, 2000), resulting in unacceptable differential settlement of track and pumping of underlying soft subgrade soils which eventually necessitate frequent costly track maintenance (Indraratna and Nimbalkar, 2011). Accordingly, to study the influence of particle breakage on the strength and deformation characteristics of granular soils, numerous experimental (Shenton, 1975; Lackenby et al., 2007; Aursudkij et al., 2009; Indraratna et al., 2010b), and numerical-discrete element (DE) (Lim and McDowell, 2005; LoboGuerrero and Vallejo, 2006; Lu and McDowell, 2010) studies have been conducted. Numerical studies have been frequently used for the ballast breakage studies as they give a complete picture of where the degradation starts, how it evolves, and how it affects the permanent deformation of railroad tracks (LoboGuerrero and Vallejo, 2006).

Most three-dimensional DE codes use spheres to represent particle shapes due to their ease in the detection of inter-particle contact and force calculation. Models to describe the physics of contacts between spheres with walls, and with each other are relatively well-known (Di Maio and Renzo, 2004; Kruggel-Emden et al., 2007). Spheres can be used to approximate the qualitative behavior of the granular materials, but cannot give the information on the quantitative measures and therefore show deviant mechanical behavior than the actual material in consideration (Latham and Munjiza, 2004; Kruggel-Emden et al., 2009; Ferellec and McDowell, 2010). Particle shape, size, and minerology have been found to have significant impacts on the mechanical behavior of the granular materials (Azéma et al., 2013; Höhner et al., 2013). However, representation of the particle shape has been one of the main challenges of DE simulations (Markauskas et al., 2010). There are Discrete Element Method (DEM) tools such as like BLOKS3D that can simulate actual interaction between the polyhedral shapes (Zhao et al., 2006) but are not able to simulate breakages. Hence, using PFC3D this research study uses the concept of agglomerates where spheres are bonded together to represent the complex and irregular polyhedral particle shapes as they can make the contact detection and force calculation easy. Similar approaches have been used by (McDowell and Harireche, 2002; Thornton and Liu, 2004; Alshibli and Cil, 2014) to model the particle breakage. It should be noted that in the remainder of this manuscript, polyhedral and ellipsoid shapes were actually approximated by using agglomerates of bonded spheres. To understand the significance of ballast breakage considerations on permanent deformation behavior of railroad ballast layer under repeated loading, the authors initially conducted a study using simple ellipsoid shapes as representative ballast shapes and reported the findings in an earlier manuscript (Dahal et al., 2018). However, since the ballast particles in Dahal et al. (2018) were simulated as ellipsoids, the results were arguably different from real-world scenarios where the ballasts are polyhedral in nature. Besides, the bond strength values (primarily governing the particle breakage) assigned to the ballast agglomerates were modified from the values obtained from the literature (Indraratna et al., 2010b). Accordingly, subsequent research tasks focused on modeling railroad ballast behavior under repeated loading by considering realistic particle shapes as well as breakage potential. This manuscript presents major findings from these efforts.

# RESEARCH OBJECTIVES

The primary objective of this research effort was to study the importance of considering realistic particle shapes during discrete element modeling while studying the permanent deformation and breakage behavior of railroad ballast under repeated loading. To achieve this objective, one of the tasks involved conducting laboratory tests to quantify the crushing strengths for actual ballast particles, and calibrating DE models of the crushing strength test to establish relevant model parameters to be used in subsequent modeling tasks. The final objective of this manuscript was to compare the results obtained for complexshaped ballast particles (simulated using the agglomerates of spheres) to those obtained for simplified ellipsoid shapes; this would help identify the importance of considering realistic particle shapes during discrete element modeling.

# ACQUISITION OF REAL BALLAST SHAPES

Several researchers have used expensive and complex image analysis approaches to create polyhedral particles for use in DE simulations (Paixão et al., 2001; Le Pen et al., 2013; Qian et al., 2013; Sun et al., 2014). However, such expensive imaging equipment to obtain the scanned 3D image of the ballast shapes are not readily accessible. In such instances, exploration of alternative inexpensive alternatives to digitize realistic ballast particle shapes becomes imperative. With the advent of modern smartphones, it has become increasingly easier to capture highresolution images that can be easily digitized. This research effort utilized such an approach to digitize actual ballast particle shapes that were subsequently imported into the DE model.

Three different ballast shapes were selected to demonstrate this approach. To capture the images, each ballast was placed on a flat surface such that the flattest part of the particle rested on the platform. More than twenty (20) sequential images of the ballast particles were taken by rotating the particles from 0◦ to 360◦ on a horizontal plane, and the images were then imported into the Autodesk <sup>R</sup> Recap Software (Autodesk Inc., 2018). The flat surface was closed using a built-in feature of the software, and was finally exported in ".stl" format to be used in PFC3D (DE software). **Figures 1A–F** show the chosen ballast shapes and exported images of ballasts from the software for Ballast 1 (B1), Ballast 2 (B2), and Ballast 3 (B3), respectively. These shapes were approximated within PFC3D using the agglomerates of spheres to be used in the subsequent simulations.

# LABORATORY SPECIFIC GRAVITY TEST

The specific gravity test was conducted to accurately simulate the ballast particles in the DE model. Specific gravity of the selected ballasts was calculated using the following procedure. The test apparatus was set up and the temperature of water was

set to 25◦C. Each ballast particle was taken separately and was weighed to calculate its dry weight (W1). The particle was then put in water and weighed (W2). Finally, the particle was removed from water, surface dried with a towel by rolling inside a towel for ∼5 s, and weighed to obtain the surface dry weight (W3). The specific unit weight or specific gravity was calculated using the equation below.

$$\text{Specific Gravity} = \frac{\text{W}\_1}{\text{W}\_3 - \text{W}\_2}$$

The specific gravity values for B1, B2, and B3 were obtained to be 2,741, 2,555, 2,534 kg/m<sup>3</sup> , respectively.

# SINGLE PARTICLE CRUSHING TEST

Several different factors such as angularity, uniformity of gradation, lower particle strength, coarseness, stress level, and anisotropy govern the crushing of railroad ballast particles (Bohac et al., 2001). The most important of all these factors is the crushing strength of ballast. Appropriate particle crushing strength values can be established by performing the Single Particle Crushing Test (SPCT) in the laboratory. The SPCT is an indirect tensile test that is conducted by compressing individual particles between two flat platens to induce tensile stresses in the particles (Lim, 2004). A typical result of this crushing test is a plot of force against deformation. The maximum peak load is the point at which major fracture occurs along the loading direction whereby the particle splits into two or more pieces.

# Laboratory Crushing Test

A Universal Testing Machine (UTM) was used to perform the laboratory tests on the three ballast particles. Since the shapes of the particles were irregular and coarse, they were aligned on the bottom platen in such a way that each particle rested with its flattest portion on bottom platen of the machine. **Figures 2A,B** show the test set up for SPCT and B1 particle in the test setup for crushing, respectively. Load was then applied by moving the top platen downwards at 1.9 mm/s. Force and displacement data were continuously collected throughout the test. **Figure 2C** shows the force-displacement plot obtained for particles B1, B2, and B3. From the figure, the load-displacement plot went through several undulations owing to shifting and reorientation of the ballast particle during the loading process. These displacements were not included in the final force-displacement calculations. Similar results were obtained for ballast particles B2 and B3. **Figures 2D–F** show photographs of the ballast particles after they have been crushed in the laboratory. After the crushing strength for each ballast particle was established, the next step involved DE modeling and calibration of this crushing test.

# Discrete Element Simulation of the Particle Crushing Test

The material modeling support package of PFC was used to simulate the Single Particle Crushing Test (SPCT) of individual ballast particles (Itasca Consulting Group, Inc., 2018). This support package can be used to create an assembly of packed particles under a specified material pressure and porosity, and assign the contact model of one's choice, like the linear contact model, bonded particle model, etc. On these packed assemblies, different tests such as the tension test, compression test, and diametral compression test can be simulated. Detailed information about the material modeling support package can be found in the Itasca technical memorandum by Potyondy D. (2018). For simulating the SPCT on irregular ballasts shapes, the material modeling support package could not be used directly. Accordingly, some modifications were implemented, and new algorithms were developed to ensure the simulations accurately represented the particles tested in the laboratory.

#### Creation of Irregular Ballast Shapes

As already mentioned, individual particles in PFC3D (version 5) are represented as spheres. It should be noted that PFC3D version 6 does have the ability to model truly polyhedral-shaped particles; however, the research work reported in the current manuscript were carried out using PFC3D version 5. Therefore, to model non-spherical particles, special care needs to be taken to combine spheres of different sizes in certain fashions so that the resulting "agglomerate" represents the irregular shape being targeted. Once the three different ballast particles were digitized from the photographs, and their corresponding ".stl" files were imported into PFC3D, the next step involved recreating these particle shapes using agglomerates of spheres. This was accomplished through the following steps: First, cuboid boxes of dimensions larger than corresponding ballast particles were created, and filled with spheres of diameter 5–6 mm, and compacted to achieve a porosity value of 20% and a material pressure of 150 kPa. Note that to check whether the packing pressure affects the ultimate strength of the specimen, the pressure was reduced from 150 to 5 kPa; no change in the crushing strength was observed. Once the packed material was created, the spheres were bonded with linear parallel bonds with certain bond strength and modulus values. The detailed description of the linear parallel bond can be found elsewhere (Potyondy and Cundall, 2004; Potyondy D. O., 2018). The elastic modulus of the spheres was chosen such that the stiffness of the spheres matched the stiffness of respective ballast particles. The moduli and ratios of normal to shear stiffness values for spheres as well as cement (bond) of linear parallel bond were set equal to reduce the number of free parameters (Potyondy and Cundall, 2004; Potyondy D. O., 2018). The radius multiplier term, which refers to the radius of the parallel bond such that the bond radius equals the radius multiplier times the radius of the smaller sphere in contact, was set to one. Once the packed assembly was created, the geometry of the corresponding ballasts (see **Figure 1**) was imported into the model. An algorithm was developed to detect all spheres that had centroids within the geometry of corresponding shape being generated. All particles with centroids outside the imported geometry were deleted. **Figure 3** shows clusters of bonded spheres that were within the geometry of corresponding particle being generated. This cluster of bonded spheres now represented the irregular ballast shapes B1, B2, and B3 which were composed of 650, 216, and 254 individual spheres, respectively. Note that due to its complex shape, particle B1 required significantly higher number of spheres to generate a representative shape.

## Simulating the Diametral Compression (DC) Test

To subject the individual ballast particles to diametral compression test, the top and bottom walls or platens were created ensuring that they were placed exactly at the top- and bottom-most points of the representative particles. Each of the ballast particles were oriented such that they resembled their orientation in the lab test. Tilting of the loading platens as well as that of the ballast particle were not allowed during the simulated crushing tests. The system was brought to equilibrium, displacements of the spheres were set to zero, and the loading was applied by moving the top and bottom walls at a rate of 0.2 mm/s. Note that the simulated loading rate was intentionally set to a value significantly lower than that in the laboratory; this was necessary to ensure quasistatic conditions throughout the simulation. The frictionless top and bottom plates were moved toward the bonded assembly representation of ballasts to simulate the crushing. Microstructural monitoring to visualize the fragments caused during the test was performed whereby small red discs were created upon bond breakage at the breakage points. The test stopped when the load carrying ability of the ballast was compromised, and the axial force magnitude fell below a specified fraction of its peak value. The axial force was taken as the average force of the opposing top and bottom walls. The axial strain was calculated based on the change of distance between the opposing walls. **Figures 4A–I** represent, the PFCgenerated irregular ballast particles, DC test setups for the individual particles, and the broken ballast particles under peak loading, respectively.

**Table 1** lists the important parameters used for calibration of the DC tests. **Table 2** lists the peak force and displacement values from the laboratory and the DC crushing tests for the three ballast particles. Even though the peak force and displacement values did not match exactly, the calibration was considered complete, as the primary intention of this research effort was to get representative ballast calibrated strength model parameters which could be used to study the effect of real polyhedral ballast layer response on ballast breakage and PD.

# CYCLIC LOADING OF BALLAST LAYER COMPRISING BREAKABLE COMPLEX-SHAPED PARTICLES

# Specimen Preparation

First, a specimen box of size 600 × 600 × 300 mm was created using PFC3D. Targeting a porosity value of 38%, a total

3,235 ballast particles were then distributed in the specimen box using the particle size distribution shown in **Figure 5A**. Each size distribution range (varying from 15 to 50 mm) as shown in the inset table of **Figure 5A**, were composed of equal proportion of ballast particles conforming to the three different ballast shapes: B1, B2, and B3. The ballast particle templates used in the simulation are shown in **Figure 5B**. Each of the ballast particles B1, B2, and B3 comprised 12, 7, and 12 spheres, respectively. Note that each particle shape was represented by combining smaller number of spheres (12, 7, and 12, respectively) as the resulting reduction in computational times was necessary to complete all research tasks on time. For particles represented by agglomerates of larger number of spheres, the computational time requirements were unreasonably high (the simulation did not finish after running continuously for 75 days on a desktop computer with high processing power). While reducing the number of spheres in each agglomerate, special care was taken to retain the original shape of the ballast particles.

Once generated, the ballast particles were allowed to settle under gravity. After this step, a material pressure of 20 kPa was applied to the specimen by moving the surrounding walls via a servo mechanism to account for the confinement provided by the crossties, shoulder ballast, and crib ballasts. Note that researchers in the past have reported typical confining stress values in the range of 10–60 kPa (Indraratna et al., 2010a). **Figure 6A** shows the packed ballast assembly after this stage represented in the form of spheres of ballast templates shown in **Figure 5B**. **Figure 6B**, on the other hand, shows the packed assembly represented in the form of geometry of the corresponding particles in the simulation. The porosity at the end of this isotropic stress stage was 38.29%. All the model parameters used in the simulation have been listed in **Table 3**.

The value of effective modulus used in the simulation was taken as the average of the modulus values for the three ballast particles tested in the laboratory. The friction coefficient for ballast-ballast interaction was taken as 0.6 (tan 31◦ ) after Kwan (2006). The ballast layer was then applied with a cyclic loading around a mean stress of 232 kPa via the top wall of the specimen box.

# Contact Model for Ballast-Ballast Interaction

Most of the researches conducted on the DE modeling of railroad ballast use the linear-elastic contact model to represent the inter-particle contacts. However, the linear-elastic approach can only give a crude estimate of the ballast layer response. The behavior of ballast under loading has been established to be nonlinear (Lekarp et al., 2000). For this reason, in this research study, the ballast-ballast contact was modeled using a contact model recently introduced by Itasca: the Hill Contact Model (Potyondy, 2016). This contact model simulates the behavior of an infinitesimal, non-linear elastic (no tension) and frictional interface that carries a compressive surface interaction force and may carry a tensile moisture force. The surface interaction force model is based on the Hertz-Mindlin contact theory along with a damping mechanism and Coulomb sliding friction TABLE 1 | Important parameters used for DC test calibration.


TABLE 2 | Final calibrated results for single particle crushing test.


(Tsuji et al., 1992). The surface interaction force consists of Hertzian and dashpot components with Hertz-Mindlin springs providing the non-linear force-displacement response. The Hill Contact Model simulates the contact behavior between two locally elastic spheres that may have a liquid bridge and the material behaves like an unsaturated granular material. The liquid bridge is present if the moisture state is wet, and absent if the moisture state is dry. Detailed information on the Hill Contact Model can be found in the Itasca Technical Memorandum (Potyondy, 2016). In this simulation, the liquid bridge was not modeled, and hence the ballast-ballast contact model acted as a non-linear contact model; this can be said to be a better representation of the nature of interaction at interparticle contacts compared to linear-elastic assumptions. For the ballast-specimen box interaction, the linear contact behavior was assumed.

# Simulation of Breakable Ballast Particles

Once the mean stress of 232 kPa was applied to the system, the non-breakable ballast agglomerates were freed, and replaced with clusters of spheres with same radii at the same positions, and bonded together using linear-parallel bonds. For simplicity, the average of bond strengths and stiffness values obtained from the DC calibration tests were used in this simulation. Breakage of the bonds within the spheres in each cluster was considered as particle breakage. Upon breakage of the linearparallel bonds under the load, the Hill contact model was assigned to each new ballast-ballast contact. **Table 4** shows the bond strength parameter values for the linear parallel bonds used in this study.

After this step, cyclic loading was applied to the ballast layer via the top wall of the specimen box. Detailed description of the cyclic load applied in this simulation can be found in Dahal et al. (2018).

# RESULTS AND DISCUSSIONS

The ballast layer was subjected to 50 load cycles, and the stressstrain responses were recorded. Strain was calculated in terms of percentage, and was calculated from the distance between top and bottom wall before and after loading. Stress was calculated as the ratio of the load experienced by the top wall to the area of the top wall.

# Study of Ballast Breakage

**Figures 7A,B** show plan views of the ballast box before and after application of the cyclic loading, respectively. The ballast particles can be seen to have undergone significant breakage, and the initial and final ballast particle orientations cannot be distinguished due to the rearrangement of the broken ballast particles. These broken ballast particles eventually move to the void space in the assembly, and result in permanent deformation accumulation. Therefore, particle breakage and subsequent rearrangement of the broken particles act as governing factors for cyclic densification of the assembly.

To quantify the degree of particle degradation in terms of bond breakage in the cluster of spheres forming ballast particles, the number of linear parallel bonds before and after the cyclic load application were tracked, and were found to be 56,652 and 1,369, respectively; a total of 55,283 bonds were broken under cyclic loading. Most of these bonds were found to be broken during the initial load cycles leading to sudden increase in permanent deformation accumulation when loading was initiated. Once the bond breakage

#### TABLE 3 | Parameters used in the simulation.


TABLE 4 | Bond strength parameters for the linear-parallel bonds.


deformation accumulation under loading. **Figures 7C,D** show the linear parallel bonds in the model before and after cyclic load application, respectively. The extent of bond breakage is evident from the scarcity of intact bonds in **Figure 7D**.

The Contact Force (CF) chains between ballast particles before and after cyclic loading have been plotted in **Figures 7E,F**. As seen from these figures, the contact force chains are more

ceased, no significant increase in permanent deformation accumulation was observed. This clearly highlights the fact that particle breakage is the major source of permanent

uniform and widely distributed in the after-cyclic loading configuration, indicating good contacts between the ballast particles. This is because under loading the ballast particles undergo breakage resulting in the rearrangement of broken ballast particles in the void spaces present in the assembly. Note that the contact force chains in the plots are scaled by force, and the thickness of the lines correspond to the force magnitudes. Observing the contact force distribution after cyclic loading, it can be concluded that particle breakage plays a critical role in controlling the development of CF chains and governing the load distribution pattern in granular media under cyclic loading.

# Study of Ballast Permanent Deformation

**Figure 8A** shows the evolution of vertical strain with the number of load cycles under cyclic loading. As seen from the figure, there was an initial bulging (dilation) of ∼8.39% at the very beginning of the load application stage. This can be attributed to the phenomenon in DE modeling, where breakage of a bond between two overlapping spheres results in elastic rebound of the resulting two particles away from each other. Once this elastic energy has been released, there is a sudden increase in the vertical axial strain with the initial load cycles due to rearrangement and packing of the freshly generated broken particles. After the initial cycle, the axial strain is found to

increase in a steady manner and toward the end of 50 load cycles, the cyclic part of the axial strain is constant, and the permanent deformation accumulation is found to have stabilized. This can be explained by the fact that after the initial breakage, the broken ballast particles rearrange themselves in the void spaces, thereby densifying the ballast. The total vertical permanent strain accumulation was calculated as 14.14 mm (5.09%) after 50 load cycles.

**Figure 8B** shows the applied cyclic stress plotted against the vertical strain for 50 load cycles. Initially, the axial strain is negative in the figure which corresponds to bulging of the specimen. After the bulging is over, the axial strain accumulation after first few load cycles is significantly higher than under successive load cycles (about 5% permanent deformation accumulation). In this plot, the spikes in the stress-strain response during the first few load cycles can be attributed to particle breakages. However, once the initial load cycles have been completed, the ballast stress-strain response is more or less elastic showing a near-complete hysteretic behavior toward the end of 50 load cycles. This demonstrates the resilient behavior of ballast layer once it has been subjected to a large number of load cycles. Similar trends were reported by Indraratna et al. (2010b). The porosity at the end of 50 load cycles was found to be 33.44%, which was a reduction of 4.85% compared to the initial packing porosity. This can be attributed to the ballast breakage and the eventual rearrangement of the broken ballast particles under cyclic loading.

# Comparison of Results for Ballast Layers Formed With Only Single Ballast Shapes (Only B1, Only B2, and Only B3)

All the simulations reported above were also repeated on ballast assemblies comprising one ballast shape (B1, B2, or B3) only; the objective was to study the effect of each ballast shape on permanent deformation and ballast breakage under cyclic loading. The same specimen size and particle size distribution curve as shown in **Figure 1** was used for the study. The properties of the ballast layers were taken from the corresponding properties in **Table 1**. For direct comparison, the ballast layers were generated such that they had similar initial porosity (after application of the isotropic stress state). The initial porosities for ballast layers comprising only B1, only B2, and only B3 ballast particles were 38.30, 37.38, and 38.35%, respectively; it should be noted that matching the porosities of two models comprising differently-shaped particles in DE modeling is virtually impractical.

**Table 5** lists the final porosities, initial bulging, and permanent strain accumulation and change in no. of linear parallel bonds after the application of 50 load cycles for the four different models (one model with equal proportion of the three ballast shapes, and the three models comprising each of the different ballast shapes). **Figure 9A** shows the permanent strain accumulation with the number of load cycles for the four different ballast assemblies. As seen earlier, after the energy has dissipated from the release of overlaps between the spheres forming ballasts shapes, there is a sudden increase in the axial strain accumulation for all the ballast layers. The permanent strain (PS) curve has slowly stabilized and remained relatively constant toward the end of 50 load cycles.

From **Figure 9A** and **Table 5**, the ballast layer formed by only B2 shapes has undergone the highest PS, even though the peak force taken by B2 ballast is intermediate to B1 (highest) and B3 (lowest) (based on laboratory crushing test). This could be because of the improper packing of the B2 ballast particles owing to its more irregular and angular shape. Besides, the PS for B3 particles is less than that for the B1 particles, which could be attributed to more rounded shape of the B3 ballast compared to B1 ballast. B3 ballast is also shown to have the maximum densification toward the end of 50 load


TABLE 5 | Final porosity, initial bulging, permanent strain and no. of linear parallel bonds after 50 load cycles for different ballast layers.

cycles, as the porosity is the lowest of all three single ballast layers. Moreover, comparing the PS for ballast layers comprising single ballast shapes with that of the ballast layer comprising equally proportioned B1, B2, and B3 ballast shapes; the PS is observed to be in between lower (B1, B3) and higher (B2) range. This could be attributed to the dense packing created by the mixture of the elongated (B2) and more rounded (B1 and B3) ballast shapes.

**Figure 9B** and right-most three columns of **Table 5** show the shift in particle size distribution curve and no. of linear parallel bonds before and after cyclic loading for ballast layers formed with all ballast shapes (equally distributed B1, B2, and B3), only B1, only B2, and only B3, respectively. The area of particle breakage (area between the particle size distribution curves before and after cyclic loading) as found from **Figure 9B** are 1,643.31, 1,740.99, 1,368.66, and 1,748.27 mm<sup>2</sup> for all ballast shapes (B1, B2, B3), only B1, only B2, and only B3 ballasts layers, respectively. From these numbers and also from **Table 5** it can be seen that the model with B3 particles underwent the highest breakage (indicated by shift in the gradation curve after loading), which is consistent with the data presented in **Table 2** (B3 exhibited the lowest particle crushing strength). It is important to note that even though the model with B3 particles exhibited the highest particle breakage, the particles were immediately able to attain a "dense" configuration, thereby resulting in low permanent deformation accumulations.

# Importance of Accurate Representation of Ballast Particle Shape During DE Simulation

To assess the importance of accurate representation of ballast particle shapes during DEM simulations, the exact same simulations steps described above, were repeated for ellipsoid ballast particles. **Figures 10A,B** show the plots for permanent strain accumulation with number of load cycles and change in particle size distribution before and after cyclic load for the complex-shaped as well as ellipsoid ballast particles. Note that in the following discussion, the particles B1, B2, and B3, have been referred to as complex-shaped particles. It is important to note that this research effort did not use actual polyhedral particles in the simulation; rather, agglomerates of spheres were used to represent different particle shapes.

The initial bulging of the specimen for the model with complex-shaped particles was 8.39%, while the same for the model with ellipsoid particles was 3.59%. The accumulated permanent strain at the end of 50 load cycles was 5.09% for the model with complex-shaped particles, and 3.36% for the ellipsoid particles. The accumulated permanent strain for the model with the complex-shaped particles was approximately 51.5% higher than that with ellipsoid particles. Moreover, comparing the shift in particle size distribution before and after the cyclic loading for the two different particles shapes, the

complex-shaped particles undergo significantly higher degree of breakage compared to ellipsoid particles. The areas between the particle size distribution curves before and after cyclic loading were calculated to be 884.75 and 1,643.31 mm<sup>2</sup> for ellipsoid and complex-shaped particles, respectively. This shows that complex-shaped ballast particles underwent significantly higher breakage (about 85% more) compared to the ellipsoid particles under cyclic loading. From these results it can be concluded that accurate representation of ballast particles shapes is of utmost importance during DEM simulations to accurately capture particle breakage and permanent strain accumulation trends under cyclic train loading. At this point, it should be noted that differences in particle shape did not have a significant effect in the total permanent deformation accumulation for non-breakable particles, as reported in Dahal et al. (2018). However, when particle breakage is taken into consideration, particle shape plays a significant role in governing the total permanent deformation accumulation under loading.

# SUMMARY AND CONCLUSION


# LIMITATIONS AND SCOPE OF FURTHER RESEARCH

The modeling reported in this manuscript has the following limitations. Current and future research tasks by the authors aim to address some of these limitations.


3. The complex-shaped as well as ellipsoid particles used in this research effort were approximated using the agglomerates of bonded spheres. As an extension to this research, the model is currently being recreated using actual polyhedral particles (PFC Version 6.0). Findings from this modeling effort will be published in future manuscripts.

# DATA AVAILABILITY STATEMENT

The datasets generated for this study are available on request to the corresponding author.

# REFERENCES


# AUTHOR CONTRIBUTIONS

DM and BD: study concept. BD: numerical modeling and manuscript draft. DM: manuscript review and finalization.

# ACKNOWLEDGMENTS

The authors greatly acknowledge all the guidance and help provided by Dr. David Potyondy of Itasca Consulting Group. Mr. S. M. Naziur Mahmud, former graduate student at Boise State University helped with the initial modeling tasks.


**Conflict of Interest:** BD is currently employed by the company Bailey Engineering, Inc. She was a graduate research assistant at Boise State University during the performance period of this research work. No part of the research work was performed when BD was an employee of Bailey Engineering, Inc.

The remaining author declares that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2020 Dahal and Mishra. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Numerical Simulation of the Dynamic Response of Ballasted Track Overlying a Tire-Reinforced Capping Layer

Qideng Sun<sup>1</sup> \* † , Buddhima Indraratna2‡ and Jim Grant <sup>3</sup>

*<sup>1</sup> Coffey Services Australia Pty Ltd, Chatswood, NSW, Australia, <sup>2</sup> Transport Research Centre, University of Technology, Sydney, NSW, Australia, <sup>3</sup> Ecoflex International Pty Ltd, Avoca Beach, NSW, Australia*

#### Edited by:

*Sanjay Shrawan Nimbalkar, University of Technology Sydney, Australia*

#### Reviewed by:

*Cheng Chen, Wuhan University of Technology, China Yang Gui, Hohai University, China*

#### \*Correspondence:

*Qideng Sun Peter.Sun@Coffey.com*

*†Formerly worked as an Associate Research Fellow at University of Wollongong, Australia*

*‡Formerly worked at ARC Industrial Transformation Training Centre, ITTC-Rail, c/o University of Wollongong, Wollongong, NSW, Australia*

#### Specialty section:

*This article was submitted to Transportation and Transit Systems, a section of the journal Frontiers in Built Environment*

> Received: *23 August 2019* Accepted: *14 January 2020* Published: *05 February 2020*

#### Citation:

*Sun Q, Indraratna B and Grant J (2020) Numerical Simulation of the Dynamic Response of Ballasted Track Overlying a Tire-Reinforced Capping Layer. Front. Built Environ. 6:6. doi: 10.3389/fbuil.2020.00006*

This paper describes a 3D finite element (FE) model developed to understand the dynamic response of a ballasted track in which the underlying capping layer is reinforced using recycled rubber tires. Track deflection, the lateral spreading of ballast and vertical stress transmitted from the capping layer to the subgrade are discussed by considering the effect of reinforcement provided by these infilled tires. In this respect, the capping layer is confined and has improved damping properties. The cellular structure of the rubber tire assembly can radially confine the infilled materials, and thus reduce excessive lateral spreading and vertical displacement that would otherwise occur in a conventional track. At the same time the tire and gravel composite layer acts like a stiff but flexible "mattress" that controls the stress transmitted to the underlying subgrade while making it more uniform. Typical soft and stiff subgrade materials were used to investigate the dynamic response of track, and the stress paths of subgrade at different depths have been studied. It is noted that the effect of the tire assembly on the stress distribution within the subgrade decreases with depth, and the tire-reinforced track deflects less than its unreinforced counterpart at any given train speed.

Keywords: track dynamics, finite-element modeling, vibration, reinforced soils, capping layer, scrap tire

# INTRODUCTION

The need for long-term performance of rail infrastructure has been increasingly higher in the past decades as rail tracks systems supporting the transport network are expected to withstand higher speeds and larger loads. Heavier and faster trains could exert higher dynamic wheel loads on track, and therefore a soft subgrade may experience higher repeated stresses which may lead to excessive deformation and progressive shear failure. Moreover, under larger dynamic loads (e.g., rail freight traffic), an existing track may degrade further and faster due to excessive track deformation and the lateral displacement of ballast requiring more frequent maintenance. Geocells offer a plastic cellular confinement system with a honeycomb-like structure that can be filled with an appropriate granular material, and has been successfully used in ground improvement applications including transport infrastructure (Leshchinsky and Ling, 2012, 2013; Indraratna et al., 2015). During loading, additional confinement is provided by the geocells as they mitigate the granular mass of subballast from spreading laterally, increase the rigidity of the infill, and improve the loadcarrying capacity (Zhang et al., 2010; Leshchinsky and Ling, 2012; Sitharam and Hegde, 2013), and thereby improving the overall track performance. Similarly to geocells, scrap tires can provide additional cellular confinement to the infill materials and also improve the strength and stiffness of a railway track. Indraratna et al. (2017) proposed a rubber tire based capping layer for railway track where one sidewall is removed and the tire is then filled with gravel; a geotextile can be placed between the rubber tire and the ground as a separator. Plate load tests (Indraratna et al., 2017) and prototype process simulation carried out by Indraratna et al. (2018) revealed that a tire cell has three primary engineering benefits: (i) the confinement provided by the cellular assembly can increase the stiffness of the contained aggregate, which then reduces lateral spreading and vertical deformation within the capping and ballast layers; (ii) the tires and gravel composites allow a reduced and more uniform stress distribution to the subgrade; and (iii) tire cells can enhance the damping properties of the system, and thus its ability to attenuate dynamic forces imposed by rail traffic. However, since the dimensions are limited by the size of the experimental facility (i.e., 800 × 600 × 600 mm), a light vehicle tire having a diameter of

560 mm was tested. In reality, heavy vehicle tires with a diameter of about 1 m are preferred because they have much larger stiffness compared with car tires. The tire-reinforcement is very similar to geocell reinforcement. The benefits of geotextile and geogrid reinforcement are: (i) geotextile reinforcement restricts the lateral deformation of the soil; (ii) strength improvement of the reinforced soil is attributed to openings in the geogrid that causes interlocking; and (iii) reinforcement effect in undrained condition is an increase in cohesion and the effect in drained condition is an increase in internal friction.

Theoretical and numerical models have been developed by numerous studies to better understand the amplitude of vibration of the track-ground system during train passage. The approaches adopted including: (i) analytical models (Kenney, 1954; Krylov, 1995; Sheng et al., 1999; Kaynia et al., 2000; Madshus and Kaynia, 2000; Kargarnovin and Younesian, 2004; Kargarnovin et al., 2005; Picoux and Le Houedec, 2005; Takemiya and Bian, 2005; Karlstrom and Bostrom, 2006); (ii) numerical models based on the finite element (FE) and boundary element (BE) methods (Hall, 2003; Sheng et al., 2003, 2005; Yang et al., 2003; Kouroussis et al., 2009; Lombaert and Degrande, 2009; Costa et al., 2010, 2012; Ju et al., 2010; El Kacimi et al., 2013; Shih et al., 2016, 2017). In design of railway embankments, simplified theoretical and empirical methods are typically used by assuming a homogeneous half-space for all the track layers, without considering the individual layer properties (e.g., Okabe, 1961; Heath et al., 1972; AREA, 1996). Multilayer track models, such as ILLITRACK (Robnett et al., 1975), GEOTRACK (Chang et al., 1980), and KENTRACK (Huang et al., 1986) were developed to analyze stresses in the track and subgrade. These models were based on the assumption that substructure materials are purely elastic and this will lead to inaccurate prediction. Recently, a finite element 3D model has been developed to investigate the load transfer mechanism between tires and infill gravels under static loading (Indraratna et al., 2017). It was reported that the confining effect causes the tire and gravel infill composite to act as a stiff but flexible "mattress" that can effectively reduce the amplitude of stresses transmitted to the subgrade.

In this paper, a 3D FE dynamic model is developed to simulate both the track super and substructure and rail traffic subjected to moving point loads that apply onto the beam elements representing the rail system. A realistic scenario to validate the 3D FE model is the well-documented case study of ground vibrations generated by X2000 trains at the Ledsgard site in Sweden (Hall, 2003). Then, for simplicity and ease of simulation, another track with a simpler substructure is adopted, and a typical three-wagon Australian freight train is simulated accordingly. The effect of tires on infilled gravels is simulated by changing material properties of capping, namely, Young's modulus (Ec), apparent drained cohesion (Cc), and damping ratio (ξ ). The friction angle of the granular mass was kept constant.

# MODELING THE MOVING VEHICLE PROBLEM

To model the dynamic response induced by a moving train, separate models may be used to represent the vehicle and the track/ground system; elements with rigid bodies can be used for vehicle components, such as car bodies, bogies, and wheelsets, while the primary and secondary suspensions can be modeled by the springs and dampers that connect to the rigid bodies. With a moving vehicle problem, the loads depend on the dynamic response of the vehicle and the track system, so a user-defined subroutine is often used to combine both programs despite having significant impact on the computational time. However, two other approaches can be adopted in ABAQUS to integrate both systems one is by using the large (finite) sliding contact model in ABAQUS (Hibbitt, Karlsson & Sorensen, Inc., 2014) to connect the vehicle and track/ground system, and the other is calculating the moving load by specifying the nodal forces that correspond to the moving load as a function of time history. While the latter approach omits the vehicle dynamics, it can greatly reduce the associated computational time and thus adopted in this study as shown in **Figure 1A**.

In this approach, the rail is essentially built up from nodes and beam elements. In this paper the nodes representing the rail are referred to as loading nodes, and every fifth loading node is connected to the elements forming the sleepers (ties) by spring elements as shown in **Figure 1A**. Since the spacing between the sleepers is 0.6 m, the spacing between the loading nodes is 0.12 m. Point loads are applied at the loading nodes and for a given speed (V), the loads can be thought of as triangular pulses distributed between three nodes, that can be moved from node to node by a time step (1t) given by 1t = 1x/V, where 1x is the node spacing as shown in **Figure 1A**.

This loading model is similar to the procedure described earlier by Hall (2003). It should be noted that the "rail" is only connected to the rest of the finite element model by the loading nodes in the "sleepers." The load on the "rail" is therefore only transferred to the rest of the model through the "sleepers." This is very similar to reality. It is acknowledged that numerical models have been developed where no beam elements were used. Load distributions was calculated for application to the sleepers (e.g., Feng et al., 2019). Comparison for a one-point load was

made between two load models and it was indicated that the two loading models gave approximately the same results (Hall, 2003). This approach has been tested on a simple moving vehiclebridge interaction problem (**Figure 1B**). The large (finite) sliding contact model that simulates the moving oscillator in ABAQUS (Hibbitt, Karlsson & Sorensen, Inc., 2014) is also adopted for comparison. The material parameters for the beam are: Young's modulus E = 2.87 GPa, Poisson's ratio ν = 0.2, second moment of area I = 2.94 m<sup>4</sup> , mass per unit length ρA = 2,303 kg/m. The suspension stiffness (kv) is 1,595 kN/m, and the mass (Mv) is 5,750 kg. The vehicle speed (V) is 100 km/h and the length of the beam L = 25 m.

The dynamic response of the moving load and moving oscillator problems obtained using ABAQUS is compared to the numerical solution obtained by Yang and Yau (1997) in **Figure 1C**. The displacements obtained for the beam agree with the analytical result, and the difference in the beam response introduced by the dynamics of the spring mass is insignificant compared to the result for a moving load. Hence, in this study the moving load procedure has been considered.

# CASE STUDY-3D FE MODEL FOR LEDSGARD TRACK

Having established a method for simulating a moving vehicle in the FE model, a 3D model of the vehicle-track-ground system has been developed using ABAQUS 6.14-2. To ensure the modeling process can provide reliable outcomes, a case study that is well-documented in the literature (Hall, 2003) has been used to ascertain the reliability and accuracy of the model. The FE mesh of the track and ground at the Ledsgard site is shown in **Figure 2A**. This FE model makes use of symmetry in the x-y plane and fixed boundaries are used. The 3D model measures 104 m long × 40 m wide × 30 m high and consists of 242,962 elements. The length of the model required to achieve convergence depends on the load speed, so longer models are needed if the load speed approaches or exceeds the critical speed. In this instance the model is 104 m long, which is sufficient for convergence (Shih et al., 2016). If the model is wide enough there is little benefit in using absorbing elements at the boundaries. This is because a damping model with a large enough massproportional term allows the energy to dissipate and avoids any reflections interfering with the results (Shih et al., 2016). Hence, 30 m deep by 40 m wide FE mesh with Rayleigh damping model is used in this study. The damping ratio in the Rayleigh damping model is frequency dependent, so for the parameters adopted in the damping model, a damping ratio of about 4% (**Table 1**) in the frequency range of interest is obtained, as suggested by Hall (2003).

The mesh of the embankment and ground are mainly composed of hexahedral 8-noded elements (C3D8). The rail is modeled with Timoshenko beam elements (B31) that include the effects of shear deformation and rotational inertia. While both of these effects are ignored in standard Euler-Bernoulli beam theory, they are better suited to modeling thick or slender beams. The beam is given the geometry and properties of regular railway rails (UCI60), and to represent the rail pads, linear springs with a vertical stiffness 4.8 × 10<sup>8</sup> N/m connect the rail to each sleeper. Discrete sleepers with a spacing of 0.6 m are included; they have a half-length of 1.25 m, a height of 0.2 m, a mass density of 2,400 kg/m<sup>3</sup> , a Young's modulus of 3 × 10<sup>10</sup> N/m<sup>2</sup> and Poisson's ratio of 0.15. The ballast is embedded into the upper ground layer to a depth of 0.3 m. A 0.3 m thick interface layer is laid underneath the ballast. The smallest element is 0.25 m near the track, but it gradually increases in size with a stretch factor of 1.2 in the horizontal direction outside the width of the track, and in the vertical direction for the clay layer. The total number of degrees of freedom in the model is 806,913. The soil properties for the case study are available from various in-situ tests, such as seismic core penetration testing, cross-hole tests, and spectral analysis of

surface waves (Hall, 2003; Costa et al., 2010). There are five soil layers, namely the ballast, interface, dry crust, organic clay, and clay. The material properties used for the various layers are listed in **Table 1** (Hall, 2003; Sayeed and Shahin, 2016). The results for a set of moving axle loads are presented, but the surface roughness excitation is omitted. The geometry of the high speed train X2000, is reported in more detail by Kouroussis et al. (2014) and represented in **Figure 2B** and **Table 2**, but a full train set consisting of a driving trailer vehicle, three passenger carriages, and a locomotive is considered for the current analysis.

The time-history response of track displacement during the passage of train at speeds of 70 and 200 km/h are calculated at the track center, and then the results are compared to the corresponding field measurements, as shown in **Figures 3A,B** (downward negative). Note that at 70 km/h, only quasi-static deflection appears when the load moves over the point of concern, but an oscillatory response occurs at a higher speed of 200 km/h. The FE predictions agree reasonably well with the field measurements, which proves that the FE modeling process in this study is reliable and can be used with confidence to predict the


TABLE 1 | Track and ground properties adopted for FE model for the site at Ledsgard (Hall, 2003; Sayeed and Shahin, 2016).

TABLE 2 | Geometry and axle loads of the X2000 high speed train (Kouroussis et al., 2014).


track behavior under a moving train load. To investigate the effect that material elasto-plasticity has on the dynamic track response, the traditional Mohr-Coulomb model is adopted to simulate ballast with a friction angle φ = 50◦ . The model predictions in terms of (a) track deflection at the midpoint of track, (b) vertical stress transmitted from ballast to the interface, and (c) lateral displacement of ballast, are compared to the predictions made with an elastic model at a speed of 70 km/h. **Figures 3C–E** shows the differences in predictions based on the elastic and the Mohr-Coulomb model. The latter predicts greater track deflection, larger vertical stress and greater ballast spreading. This occurs because the elasto-plastic analysis allows plastic deformation to develop. So in the following analysis the Mohr-Coulomb elastoplastic model is considered to model the track materials.

# MODELING TIRE-REINFORCED CAPPING LAYER

To study the effects induced by rubber tires in the capping layer, an FE model with a simplified profile is used (**Figure 4A**). Traditional capping layer material (i.e., crushed basalt) is modeled with Young's modulus equals to 140 MPa, frictional angle equals to 38◦ and dilation angle equals to 5◦ . The model consists of a 0.35 m thick ballast layer, and a 0.15 m thick capping layer founded on 30 m thick ground. A train with three typical freight wagons has been modeled in this analysis to mimic Australian rail freight traffic (ARTC, 2011). Details of the freight train are shown in **Figure 2C** and **Table 3**. A 25 t axle load is applied in the model, and the train runs along the rail in the positive direction of the z axis at various speeds (V).

Laboratory tests revealed that the cylindrical structure of a tire confines the infill material and reduces its lateral displacement (Indraratna et al., 2017), whereas the vertical load could produce circumferential tensile strain in the tire that would induce more confining stress 1σ<sup>3</sup> ′ from the rubber tire to the capping gravels. By assuming that the internal friction angle of gravel for the reinforced sample remains constant (φ = 38◦ ), the apparent cohesion C<sup>c</sup> can be used to account for the increasing strength of the tire cell (Bathurst and Karpurapu, 1993) and can be readily incorporated in ABAQUS:

$$C\_c = \frac{\sigma\_3 r}{2} \cdot \tan\left(\frac{\pi}{4} + \frac{\phi}{2}\right) \tag{1}$$

where φ is the friction angle of the infill material. The increased confining stress 1σ<sup>3</sup> ′ , can be estimated by the following equation (Bathurst and Karpurapu, 1993).

$$
\Delta \sigma\_3^{'} = \frac{2M\_l}{D} \left[ \frac{1 - \sqrt{1 - \varepsilon\_a}}{1 - \varepsilon\_a} \right] \tag{2}
$$

where M<sup>t</sup> (=45 MN/m) is the tensile stiffness of the tire, D (=1 m for a truck tire) is the initial diameter of the tire, and ε<sup>a</sup> is the axle strain of the sample. Bathurst and Karpurapu (1993) used a similar method to account for the increasing strength of geocell reinforcement. In the 3D FE model, two values of C<sup>c</sup> (i.e., 5 and 46 kPa corresponding to 1σ<sup>3</sup> ′ of 4.5 and 45 kPa at ε<sup>a</sup> = 0.01 and 0.1%, respectively) are used to study the increase of additional confinement.

According to Indraratna et al. (2018), there is ∼15% improvement in the damping ratio of the sample with a tire compared to the sample without a tire. A viscous damping model based on Rayleigh damping, is used in the 3D FE model. Rayleigh damping is based on the two parameters α and β, both of which allow the damping matrix C to be determined from the mass and stiffness matrix M and K (Hibbitt, Karlsson & Sorensen, Inc., 2014):

$$C = \alpha M + \beta K \tag{3}$$

FIGURE 3 | (A) Comparison between FE simulation and field measurement of track response for X2000 train with a speed of 70 km/h; (B) with a train speed of 200 km/h; (C) comparison of predictions with different ballast models for track deflection at the midpoint of the track; (D) vertical stress transmitted from ballast to interface; and (E) lateral displacement of ballast.

This allows the equivalent loss factor η or damping ratio ξ to be obtained as a function of frequency as (Hibbitt, Karlsson & Sorensen, Inc., 2014):

$$\frac{\eta}{2} = \xi = \frac{\alpha}{2\alpha} + \frac{\beta\alpha}{2} \tag{4}$$

where ω is the circular frequency where the loss factor η applies. The circular frequency is usually selected to correspond to the resonance frequency in the system because the damping resonances are important. Rayleigh damping (stiffnessproportional) that is equivalent to Kelvin-Voigt damping can be obtained by setting α equal to zero (Connolly et al., 2013), while mass-proportional Rayleigh damping can be obtained by setting β equal to zero (El Kacimi et al., 2013). Combined Rayleigh damping can be obtained by fitting the damping ratio to that measured at the site, as proposed by Hall (2003). In the 3D FE model, damping ratios ξ = 0.09 and 0.18 are used with and without tires. The effect that rubber tires have on the stiffness of the capping layer is reflected by Young's modulus variation, i.e., 220 and 500 MPa were used in the 3D model to study the increase in Young's modulus. The values of the model parameters are shown in **Table 4** (Indraratna et al., 2017, 2018).

TABLE 3 | Geometry and axle loads of the idealized three-wagon Australian freight train (ARTC, 2011).


# MODEL PREDICTIONS

# Dynamic Track Response

A series of simulations have been carried out to determine how the tires improve track stiffness, and how they affect the increment of additional confinement and damping ratio. To compare the results, each scenario has been analyzed with and without tire reinforcement for comparison. During this simulation, track deflection, lateral displacement of the ballast layer and the vertical stress transmitted from the capping layer to the subgrade layer were observed.

Deflections at the midpoint of the track are plotted in **Figure 5**, and it shows that with tire reinforcement the track experiences less deflection than that without tires. **Figure 5A** shows that the increment of Young's modulus in the capping layer only has a marginal effect on track deflection, but as drained cohesion increases, the track deflection decreases (**Figure 5B**). Moreover, incremental increases in damping in the capping layer provided by the inclusion of tires seems to have minimal influence in the track deflection even though the damping ratio increases from 0.09 to 0.18 (**Figure 5C**).

**Figure 6** shows how tire reinforcement influences the lateral displacement of ballast particularly at the shoulder location as shown in **Figure 6A**. The results indicate that rubber tires could help to reduce the lateral displacement of ballast, and with the increase of Young's modulus and cohesion in the capping layer, the ballast layer experiences less lateral displacement, as shown in **Figures 6B,C**. For instance, at the crest of the ballast layer a reduction in 0.6% is observed with the inclusion of tires (**Figure 6B**). This improvement in the stiffness of the capping layer by additional confinement by tires also reduces lateral displacement of the capping layer and ballast layer and reduces the vertical displacement.

Rayleigh damping is based on the parameters α and β, which allow the damping matrix to be determined from the mass and stiffness matrices. In this study a constant damping ratio (ξ )

TABLE 4 | Track and ground properties used to investigate the effect of rubber tire (Indraratna et al., 2017, 2018).


of 0.09 is assumed in the model without reinforcement; the corresponding parameters α = 0.464 and β = 5.99 × 10−<sup>4</sup> at 10 Hz. The damping ratio ξ , of 0.18 is used for the model with reinforcement. The stiffness-proportional Rayleigh damping is studied by allowing α to equal 0.464. Hence, the increment of ξ from 0.09 to 0.18 is reflected by an increment of β from 5.99 × 10−<sup>4</sup> to 1.315 × 10−<sup>3</sup> . Conversely, mass-proportional Rayleigh damping can be obtained by setting β equal to 5.99 × 10−<sup>4</sup> , this results in α = 3.2896 and 0.464, which corresponds to ξ = 0.18 and 0.09, respectively. **Figure 6D** shows that the value of β that corresponds to the stiffness-proportional Rayleigh damping has an apparent influence on the lateral displacement of ballast.

**Figures 7A–C** shows that rubber tires can reduce the pressure transmitted from the capping layer to the subgrade layer, but the variation in Young's modulus and the damping ratio have negligible effect on the vertical stresses. This reduction in vertical stress in the capping layer incorporating rubber tires is due to the additional confinement provided by the rubber membrane which also reduces the vertical stress, as shown in **Figure 7B**. Intuitively, the confining effect causes the tires and gravel infill composite to act as a stiffer, flexible "mattress" which results in a reduced and more uniform stress distribution.

# Effect of Different Types of Subgrade

The dynamic response of the track-ground system has been investigated considering two typical subgrades, one with soft soil (i.e., fat clay of the Monroe dam, E = 6.9 MPa; Duncan et al., 1980) and the other with stiff soil (i.e., low density sand, E = 58.5 MPa; Al-Shayea et al., 2003). The properties of the subgrades are summarized in **Table 4**. Two different scenarios (i.e., with and without reinforcement) are modeled with each type of subgrade. Model properties for reinforced capping layer are E<sup>c</sup> = 220 MPa, C<sup>c</sup> = 5 kPa and ξ = 0.18 with α = 0.464 and β = 1.315 × 10−<sup>3</sup> . **Figure 8A** shows deflection obtained at the track midpoint with different subgrades. As expected, the track experiences much less deflection on a stiff subgrade than on a softer subgrade. For example, at a train speed of 216 km/h the maximum deflection of the unreinforced track with soft and stiff subgrade is 0.024 and 0.003 m, respectively, which is because track built on soft subgrade usually yields high ground vibrations at low train speed than those founded on stiff subgrade. The reduction of settlement for a track substructure reinforced with tires is more obvious for the cases where the track is supported by a soft subgrade material. For instance, while a reduction of 5.3% is observed for soft subgrade case, no discernible difference was observed for the stiff subgrade (**Figure 8A**).

One obvious advantage of the tire is its ability to redistribute the stress and reduce the magnitude of the subgrade stresses as shown in **Figure 8C**. Unlike the unreinforced track, the peak stress of reinforced track decreased by ∼30.9 and 10.6% with soft and stiff subgrade, respectively. It is noteworthy that the vertical stress on a stiff subgrade is higher than on a soft subgrade. This is probably because the lower damping ratio of the stiffer material leads to a higher resistance. Intuitively, lowering the vertical stress on the subgrade reduces the vertical and lateral displacement of the track. Using rubber tires in the capping layer decreases lateral spreading by 47.3 and 83.9% with soft and stiff subgrade, respectively (**Figure 8B**).

# Stress Paths in Subgrade

drained cohesion, and (D) damping ratio.

The dynamic stress induced by train passage in the track substructure is one of the most important parameters in railway design and maintenance (Bian et al., 2014). Determining the intensity of dynamic stress at the ground surface and its attenuation in the subgrade is crucial to proper design. The dynamic stress-time history curves and stress paths of the models are compared with and without rubber tires. There are three observation points in the mid cross-section of the model at depths of 0.4, 2.2, and 5.8 m, these points are V1, V2, and V3, as shown in **Figure 4B**. **Figure 9** shows a comparison of the dynamic soil stress-time history curves under a load of three moving carriages; they consist of the vertical stress Syy and shear stress Syz at points V1, V2, and V3. It is important to note that the vertical stress and shear stress in the subgrade decrease with depth. It can be observed that the tire reinforcement can reduce the vertical and shear stress, but the reduction is more evident for shallow depths. In fact, for V3 point located at 5.8 m the reduction is marginal. **Figure 9A** also shows that the train axles are not visible in the dynamic stress simulated at the subgrade level, and only the distribution of adjacent bogies can be clearly seen at V1. In contrast, for the patterns of vertical stresses simulated at V3, the dominant loads are associated with train carriages. This difference indicates that there are different dominant frequencies across the subgrade layer determined by the different depths.

**Figure 10** shows the stress paths (i.e., vertical stress Syy against shear stress Syz) followed by the soil element at V1 during the passage of the first bogie (i.e., 1-2-3-4-5) and second bogie (i.e.,

5-6-7-8) of the first carriage. In stage 1–2, the vertical stress and shear stress increase as the first bogie approaches, while in stage 2–3, the shear stress decreases while the vertical stress continues to increase until it peaks when the first bogie arrives at position 3.

After, the vertical stress begins to decrease, while the shear stress reaches the opposite peak at position 4. Finally, the vertical stress and shear stress both decrease as the first bogie arrives at position 5. As the second bogie approaches, the soil element experiences a different stress path because of the stress superposition caused by two bogies between adjacent carriages. For example, when the second bogie arrives at position 8, the second loading cycle for the soil element is complete, but the vertical stress at position 8

is much higher than at position 5 at the end of first loading cycle. At a lower depth (V2) beneath the track, one bogie corresponds to one loading cycle, but as the depth of soil increases the vertical stress at the end of stage 7–8 moves toward its maximum value, and this may result in two bogies corresponding to one loading cycle. Not surprisingly, the soil near the track experiences more and larger stress cycles than soil at greater depths. This increase in the number of cycles and magnitudes of train loads is likely to result in larger permanent deformation, so the soil at shallower depths is more likely to experience increased displacement and possible shear failure than soil in deeper ground. The use of rubber tires can effectively reduce the magnitude of stress.

# Track Deflection at Critical Speed

**Figure 11A** shows the predicted track deflection vs. train speed for tracks with and without tires. It can be observed that deflection generally increases as the train speed increases

and reaches a maximum deflection at the critical speed, and then decreases as the train speed increases further. Moreover, deflection at high speed with tires is less than without tires. However, the results in **Figure 11A** indicate that the critical speed for track is not sensitive to the inclusion of tires in the capping layer. For the purpose of illustrating the impact of train speed on the track, the contour plots of vertical deflection along the track with tire reinforcement are shown in **Figures 11B–D** at three different speeds (i.e., 36, 144, and 360 km/h). **Figure 11B** shows that at a train speed of 36 km/h, which is considerably lower than the critical speed, vertical deflection is mainly induced near the axles and as expected, there is a slight propagation of wave to the surrounding ground. However, **Figure 11C** shows that at a critical speed of 144 km/h, vertical deflection is induced near the axles and in the surrounding ground. It is noteworthy that the series of wave fronts radiating from the loading positions appear as a shockwave in the ground that is known as the Mach cone (Krylov, 2011). **Figure 11D** shows that the vertical deflection of the ground at a train speed 360 km/h which is greater than the critical speed. For this case, the loading speed is greater than the wave speeds and the source passes through wave fronts.

# CONCLUSIONS

A 3D track/ground FE model was developed to investigate the dynamic response of ballasted railway track with a rubber tirereinforced capping layer. The moving load was considered by specifying the nodal forces corresponding to the moving load as a function of time. The 3D model was verified and good agreement was found when compared to the field measurements taken at the soft soil site in Ledsgard, Sweden. Various train speeds (i.e., 36–360 km/h) have been simulated to identify the critical train speed for an idealized ground condition. It was found that the critical train speed for the ground condition in this study was 144 km/h. The impacts of break and acceleration of the trains

on the railway track degradation and deformation have attracted research interests recently. Two typical subgrades (i.e., fat clay of the Monroe dam representing soft soil and low density sand representing stiff soil) have been considered to investigate the effect of tire-reinforcement on the dynamic response of trackground system. One of the purposes of this research was to inspire engineer when they design railway tracks. The author acknowledges further detailed research needs to be done on the above mentioned points in the future.

Numerical predictions indicate that track with tire reinforcement experienced less deflection than tracks without tires; this is mainly due to the additional confinement provided to the infill materials that reduces lateral spreading, both of which help to reduce the vertical and lateral displacement of ballast. However, it is acknowledged that further experimental and field test results are needed to validate the conclusion obtained from the numerical simulation. For modeling, an increase in apparent cohesion could be considered to simulate the tire-fill composite. Moreover, the increasing damping property in the capping layer influences the lateral displacement of ballast, an effect that is dominated by the stiffness-proportional Rayleigh damping. As expected, the soil close to the track experiences more and larger stress cycles than soil at greater depths. This increase in the number of cycles and magnitude of train loads will lead to larger permanent deformation, so the soil at shallower depths is more likely to experience increased displacement and possible shear failure than the soil located in deeper ground. The additional confinement provided by tires could reduce the vertical and shear stress, but the improvement is marginal for greater depths. Tire reinforcement is very useful when the substructure overlies a soft foundation; for example, rubber tires in the capping layer reduced lateral spreading by 47.3%. Moreover, track deflection at critical speed with tires is less than without tires, but the critical speed for track is relatively insensitive to the presence of tire reinforcement.

# DATA AVAILABILITY STATEMENT

All datasets generated for this study are included in the article/supplementary material.

# AUTHOR CONTRIBUTIONS

BI had provided very critical and constructive review comments at different stages of this work, and was the project leader of this sponsored industry-based research. JG had participated in discussions and provided inspirational comments at various stages of this work, and his company EcoFlex has already put to practice the rubber-tyer assembly concept in numerous Australian road projects. QS approved the article for publication.

## ACKNOWLEDGMENTS

The first author formerly worked as an Associate Research Fellow at University of Wollongong. The authors acknowledge the financial assistance provided by the NSW Environmental Trust, TSA, Ecoflex International Pty Ltd and University of Wollongong. The original concept of tyer assembly for roads was an invention by EcoFlex. The authors are

## REFERENCES


thankful to Dr. Trung Ngo and Dr. Ana Heitor for their valuable review comments to this paper. Assistance of A/Prof. Cholachat Rujikiatkamjorn at various times during this project is also appreciated. Contributions from Dr. Mahdi Biabani through similar conceptual research based on geocell reinforced tracks under the guidance of Distinguished Professor Indraratna supported by EcoFlex also inspired this particular study.

Indraratna, B., Sun, Q., Heitor, A., and Grant, J. (2018). Performance of rubber tyre-confined capping layer under cyclic loading for railroad conditions. J. Mater. Civil Eng. 30:06017021. doi: 10.1061/(ASCE)MT.1943-5533.0002199


**Conflict of Interest:** QS was employed by the company Coffey Services Australia Pty Ltd. JG was employed by the company Ecoflex International Pty Ltd.

The remaining author declares that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2020 Sun, Indraratna and Grant. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# LIST OF SYMBOLS


# Corrigendum: Numerical Simulation of the Dynamic Response of Ballasted Track Overlying a Tire-Reinforced Capping Layer

Qideng Sun<sup>1</sup> \* † , Buddhima Indraratna2‡ and Jim Grant <sup>3</sup>

*<sup>1</sup> Coffey Services Australia Pty Ltd, Chatswood, NSW, Australia, <sup>2</sup> Transport Research Centre, University of Technology, Sydney, NSW, Australia, <sup>3</sup> Ecoflex International Pty Ltd, Avoca Beach, NSW, Australia*

Keywords: track dynamics, finite-element modeling, vibration, reinforced soils, capping layer, scrap tire

#### **A Corrigendum on**

#### Approved by:

*Frontiers Editorial Office, Frontiers Media SA, Switzerland*

# \*Correspondence:

*Qideng Sun Peter.Sun@Coffey.com*

*†Formerly worked as an Associate Research Fellow at University of Wollongong, Australia*

*‡Formerly worked at ARC Industrial Transformation Training Centre, ITTC-Rail, c/o University of Wollongong, Wollongong, NSW, Australia*

#### Specialty section:

*This article was submitted to Transportation and Transit Systems, a section of the journal Frontiers in Built Environment*

> Received: *10 June 2020* Accepted: *22 June 2020* Published: *07 August 2020*

#### Citation:

*Sun Q, Indraratna B and Grant J (2020) Corrigendum: Numerical Simulation of the Dynamic Response of Ballasted Track Overlying a Tire-Reinforced Capping Layer. Front. Built Environ. 6:115. doi: 10.3389/fbuil.2020.00115* **Numerical Simulation of the Dynamic Response of Ballasted Track Overlying a Tire-Reinforced Capping Layer**

by Sun, Q. (2020). Front. Built Environ. 6:6. doi: 10.3389/fbuil.2020.00006

Buddhima Indraratna and Jim Grant were inadvertently omitted as authors in the published article. The corrected Author Contributions Statement and Affiliations appear below:

The following affiliation for the 2nd author Buddhima Indraratna was added: "Transport Research Centre, University of Technology, Sydney, NSW, Australia. Formerly worked at ARC Industrial Transformation Training Centre, ITTC-Rail, c/o University of Wollongong, Wollongong, NSW, Australia."

The following affiliation for the 3rd author Jim Grant were added as "Ecoflex International Pty Ltd, Avoca Beach, NSW, Australia."

The corrected Author Contributions Statement appears below:

BI had provided very critical and constructive review comments at different stages of this work, and was the project leader of this sponsored industry-based research. JG had participated in discussions and provided inspirational comments at various stages of this work, and his company EcoFlex has already put to practice the rubber-tyer assembly concept in numerous Australian road projects. QS approved the article for publication.

Furthermore, in the original article, there was an error in the acknowledgments section. Funding and an acknowledgments of contributions of additional researchers have been added. The following correction has been made:

"The first author formerly worked as an Associate Research Fellow at University of Wollongong. The authors acknowledge the financial assistance provided by the NSW Environmental Trust, TSA, Ecoflex International Pty Ltd and University of Wollongong. The original concept of tyer assembly for roads was an invention by EcoFlex. The authors are thankful to Dr. Trung Ngo and Dr. Ana Heitor for their valuable review comments to this paper. Assistance of A/Prof. Cholachat Rujikiatkamjorn at various times during this project is also appreciated. Contributions from Dr. Mahdi Biabani through similar conceptual research based on geocell reinforced tracks under the guidance of Distinguished Professor Indraratna supported by EcoFlex also inspired this particular study."

The Conflict of Interest statement has also been updated to reflect the addition affiliation of Ecoflex International Pty Ltd. The new Conflict of Interest statement reads as follows: "QS was employed by the company Coffey Services Australia Pty Ltd. JG was employed by the company Ecoflex International Pty Ltd.

The remaining author declares that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest."

The authors apologize for these errors and state that this does not change the scientific conclusions of the article in any way. The original article has been updated.

Copyright © 2020 Sun, Indraratna and Grant. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Field Assessment of Gravel Loss on Unsealed Roads in Australia

#### Vasantsingh Pardeshi\*, Sanjay Nimbalkar and Hadi Khabbaz

*School of Civil and Environmental Engineering, University of Technology Sydney, Ultimo, NSW, Australia*

The gravel loss is a major limitation for unsealed roads and it needs major maintenance annually. The continual process of gravel loss leads to the unsustainability of these roads. The unsealed road management faces several issues, viz., difficulty to forecast behavior, huge data collection needs, and a vulnerability in the service and maintenance practices. The quality of gravel material also plays a major role in the process of gravel loss. In view of the aforementioned, appropriate revisions to ARRB material specifications are proposed in this study. The gravel material as per modified ARRB specifications is used on the unsealed road network in the Scenic Rim Regional Council in the state of Queensland. Gravel loss monitoring stations were established over the entire region in order to assess the gravel loss and the implication of using a better quality of gravel material. This study discusses the gravel loss monitoring approaches, data analyses, and improved material specification for gravel. It is found that the modified gravel used on unsealed road performs better than conventionally used gravel.

#### Edited by:

*Sakdirat Kaewunruen, University of Birmingham, United Kingdom*

#### Reviewed by:

*Qideng Sun, Coffey International, Australia Shunxiang Song, Hohai University, China*

#### \*Correspondence:

*Vasantsingh Pardeshi vasantsinghbhimsingh.pardeshi@ student.uts.edu.au*

#### Specialty section:

*This article was submitted to Transportation and Transit Systems, a section of the journal Frontiers in Built Environment*

> Received: *26 October 2019* Accepted: *10 January 2020* Published: *13 February 2020*

#### Citation:

*Pardeshi V, Nimbalkar S and Khabbaz H (2020) Field Assessment of Gravel Loss on Unsealed Roads in Australia. Front. Built Environ. 6:3. doi: 10.3389/fbuil.2020.00003* Keywords: gravel loss, gradation, prediction model, modified gravel, unsealed road

# INTRODUCTION

Unsealed roads are the economic backbone of Australia that depends on mining, farming, and forestry. The significant concerns across many councils about unsealed roads are related to dust, potholes, and either slippery road after wet weather or loose stones after a long period of dry weather. For about 500,000 km length of unsealed roads in Australia, approximate maintenance cost per year is 1B AUD (de Percy, 2018). Scenic Rim Regional Council's annual maintenance budget on the unsealed road is about 2M AUD per year [Scenic Rim Regional Council (SRRC), 2018].

The problems associated with these roads are gravel loss, shape loss, and rideability. These issues are the result of deterioration factors: inadequate drainage capacity, dust, corrugations, potholes, ruts, loose gravel, and frost damage. Alzubaidi and Rolf (2002) reported the systems for rating unsealed roads in Sweden, Finland, Canada, USA, New Zealand, and Australia. Linard (2010) used the System Dynamics Modeling (SDM)-based pavement management model, which is using the Australian Road Research Board (ARRB) research for pavement deterioration and rehabilitation of unsealed roads. Henning et al. (2008, 2015) developed a framework for the maintenance of unsealed roads. They point out that the authorities did not commonly use prevailing asset management systems due to resource-demanding and misleading of the actual conditions. It was demonstrated that dependence on condition data collection could be limited by the use of the material properties of surface aggregates to predict surface performance. The study used historical maintenance records and visual inspection to predict the rate and mode of failure. In addition, Henning et al. (2008) proposed the inclusion of new performance indicators, shape loss, and slope loss, to predict the gravel loss model. Van Zyl (2011) identified a large number of variables influencing the deterioration of unsealed roads. These variables pose difficulty in the development of a reliable model based on statistics of the available data.

The Paige-Green (1989) concept is widely used in a number of design guidelines including the Austroads Part 6 and the TRH 20 Unsealed Roads Design Construction and Maintenance. Ellis and Andrews (2013) concluded the optimized gravel selection on a network basis. Tim and Choummanivong (2016) established the first countrywide Australian local road deterioration model with collaborative effect of 236 road agencies. Those models were developed by analyzing the results of a long-term monitoring program covering 500 sealed and 100 unsealed local road sites in various traffic and climatic environment in Australia. Uys (2011) selected various internationally recognized gravel loss deterioration models to compare the predicted gravel loss with the actual measurements. Tim and Choummanivong (2016) proposed Road Deterioration (RD) and Work Effects (WE) models. The Local Road Deterioration Study (LRDS) commenced during 2000 has produced RD models that can be readily incorporated into Pavement Management System (PMS). Tim and Choummanivong (2016) considered RD models catering for roughness, gravel loss of surface, and shape loss of pavement as well as interim WE model considering the impact of light and medium grader blading and granular resheeting.

In Australia and New Zealand, the outcome from models and research so far is the development of the "Unsealed Roads Manual, Guidelines to Good Practice," 3rd edition (Giummarra, 2009), which covers management procedures and practices for both the Australian and New Zealand unsealed road network. ARRB has clearly indicated that unsealed roads are presently managed with a few technical input. Due to lack of technical input, the full benefit is not yet realized from the available funding. The "Unsealed Roads Tactical Asset Management Manual" (Henning, 2015) developed in New Zealand is another valuable document. Ellis and Andrews (2013) used the "Unsealed Roads Tactical Asset Management Manual" as a beginning argument and developed a comprehensive strategic forecasting tool for the management of unsealed road networks. This manual emphasizes substantial dependence on rigorous data collection.

Most of the studies reported the complexity in accurate prediction of deterioration models for unsealed roads due to a large number of variables involved. Managing unsealed roads frequently encompasses operational issues, because unsealed roads change quickly and when defects are evident, they must be resolved within a short response period. Due to this, the most routine and cyclic maintenance are programmed according to routine inspections and experience from road operators. However, long-term maintenance activities such as regraveling and surfacing of these roads need a more sophisticated process that includes predictive models (Paige-Green, 2007). A most important consideration during these analyses involves the economic assessment of diverse maintenance options and timings of intervention. A typical Unsealed Road Management issue is identified as (i) difficulty to forecast behavior, (ii) significant data collection needs, and (iii) variability in the level of service and maintenance practices. Based on the review of available literature, research gaps are identified as follows: (i) requirement to develop better material specification by considering local conditions, (ii) developing Martin's model further for localized condition, (iii) including effects of proper maintenance and effects of blading after the use of better specified material, and (iv) establishing better correlation between gravel loss and rough meter based on modified gravel specification.

# STUDY AREA AND FACTORS INFLUENCING GL

Scenic Rim Regional Council (SRRC) maintains an extensive road network of sealed and unsealed roads. The Council provides a road network of 1,816 km, which consists of 53% of the unsealed road networks [(Scenic Rim Regional Council (SRRC), 2018)]. On 30 and 31st March 2017, the rainfall produced by Cyclone Debbie and the cold front meeting over the SRRC ranged from 350 mm in the West of the Scenic Rim region to 800 mm in the East of the Scenic Rim region in a 24 h period. The annual average rainfall for Scenic Rim was 892 mm (Cryna weather station). The 24 h flood event was approximately equal to the annual rainfall. Road drainage is not designed to withstand this level of rainfall in such a short period of time. **Figure 1** shows the Scenic Rim Regional Council and adjoining council area.

SRRC has a planned maintenance schedule resulting in a fairly well-maintained gravel road network. Prior to the cyclone Debbie event, SRRC contracted a private company, IMG to rate all the roads in the region. The majority of the unsealed roads have been rated as three (better) on a 1–5 scale, where 1 is excellent. In many instances, the damage is not immediately apparent as there is no evidence of destructive damage (washes, slips, major erosion). The damage is in the loss of surface material across the entire road surface, loss of shape, loss of fines, and washes in wheel tracks. The volume of water on the roads in a short period has resulted in surface erosion of almost all unsealed roads to some extent. By using the recent information prior to the flood and the information from data collected, the range of restoration treatments proposed are a minimum of a medium grade, a heavy grade with a 50 or 75 mm top-up, and a heavy grade with a 100 mm top or 100 mm resheet and a full 150 mm resheet.

There are 986 unsealed road sections, ∼861 km in length, in the Scenic Rim region. Due to a large amount of gravel road involved and many roads now getting resheeted, the SRRC has initiated gravel road-related research project. **Figure 2** represent locations of Gravel Loss monitoring stations. The aim of this project is to enhance the existing gravel material specification, measure gravel loss, and calibrate the existing gravel loss models. Based on the developed gravel road maintenance strategy, it is clear that the improved gravel loss models are suitable to the SRRC area.

Various past studies on gravel road performance have argued that most existing international gravel loss prediction models are unable to capture the intangible characteristics of the locality. Also, the proper handling with inconsistency in gravel materials and climate is not captured. Uys (2011) reported that gravel loss prediction models of South Africa were neither accurate nor transferable due to their local influences. The international

FIGURE 1 | Scenic Rim and surrounding Council in South East Council, Queensland (map reproduced from the Department of Local Government, Queensland website).

prediction models assessed inaccurate gravel loss, re-graveling quantity, and design inputs for estimation of gravel layer thickness, whereas the pavement management theory is based on the precise prediction rate of roadway deterioration.

Various factors such as traffic, weather, material property, the geometry of road, and maintenance practices affect gravel loss. In this study, straight and flat roads are considered (i.e., exclude roads with varying vertical geometry). Three factors, viz., material properties, traffic, and weather condition, are considered and their impacts will also be assessed.

# Gravel Loss Prediction Models

The primary maintenance activities of unsealed roads are grading and resheeting. The frequencies at which grading and resheeting should be applied depending on the economic trade-off between the costs of the grading/resheeting and the benefit to be gained from reducing road-user costs (Paterson, 1991). The loss of gravel material is caused by natural weathering and friction/whip-off from vehicles (Paige-Green, 1989). Most of the gravel loss forecast models enumerated that the rate of gravel loss is a variable of the traffic, road geometry, material properties, and weather. The forecast of the expected gravel loss from a gravel road is significantly important for unsealed road design and maintenance planning, as graveling and re-graveling operations are quite expansive procedures (Paige-Green, 1989). The selection of gravel loss prediction model is based on the optimal re-graveling schedule for effective maintenance management of unsealed roads.

For a management model, the life cycle of unsealed road deterioration and maintenance can be represented by the extents of material loss over time (Paterson, 1991). The mean gravel surfacing material design thickness tends to be reduced progressively under the action of traffic and climate-related deterioration. This reduction in thickness, by keeping other variables constant, may vary with material characteristics and the effectiveness of maintenance practice, until it reaches the minimum depth that shall trigger re-graveling, and then the cycle repeats itself (Paterson, 1991).

Although it has limitations, the gravel loss measurement is among the tools available for implementing maintenance management goals in terms of depletion rate of gravel loss. Therefore, knowledge regarding rate of gravel loss is crucial to an accurate and ideal decision on its re-graveling cycle (Mwaipungu and Allopi, 2014). Bannour et al. (2019) concluded that HDM-4 pavement deterioration models are insufficient to develop an optimal maintenance management strategy for road networks in Morocco.

Pardeshi et al. Gravel Loss on Unsealed Roads

The predictions from four different gravel loss models are compared with actual measured gravel loss in the field (see **Figure 3**). These models include (i) the HDM-4 gravel loss deterioration model, (ii) the TRH 20 gravel loss deterioration model, (iii) the Australian GL model (i.e., ARRB model) proposed by Martin et al. (2013), and (iv) the Brazilian gravel loss deterioration model. These four models are more relevant to this study due to the similarity in variables and resemblance to SRRC conditions. To obtain the predicted GL, the input parameters for traffic data, rainfall, and material quality (plasticity and grading values) used are the same for all the models. Percentage differences were calculated for selected stations due to the availability of level data for those stations. The median percentage difference in GL for selected stations was 2% between the actual and ARRB models. There were significant differences between the other three models and actual GL median values. The ARRB model is used for further analysis as it is with least differences with actual GL and is more suitable for Australian conditions. More details of these models are listed below.

#### HDM-4 Gravel Loss Deterioration Model

$$\begin{array}{rcl} \text{AGL } = & k\_{\text{g}l} \times 3.65 \left( 3.46 + 2.46 \times \text{MMP} \times \text{RF} \times 10^{-4} \right) \\ &+ \text{KT} \times \text{AADT} \end{array} \tag{1}$$

where AGL is the predicted annual material loss (mm/year), RF is average rise and fall of the road (m/km), MMP is mean monthly precipitation (mm/month), AADT is annual average daily traffic (veh/day), KT is the traffic-induced material whiff-off coefficient, and kgl is gravel material loss calibration factor.

To derive the GL values using this model kgl, KT was derived based on actual traffic volume, material quality, and weather data. Also, roads with no curvature were selected for comparison. The gravel loss is calculated using HDM-4 gravel loss deterioration model as:

$$\mathrm{GL} = k\_{\mathrm{gl}} \mathrm{D} \left( 3.46 + 2.46 \times \mathrm{MMP} \times \mathrm{RF} \times 10^{-4} + \mathrm{KT} \times \mathrm{AADT} \right) \tag{2}$$

where GL is the gravel loss (mm) and D is time period under consideration in hundreds of days (days/100).

#### TRH 20 Model

$$\text{AGL} = \text{3.65(ADT} \times \text{(0.059 + 0.0027N - 00006P\_{26}) - 0.00474P\_{26})} \tag{3}$$

where ADT is the average daily traffic (total number of vehicles/day), N is Weinert N value, P<sup>26</sup> is the percentage by mass of material passing a 26.5 mm sieve, and PF is the product of plastic limit and percentage passing a 0.075 mm sieve.

To derive the GL values using this model, N was assumed 4 based on Australian Climatic conditions. The value of gravel loss is calculated from the following equation:

$$\text{GL} = \text{D(ADT} \times \text{(0.059 + 0.0027N - 00006P\_{26}) - 0.00474P\_{26})} \tag{4}$$

## ARRB Model

$$\text{GL} = k\_{\text{jl}} \text{D} (-0.00985 \text{ADT} - 0.02991 \text{MMP} - 0.00583 \text{PF}) \text{ (5)}$$

where MMP is the mean monthly precipitation (mm/month), PF is the plasticity factor (= PI × P0.075), P0.075 is the percentage by mass of material passing a 26.5 mm sieve, and PI is the plasticity index. kgl was derived based on the measured actual gravel loss and the median value was used, even though the ARRB model uses negative coefficients in the equation, which returns net GL as a negative number. In this study, while deriving the GL numbers, kgl is assumed negative so that the GL values become positive. It can be compared with other GL model values.

#### Brazilian Gravel Loss Deterioration Model

$$\text{GL} = \text{D} \left[ 1.58 + 0.366 \text{(G)} + 0.083 \text{(SV)} - 0.210 \text{(PI)} \right]$$

$$+ 0.0132 \text{(NC)} + 0.008 \text{(NT)} + \left\langle \frac{420.45}{R} \right\rangle \text{\(} \tag{6}$$

where G is the absolute value of grade (percent), SV is the percentage of the surfacing material passing the 0.075 mm sieve, NC is an average daily car and pick up traffic in both directions, NT is average daily truck traffic in both directions, and R is the radius of horizontal curvature (m). To derive the GL values using this model, a road with no curvature was selected for comparison.

There are significant differences with actual GL measurement and the predicted GL. Model values are presented in **Figure 3** for the selected stations.

McManus (1994) concluded that deterioration models were produced from substantial studies conducted in other countries. Each individual model reflected the characteristics of the particular country for which they were developed. Before the development of the Australian individual GL model, the models were based on the performance of overseas pavements (Giummarra et al., 2007). Therefore, they were not able to predict pavement performance precisely under Australian conditions. These caused local government authorities across Australia to have difficulties with their Pavement Management Session (PMS) package. In this regard, it is emphasized for the importance of having gravel road performance prediction models that reflect the local characteristics. The performance of prediction models must imitate the local conditions and should be developed from local data or substantiate based on these data (Mwaipungu, 2015).

As explained earlier, there are large differences between the actual measured GL and calculated by the ARRB model. Thus, further study is needed considering three parameters and their effects with maintenance activities. The Scenic Rim region has different climatic conditions, varying subgrade soils, terrain, and marginal gravel materials that affect gravel road performance. As a result, this requires the parameters identified above to be reflected in any gravel loss prediction model for it to be applicable in the region. The following sections describe the effects of GL model variables based on actual GL measurement carried out at the Scenic Rim region in the state of Queensland during 2018.

This study intends to use traffic survey data and analyze traffic variables in further detail. Traffic variables such as the percentage of heavy vehicles and speed need to be considered. The speed environment also affects the distribution of loose gravel aggregate. In this study, the Thornthwaite Moisture Index (TMI) is proposed to account for weather-related variables. TMI considers temperature and rainfall, which are proposed in this study. The Scenic Rim region has three major types of areas based on weather characteristics: dry, arid, and dry temperate. The unsealed roads behave differently according to local weather characteristics. There is also an extreme weather pattern developing all over Australia. The Scenic Rim region has experienced very wet weather during 2017 and very dry weather during 2018. The extreme weather patterns need to be accounted for and the GL model needs to be dynamic enough to account for those extreme weather patterns. Provided those weather patterns are considered in the model, the gravel maintenance program will be adjusted based on GL prediction. The inclusion of TMI becomes a necessity.

The material variability can be minimized by specifying material quality. In this study, it was prioritized to modify material specifications. Based on the field experience, some changes are proposed for the gravel, which is used for the regraveling and resheeting work. The results are looking promising and discussed in the following section.

# PROPOSED MODIFICATION TO GRAVEL MATERIAL SPECIFICATION

# Background on Gravel Road Material

The ARRB manual (Unsealed Roads Manual: Guidelines to Good Practice, 3rd edition, March 2009) is an excellent source of information for the construction and maintenance of unsealed roads. It contains a number of specifications for gravel to be used on unsealed roads (section proposed modification to gravel material specification). Many of these specifications have been developed over many years in different parts of the world and proven to provide good gravels.

The damage to unsealed roads caused by Ex-Tropical Cyclone Debbie in the SRRC area amounted to 500,000 tons. SRRC decided to rebuild the gravel road network using the ARRB gravel specifications. The specification used is the combined wearing course/base course material described in section 3.5.2 of the ARRB Unsealed Roads Manual (Giummarra, 2009).

Cocks et al. (2015) used the Paige–Green grading coefficient and the shrinkage product concept to improve limits on mine haul roads in Western Australia. It is important to understand the intent of the specification that is described in the preamble of the ARRB specification. The intent is to create a practical specification not constrained by tight limits such as a narrow range of PI and a tight grading curve. It is recommended that the specification is adapted to the available material and to the circumstances. The basis of the ARRB specification is that the road needs structure in terms of the grading coefficient and binder in terms of the shrinkage limit. The goal is to find a balance between these two to provide a material that performs in all weather and is easily constructible. The ARRB nomograph provides a means to obtain this balance.

The SRRC experience is detailed below and hopefully provides a deeper understanding of the ARRB specification. For ease of reference, we have called the unsealed gravel material "Type 4.5."

# The Experience on Gravel Material and Proposed Modifications for Gravel Material Specification

Initially, SRRC worked with a single local quarry. The quarry produced the material giving results within the envelope on the ARRB nomograph. The SRRC decided to allow the full range of the shrinkage product from 100 to 350. The Shrinkage Product was within the desired range (around 180). The grading coefficient was at the midpoint of allowable range and the CBR averaged as 30. A production run was made following the testing. The test results confirmed that the material can be used for road construction. This material was trialed on three roads and performed beyond expectations. It had good constructability and excellent wearing ability. The wearing surface provided a smooth ride, required little maintenance, and created little dust. Tenders were called for to create a panel of gravel suppliers using the ARRB specification.

The next production run from the same quarry tested within specification but much higher on the shrinkage product, reaching 300, although, with the same specification, the material was not easy to use. The slightest rain caused a very slippery surface and vehicle bogging. Over the next few weeks, the in-specification material continued to be hit and miss. Many tests were conducted to overcome this issue. However, there is no clear indicator about the inconsistency.

To fill this gap, it is crucial to understand the nature of the gravel. The material has two main parts: (i) the grading must be such that there is a good structure providing interlocking and support, and (ii) the material fines (percentage passing 0.425) must have the ability to bind the structure without being overly plastic. The bonding with these parts is sufficient to create the good gravel material.

Approximately 80 material quality tests were conducted from numerous stockpiles. All the test results were listed and the material was rated by constructability and service. There was no immediate pattern as to why some materials were good and others were poor. In particular, there was a tendency to produce material that reduced to mud at the slightest wetting. From visual inspections, supervisor comments, and the tests, it was realized that the muddy material had little structure, although it was within the limits of the ARRB grading coefficient specification.

After some analysis, it revealed that the better-performing materials have a good structure defined by the percentage passing 19 mm minus the percentage passing 4.75 mm (i.e., P19.0 – P4.75). The value of better material is between 35 and 50, and 40 is performed as best. This "structure" factor is not the same as the grading coefficient. However, this still does not fully coincide with the experienced performance. It was then decided to investigate the binder or clay portion. After some analysis, a correlation emerged showing that the term (PI × P0.075) was significant. The plasticity index (PI) is the size of the range of water contents where the soil exhibits plastic properties. The PI is the difference between the liquid limit and the plastic limit. It represents the amount of plasticity in the material. By dividing this by the structure calculation above, a definite correlation was found for wet performance. If this value is <5, the material was easy to construct and supported the traffic without rutting or tracking. Between 5 and 10, the material is sticky, more difficult to work with. Above 10, the material shows little support and rutting is evident immediately. This factor is termed as the Soak Performance Indicator (SPI).

$$SPI = PI \times \frac{P\_{0.075}}{P\_{19.0} - P\_{4.75}} \tag{7}$$

where P0.075, P4.75, and P19.0 are the percentages passing through 0.075, 4.75, and 19 mm sieves, respectively. Note that the WtPI is an indicator developed by SRRC as an additional indicator for the suitability of material with high shrinkage product. There is not sufficient evidence to support that it will provide accurate indications for all materials.

Using this factor, a clear correlation with the performance was apparent. The WtPI ascertains the importance between the ratio of binder and the structure. The ARRB specification allows a maximum shrinkage product of 350, and prefers 240. The results showed that 240 is the limit and even that is not suitable in continuously moist areas such as roads under trees on the south side of a hill. The relationship between the binder and the structure is examined in greater detail. An adequate binder is achieved with low PI material. The low linear shrinkage (normally 50% of the PI) and more fines (0.425 and less) are required to create a material that provides an adequate binder. However, if the percentage fines become too high, the structure deteriorates and the material is too sandy. From field experience, this starts to happen when the percentage passing 0.425 mm approaches 30%. This roughly equates to a PI of 7 (LS of 3.5) to get the minimum 100 shrinkage product. From this, our specification starts with the PI at 8.

In contrast, when the PI is high, the amount of fines must be reduced. There is a limit to how little fines are required; if the fines are not enough, the binding action does not take place. For the modified ARRB nomograph, a maximum of 220 for the shrinkage product is used. The lower limit of workability for 0.425 mm is 15%. Below this, the material is too bony. This equates to a shrinkage limit of around 15 or a PI of around 30. Working with a PI of 30 is not easy because as the material sticks to the machinery, it tends to lump the fines together, and the clay prevents the water from penetrating. From experience on site, it seems that a PI of around 25 is the limit for workability. This equates to a shrinkage limit of around 12. To remain within specification, the maximum fines (0.425 mm) are around 18%. As 15% is the minimum for workability, this small range is not possible to work with. It is almost impossible for quarries to produce this fine tolerance. For this reason, the PI is limited to 20. The shrinkage limit (SL) is around 10 and passing percentage through a sieve size of 0.425 is 22%. This leads to the limits of the size 0.425 mm being 15–22%. We have permitted extra finer material by shifting 15–25%. Further investigations are underway to analyze material if it is suitable outside this range.

As the PI is fixed to the material available, it is important for the quarry to vary the fines and the grading curve to match the PI. Although the other factors like the liquid limit (LL) and SL are important, they are less sensitive for this material, provided the PI, percentage fines, and structure (19–4.75) are within the required limits.

The fines ratio should be around 0.65, but there is no evidence that departure from this causes much of a problem. Note that the fines ratio for type 4.5 is typically higher than most other materials. On the ARRB nomograph, the climate must also be taken into account. In wetter areas, it is advisable to aim lower toward 100, and in drier areas, it is recommended to aim toward 220. Refer **Figure 4** which shows relationship between modified shrinkage product limits and grading coefficient. We found that due to variations of material between quarries and

Frontiers in Built Environment | www.frontiersin.org

also variations of materials within quarries, the specification needs to be adaptable to the specific material. **Figure 5** shows the relationship between shrinkage product (SP) and grading coefficient (GC) on the performance of the base/wearing course.

$$SP = SL \times P\_{0.425} \tag{8}$$

where SL is shrinkage limit and P0.425 is the % passing the 0.425-mm sieve.

$$GC = \frac{(P\_{26.0} - P\_{2.0}) \times P\_{4.75}}{100} \tag{9}$$

where P26.0, P4.75, and P2.0 are the percentages passing 26, 4.75, and 2 mm, respectively. The blue dots in **Figure 5** indicate that the gravel used for reconstructing unsealed roads within the SRRC area is conforming material. A total of 696 number of test results are plotted in this graph and the majority are within "Good" area E of this graph. The original "shrinkage product and grading coefficient" monograph is modified. The zone V, which is noted as "Good," is split into two zones, Zone V and Zone VI. Zone V is noted as "Good" and Zone VI is noted as "Acceptance in certain conditions." Zone V and VI are split by using a shrinkage product line of 220. Material with shrinkage product higher than 220 falls under Zone VI and below 220.

# Field GL Measurement Outcome on Existing Gravel Material and Proposed Modified Gravel Material

Gravel loss monitoring stations are established across the region to monitor and compare gravel loss. **Figure 6** is a photo of a typical GL monitoring station installed adjacent to a road section. This station is used to record levels over period of time. The permanent control mark was established by driving 1.8 m star picket deep into the ground. Two-star pickets were driven on both sides of the road. Precast concrete surround was placed to make the top surface even with surrounding ground level. Due to safety reasons, a plastic cap was placed on the concrete surface.

Twenty road stations with existing material (EM) plus modified material (MM) were compared with 10 stations of EM to compare gravel loss across a section of road. For all stations, ADT varies from 80 to 1,000 and MMP is varied from 50 to 117 mm. Material property noted as PF is the product of plasticity index and material passing through a 26.5 mm sieve. PF of new material varied from 56 to 268. The new material used is as per modified specification and referred to as modified material discussed.

**Figure 7** presents the variation in gravel loss (GL) in mm over a number of days (t). On the existing gravel material road sections, the gravel loss was at a higher rate for a shorter duration. For a fair comparison, GL at each station is converted for 365 days (1 year) by extrapolating the recorded gravel loss over specific measured days. On EM, maximum gravel loss was 214.13 mm. The minimum gravel loss was 7.64 mm within 239 days on one location, which is an exception. On the new gravel material road section, the GL was at a lower rate than old gravel and for a longer duration. The maximum GL was 44.37 mm during 365 days and even a marginal GL of 1 mm was recorded over 230 days. The median value of GL is 48.66 mm on EM for 1 year and 7.93 mm on MM.

The GL rate is significantly high for old material, which is not modified material. However, no significant maintenance was

FIGURE 6 | A typical GL monitoring station used on the field.

required for the improved gravel road. There was a lesser need for regrading the roads due to this improved gravel road due to the use of modified gravel material. The field crew confirmed that a gravel road with improved material performed well. Generally, these roads need to be graded within 1 year. Due to the modified material, there is no need to grade for almost 18 months. This demonstrates the benefits of good gravel material. The modified material as per **Figure 7** shows consistent results of reduced GL.

# CONCLUSIONS

There is a huge potential for cost savings by reducing gravel loss on unsealed roads. The majority of the gravel loss rate is between 6 and 10 mm per year. Many studies conclude difficulty predicting deterioration models very accurately for unsealed roads due to a number of variables involved and the requirements of time and resources. Even though some of the models have predicted gravel loss, there is a need to develop a gravel loss model suitable for local conditions.

In this study, a number of gravel loss monitoring stations are established and actual gravel loss is measured. Actual measured gravel loss is compared with the predicted gravel loss from existing models. Based on the field data, the following specific conclusions could be drawn:

a) The current models are not able to capture gravel loss accurately and further data are needed to calibrate parameters

# REFERENCES


resembling localized conditions. While completing the resheeting and re-graveling, the gravel material specification is amended and improved.


# DATA AVAILABILITY STATEMENT

The datasets generated for this study are available on request to the corresponding author.

# AUTHOR CONTRIBUTIONS

VP prepared the draft. SN and HK edited the manuscript.

# ACKNOWLEDGMENTS

The authors would like to acknowledge Chris Gray, General Manager, Infrastructure Services from the Scenic Rim Regional Council, for providing support for this research project.


**Conflict of Interest:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2020 Pardeshi, Nimbalkar and Khabbaz. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Liquefaction Characteristics of Sand With Polypropylene Fiber at Low Confining Stress

Sudhir K. Tripathi, Prabir K. Kolay\*, Vijay K. Puri and Sanjeev Kumar

*Department of Civil and Environmental Engineering, Southern Illinois University Carbondale, Carbondale, IL, United States*

The present study evaluates the liquefaction characteristics of Ottawa sand with polypropylene fiber using cyclic triaxial test. A series of stress-controlled cyclic triaxial tests were performed at 34.47 kPa (5 psi) effective confining pressure. Specimens of clean Ottawa sand, and sand containing 0.05, 0.075, 0.1, and 0.3% polypropylene fiber by dry weight of sand were tested at 0.2, 0.25, 0.3, and 0.4 Cyclic Stress Ratio (CSR). Results show a significant improvement in liquefaction resistance when the polypropylene fiber content exceeded beyond 0.075% at 34.47 kPa effective confining stress.

Keywords: sand, polypropylene fiber, cyclic triaxial, liquefaction, cyclic stress ratio

#### Edited by:

*Sakdirat Kaewunruen, University of Birmingham, United Kingdom*

#### Reviewed by:

*Sujit Kumar Dash, Indian Institute of Technology Kharagpur, India Nhlanganiso Mlilo, Coventry University, United Kingdom*

> \*Correspondence: *Prabir K. Kolay pkolay@siu.edu*

#### Specialty section:

*This article was submitted to Transportation and Transit Systems, a section of the journal Frontiers in Built Environment*

> Received: *14 June 2019* Accepted: *28 February 2020* Published: *17 March 2020*

#### Citation:

*Tripathi SK, Kolay PK, Puri VK and Kumar S (2020) Liquefaction Characteristics of Sand With Polypropylene Fiber at Low Confining Stress. Front. Built Environ. 6:28. doi: 10.3389/fbuil.2020.00028*

# INTRODUCTION

Generally, earthquake damages are depicted based on structural and non-structural failures of superstructures and substructures; but the underlying ground on which a structure is founded has a fundamental role in extents of damage. One of the effects of earthquake to the underlying ground is liquefaction of soil. Liquefaction of saturated soil deposit is an important phenomenon, which governs the behavior of ground and the structure overlying on it. Liquefaction occurs when a saturated soil deposit undergoes significant reduction in strength and bearing capacity when pore pressure rises and equals effective confining pressure due to undrained loading condition created by transient strong shaking during earthquake. Liquefaction phenomenon received major attention after 9.2 Richter scale magnitude earthquake in Nigata, Japan (1964) and 7.5 Richter scale magnitude earthquake in Anchorage, Alaska (1964). In both earthquakes, liquefaction induced damages such as slope failures, bridge and building foundation failures, and floatation of structures were observed (Dobry et al., 1982). In addition, Kramer (1996) also observed liquefaction induced excessive settlement and tilting of buildings. The initial laboratory studies during 1960–80 mainly focused liquefaction of clean sands. However, later, Standard Penetration Tests (SPT), Cone Penetration Test (CPT) to evaluate liquefaction potential revealed that liquefaction may also occur in sands containing fines. This fact was supported by the liquefaction damages that occurred during Central Chile earthquake (1985) (de Alba et al., 1988); Saguenay earthquake (1988) (Tuttle et al., 1990), Superstition hills earthquake, California (1987) (Holzer et al., 1989).

The liquefaction potential is influenced by many parameters, which directly or indirectly affect the behavior of soil, viz. degree of saturation, relative density, earthquake magnitude, ground motion characteristics, fines content, and effective overburden pressure. Several researchers e.g., Seed and Lee (1966), Castro (1969), Youd and Hoose (1977), Ishihara (1984), and Evans and Seed (1987) have investigated the effect of those and several other parameters. Lee and Aurelio (1974), Tatsuoka et al. (1984), Tokimatsu and Seed (1987) have studied earthquake induced settlements and found that significant hazard may result from the dissipation of developed pore pressure after liquefaction. Meanwhile, several research works are focused on measures to prevent or reduce liquefaction potential. One of the techniques that has gained significant interest is the use of geofibers to improve potentially liquefiable soils. While some research suggests that the use of geofiber in soil improves the static characteristics of sands, and sands with fines, some research works suggest otherwise. On the dynamic characteristics part, Noorany and Uzdavines (1989), and Maher and Woods (1990) observed an increase in soil shear strength and cyclic resistance with the addition of geofibers.

The main purpose of this research is to study the effect of polypropylene fiber on liquefaction characteristics of sand at low confining stress. This study presents the results obtained from cyclic triaxial tests carried out on mixtures of sand with various percentages of polypropylene fibers. The sand-fiber specimens contained 0.05, 0.075, 0.1, and 0.3% of fiber by dry weight of sand. The results from cyclic triaxial tests are evaluated in terms of deformation and pore pressure generation characteristics with respect to application of cyclic loads. The Cyclic Stress Ratio (CSR) for the tests varied from 0.2 to 0.4.

# MATERIALS AND METHODS

# Materials

The materials used for the cyclic triaxial test specimens are Ottawa 20–30 ASTM C 778 sand and Monofilament Polypropylene fibers. Ottawa sand is naturally occurring silica sand passing through sieve no. 20 (0.84 mm) and retained on sieve no. 30 (0.60 mm) supplied by Humboldt. Quanta chrome gas pycnometer was used to determine the specific gravity of sand tested in accordance with ASTM D5550-14 (2014) and the value was found to be 2.69. The properties of sand are presented in **Table 1.**

The fiber used in this research is Virgin Homopolymer polypropylene monofilament supplied by FORTA corporation, USA. Each fiber is 6 mm in length. These fibers are surface treated with special chemical, which resist moisture and aging. Because of their durability, polypropylene fibers have gained popularity in geotechnical application. The properties of fiber are listed in **Table 2**.

# Methods

The maximum and minimum void ratio were determined in accordance with ASTM D4254-16 (2016) and ASTM D4253-16e1


(2016), respectively. The maximum and minimum void ratio of sand was found to be 0.706 and 0.513, respectively. ASTM method 1a was used to perform the test. ASTM D422-63 (2007)e2 (2007) was used to perform sieve analysis of Ottawa sand; and, the grain size distribution was plotted. The sand is classified as poorly graded sand (SP) under Unified Soil Classification System (USCS). The values of D10, D30, and D<sup>60</sup> were obtained from the particle size analysis graph as 0.48, 0.50, 0.54 mm, respectively. The value of coefficient of uniformity, Cu, and coefficient of curvature C<sup>c</sup> were found to be 1.125 and 0.965, respectively.

# Sample Preparation

The sand-fiber mixture was prepared by selecting various percentages of fiber by dry weight of sand. The percentage of fibers that were tested in this study were 0, 0.05, 0.075, 0.1, and 0.3%. The fibers were mixed in sand by dry mixing. The fiber strands were separated manually by hand and mixed in sand in small proportions to ensure uniform distribution of fibers within sand. Proper care was taken to avoid segregation of fibers in sand-fiber mixture.

Quanta chrome gas pycnometer was used to determine specific gravity of sand fiber mixture based on ASTM D5550- 14 (2014) standard. The tests were repeated until the specific gravity was within standard deviation of 0.005%. The decrease in specific gravity of sand-fiber mix was observed with the increase in percentage of fiber. **Table 3** presents the specific gravity of various percentages of fibers mixed with sand.

Cylindrical samples for cyclic triaxial test had length and diameter of 5.7′′ and 2.85′′, respectively. The samples were prepared at the relative density of 30% and degree of saturation of 30%. Dry and wet unit weight of samples were calculated by using void ratio, specific gravity, water content, and degree of saturation. The sample was compacted in six layers. Method of undercompaction suggested by Ladd (1978) was used to prepare specimens at constant relative density of 30%.

TABLE 2 | Properties of polypropylene fibers.


TABLE 3 | Specific gravity of sand with various percentages of fiber mixtures.


The height of each layer was then calculated based on the percentage of undercompaction that was needed. The percent of undercompaction for each layer was as follows:

Percent under compaction for nth layer can be calculated using Equation (1) (Ladd, 1978).

$$U\_n = U\_{ni} - \frac{U\_{ni} - U\_{nt}}{n\_t - 1} \* (n - 1) \tag{1}$$

Where, Uni = Percent under compaction for the first layer

n<sup>t</sup> = total number of layer

n = number of layer being considered.

The height of compacted layer is calculated using Equation (2) (Ladd, 1978). The compacted height for bottom most layer was 26.16 and 23.62 mm for all other layers

$$h\_n = \frac{h\_\ell}{n\_\ell}(n-1) + \left(1 + \frac{U\_n}{100}\right) \tag{2}$$

Where, h<sup>t</sup> = total height of sample

h<sup>n</sup> = height of specimen at the top of the nth layer.

# TEST RESULTS

The results and findings of cyclic triaxial test with fiber reinforced sand specimens are presented in this section. First, the details of tests conducted during study are presented in the form of table. Finally, the results of cyclic triaxial test on clean sands and sands with various percentages of fiber i.e., 0.05, 0.075, 0.1, and 0.3% at 34.48 kPa confining pressure have been discussed.

# Test Details

In this study, numerous cyclic triaxial tests were carried out in clean sand, and sand-fiber mixtures containing 0.05, 0.075, 0.1, and 0.3% fiber by dry weight of sand. All the samples were prepared at 30% relative density. The samples were tested at 34.48 kPa confining pressure and CSR was varied from 0.2 to 0.4. **Table 4** presents all the tests that had been conducted in this study.

# Effect of Fiber on Liquefaction Resistance

To evaluate the effect of inclusion of fiber in Ottawa sand with respect to liquefaction, several tests were conducted on clean sand, and sand containing 0.05, 0.075, 0.1, and 0.3% fiber, respectively. The tests were conducted at 34.48 kPa confining pressure. The initial liquefaction is based on the pore pressure generation criteria i.e., initial liquefaction is considered to have taken place when pore pressure equals the effective confining pressure. In terms of deformation, number of cycles required to trigger 2.5 and 5% of DA strain have been taken into account to compare resistance against liquefaction.

**Figure 1** presents a plot of cyclic stress ratio against loading cycles required for initial liquefaction for clean sand and sand containing various percentages of fiber. It can be observed that the initial liquefaction resistance tends to increase and then decrease with the addition of 0.05 and 0.075% fiber, respectively; beyond 0.075% fiber content, initial liquefaction resistance tends to increase with further increase TABLE 4 | Cyclic triaxial tests on sand-fiber mixtures at initial effective confining pressure of 34.48 kPa.


of fiber content. This may be ascribed to the fact that the voids in the sand-fiber matrix decrease with the increase of fibers which leads to some densification of sand-fiber matrix. Numerically, for initial liquefaction of clean sand, sand containing 0.05, 0.075, 0.1, and 0.3% fiber at 50 loading cycles, it would require CSR values of 0.31, 0.33, 0.30, 0.323, and 0.385, respectively.

Similarly, **Figure 2** presents the graphs of CSR against loading cycles for 2.5% strain. It can be observed that when the fiber content increases beyond 0.075%, the 2.5% DA lines plot higher than 2.5% DA line at 0.05% fiber content. Numerically, CSR values of 0.3123, 0.33, 0.3267, 0.323, and 0.385 would be required

to trigger 2.5% DA strain at 50 loading cycles in clean sands, and sands containing 0.05, 0.075, 0.1, and 0.3% fiber.

Finally, **Figure 3** presents the plot of cyclic stress ratio (CSR) against loading cycles for 5% DA strain. From the **Figure 3**, it can be observed that CSR required to cause 5% DA strain tends to increase and then decrease with addition of fiber up to 0.075% fiber, and consequently, required CSR increases again beyond fiber content of 0.075%. Numerically, it requires CSR values of 0.312, 0.33, 0.3, 0.324, and 0.385, to trigger 5% DA strain 50 cycles for clean sands and sand containing 0.05, 0.075, 0.1, and 0.3% fiber. The addition of fiber may provide reinforcing effect on sand which is the main reason for increase liquefaction resistance.

# Performance of Fiber

**Figure 4** presents the plot of CSR against percentage fiber content for initial liquefaction at 10 cycles at 34.48 kPa confining pressure. It can be observed that for 34.48 kPa effective confining pressure the CSR required for initial liquefaction

at 10 cycles increases with increment in percentage fiber by dry weight of sand. With further increment of fiber content, the CSR required for initial liquefaction at 0.1 and 0.3% fiber content is 0.4 and 0.458 at 34.48 kPa effective confining pressure.

**Figure 5** presents the plot of CSR against percentage fiber content for 2.5% DA strain at 10 cycles at 34.48 kPa confining pressure. It can be observed that for 34.48 kPa effective confining pressure the CSR required for 2.5% DA strain at 10 cycles increases with increase in fiber percentage by dry weight of sand. With further increment in fiber content, the CSR required to get 2.5% DA strain at 0.1 and 0.3% fiber content is 0.4 and 0.457 at 34.48 kPa effective confining pressure.

**Figure 6** presents the plot of CSR against percentage fiber content for 5% DA strain at 10 cycles at 34.48 kPa confining pressure. It can be observed that for 34.48 kPa effective confining pressure the CSR required for 5% DA strain at 10 cycles increases with increment in percentage of fiber by dry weight of sand. With further increment in fiber content, the CSR required to get 5% DA strain at 0.1 and 0.3% fiber content is 0.39 and 0.458 at 34.48 kPa effective confining pressure.

## Summary

In this study, results of cyclic triaxial test performed on sand and sand containing various percentages of fiber were analyzed and

discussed. The specimens comprised clean sand, and sand mixed with 0.05, 0.075, 0.1, and 0.3% fiber by dry weight of sand. To study the effect of cyclic stress ratio (CSR), the specimens were tested at 0.2, 0.25, 0.3, and 0.4 CSR. Also, to identify the effect of confining pressure, specimens were tested at 34.48 kPa initial effective confining pressure. The results for clean sand based on pore pressure generation, 2.5 and 5.0% DA strain were presented first, followed by results of sand specimen containing 0.05, 0.075, 0.1, and 0.3%, respectively. Next, the results discussing the effect of fiber content were presented based on confining stress and cyclic stress criteria.

The results revealed that the liquefaction susceptibility increases with increase in cyclic stress ratio (CSR). For instance, it required fewer loading cycles for all specimens to liquefy as the cyclic stress ratio increased. Next, at 34.48 kPa initial effective confining pressure, the addition of fiber seemed to improve the liquefaction resistance at fiber contents beyond 0.075%. This is evident from the **Figures 1**–**3**, which shows that the liquefaction

## REFERENCES


resistance initially decreases with addition of fiber up to 0.075% fiber content beyond which it significantly improves liquefaction resistance. This finding is consistent with the observation made by Noorzad and Amini (2014) on Babolsar sand.

# CONCLUSIONS

The present study summarizes the results obtained from number of cyclic triaxial tests conducted on 20–30 Ottawa sand with various percentages of monofilament polypropylene fibers. The tests were undertaken with the aim to investigate the liquefaction characteristics of Ottawa sand with the addition of monofilament polypropylene fibers. Based on the series of cyclic triaxial tests on clean sand and sand containing various percentages of fiber the following conclusions can be drawn:


# DATA AVAILABILITY STATEMENT

The datasets generated for this study are available on request to the corresponding author.

# AUTHOR CONTRIBUTIONS

All authors listed have made a substantial, direct and intellectual contribution to the work, and approved it for publication.

Bureau of Standards, US Department of Commerce, US Governmental Printing Office, Washington, DC.


Geoenviron. 116, 1116–1131. doi: 10.1061/(ASCE)0733-9410(1990)116: 7(1116)


**Conflict of Interest:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2020 Tripathi, Kolay, Puri and Kumar. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

digital media

of impactful research

article's readership