NUMERICAL MODELLING OF DISPLACEMENT PILE RESISTANCE IN SAND GROUND. PART 2: DISCRETE MODEL, AT REST STAGE, LOAD TEST

. The article presents a continuation of research results on numerical modelling of displacement pile prototype ground (sand) resistance. Numerical modelling of ground resistance consists of three stages: restitution of initial soil stress state, the restoration of stress state developed during pile installation and the modelling of pile resistance against applied load (pile load test). Restitution of initial ground state induced in artificially created sand deposit after compaction procedures was realized by creating stress history and defining experimentally determined over consolidation stresses. Installation effects were restored using a well-known empirical approach based on relation between CPT test data and radial stress increase. Hardening soil model and its parameters were employed for modelling pile test substantiated in the first paper (Part 1) of the research. The discrete model of pile, soil layers, and pile-soil contact was created. Modelling


Introduction
This paper is continuation of Part 1 (Martinkus et al., 2021) of research on numerical modelling of displacement pile (DP) resistance in sand ground.The reliability of modelling DP foundation resistance evolution, including stress and state corresponding ultimate load, depends on two groups of factors: qualitative evaluation of deformed foundation resistance measures and relevant discrete and mathematical models, both created and calibrated via back analysis of soil resistance within predicted stress and strain ranges (different depths of corresponding vicinities).Any proposed method requires its validation by specified tests (clearly identified physical and mechanical properties with full scale or larger prototypes, reducing the scale factor effect aiming at minimizing the side factors).
Suitability for practical application of proposed method for design practices is expected to be not too complex, but reflects the main resistance peculiarities at certain loadings ranges, also evaluates larger variability of soil profile and properties, not completely investigated compared with comprehensive soil testing, specified profile preparation and testing.Therefore, for major cases practical design approach is of empirical nature.Design procedures are related either to LRDF (load resistance design factor) or global factor (reliability based design) in respect of ultimate state implementation.In normative geotechnical design documents or recommendations, an extra overdesign factor (model coefficient) is allowed to be employed.Large scatter resulting overestimation as well as underestimation of pile resistance predictions versus experimental, when employing various proposed methods and approaches (in some cases even being in conflict with similar foundation conditions), can be met (see, e.g., comparative analysis for pile group in Norkus & Martinkus (2019).
Numerical modelling techniques, available for more complex numerical simulation tools and convenient for practitioners, are under permanent demand.The proposed method or techniques should evaluate governing resistance measures and peculiarities of pile foundation resistance evolution, compatible to experimental results (some underestimation or overestimation is obviously met).General approach is to use methods based on complex proper stress and strain evolution analysis, reducing "adjusting factor" to compensate modelling versus experimental results, e.g., pile load-settlement relation.The main aim of the current research is to determine governing resistance evolution simulation paths applying modelling techniques.The obtained results for practical design can be used by introducing either smaller partial resistance factors or global safety factor, ensuring design reliability according to normative documents.
Investigations of piled foundations with the aim for practical applications employ different numerical and soil models.Comodromos et al. (2003) used numerical simulation to create loadsettlement relation for bored single piles and pile groups.Experimental results were used in numerical back analyses to calibrate and verify the soil mechanical properties with the aim to establish numerically the unaffected behaviour by the interaction single pile via load-settlement graph.The elastic perfectly plastic Mohr-Coulomb model in compliance with the non-associated flow rule was employed to model deformed soil behaviours of layers.Despite more complex soil models (e.g., HSM and other) that allow evaluating soil dilation and compaction processes, the application of simple Mohr-Coulomb model was sufficient considering that the anticipated stress paths dominated by shear failure.The numerical simulation of the pile load test was performed through the finite difference code FLAC (FLAC 3D, 2000).Mesh generation was parametrically defined for geometrical variations.Han et al. (2017) performed a 3D FEM analysis applying advanced, two-surface-plasticity, constitutive sand model to model resistance evolution of centrically loaded bored piles.The authors developed a mathematical model and subroutine (VUMAT) compatible to be implemented in ABAQUS (ABAQUS, 2012).A semi-implicit backward-Euler (cutting-plane) algorithm was adapted via the sub-incrementation and error control, employed to update stress and user-defined state variables from strain increments at each Gauss point.It was concluded that realistic FEM (finite element model) analysis required formation and evolution of shear bands in zones close to pile depending on stress and strain evolutions and boundary conditions (confining stresses), as well as careful FEM generation.Imseeh and Alshibli (2018) proposed the FEM within framework of discrete element model (DEM) for simulation of load transmitting.It was stated that in terms of evolution of force chains, fracture modes, and stress-strain relationships resulted in good agreement with testing results presented by Cil et al. (2017).It was stated that both discrete and continuum models successfully provided effective stress-strain simulations for granular materials.It was obvious that continuum model for soil ignored particular nature of disperse soils.DEM assumed that soil particles were absolutely rigid and particular material stressstrain relations were calculated via model contact forces.Therefore, taking into account computational resources required for foundation resistance analysis (compared with simulation of laboratory tests), the mixed finite element method (FEM) was combined with DEM.DEM was employed only for certain contact zones of structure and soil contact.Thus, a combined approach that addresses both the continuum and discrete behaviour of granular material has been established in recent years, which is known as a finite element-discrete element method (FEM-DEM), e.g., employed for compression test (Druckrey & Alshibli (2016)) within code ABAQUS (ABAQUS, 2012).It adopts FEM to simulate the constitutive behaviour of discrete particles, in which each discrete particle is meshed into deformable finite elements.
The simplified piled foundation resistance prediction methods were also reported.Matsumoto et al. (2010) developed a simplified 3D deformation analysis method for eventual application in practical design.Interactions of foundation structural element -soil -was based on Mindlin solutions (Mindlin, 1936).Nonlinear foundation response to loading was calculated by employing the simplified perfectly elastic model of soil springs in compliance with the simple analytical (empirical, semi-empirical) relations (e.g., ultimate base, skin resistances, other) for limit states of foundations proposed by other researchers.Validation with tests was performed with the aim for applicability in design.
The study aimed at evaluating the governing resistance measures and peculiarities of displacement pile foundation resistance evolution on the stress and strain evolution analysis applying modelling techniques, validating the obtained results with the experimental ones.
Reliability of applied techniques for numerical modelling of the displacement pile ground resistance depends on relevant representation of three main stages: (1) initial stress state before pile installation; (2) recreating ground stress induced during pile installation (initial zero stage prior to pile loading) and (3) modelling pile resistance versus loading.The second stage includes evolution ground stress and strain stage from initial "at rest" state to "zero state", developed by penetration of displacement pile to design depth by installing load and having removed installing load.The third stage includes evolution of pile ground resistance via loading-unloading starting from "zero state".At the beginning of this stage, i.e., prior to loading DP, pile weight load is in equilibrium with shaft and base resisting forces, resultants of contact stress distribution at pile side and base surfaces.Determination of the main physical properties and implementation of DP prototype testing program in artificially created sand deposit with monitored properties per depth are described in Part 1 of the paper (Martinkus et al., 2021).Pile numerical modelling of ground resistance against loading includes three main phases: 1. Evaluation of geostatic "at rest" state in artificially created sand deposit; 2. Restoring of installation (pile penetration) effects, resulting in zero ground stress and strain state, developed prior to DP loading; 3. Incremental loading and unloading of DP.An explicit modelling of pile installation effects, similar to modelling of cone penetration test, is a complex problem of large deformations.In the current research, a simpler practical approach was employed.The approach allows using the FEM software Plaxis for all aforementioned modelling phases.Substantiation of applying the hardening soil physical model (HSM), selection for modelling FEM software Plaxis, the determined and calibrated HSM parameters are given in Table 2 of Martinkus et al. (2021).
Creation of discrete model for numerical simulation of DP ground resistance and pile evolution modelling by three aforementioned phases is analysed.

Discrete model and boundary conditions
Ground volume is required for modelling.Taking into account pileground interaction symmetry, 3D problem can be transformed to 2D the axisymmetric (plane) problem.The use of this transformation for modelling of DP resistance results in the significant savings of computational resources.
Discrete model for single steel displacement pile (SDP) of length L = 1.445 m and diameter D = 0.219 m was created by applying the widely used "floating pile" concept (Said et al., 2009;Mascarucci et al., 2013;Mascarucci et al., 2016).Applying the concept, the necessary deformed ground volume was created.Regular rectangular ground volume of prism form was constrained in the following way: bottom surface in vertical and horizontal directions, side surfaces in lateral directions.According to Randolph and Wroth (1978), ground deformations at distance L are vanishingly small, the ground lateral constrains are introduced at distance 7.5 D = 1.14 L from pile centre.The bottom constraints are introduced at depth L form pile base.Sufficiency of L depth was also confirmed by Mascarucci et al. (2013), (2016) and by analysing ground behaviour results obtained by implementing the SDP test program at Geotechnical Laboratory of Vilnius Gediminas Technical University (testing program described in Martinkus et al. (2021)).The ground volume along 2L depth from surface depth is divided into separate layers (zones).Layers are introduced in accordance with physical properties and stress state.Different mechanical parameters (peak friction angle j′ p , dilatation angle ψ′ p and friction angle at critical state j′ cs ) for layers are assigned.Following Bolton (1986), the relation is used for this case.The magnitudes of assigned aforementioned parameters depend on stress state, which is different for deformed ground depths of loaded pile.Based on ground average stress levels per depth, the five design layers with relevant angles j′ cs and ψ′ are introduced.The created 2D axisymmetric discrete model of layered ground is given in Figure 1a.
Pile-soil and separate soil-soil interaction is described by contact parameters.According to the applied FEM Software Plaxis 2D (PLAXIS, 2016) and hardening soil model (HSM), the contact is described by friction angle of interface δ inter in relation with soil angle of internal friction.Larger pile load values result in the development of significant magnitude radial (horizontal) normal stresses s′ r , close to pile side surface.Subsequently, the increment of radial stresses (the components of shear stress τ in sand soil), leads to development of shear bands (Lings & Dietz, 2005;Loukidis & Salgado, 2008;Yang et al., 2010).Outer surface of band (see Figure 1b) is characterised by distance from pile side surface (thickness t s ), forms slip zone, where limit shear stresses τ u develop.The thickness t s mainly depends on the average size of soil particle d 50 (Uesugi & Kishida, 1986;Viggianni et al., 2001;Frost et al., 2004) and vary within bounds of (5-20) d 50 (Uesugi et al., 1988;Nemat-Nasser & Okada, 2001).Actual thickness t s = 9 mm was measured having pulled DP of ground after DP test.It should be noted that measured shear band thickness was out of bounds, namely t s = 27, d 50 = 9 mm, i.e., it was larger than the values stated in the aforementioned references.On the basis of measured t s , it was concluded that developed ground slip zone along DP side surface corresponded to soil-soil (versus pile-soil) interaction and the pile-soil contact friction coefficient was assigned to δ inter = j′ cs.= 30.5°.
The discrete model ground and pile-soil contact properties are given in Table 1.Creation of discrete model (Figure 1a) and assignment of design parameters for discrete model (Table 1) are related with the use of comprehensive ground laboratory test results, taking into account soil responses within loading ranges.Selection of parameters for HSM implemented by PLAXIS (PLAXIS 2D, 2016) was explained in the first paper of investigation (Martinkus et al., 2021).
Depending on the state of sand prior to pile loading (relative density, confining stress, compressibility index, reported in Norkus & Martinkus, 2019), it may either dilate (denser soil) or contract (looser soil) during shearing before it reaches critical state.For critical state adopted, no longer any change in volume or stress and strain softening after the maximum rate of dilation and peak shear strength is observed (Han et al., 2017;Comodromos et al., 2003).A mathematical model requires to capture the development of peak and critical-state strengths during simulation of pile foundation response in the loading process.
Assigning parameters for layers 1-5 (5 and 6 rows of Table 1) was based on dilatation effect evaluation (depending on positive volumetric strain magnitude) of soil element under stress states (measured during pile test (in vicinity of pile base and below and above pile base by surcharge load calculations) corresponding to ranges of foundation states of soil layers.24 direct shear tests (under the constant normal load) were conducted for 8 different magnitudes (3 tests for each magnitude) of normal stresses of 10 kPa, 20 kPa, 30 kPa, 100 kPa, 200 kPa, 300 kPa, 400 kPa and 500 kPa.Volumetric positive strains appeared only for normal stresses of 10-100 kPa (dilation identified); for remaining rest normal stresses, the negative volumetric stresses appeared.Friction angle at critical state j′ cs = 30.5°derived from the shear failure envelope taking into account only failure points where no volumetric strains appeared and setting cohesion was zero value.Dilation effect was also confirmed by three consolidated triaxial tests performed.
The ground discrete model finite element size was chosen based on sensitivity analysis, resulting in the average finite element size of 1.35 D and introducing the smaller size of 0.4 D in finite elements in vicinity of pile base and side surfaces (where larger ground deformations were expected).FEM analysis via mesh generation to capture local behaviour of sand bear contact zones was also reported by Han et al. (2017) and Comodromos et al. (2003).

Pile ground modelling phases
As mentioned in the introduction, the procedure for modelling of displacement pile (DP) resistance evolution is divided to three main stages of analysis.The first stage corresponds to restoring "at rest" stage; the second stage deals with the ground stress state after DP installation before applying the load (determining the actual DP ground zero state).The third one corresponds to modelling ground stress and strain evolution during DP loading -unloading stages.

Evaluation of "at rest" stage
Ground "at rest" or geostatic stage, existing before pile penetration into ground, is described by relation of effective vertical c V v0 and horizontal c c V V h r 0 0 (as in the considered case it corresponds to radial, perpendicular to radius r circle pile side surface) normal stresses: where K 0 is a coefficient of ground "at rest" state.
The vertical c V v0 stresses in ground are calculated as follows: where ρ = 1.64 g/cm 3 is soil density, g is acceleration of free fall, h is depth of level i below ground surface.The distribution of vertical stresses in laboratory pit of sand deposit is given in Figure 2. Coefficient K 0 depends on soil friction angle at critical state j′ cs and soil deposit formation (consolidation) history, evaluated by coefficient of overconsolidation ratio OCR, the ratio of past (preconsolidation) and of present day (overbunden) stresses.For primary or virgin (no preconsolidation) ground formation path this coefficient is considered as the one for normally consolidated (NC) soil deposit.For this case, coefficient K 0 can be determined as follows (Jaky, 1948): This expression is introduced in many design codes.
When additional external stressing is applied and removed (preconsolidation OCR applies), K 0 corresponds to overconsolidated stage of ground.In this case, the coefficient of ground "at rest" state K 0 can be determined as follows (Mayne & Kulhawy, 1982): It was proved (Broms, 1971;Duncan & Seed, 1986) that compaction of soil induces complimentary residual stresses and creates stress history, i.e., overconsolidation.The artificial sand deposit was formed stage by stage compacting the sand layers (Martinkus et al., 2021).Soil before filling was normally consolidated, but after filling and compaction, it became overconsolidated.In geotechnical engineering practice adopted, overconsalidation is more typical of fine soils.In our case (compacted by layers sand deposit), overconsolidation phenomenon can be recognized by analysing cone penetration test (CPT) graphs, yielding the nonlinear relation of cone penetration resistance (q c ) versus depth (see Figure 1 in Martinkus et al. (2021)).For NC soils q c versus depth is linear (Krasinski & Kusio, 2014).
Evaluation of preconsolidation pressure (stress) is provided below.Aiming to ensure stable properties for artificial sand deposit, the constant initial void ratio e 0 = 0.689 has to be reached per deposit depth.For this purpose, each filled soil layer should be compressed by certain value static stress or vibrating tool (in our case the single direction vibrating plate employed), inducing the certain value vertical stress magnitudes.This stress theoretically is equal to overconsolidation (compaction) stress: (5) The odometer tests, performed with samples (undisturbed, taken from laboratory pit deposit sand) and disturbed (disturbed analysed sand, compacted to reach e 0 = 0.689) resulted in the same behaviour (see average results in Figure 5).Therefore, an assumption that sand deposit consolidation history is analogous to the one of odometer test was applied.Based on this assumption, the overconsolidation stress c V p = 50 kPa was determined using the Casagrande (Casagrande, 1936) method (see c V p in Figure 3).Combining overconsolidation c V p = 50 kPa and vertical c V v0 stress (Equation (2), Figure 2) one can find relation of OCR versus depth (see Figure 4).
Horizontal stresses c V v0 of soil deposit are calculated applying Equation (1) to include coefficient for ground "at rest" state K 0 .Combining relation for horizontal stresses (1) with relevant relation of K 0 for NC (filled, not compacted, Equation ( 3)) and overconsolidated (filled and compacted, i.e., OCR > 1, Equation ( 4)), one can evaluate the distribution of radial stresses per pile depth "at rest" state of compacted sand deposit (see Figure 5, path for c V r0 ) and evaluate qualitatively the influence of compaction (Figure 5, path for c V r0 against path for c V r0 NC , ) as well.Shaft stresses.Pile penetration causes limit shaft stresses.Thus, restoring installation effect, the method for prediction limit shaft stresses can be employed for evaluation of residual stress and strain state referred to DP installation effect.Comparing measurements of ground stress and stress state at state after removal of installation load (contact stresses along shaft and base are in equilibrium with pile weight load), it was found that the best match of the measured radial stress corresponded the ones determined by the semi-empirical Jardine method (Jardine et al., 2005).Applying the method, distribution of horizontal stresses c V rc along depth h is as follows: where p a = 100 kPa is atmosphere pressure, h is depth from soil surface.Using the Jardine method, the distribution of horizontal stresses c V rc was restored in iterative way, by introducing different magnitudes of average volume strains for layers of discrete model.An evaluation of c V rc distribution per pile depth is given in Figure 6.It should be noted that model pile weight of 75 kg gives load 0.75 kN much less the ultimate load F u = 154 kN corresponding to limit state at conditional settlement s u = 0.1 D = 22.5 mm, so in the analysed case it is insignificant for prediction of actual zero state prior to loading.Base stresses.Pile penetration effect to stress state under base was restored via simulation of pile primary loading during the installation process.For case the contact stresses equal to yield stress of magnitude 2225 kPa were assigned on the basis of pile load test graph (see test graphs in Figure 7).Yield stresses below base were mobilized for total pile load of F = 120 kN (being in equilibrium to pile shaft and base contact stress resultant forces F s and F b ).As model pile weight load was insignificant, the pile penetration effect on DP base resistance against loading was generally conditioned by increased stiffness, i.e., due to preconsolidation influence of compacted sand.
When pile ground stress and state prior to DP loading, i.e., the zero state were restored, the incremental DP loading-unloading test could be simulated.

Modelling of displacement pile loading and reloading stage and analysis of the results
Results of modelling the DP ground resistance evolution during loading via graphs of point forces (total load F, pile shaft resisting force F s , pile base resisting force) against pile settlement s are given in Figure 7.
One can find that simulated F s at limit state results in the magnitude that is larger by 12% compared to the one obtained in the pile test.one of s = 0.02 D were obtained experimentally.Parametric analysis results show that within zone of elastic ground deformations the main parameters characterising relation of radial stresses c V r against pile settlement s are soil unloading-reloading stiffness E ur and shear band thickness t s .
The measurements of pull-out DP surface proved that pile-soil interaction corresponds to soil-soil type, resulting in contact friction coefficient δ inter = c M cs = 30.5°.According to Terzaghi et al. (1996) and Craig (2004), the c M cs for similar sand deposit mineral composition analysed varies insignificantly (within bounds of 30-32); therefore, it is reasonable to use constant magnitude of δ inter = c M cs = 30.5°,when determining c M cs = 30.5°via simple direct shear tests.One can find that mobilized limit shaft stresses (developed at s = 0.01 D) depend mainly on residual radial stresses , induced by pile penetration.
Comparing the graphs of F s versus s obtained via numerical simulation and test results, one can state that restoring distribution of radial stresses by applying the semi-empirical Jardine method, one obtains sufficiently exact result for geotechnical design procedures.
One can find that simulated F b corresponding to pile limit state results in the magnitude that is larger by 13% than the one determined experimentally.According to parametric analysis, ground behaviour below pile base mainly depends on overconsolidation stress induced during pile installation, resulting in mobilized value unloading-reloading stiffness E ur (see Figure 8).This parameter prescribes pile base resistance (F b versus s) for short displacement pile, prior to development of large plastic deformation of ground.One must note that definition of F b versus s relation is the most important for geotechnical practice to determine actual pile settlement s when checking its serviceability limit condition.

Conclusions
Numerical modelling techniques for piled (non-displacement and displacement) foundations are under permanent demand.Complexity of characterising initial conditions prior to loading for displacement pile foundation resistance and limit state characterization is mainly based on employment of empiric or semi-empiric relations and relevant adjustment of modelling results with test data.According to normative requirements, a set of partial factors or global safety factor (also extra model factor) is introduced to ensure design reliability as a whole.More exact modelling techniques, based on piled foundation stress and strain evolution analysis in compliance with validation of specified test results, allow minimizing overdesign result, thus correcting the aforementioned factors.The performed numerical modelling of displacement pile resistance analysis results, validated by the analysis of specified load tests, resulted in the following conclusions: 1.The governing resistance measures and peculiarities of displacement pile foundation resistance evolution were evaluated based on the experimental results by applying the proposed method and FEM modelling techniques.Modelling of DP pile resistance evaluation in terms of point forces of external load versus (shaft plus base resistances) was in good agreement with numerical load tests.2. The proposed techniques for qualitative evaluation of stress and strain of ground restitution prior to pile loading and DP ground resistance evolution predictions with the aim of checking ultimate limit condition and tolerable settlement can be employed in geotechnical engineering practice.Hardening soil model applied via PLAXIS serves as an appropriate solution.3. Creating the discrete (layered soil) model for DP resistance analysis, the proper deformation measures of ground along deposit depth should be evaluated, by assigning relevant magnitudes of soil mechanical parameters for soil layers.Pilesoil interaction should be carefully evaluated, as pile ground resistance depends on shear band evolution and its thickness.4. Dilation angle values for soil layers should be set on the basis of stress states (according to the predicted ranges of stress states developed as a solution in pile loading).5. Initial stress state of pile foundation prior to loading should be evaluated during installation of mechanical properties and residual state in pile vicinity accounting for pile self-weight, if it is significant compared to ultimate load.6. Radial stress evaluation using semi-empirical (Jardine) method, vertical "at rest" stresses and actual OCR coefficient restitution give the good agreement with DP test results.7. Base resistance mainly depends on the stresses, induced during pile installation and soil unloading-reloading stiffness E ur .

Figure
Figure 2. Vertical stress versus ground depth

Figure 5 .
Figure 5. Horizontal stress versus ground depth for NC and overconsolidated soils

Figure 6 .
Figure 6.Distribution of radial stresses after pile penetration Figure 7.Comparison of DP ground resistance evolution during loading and unloading stages (load test and numerical simulation)

Figure 8 .
Figure 8. Hyperbolic deviatoric stress vs. axial strain relation in primary loading for a standard drained triaxial test (after PLAXIS, 2016)

Table 1 .
HSM Parameters for discrete model