3D Numerical Modeling applied to analysis of piled foundations

Applying numerical modeling to geotechnical engineering problems has become a powerful tool to interpret the behavior of foundations and earthworks. However, not all numerical analyses follow the requirements for the model to be representative of the experimental condition. This work describes the steps required during modeling to determine the appropriate dimensions of the half-space and the respective boundary conditions, as well as the tests appropriate to perform this determination. On the other hand, the results of comparisons between the numerical and the experimental model for a piled raft are shown.


INTRODUCTION
Foundations are designed to convey the load of a building to the soil in a safe and economic way, thus ensuring reliability and easy maintenance of the structure [1].According to [2], over the last years a growing number of structures (particularly high-rise buildings) have been executed on combined systems instead of being considered as an "alternative" to rafts.Therefore, piled rafts are considered as a foundation element with multiple interfaces, among which the piles under the raft are interrelated.According to [1], most works about foundations on piled rafts have focused on the load vs. displacement capacity or behavior, and few were dedicated to examining the influence of the factors that govern analyses via numeric tools and their parameters on the interactions within the system, and, as a consequence, on the load sharing mechanism.
The methods to analyze the behavior of piled rafts are complex due to the large number of factors involved in the raft-soil-pile interaction [6][7][8].As a consequence, it has become common practice to use computer tools to enable numerical simulation of complex piled raft foundations.The key methods are: Finite Element Method (FEM), Boundary Element Method (BEM), Finite Layer Method (FLM), Finite Difference Method (MDF) or a combination of two or more of these methods.The finite element method can be considered as a tool that can be used to analyze the behavior of foundations.This method takes into account the effect of interaction factors on the process of analysis, such as: pile-pile, pile-raft, raft-soil and pile-soil interactions [9].The numerical methods depend on the number of elements that make up the finite element mesh, aiming at the adequate representation of the analyzed model, consequently requiring high speed processors and high RAM [10].
Despite the high capacity of numerical tools applied to modeling of problems of foundation engineering, some principles should be followed in order to provide reliability to numerical analyses.Three key steps stand out: • Tests involving comparison to the results from the literature obtained by other authors; • Convergence test to delimit the border; • Calibration of the model with experimental results.

NUMERICAL MODELING
The processing time in numerical analyses varies greatly as a function of the dimensions of the model and composition in terms of number of elements.In this sense, using the tool may be impossible for practical purposes due to long processing time.A reduction may be viable by analyzing only part of the problem, due to its axisymmetry, i.e., symmetry from an axis.An example is shown in Figure 1: a typical case of axisymmetry that produced a 75% reduction of the model under analysis, i.e., the modeling was applied only to a quarter of the problem.The dimensions of the numerical model were attributed based on tests performed to ensure that the boundary conditions at the corners of the problems, which might be considered as not displaceable or else which might have displacements of minor magnitude, would not impact the results of the analyses (Figure 2).According to [11], the limits of the model should be defined in sufficient spacing so as to minimize the influence of strain on the boundaries of the foundation.The nodes in both side limits of the model are fixed against horizontal movement (δx = 0), but they are free to move vertically.Meanwhile, the nodes in the lower limit of the model are fixed against vertical and horizontal movements (δx = δy = 0), while the upper limit is free to move in both directions.
Tests should be carried out at the initial phase of the analyses to validate the dimensions of the half-space.
The geometry of the model will be determined via convergence tests, which involve checking whether the results of the boundary conditions are in conformance with the definitions set in the pre-processing step.
In this work, the "calibration" of the model of finite elements was carried out by comparing the results of a load test executed on a single excavated pile (L=5 m and φ=0.25 m), as shown in Figure 3a.This pile was compression-tested by [14] at the Experimental Site of the Soils Mechanics and Foundations program of the School of Civil Engineering, Architecture and Urbanism of the University of Campinas (Unicamp), and then compared to the results obtained via numerical analysis of this very same pile to the soil parameters obtained by [15].This step helped numerical analyses become more reliable in terms of applicability and interpretation of a piled raft composed of one pile, as shown in Figure 3b.

CONSTITUTIVE MODEL
The Mohr-Coulomb model, which is perfectly elastoplastic, was used to simulate the non-linear behavior of the soil in terms of stress vs. strain.Soil properties assigned to the different soil layers followed adopted failure criteria and are given by the following parameters showed in Figure 8: unit weight (γ); cohesion (c); friction angle (φ); strain modulus (E) and Poisson's ratio (ν).These parameters had been previously obtained by [15], except for Poisson's ratio, which was adopted because of the soil behavior as appraised by the tests.For materials with a brittle behavior (Parabolic Model) such as concrete from piles and raft, laboratory-determined values obtained by [16] were attributed for compression strength (Rc), E c and ν c , and values for tensile strength (R t = 10%•R c ) and γ c were adopted.
The mesh of finite elements that was used included triangular shaped elements of quadratic interpolation, which were extruded every meter in depth to produce a volumetric element of the pentaedric type.Numerical analyses were performed using the LCPC-CESAR v.This software is a 3D FEM (Finite Element Method) tool.To adequately balance the computational effort and the convergence of the obtained results, a quadratic pentahedral element consisting of 15 nodes was used."

CONVERGENCE TEST
Prior to the numerical simulation of the result of a load test carried out in real scale at the site [14], the dimensions of the model and the respective quantities of nodes and elements were checked  Due to the symmetry of the problems under analysis, the numerical tool made it possible to carry out the analyses with a fraction of the problem.This way, the convergence tests resulted in a "block" measuring approximately 25 x 25 m and 10.8 m deep, representing a quarter of the problem, i.e., the total dimension of this test would be equal to 50 m x 50 m on the plane.Note that, from the depth of 10.8 m down, the material is impenetrable.Another key aspect of this analysis was the use of the resource of symmetry, which provided a reduction of the dimensions of the problem, and, as a consequence, the mesh of finite elements, leading to a dramatic reduction of the number of elements and nodes.This way, it became possible to simulate a model with thousands of elements and nodes within a short time.An elastoplastic model was used, which varies as a function of the stresses applied, following a non-linear behavior.
The mesh of finite elements included triangular shaped elements of quadratic interpolation, which were extruded every meter in depth to produce a three-dimensional pentaedric element.The composition of the problem resulted in a mesh of finite elements comprising 6,415 elements and 18,108 nodes (Table 1 and Figure 5).In order to supplement the validation test, refinement tests were made to the mesh of finite elements to check which changes were occurring at the maximum displacement obtained in comparison to other conditions of refinement, with less density in terms of elements and nodes.
This test increased the degree of reliability of the results obtained in the post-processing phase.Table 1 demonstrates that the results of the refinement tests of the mesh of finite elements determined the optimal quantity of nodes and elements, beyond which an increase would not translate into better quality results.
Figure 5 demonstrate that, for the 25 x 25 m mesh assessed in the convergence test, the displacement obtained in the mesh refinement by increasing the quantity of elements caused the displacement initially obtained (10.83 mm) to display an elevation to the 1st and 2nd stages of refinement and remain practically stable for the 3rd stage.Therefore the option was to keep the mesh with the same quantity of elements and nodes as the previous phase (2nd stage), since there would be no improvement with greater refinement of this mesh.
According to the tests conducted before, the geometry of the problem and the respective dimensions appropriate to the boundary conditions were obtained, as shown in Figure 6.Finally, the input parameters to represent the soil were checked and "calibrated" via comparison to the results obtained in a load test on a single pile conducted by [14], at the Experimental Site of the Unicamp.

IN SITU TESTS
The experimental step of this work was developed at the Experimental Site of the Unicamp, located in the city of Campinas, State of São Paulo, Brazil (Figure 7).The region is formed by basic intrusive rocks of the Serra Geral Formation (Diabase), of the São Bento Group.The local subsoil is composed by porous soils under unsaturated conditions, with collapsible characteristics.Figure 8 shows the streamlined geological / geotechnical profile of the experimental site, obtained from laboratory tests carried out by [15] and in-situ tests carried out by [17].
The pile was executed by perforating the soil with a tool composed of a 0.25 m diameter helical auger to excavate the soil down to the depth of 5 m.Reaction piles were executed with the same equipment for the load test, but with a 0.6 m diameter auger to perforate the soil down to the depth of 9 m.The Characteristic strength (f ck ) of the concrete used to fill the piles was equal to 25 MPa (28 days).
Figure 8.Average parameters of the geological profile at the experimental site [16].
In order to get information relating to load transfer, steel bars instrumented by strain gages were installed.This system was previously checked at a laboratory before being installed in the pile shaft.Owing to the short length of the piles, the option was to install the instrumentation at the top and at the tip of each test pile.Figure 9 shows the position of the instrumentation in the piled rafts and the respective properties of the materials.
Instrumentation with electric resistance extensometers or strain-gages was used as an indirect manner to obtain strains.The extensometers used were of the double rosette type with complete bridging.The electric extensometers were installed in CA50 steel bars (φ=12.5 mm, L=0.60 m), sheltered against humidity and mechanical shock by application of an appropriate resin and coating for protection against impact and humidity, laboratorycalibrated, and bonded on the field via sleeves to form a continuous bar.
Figure 9. Position of the instrumentation along the pile shaft and properties of the materials of the piles and raft [16].
A slow maintained load (SML) test was carried out per prescriptions of [18] and the sketch shown in Figure 10.The loadings were made in equal and successive stages, not above 20% of the work load forecasted for each piled raft that was tested.At each stage, the load was maintained until the displacements stabilized or at least for 30 min.All load readings, axial displacement of the pile and strains of the instrumentation were obtained through a data acquisition system, and managed by a reading transduction software program.
Figure 10.Sketch (front view) of the Main Reaction System [16].
By means of tests of comparison to the results obtained in a load test performed in a single pile carried out by [14] and back-analysis, it was possible to get the soil parameters suitable to the numerical representation, by comparing the experimental results of load vs. displacement.The pile analyzed by [14] refers to a load test performed on a single pile excavated by a mechanical auger (bored pile), with 0.25 m diameter and length measuring 5.0 m executed at the Experimental Site of the Unicamp.

VALIDATION OF THE NUMERICAL MODEL
The results of the load test carried out by [14] and the respective numerical analyses of this isolated pile conducted with the soil parameters originally obtained by [15] demonstrated substantial agreement among the back-analyzed experimental and numerical results, so the parameters for soil resistance calibrated in this analysis (with no contact) could be extended to the case of piled rafts (with contact).Note that the maximum test load obtained by [14] was 180 kN.However, the ultimate load was stipulated to a displacement of 10 % of the rated diameter of the pile, i.e., 25 mm.This way, the ultimate load results in 174 kN for the experimental load vs. displacement curve and 171 kN for the load numerically obtained, i.e., the load obtained numerically is 1.7% smaller than the experimental ultimate load.
The distributions of load in depth, as obtained by means of numerical analysis and experimentally, were compared for the 1st, 5th and 10th loading stages.The observation was, numerically, 6% of participation of the pile tip in the maximum test load and no participation recorded by [14] in the load test.In the previous loading stages, the observation was that pile resistance was basically the result of the resistance by skin friction.

EXPERIMENTAL AND NUMERICAL ANALYSIS OF A PILED RAFT COMPOSED BY ONE PILE
The results of the load vs. displacement curve obtained by numerical analysis for the case under analysis tend to display greater displacements in comparison to the experimental results, since the pile vs. soil sliding is different in the numerical model and in the experimental model.In this sense, Figure 11 demonstrates that the numerical model for the piled raft displayed greater displacements in comparison to the experimental piled raft.
Despite the forecasted inflection in terms of displacements in the curve or the numerical piled raft, the observation is that up to the 5th stage (section close to the workload), understood as a linear elastic section, the differences in behavior are minor.It can be stated that they have good agreement in terms of predictability of behavior.The same cannot be said of the section comprised between the 6th and the 10th loading stages, in which the numerical model displays high displacement when compared to the experimental pile raft.This section can be understood as one of elastoplastic behavior, i.e., the strains that occurred in the soil massif are permanent.
Figure 11.Load vs. displacement curves of the experimental and numerical piled rafts.
The differences and agreements in behavior between the numerical and the experimental models can be seen in Figure 11.This figure shows that only in the 1st loading stage there is good agreement between the tip portions, skin friction and raft-soil contact.Up to the 4th stage (100 kN -½⋅Q máx ) the behavior of tip is the same for the experimental and the numerical models.
Figure 12 also shows that, from the 5th stage on, the tip portion of the numerical model tends strongly to growth, and attains 11% of tip participation at the end stage, in comparison to 6% of the experimental model.This larger tip participation can be attributed to the difficulty to provide the numerical model with the in-situ conditions and the executive process of the type of pile used.The resistances by skin friction and raft-soil contact are divergent since the 2nd loading stage, as shown in Figure 12.At the end stage, the portions of resistance by skin friction and raft-soil contact are discrepant in the experimental piled raft, but coherent with the numerical model.Therefore, it can be noted that the performance of the foundation element requires further evaluation and understanding.The greater participation of the raft-soil contact and the pile tip for the group of piles may be associated to the fact that the numerical model has more pronounced displacements, which may favor exhaustion of skin friction and increase the other portions of resistance (tip and raft-soil contact).
Figure 12.Load distribution between experimental and numerical models.
Figure 13 shows the load transfer curves for the piles of the experimental and numerical piled rafts.In the 1st stage there is good agreement of behavior between the experimental and the numerical models.This evidence cannot be seen in the 5th and 10th loading stages.
Figure 13.Load transfer in the pile of experimental and numerical piled rafts.
The divergences found in the 5th stage of transfer are more influenced in the numerical model in the first 2 m of the pile length.After this section, there is good agreement with the experimental results up to the tip of the pile.For the maximum test load (Q máx -10th stage), the observation was that the numerical model absorbs less load than the experimental model; therefore it does not display good agreement in behavior, although high discrepancy is not noticed.The variation in the Q/Q max ratio reveals how each piled raft, both experimental and numerical, behaved as a function of the loads applied.Figure 14 shows the results for Q/Q max and the respective displacement for each loading stage.Figure 14 shows how the loadings influenced the reduction in the Q/Q max ratio and in the magnitude of displacement.
Figure 14.Q/Qmax ratio as a function of displacement.
Up to the 7th loading stage, the numerical model displayed Q/Q max greater than in the experimental result.In the 8th and 9th loading stages, the numerical model did not display variation in Q/Q max .This ratio was seen to be smaller than one (1) only at the final stage for the numerical model, which indicates exhaustion of this foundation when the maximum test load is applied, as shown in Figure 15.It must be pointed out that the reduction in Q/Q max of the numerical model is higher or equal to that of the experimental model, even if the numerical model produced displacements higher than those obtained in the experimental piled raft.Therefore, it is not appropriate to relate the accentuated displacement to a smaller load capacity.
Taking into account the considerations above, the displacements of the experimental and numerical piled rafts were normalized and compared to the Q/Q max ratio (Figure 15) and to the also normalized load (Figure 16).
Figure 16 shows that the reduction in the Q/Q max ratio of the experimental raft is sudden and asymptotic only when it is close to one (1), unlike what occurs with the curve of the numerical model, which starts to tend to converge when the Q/Q max ratio is smaller than 2.0.The graphs above show that the behavior of the experimental piled raft is influenced by the Q/Q max ratio when the loading stage is close to the maximum test load.The Q/Q max ratio of the piled raft under numerical analysis is influenced since loading up to half of the breaking load, equaling the behavior of the experimental piled raft only at the normalized displacement equal to 7.5%.
The behavior in terms of normalized load and displacement is shown in Figure 16.It corroborates previous analyses, i.e., only from the 5th stage on (0.5 of the normalized load) more accentuated displacements start in the numerical model, which is limited to simulating behavior of the experimental raft, which shows a curve less susceptible to loadings.The displacements of the experimental piled raft tend to convergence to "rupture" from the 6th to the 7th stage.
The ultimate load of the experimental piled raft is 200 kN and 184 kN of the numerical model, i.e., 8.7% larger, to point out that the normalized load vs. displacement curve in the numerical model is down between 0.5 and 1.0 of the normalized load (Figure 16).Therefore the numerical model has less load capacity (Q normalized = 0.89) for the same displacement imposed to the experimental piled raft (Q normalized = 0.96) for 10% of the pile diameter.
The analyses performed so far demonstrate the differences in behavior between experimental and numerical piled rafts.However, it is complex to understand the contribution of the raft-soil contact.Additional analyses are required to determine changes in behavior and the load capacity of piled rafts, in case they do not touch the soil.

CONCLUSIONS
For an appropriate numerical modeling, it is critical to perform an in-depth study on the boundary conditions of the problem.Convergence tests have shown that the optimal spacing between the mesh and the axle of the pile should be 100 times the diameter of the pile.As to analysis of depth, the problem was limited to a depth twice the length of the pile, since the impenetrable soil is found precisely at this level.
The optimal number of elements and nodes should be correctly appraised by means of convergence tests and comparison to the results from the literature before being applied.On the other hand, there is an optimal number of elements and nodes beyond which an increase in precision is negligible in comparison to the increase in processing time.
The validation of the numerical model and the respective constitutive method by means of back-analysis with experimental data ensure consistency between the numerical response and the experimental results.
The resource of symmetry proved to be appropriate to solve the problem.It not only produces coherent results, but it also makes it possible to use a model with thousands of elements and nodes within a short time, thus optimizing the time for analyses.
The numerical model produced satisfactory results as to the estimation of load transfer of the piles that compose the raft.However, it was noted that the largest divergences are in the value of the top load, precisely at the spot of the splitting of the participation of raft-soil contact and pile.Applying greater refinement of the mesh of finite elements is recommended to get more reliable results.
Analyses performed with the numerical and experimental models in terms of normalized load and displacements have demonstrated that, in absolute values, the maximum load and displacements were close.However, the differences detected are maintained as to the behavior of load transfer.
The 3D numerical modeling via finite elements is a powerful tool that enables correct, accurate analyses of geotechnical problems.Nevertheless, some important steps to understand the behavior and the responses of this tool are required.

Figure 1 .
Figure 1.Sketch of the numerical model.

Figure 2 .
Figure 2. Boundary conditions and vertical displacement of an axisymmetric 3D model.
. The geometries assessed were: 10 x 10 m, 15 x15 m, 20 x 20 m and 25 x 25 m.The boundary conditions imposed to the problem were checked at each simulation of the validation process, resulting in the graph in Figure4, in which the boundary conditions for the 25 x 25 m mesh displayed maximum strains smaller than one tenth of a millimeter, which is considered negligible for the type of problem under analysis.

Figure 4 .
Figure 4. Displacements at the corner of the mesh of finite elements.

Figure 5 .
Figure 5. Variation in displacement as a function of the number of nodes.

Figure 6 .
Figure 6.Model of finite elements for piled raft.

Figure 15 .
Figure 15.Mean normalized displacement as a function of Q/Q max .

Figure 16 .
Figure 16.Mean normalized displacement and normalized load of the piled rafts.
4.07 software, developed at the Laboratoire Central des Ponts et Chaussées (Road and Public Works Research Institute).

Table 1 .
Results of the refinement general and local tests of the mesh of finite elements.