Simple Numerical Model to Simulate Penetration Testing in Unsaturated Soils

Cone penetration test in unsaturated sand is modelled numerically using Finite Element Method. Simple elastic-perfectly plastic Mohr-Coulomb constitutive model is modified with an apparent cohesion to incorporate the effect of suction on cone resistance. The Arbitrary Lagrangian-Eulerian (ALE) remeshing algorithm is also implemented to avoid mesh distortion problem due to the large deformation in the soil around the cone tip. The simulated models indicate that the cone resistance was increased consistently under higher suction or lower degree of saturation. Sensitivity analysis investigating the effect of input soil parameters on the cone tip resistance shows that unsaturated soil condition can be adequately modelled by incorporating the apparent cohesion concept. However, updating the soil stiffness by including a suction-dependent effective stress formula in Mohr-Coulomb material model does not influence the cone resistance significantly.


Introduction
Unsaturated soils are available in surficial parts of the earth and above the water table. Dealing with unsaturated soils is inevitable in many geotechnical projects including shallow foundations, fills, embankment, dams, pavements, slopes, and landfills. The behaviour of unsaturated soils is complex and influenced by many factors including stress state, soil type, soil fabrics, mineralogy, density, and most importantly, inter-particle suction stresses within the pore space [1]. For instance, suction stress alters strength and stiffness of geomaterial by increasing the effective stress acting on soil [2,3]. Thus, a suction-dependent soil characterization approach should be implemented in order to interpret the results of laboratory and field investigations in unsaturated soils.
Cone Penetration Testing (CPT) is one of the popular field investigation methods to evaluate soil properties. The measured cone resistance (q c ) and sleeve friction (f s ) can be linked to the physical, mechanical, and hydraulic properties of soils using empirical and analytical relations [4][5][6]. CPT results (q c , f s ) are in fact controlled by soil shear strength, stiffness, compressibility, and permeability [7].
Suction stresses can be significant in unsaturated soil and influence CPT results directly or indirectly by changing the effective stress, and in turn, affecting the cone resistance. Therefore, using correlations that have been developed for dry or saturated soil conditions without considering the effect of degree of saturation may lead to misrepresentation of soil properties and requires further scrutiny.
Experimental studies by Hryciw & Dowding [8] and Lehane et al. [9] showed that cone resistance (q c ) in an unsaturated soil condition can be increased up to two times compared to that of dry and saturated soil conditions under certain circumstances. Also, results of CPTs conducted by Vanapoli & Fathi [10] in a specially designed tank showed that cone resistance and bearing capacity of sand can be influenced by suction. Russell and Khalili [11] also showed that the cone penetration resistance is comparable to the required pressure for expansion of a spherical cavity in an infinite sand, and it could be increased due to presence of suction in an unsaturated condition. Further, CPT results on saturated and unsaturated sands in a suction-controlled calibration chamber conducted by Pournagiazar et al. [12] showed that the effect of suction is more significant for low chamber confining stresses and low relative density values. This emphasizes the importance of suction in cone penetration into the shallow ground where soils are mostly unsaturated.
Numerical modelling of cone penetration process is complex due to large deformations in the model and sever mesh distortion. Several attempts have been made by different researchers to overcome the challenges involved in solution schemes and mesh generation [13][14][15][16]. Recently, Nazem et al., Sheng et al, and Yi et al. [17][18][19] used ALE method to avoid mesh distortion problem in soil. Thus far, this technique combined with adaptive remeshing showed the most stable responses.
Implementing these numerical modelling techniques to simulate cone penetration in unsaturated soil requires further examination as previous works mainly focused on dry or saturated soil conditions. In this paper, a simple numerical model for simulating CPT in unsaturated soil using Abaqus/Explicit is presented. Results of a set of CPT in unsaturated sand inside a model calibration chamber are discussed. In addition, sensitivity analyses were performed to investigate the effect of various soil parameters on the penetration response.

Model Specification
Penetration of a 500 mm long cone with a 16 mm diameter into a sandy soil inside a calibration chamber was modelled using axisymmetric Finite Element (FE) in Abaqus/Explicit. Model geometry, boundary conditions, and mesh configuration are shown in Figure 1. Soil domain is 160 mm-wide and 400 mm-high. The side boundary was located at radial distance of 20 times the cone radius to reduce the boundary effects [20]. The bottom boundary was completely fixed while a constant normal stress boundary, equal to the confining pressure in a calibration chamber, was considered for the top and sides of the model. Two confining pressures of 50 and 100 kPa were considered in this study. In order to model the penetration process, a downward vertical velocity of 250 mm/s was imposed to the cone's top nodes, which were all coupled together. Penetration was considered to begin from the soil surface.
4-node bilinear, axisymmetric quadrilateral reduced integration elements (CAX4R) were used for both soil domain and cone and they were both considered as deformable bodies. To prevent mesh distortion problem in soil elements, the remeshing method, ALE, was used which keeps the mesh size fine enough during penetration and allows independent movement of mesh from material while maintaining the topology (the element and connectivity) of the mesh.

Material model
Simple elastic-perfectly plastic Mohr-Coulomb (MC) model was chosen as the soil constitutive model. The MC model parameters were estimated based on the properties of Sydney sand [21] in Tables 1 and 2. The very small value of cohesion for this sand was to avoid numerical instability in the system. Equation 1 proposed by Bolton [22] was used to update soil friction angle for different confining pressures.
where ∅ ′ is critical state friction angle, is the soil relative density, ′ is the mean effective stress before shearing and is the reference atmospheric pressure.

Interaction model
The cone-soil interaction was considered using a surfaceto-surface contact model with pure master-slave algorithm in tangential direction. Here, the cone was chosen as the master surface while the soil domain was considered to be slave node region. To avoid slave node penetration master surface and make the interaction model work properly, slave node region (the soil elements) should have coarser mesh compared to that of master surface (cone element) [23]. The penalty method with no limit for transmitted shear stress was also used as friction formulation. The friction angle at the cone-soil interface was assumed equal to δ=25º which resulted in an interface friction coefficient of µ=tan =0.43, based on suggested value of soil-pile friction angle for medium sand (API, [24]). "Finite Sliding Algorithm" was also used for tracking contact nodes in which the nodes and elements paired between master and slave surfaces are tracked and kept linked until the calculation stops.

Remeshing method
The mesh distortion due to the large deformation around the cone could be avoided by implementing very fine mesh in the soil zone near the cone and using an adaptive remeshing technique. ALE remeshing technique is an efficient method to control the mesh distortion problem when large none-recoverable deformation occurs. It combines the features of pure Lagrangian analysis and pure Eulerian analysis. The ALE method allows independent movement of mesh form material while it doesn't alter the topology (element and connectivity) of the mesh. In this study, ALE method with 'improve aspect ratio' remeshing algorithm was used for soil elements. The remeshing frequency was set to be at every time step to improve the quality of results unless it involves high computational cost.

Unsaturated Soil Modelling
In order to capture the deformation behaviour of soils in unsaturated soil condition, a simple modification was applied to Mohr-Coulomb model. In this method, suction stress was incorporated in Bishop's effective stress equation (i.e. Eq. 2).
where ' ' is the effective stress parameter, which is a function of degree of saturation and varies between zero and unity, '( − )', also known as 'ψ', is the matric suction (i.e. the difference between pore air and pore water pressure), and '( − )' is the net normal stress.
( − ), is called suction stress and is the contribution of matric suction to effective stress [25]. Then, this modified effective stress formula was used to estimate the shear strength, as in Eq. 3. where ' ′ ' is effective soil cohesion and ′ is the effective angle of internal friction.
The extra shear strength gained from the addition of suction stress called "apparent cohesion" (i.e. " = ′ + ( − ) tan ′ ) was regarded as an input cohesion in the material constitutive model. The effective stress parameters, χ, was estimated using the empirical equation where ' ' is the suction value separating saturated from unsaturated states, i.e. the air entry value. This expression was found based on laboratory strength and volume change data for Sydney quartz sand. Considering equal to 7 kPa for Sydney sand, the effective stress parameter, , and the apparent cohesion c'', and friction angle values for each suction level under different confining pressures can be calculated using Eqs. 1, 3, and 4, which are shown in table 3.

Results
In the following sections, the results of numerical model CPTs for different saturation conditions and confining pressures are presented.

Cone resistance
Cone resistance, q c , profiles for CPT on saturated and unsaturated sands are shown in Figures 2(a) and 2(b) for confining pressures of 50 and 100 kPa, respectively. All models had an initial relative density of 61% and suction values were 25, 100, and 200 kPa. The numerically estimated cone resistance in depth was increased gradually until it gets to steadily constant value, where the simulation was stopped. This is consistent with observations from calibration chamber testing. The results indicated that suction has significantly influenced and increased q c comparing with the saturated soil condition. However, the effect of suction is more noticeable in lower confining pressures, which is evident where suction contribution to the effective stress and overall resistance is more pronounced. For example, for suction equal to 25 kPa, the q c was increased by 20 % and 30 % for confining pressures of 50 kPa and 100 kPa, respectively, while for suction equal to 200 kPa, it only increased by 9% and 26 %, respectively.

Stress distribution
The distribution contours of stress components after 50, and 150 mm of penetration are shown in Figure 3. It indicates that the vertical and horizontal stresses around the cone tip were increased greatly due to penetration and were higher near the upper end of the cone than near the tip. This shifting of stresses to the upper end of the tip has likely been resulted from the frictional interface between the cone and the sand and the cone corner geometry, which led to stress built-up [14].
In order to evaluate stress distribution in different degrees of saturation, generated horizontal and vertical stresses versus distance from the cone tip during the penetration through the depth of 0.15 m for different confining pressures and saturation conditions were obtained and presented in Figure 4 and 5 respectively. For both confining pressures, there is no change in horizontal stress at and beyond a distance approximately 6 times the cone radius. This distance is the same for CPTs on unsaturated sand with different suction levels and under the two confining pressures. Presence of suction in unsaturated sand increased the induced horizontal stress inside the zone of influence while the stress values converged outside this zone. However, for vertical stress, the stress influenced zone was about 2 times the cone radius, which is less than the one for horizontal stress.

Sensitivity Analysis
The suction stress in unsaturated soil affects the effective stress in the soil. This may change other soil properties such as friction angle and soil stiffness. However, we only considered the suction effects in an apparent cohesion, which is a strength parameter in MC constitutive model. In order to investigate the significance of carrying the other parameters, numerical CPTs were performed on unsaturated sands while the elastic modulus and the friction angle were altered and the cone resistance results were compared with the case with a confining pressure of 50 kPa and suction value of 200 kPa.
The suction-dependent effective stress was calculated using Equation 4. Then, the friction angle was modified using Equation 1 by considering the updated effective stress for unsaturated sand. The elastic modulus was also modified considering a typical square root of effective stress variation profile, as follows  The modified material parameters are summarized in Table 4. The cone resistance profiles resulted from these alternate parameters were obtained and compared with the one from the reference case, shown in Figure 6. Both plots demonstrate a slight change in cone resistance that could be ignored. The small change in the friction angle is due to the very low sensitivity of the friction angle to the induced amount of increase in the effective stress, otherwise friction angle, which is a strength parameter, play a major role in the cone penetration response modeled using MC constitutive behavior. On the other hand, the negligible effect of elastic modulus on the cone resistance is attributed to the fact that the soil during the penetration process becomes plastic. Thus, for a perfectly plastic model like MC the effect of elastic modulus on the cone resistance turns out to be insignificant.

Conclusions
A numerical approach was successfully implemented to investigate the cone penetration process in unsaturated sand inside a soil calibration chamber. The application of ALE and adaptive remeshing in Abaqus/Explicit avoided the mesh distortion in this large deformation problem. Apparent cohesion concept was employed in suctiondependent strength characterization based on Mohr Coulomb constative behaviour. Suction in unsaturated soil influenced the cone resistance significantly while the apparent cohesion efficiently was able to capture this effect. However, modifying the soil elastic modulus and friction angle based on the effective stress equation in unsaturated soil did not affect the cone resistance profiles partly due to the choice of the constitutive model and partly due to the low sensitivity of the parameters to the range of the varied effective stress.
Although the simple frictional model implemented in this study was relatively successful, it can be improved to better simulate the penetration. These include but not limited to using a stress-dependent interaction model, considering the effect of soil compressibility during the penetration, or using more sophisticated elastoplastic constitutive models.