Introduction

The structural response in a rigid pavement allows to know how the system reacts to the different solicitations to which it is subjected during its work regime, therefore, its determination is the basis of the design methods. In these pavements two failure criteria are considered, which are taken into account by the design methods; the fatigue failure of the concrete slab and the erosion failure of the foundation. The analysis of these failures is based on obtaining the structural response in terms of tensions and displacements, and checking, using empirical models, if the admissible values at critical points of the slab are satisfied. The models describe the behavior of the system as consequence of the action of repeated loads during the design period.

The models that determine the response in rigid pavements have evolved over time, but always based on the physical principles of the phenomenon, changing only the way to solve it. The first known model was created by (^{Westergaard, 1926}), obtained from the solution of the differential equation that describes the flexion of a thin slab on a Winklerian medium, at that time very complicated to solve, so simplifications were made to the phenomenon, giving solution in three singular points (center, edge and corner of the slab). With the rising of new numerical tools, together with the progress of computers, it is achieved that Westergaard's analysis could be generalized to any point of the slab, as well as a fully understand the system.

The numerical methods are used to characterize the stress-strain states in solid bodies, where the most generalized is the Finite Element Method (FEM). Among the advantages, which extend this method, are, the incorporation of better elaborated calculation tools, as well as the possibility of modeling geometries and materials of different complexity, so it is appropriate for the modeling of JPCP (^{López, 2016}).

In the analysis and design of these pavements, the assumptions made, offer certain limitations in their application. As a consequence, since the beginning of the year 1970, the MEF, which has become a widely used tool for the analysis of rigid pavements, shows some differences between the studies developed by different authors, related to the way of modeling the slab, the interface with the foundation and load transfer mechanisms (dowels). Such differences can be reflected in the most recent works of the authors (^{Bayrak & Hınıslıoğlu, 2016}; ^{Cobarrubias, 2012}; ^{López & Tejeda, 2015}; ^{Roy, Reddy, & Ramachandra, 2013}).

The FE method is a method for solving problems with complex geometries, loading, and material properties numerically. The 2D finite-element modeling had several limitations, they were unable to model more than two layers above the subgrade, and only one slab laying on top of a Winkler during temperature analysis could be considered. With increasing computers speed and power, 3D modeling was adopted by many researchers, which had significant improvement over traditional design methods, including the availability of interface algorithms and thermal modules (^{Zokaei, Fakhri, & Rahiminezhad, 2017}).

Methodology

Finite elements model development

Each part of the model developed for the study of the stresses caused by the loads is described below, composed of several sub-models: geometry, loads, materials, edge conditions, and interactions.

It has been modeled in three dimensions, taking advantage of the facilities provided by the computational tool used, trying to respect the characteristics of the geometry of the rigid pavements. The investigation follows the recommendations, in terms of dimensions, of the American Association of Concrete Pavements ACPA (^{Calo, 2012}).

Geometry

Two variants of the geometric model were considered, one variant without the concrete shoulder and the other with the slab tied to a concrete shoulder by steel bars. In both cases, it is only necessary to consider one lane, since they have a symmetric configuration (Figure 1).

The slabs composed of concrete, are connected in the transversal joints by means of load transfer pins, and are supported on a foundation floor (subgrade), represented in the model as an infinite semi-space.

Dowels: Plain steel element, represented as a cylinder of 450 mm in length and 32mm of diameter, spaced at 300 mm, for a total of 11 dowels in each transverse joint.

Tie bars: Corrugated steel element, represented as a cylinder of 900 mm in length and 13mm of diameter, spaced at 700 mm, six bars are placed per slab, to be connected with the concrete shoulder in the longitudinal direction.

Concrete slab and shoulder: Structures in parallelepiped shape, with cavities in the dowels positions, according to their diameters and length. The slabs will be 4.50 meters of length and 3.50 of width, with a thickness that can vary in each case. The shoulders will be equal to the slab in thickness and length but with 1.50 m of width.

Soil foundation: Hexahedral structure that represents the soil mass, as an infinite semi-space, whose dimensions are 13.50 x 5.00 x 3.00 m.

Material model

The constitutive model used for the slab-foundation-dowels system was selected taking into account the nature of the phenomenon. The speed of application of the loads and their magnitude, added to the high stiffness of the slab, cause that the stresses originated, remain below the elastic limits of the materials, thus discarding the phenomenon of plasticization, therefore, considered the modeling of the pavement system as a linear, elastic and isotropic medium, characterized by the proportionality between stresses and unitary deformations, where the constants of proportionality are Young's Modulus and Poisson's ratio (E, ν).

Load representation

Two types of axles were contemplated, one simple and the other tandem. The single axle, with dual wheels, weighing 100 kN (single reference axle in Cuban regulations), and the tandem axle weighing 180 kN, obtained as equivalent to the simple reference axle. The load was considered static, which constitutes the most unfavorable conditions in the calculation of stresses, using the equivalent contact area in the load configuration, a practice widely used to represent pavement loads (^{Huang, 2004}).

A square footprint was used, thereby achieving a better agreement with finite element meshing. Figure 3 shows the dimensions of the contact areas and the configurations for the two case studies.

Interactions

Interactions between contact surfaces in the finite element model were established. The interactions have different characteristics, which are important for the proper convergence of the solution and are established according to the physical relationship between its components. Figure 4 explains each one of the interactions, and their relationship with the real physical model.

Slabs-soil: For this interaction, a normal contact with a "hard-contact" algorithm was used, the interface responds to vertical stresses and simulating the slab unbound condition with the soil, allowing it to separate in some areas. This condition occurs in the real physical model.

Contiguous slabs face: It is based on the condition that the joints do not have separation between them, as happens in joints built with simulated crack. A normal contact is established with a "hard contact", "surface-to-surface" algorithm, to minimize the penetration of the slave surface into the master surface, allowing the transfer of stresses across the interface.

Dowel bar-slab: The dowel transfers the shear stresses between contiguous slabs, it is characterized by having one part embedded to the concrete and the other not, this allows it to slide without generating other types of stress in the slab caused by temperature or traffic. The non-imbedded part was simulated with normal and tangential contacts, the normal was modeled in the same way as the previous surfaces, the tangential was configured as "frictionless". The embedded part was modeled with a constraint known as "embedded region", an algorithm for modeling rebar within concrete (Figure 5).

Tie bar-slab: The tie-bar by its dimension and shape is considered as fully embedded to the concrete. The same algorithm employed in the embedded part of the dowels was taken.

Boundary conditions model

An infinite semi-space through a rectangular domain geometry for the concrete slabs foundation have been represented. Boundary conditions were applied to define the peripheric and bottom constraints. All horizontal movements (x-direction and y-direction) were constrained. The vertical displacements in both sides are not restricted, allowing the tensions and displacements on this direction to propagate until the end of the continuum, without generating errors in the numerical model. At the edge of the concrete slab corresponding to the slab center, conditions of symmetry on X axis were placed on the slab and the foundation (XSIMM U1 = U2 = U3 = 0) defining only the study section (Figure 6).

Calibration of numerical models

3D brick solid elements (hexagonal C3D8R 8-node) were used for the whole model. These elements are good for linear analysis which is the case for this model. The size of the elements varied depending on the importance of the section. The sections of interest to obtain the stress and strain (i.e. loaded section) were meshed with small sized elements (fine mesh). The remaining sections had large sized elements (coarse mesh). The variation of elements sizes allows detailed analysis but on the expense of memory and computational time. On the other hand, coarse mesh did not result in a detailed analysis but used little memory and shorter time (^{Mesa & Alvarez, 2011}).

For the selection of the appropriate mesh density a factorial design was made, as show Figure 7, combining the horizontal node span (b) and vertical node span (h); the aspect ratio of the element was respected. The study was only carried out in the center of the slab (including the concrete shoulders), in the rest of the slabs and soil foundation was maintained a constant mesh (100x100x100 mm concrete slabs and shoulders, 200x200x200 mm soil foundation). The spacing between nodes was varied from 100 mm to 20 mm in the vertical and horizontal directions. The control variable measured in each combination was the tensile stress at the point located below the load footprint applied at the edge of the concrete slab, where stresses are greatest.

The model was calibrated with software EverFE v2.24, it is a freely available 3D finite element program for the analysis of jointed plain concrete pavements developed by the Maine University and used by different researchers (^{Davids, Wang, Turkiyyah, Mahoney, & Bush, 2002}; ^{Šešlija, Radovic, & Togo, 2016}; ^{Zhang & Gao, 2016}). This software simulates a concrete pavement supported on a Winklerian medium, characterized by the reaction module K. The contrast value taken was the tensile stress, at the same point where the proposed numerical model was obtained.

The model studied presents the same characteristics of the slab as in the model developed by the University of Maine, but they differ in the way of characterizing the foundation. In order to test the solutions, achieving the equivalence between both results, the subgrade reaction module (k) required in the Maine model was transformed into an elastic module, through the correlation model obtained through the software (StreetPave, 2014).

Twenty-five combinations of mesh sizes were tested, both vertically and horizontally, as shown in Figure 8, comparing the results with those of the EverFE v2.24 program (fine mesh model). In the comparison of results, the norm of absolute and relative error was used at one point, the results of which can be seen in the graph in Figure 9.

The figure shows that the horizontal mesh density curves do not experience appreciable differences, so there is very little influence on the response variable. This is not the case with the vertical mesh since as the density increases, the tensions increase. For the selection of the most suitable combination, the curves of Figure 10 were constructed, with the least errors.

As shown in Figure 10 (a), as the mesh size in the horizontal is reduced, for maximum vertical density, the error is reduced, which occurs until it is practically kept below 3% with respect to the pattern value. This analysis confirms the efficiency of this mesh density, as indicated by (^{Zienkiewicz & Taylor, 2003}), who affirm that reducing mesh density increases convergence towards a better solution, up to a certain limit where its decrease no longer generates changes, although the calculation time continues to grow, which can be seen in Figura 10 (b), where it can be seen how the computation time increases significantly with the reduction of the mesh size.

It is concluded that the most recommended density for the available hardware is the resolution 40x40x20 mm in the slab and 200x200x200 mm in the foundation, with an approximate calculation time of 3 hours and 11 minutes.

Domain simplifications

Usually, the modeling of elastic solids in ABAQUS does not generate excessive calculation times, however, the dimensions and number of parts that make up the model, in addition to the multiple interactions (dowels and tie bars), cause delays in the convergence of the solution, which can be reduced by simplifying the domain.

Depth of continuum

It was tried to find with this analysis, the depth of soil foundation that represents the domain as the infinite semi-space in order to reduce the time of calculation. If the domain is very small, the joints placed in the contour introduce reaction forces that increase the soil response, reducing the stress by confinement effect. Different depths foundation was analyzed, verifying the displacement produced in each one, in a control node close to the contour. The most unfavorable condition of the model was taken (slab thickness 150 mm and foundation modulus 20MPa, without concrete shoulder). The results of the study are shown in Figure 11.

It is observed that when the depth in the soil foundation is bigger, it is closer to the condition of infinite space, however, the number of finite elements increases, therefore, the calculation time increases too. A depth was chosen, in which the calculation time was not high and the displacement was close to zero. The obtained depth is 1.50 m, this represents a calculation time of 33 minutes for the concrete slab without shoulders case. With a shoulder, time grew approximately 2.5 times (one hour and 23 minutes).

Geometrical simplifications

When it is located on the center edge of the slab, the presence of transverse joints has very little influence on the stresses, due to the relative distance between the load and these joints. The stresses at the center edge of the slab are independent of the type of load transfer mechanism (^{Pickett & Ray, 1951}). Geometric simplifications were carried out, using a model from a single slab and the results were compared with those obtained with the three-slab model, initially designed, under the same load conditions and materials (Figure 12).

In the calculations, different slab sizes were used, with the same lane width (3.50 m) in both models. In each case the stresses below the center of the footprint located at the center edge of the slab were measured, comparing the results by means of the graphs of Figure 13, this reflects the differences and errors in both models, with various sizes of slabs.

The variation of the stresses obtained in both models is shown on the Figure 13 (a). It can be seen that for slab lengths less than 4.5 m, on the three slab model, the stresses are lower than those of single slab model, which is explained by the influence of the load transfer mechanism, which contributes to minimize the stresses caused by the slab of the center, distributing part of them to the contiguous slabs. For lengths equal to or greater than 4.5 m, the stress does not differ, which can be better seen in the error graph of Figure 13 (b). It is observed that from 3.5 m, the error that is committed using one or the other model, is less than 5%, this confirms what was proposed by the authors, so it is more rational from the numerical point of view to use the model of single slab.

Comparative analysis of the numerical model

To validate the pre-established conditions in the mathematical calibration process, three models were compared: the ABAQUS / CAE model, the EverFE software (both developed with the finite element method) and the classic Westergaard´s response model for a load applied on the slab edge eq. (2). The equivalence between the reaction module and the elastic module for the soil was applied in the Eq. (1) for to achieve physical correspondence between the three models. The thicknesses were varying taking into account the practical dimensions, having in each case the bending-tensile stress originating in the center edge of the slab (Figure 14).

When:

(: Poisson's Coefficient of concrete

*h*: Slab thickness (cm).

*P*: Applied load (kg).

*E*: Elasticity modulus of concrete (kg/cm^{2}).

*k*: Reaction modulus of foundation (kg/cm^{3}).

*a*: Radius of the footprint load (cm^{2}).

*l*: Relative stiffness ratio (kg/cm^{4})

The EverFE computer program was used to compare the stresses obtained using the ABAQUS numerical solution. This software was developed by the universities of Maine and Washington, and applies three-dimensional (3D) Finite Element (FE) analysis in the design and revision of concrete pavements. It is equipped with qualitative tools to create 3D FE models. In the creation of the concrete pavement model its dimensions are determined first and the number of layers is defined, each represented by a specific modulus of elasticity (E), Poisson ratio (ν) and layer thickness (d) (^{Šešlija et al., 2016}).

The graph shows that using the ABAQUS numerical model, lower stresses are obtained, than in the other two models, all the graphs present a decreasing tendency as the thickness of the slab increases. The Westergaard´s model presents the most conservative results due to the considerations made by the author. The ABAQUS and EverFE models, despite sharing the same numerical tool (FEM), have differences but not so significant compared to the Westergaard´s model. These differences are due to the fact that they do not share the same physical phenomenon, since the resistance in the foundation soil is not given by the same parameter, in one case the modulus of elasticity (E) is used and in the other the reaction module is used (k), the difference was reduced by using a mathematical relationship between the two parameters. The stresses differences between the model studied and that of EverFE are not higher than 10%, and therefore, the model proposed for the investigation could be accepted.

Parametric study

In order to study the influence of some factors that cause variations in the pavement structural response, the geometry, the materials properties of the system and the tire contact pressure were analyzed. In each case, the parameter studied was varied, while the rest of the factors remained constant.

Slabs geometry

The slab geometry is composed of the dimensions: length, width and thickness. The thickness of the slab is the factor of greatest influence on the load stresses, and is the most conditioning factor in the pavement design. The width usually depends on the lane width and therefore on the road category, while the length varies depending on the available technological conditions, in simple concrete pavements it does not usually exceed six meters. A factorial design was programmed where the lane width is varied from 3.0 m to 3.75 m and the length between 3. 0 m and 6.0 m, as shown in table 1. Rectangularity is the relation between the length and width of the slab, used as an analysis variable.

Figure 15 shows the behavior of the tensions as the rectangularity increases, represented in four lane widths. The tensions decrease with the increase in rectangularity, however, these differences are not significant, since the maximum difference between values does not exceed 5%, which is why it is considered as a low impact effect, which allows us to assume in the study of tensions, settled dimensions within the studied interval.

*Concrete compression Strength*

A simple compressive resistance was considered, after 28 days of curing, as a variable to be taken into account in the calculation of stresses, from which the modulus of elasticity was estimated, by means of correlation (MPa). In general, the 28-day f'c values for concrete used in rigid pavement vary between 20 MPa and 40 MPa. In the study, the resistances in that range are varied, for each thickness between 150 and 250 mm.

Figure 16 contains the stresses obtained at the edge of the slab, varying the slab thickness and the compressive resistance of the concrete. Observe how the quality of the concrete does not cause a significant variation in the stresses, with maximum differences lower than 9%. By performing a multiple regression of the parameters and applying a significant test to each one, it was obtained that the compression resistance parameter has a p-value greater than 0.05, which is why it is not statistically significant, for a reliability level of 95.0%.

Support

In the study of soil influence, the resilient modulus was considered as resistance, obtained through its correlation with the CBR (California Support Index). Soil qualities between 50 and 800MPa were selected, which is approximately equivalent to CBR values between 5 and 80% respectively. The thicknesses were varied between frequent values used by the international standards in highway and street concrete pavement design. The results are shown in Figure 17.

From the graph it is observed that, as the thickness of the slab increases, the variations of the stresses decrease, evidenced by the reduction in the curvature, which shows the stiffness contribution of the slab with the action of the load, supported every time in better foundation resistance. It is also observed that from 350 MPa, the soil tends to have less stress influence, corroborating what was expressed by authors and institutions (^{Cajka, Burkovic, Buchta, & Fojtik, 2014}; PCA, 1984; ^{Yoder & Witczak, 1975}), who point out that improvements of the support layer do not have a structural purpose, these layers are placed to reduce the effect of erosion and the rise of fine materials to the surface (pumping). However, the improvement of the properties of the foundation, contributes slightly to improve the behavior of the system, which is a factor to be taken into account in the obtaining of the work stresses, especially when the slab is thin.

Contact pressure (q)

It is of interest to evaluate the influence of the contact pressure on the stresses originated in the edge of the slab, since during the exploitation of the pavement the contact pressure can vary significantly, due to the increase of the temperature inside the tire, by the impacts produced by irregularities of the surface. In the calculations, several contact pressure values were considered; those shown in Table 3, in addition to the dimensions of the equivalent footprint on the surface. With the load placed on the edge of the slab, stresses were obtained every 20 cm in the traffic direction, these results are shown in Figure 18.

It was found that the increase in contact pressure increases the stresses at the edge of the slab, for the same load conditions, evidencing greater increases in the center of the edge, with reduction of these differences as it moves away from the center, so that the stresses at the ends of the slab are very reduced, practically equal, for any value of contact pressure, being verified what was published by (^{Pickett & Ray, 1951}) regarding the few influence of the load on the transversal joint of the slab.

Results and discussion

The use of computational modeling in mechanistic empirical analysis for the design of simple concrete pavement

In the empirical-mechanistic design for concrete pavements, the system formed by the slab and the foundation is modeled, with dimensions and resistant characteristics adopted in the design, to provoke the response of the structure, in tensions and deformations, to loads imposed on specific critical areas, so that the thickness of the slab can be obtained capable of withstanding repeated solicitations until the failure occurs. The response of the system is obtained with the help of the computational model built using finite elements, placing the loads at the critical points.

In the design, as provided by the PCA procedure, it is determined, for a given slab thickness proposed, and depending on the variables established in the project, the permissible repetitions for each single, double and triple axle load, both for the criteria of fatigue as erosion.

The computational modeling then allows obtaining the tensions originated at the critical points, in different conditions and materials of the foundation and the slab, obtained from applying reference loads, quantities that must correspond to the design regulations. These tensions are reduced to values that can be obtained at a distance of 60cm from the edge, for which the factors offered by the PCA (Portland Cement Association) are used, and whose research has shown that only 6% of trucks pass to this distance.

The PCA method bases its analysis on the verification of the two main failure modes of rigid pavements. The Fatigue criterion is the one that allows to maintain the pavement stresses, produced by the repetitive action of loads, within the safety limits and thereby prevent fatigue cracking. While the Erosion criterion is concerned with limiting the effects of deflection of the pavement on edges, joints and corners of the slabs, thus controlling the erosion of the materials of the lower layers (^{Calo, 2012}). The load positions for the analysis and design, in which the tensions that can cause the failure are generated by: the load located in the middle of the longitudinal edge of the slab, for fatigue failure, and the load in the corner for erosion failure.

The stresses that are produced by each type of load on the edge of the slab, can be used to determine the relationship of stresses between the tension produced by each load and the modulus of breakage of the concrete, according to the fatigue failure. In the same way, the stresses produced at the corners of the slab, calculated by modeling, are compared with the admissible stresses obtained from the erosion model. By means of failure models, the relationship between the expected repetitions of each axle and the admitted ones is determined, obtaining the consumption of fatigue or erosion damage for each of the planned loads. The sum of fatigue or erosion consumptions in each of the planned loads must not exceed 100%.

There is a great dispersion among the models used internationally to evaluate fatigue failure, the result of research that offers laws developed both in the laboratory with beams and in experimental sections, where the most widespread is that of the Portland Cement Association (PCA). The PCA fatigue model is based on information originating in fatigue tests on beams, developed during the 50s and 60s. ACPA improved the PCA model, including reliability as a parameter for cracking prediction in concrete pavements, depending on the importance of the road.

The heterogeneous and random nature of the parameters involved in the tensions produced by the loads have led to the establishment of procedures that allow these tensions to be determined from a model conceived in specific working conditions. Once the computational model has been calibrated and validated, its importance lies in the fact that it can be used to generalize any specific condition of solicitations, both loads and traffic, modifying the variables involved, specific to a given territory, such as the conditions of the foundation and of the subbase, as well as the characteristics of the concrete and climatic conditions.

Conclusions

The numerical method as a solution in the modeling process through the MEF, with the use of multipurpose software, is a powerful tool that allows analyzing different working conditions of a rigid pavement slab. But it should be kept in mind that the confidence in these procedures lies mainly in the selection or obtaining of appropriate behavior models for the materials, in correspondence with the environmental conditions of the place where it is applied.

After addressing the procedure and general considerations about the computational modeling process, the following conclusions can be reached:

For the sizes of slabs commonly used in simple concrete pavements (L ≥3.5 m) there are no significant differences in the stresses (at the study control point) between the modeling of a system with three slabs (with dowels) and the model of an isolated slab. The application of the isolated slab, as a definitive model, speeds up the convergence, reducing the calculation time to 25 minutes (with concrete shoulder).

With the proposed model, the stress calculated was lower than the solutions obtained by the Westergaard and EverFE software models. The main differences in the tensions are due to the use of a continuous model in the foundation and improvements in the formulation of the finite element method, however, the proposal follows the same tendency in the tensions, for each one of the factors in the models analyzed, which shows a correct modeling of the phenomenon.

In a general way in the research, studies were carried out that aim to reduce costs in experimental tests from combining computer simulation and experimentation, an approach that is integrated into a process, where after numerically simulating the phenomenon, a mathematical and physical calibration is performed, with a validation, taking into account the same conditions of the experiment.