Combined DEM-FDM modelling of encased stone column

Combination of the continuum-based numerical methods and the discrete element method (DEM) could be a powerful way of simulating complex problems. This approach benefits from the capabilities of both methods. The main feature of the discrete element method is that the soil grains are considered as individual particles without need to impose any behaviour law in modelling the medium. The limitation of this method is, however, its high computational demand. In continuum based methods, on the other hand, it is impossible to trace micro scale phenomena. According to these facts, combining continuum and discrete methods is an optimal way in approaching geotechnical problems which deal with granular soils. In this approach, the coarse grain zone (medium) is modelled using DEM and the surrounding media are modelled using the continuum methods. Stone columns that are widely used for improving and/or increasing the strength of weak soils could be modelled using this type of coupled simulation. The Coarse aggregates present in the stone column make it appropriate for the coupled modelling. In this paper, the ordinary and encased stone columns have been simulated by combining 2D DEM and finite difference method (FDM). Clump technique was employed to achieve the interlocking of aggregate particles in DEM, and the surrounding cohesive soil was modelled using FDM. The obtained results were validated by the reported experimental results in the literature, indicating that the coupled DEM-FDM method is a robust way to simulate stone columns.


Introduction
A robust way of simulating geotechnical problems is the combination of different numerical methods and utilizing the capabilities and benefits of each method. A new coupling method that has been employed in recent years by researchers for geotechnical modelling is the combination of continuum-based numerical methods -i.e., finite element method (FEM) and/or finite difference method (FDM) and discrete element method (DEM). The constitutive models in the continuum methods are based on the experimental parameters, empirical assumptions and some simplifications that while predicting valuable information of the materials behaviour in macro scale, ignore some important granular materials characteristics in micro scale such as micro rotation, size, shape, interlocking of particles. Moreover, constitutive models are based on statistics and don't make physical sense [1]. Materials in the discrete element method, on the contrary, are considered as individual particles without imposing any behaviour law. This novel discrete-nature numerical method was introduced in 1971 by Cundall [2] to analyse rock mechanics problems. Then it was developed for soil mechanics problems by Cundall and Strack in 1979 [3]. DEM is based on the simple physics laws, whereas the particle contact patterns could change arbitrarily during the simulation. DEM simulation deficiency, however, lies in the numerous particles needed for the realistic simulation of a medium. This causes a high computational demand which limits the number of particles, as DEM needs to detect the contacts and calculate the contact forces at every timestep. Most geotechnical problems are complex in nature and their modelling might be impossible using only DEM or applying it might be computationally nonoptimal. Instead, the expected complex zone in a simulation could be modelled by DEM and the remaining portion could be simulated using the continuum-based methods [4]. The stone column that is widely used to improve soil strength and stiffness is an appropriate problem for combined modelling. Furthermore, the inter-particle micromechanical behaviour and coarse aggregates of the stone columns make them suitable for this type of simulation. The soilstone system could be considered as two separate zones, providing possibility of simulating each zone in this system based on the corresponding material behaviour. The aim of this paper is to simulate ordinary and geotextile encased stone columns in the clay bed using the edge-to-edge DEM-FDM coupling method. The results were validated using the experimental study which showed the robustness of this approach in modelling the complex systems.

Numerical modelling of stone column and surrounding clayey soil
To capture a complete and realistic behaviour of the stone column there is need for three dimensional modelling. 3D DEM modelling, however, needs much more input parameters and calculation power alongside a tougher calibration process. Hence, to overcome these difficulties, 2D DEM modelling was employed in this study. On the other hand, stone column has an inherent axisymmetric behaviour, but particles in 2D DEM are defined as a 2D disk. Therefore, the stone column is actually modelled with the plain strain conditions. This approach is also used by other researchers [5][6][7]. The schematic coupled DEM-FDM model based on an experimental setup as reported by Ghazavi and Afshar [8] is depicted in Fig.1. Clayey soil in the bed is modelled using the finite difference method. The Mohr -Column failure criteria were chosen to represent the surrounding soil behaviour. The input parameters of clay are listed in Table 1. As both the numerical and experimental models are subjected to constant displacement compression loading with a fast displacement rate, therefore the undrained conditions are prevailing.  As shown in Fig.1, the stone column system consists of a stone with 60-mm diameter and 300-mm height, constructed in the clay bed with the dimensions of 1.2×1.2m in plane and 0.9-m height [8]. Based on the previous paragraph, half of this medium is modelled.
The clump technique was applied in DEM to simulate the none-circular particles in order to model the interlocking of particles in stone materials. In this method which has been adopted by many researches for modelling the irregular shaped particles [9][10], each particle has been considered as a set of rigid circular particles with diverse sizes. Grains with sizes varying between 2 and 10 mm were generated in a way that the grain size distribution of the simulated materials is similar to the experimental gradation curve as presented in Fig.2. In this research, the linear contact model has been applied for particle-particle and particle -wall contacts. This type of contact model consists of the inter-particle friction coefficient; linear (none-tension) and shear spring parallel with linear and shear dashpots. Modelling of the geotextile has been done using a row of bonded circular particles with a diameter equal to the geotextile thickness used in experimental test (i.e. 1.8 mm). The linear parallel bond logic was used to bond particles to each other.

Coupling DEM -FDM
Combining DEM and the continuum-based methods has many advantages. First attempts to couple these methods were made in the 1980's to study the mechanical-thermal behaviour of materials at molecular scales [11]. Since then, coupling strategies have been developed extensively and are used in various contexts [12]. One of these strategies is the Edge-to-Edge method that has been employed to couple DEM and FDM in this study. In this method, the common nods and vertices at the joint interface of the continuum and discrete models are defined in a particular way, so that forces and displacements could be transferred correctly between both models [12]. After creating DEM stone column and FDM clay bed models, both were cycled to reach equilibrium state. Then, to couple DEM and FDM, the right wall of the DEM model was divided into thirty segments so that each vertex of these wall segments corresponded to a node on the FDM model. Afterwards, each DEM contact force was divided into two parts based on the contact point location on the wall piece and was transferred to The timestep assumed for both models was identical to ensure the synchronization of the displacements in discrete and continuum zones and proper data transfer between DEM and FDM models. Because of the dynamic formulation of DEM, the timestep was selected small enough to ensure model stability. Due to the full dynamic formulations of DEM, the energy dissipation of interparticle friction is not enough to reach steady state in the model; therefore, extra energy dissipation is needed which is provided by the local dashpots [13]. The local damping adopted was equal to 0.1 in the DEM model for particles and the walls. Fig.3 illustrates the coupled DEM-FDM model of the encased stone column.

Validation of coupled model
Input model parameters (i.e. strength and stiffness parameters) in the continuum based numerical methods are determined by conducting experimental tests. In contrast, micro parameters (i.e. normal and shear stiffness ...) in DEM models could hardly be obtained using the experimental tests and determination of these parameters in a way that the numerical model mimics the experimental model behaviour is a challenge for the DEM simulation. The common approach to determine the micro parameters is back modelling or calibration [14]. At the first stage, by conducting a large number of simulations with different combinations of the micro parameters of the ordinary stone column and comparing the load-settlement response with the experimental one reported by Ghazavi and NazariAfshar [8], a set of micro parameters for the particles and the walls of ordinary stone column were adopted, which is listed in Table.2. Fig. 4 shows plots of the axial vertical loading versus the settlement of ordinary stone column obtained from the coupled model and the laboratory test. Obviously, the numerical and experimental results are in good agreement, and the laboratory trend of the stone column is accurately predicted by the coupled simulation. Inter-particle coefficient of friction, µ 0.4 Contact normal and shear stiffness of wall-particle, kn-wall (N/m)  Fig.5. As a natural behaviour of the ordinary stone column under vertical stress that has been reported in the literature, the lateral bulging which is developed along the column, increases with increase in the settlement of the stone. The maximum lateral deformation of the stone per 20mm settlement is 2.8 mm. Another important observation is that the maximum bulging occurs approximately at a depth equal to 2.5 times the diameter of the column from the top, which is a little deeper than the bulging location usually reported in the literature (i.e. 1.5-2 times the diameter of column). The reason is that in this study, the loading plate diameter is greater than 3 times of the column diameter. Thereby vertical loading at the surrounding clay in the FDM model provides extra confining for the upper length of the stone column and consequently results in a deeper bulging. Fig.6 illustrates the ordinary stone column particles force chains as lines at various settlements ranging from 0 to 20 mm. The contact force line thickness is scaled by the force magnitude. As expected, the particles force network develops deeper in the column by the increase of the column settlement. The maximum magnitude of the contact force is 557 N at settlement of 20 mm. The coordination number for settlements of 0, 5, 10, 15, and 20 mm is 5.34, 6.08, 7.14, 7.66, and 7.91, respectively.

Fig. 5. Deformation of ordinary stone column at various settlements
After calibration of the coupled model for ordinary stone column, the geotextile cover was inserted into the DEM model and the calibration process for geotextile micro parameters was performed. The values of the particles parameters in the encased stone column were equal to those of the ordinary stone column. The parameters selected for the geotextile are listed in Table.3. It must be noted that in order to achieve a comparatively smooth surface of geotextile, and to avoid sticking of the particles to the geotextile, the geotextile particles were placed within each other. Therefore, to maintain the stability of the model and prevent dispersion of the particles, small contact stiffness was selected for the geotextile particles. Fig.7 compares the load -settlement response of the encased stone column for the numerical and experimental results as reported by Ghazavi and NazariAfshar [8]. It was observed that the employed method to model the geotextile in the encased stone column state increased the bearing capacity of the ordinary stone column. Furthermore, the predicted response obtained from the coupled model almost fits the experimental result. A small deviation, however, was observed in the numerical response for both the ordinary and encased states that maybe due to the simplified assumptions in this study. Inter-particle coefficient of friction, µ 0.4 Contact normal and shear stiffness of wall-particle, kn-wall (N/m)

Fig. 6. Force chain of ordinary stone column at various settlements
Deformation of the encased column for various settlements values ranging from 0 to 20 mm is depicted in Fig.8. Previous experimental studies have shown that using geotextile as a cover for the stone column provides a lateral confinement, which increases the bearing capacity of the column and vanishes the lateral bulging. It can be concluded from Fig.7 and Fig.8 that these phenomena are observed by the coupled model. The maximum lateral displacement of the encased column is less than 1 mm.
In Fig.9, the scaled force chains of the encased stone column are illustrated for different loading step settlements of: a) 0 mm, b) 5 mm, c) 10 mm, d) 15 mm, e) 20 mm. The geotextile prevents particles from moving laterally, so the max magnitude of particles contact force in this case is 761 N. The coordination numbers also are; 6.08, 6.94, 7.58, 8.22, and 9.04, respectively, for settlement values ranging from 0 to 20 mm.

Fig. 8. Deformation of encased stone column at various settlements
The yield state of the surrounding soil for the ordinary and encased stone column at settlement of 20 mm is shown in Fig. 10 and Fig. 11, respectively. In the ordinary stone column at the final step of loading, the shear failure of clay occurred near the bulging zone which resulted in establishment of the shear band in clay. In the encased stone column, although bulging almost did not occur in the column, but shear failure occurred in the soil due to compression loading at top of the soil. Using the geotextile as a cover for column reduces the soil failure at the bottom of the stone column.

Conclusions
In this study, the coupled DEM-FDM numerical modelling was employed to simulate single ordinary and encased stone column in a soft clay bed. The Stone column was represented by the clumped irregular shaped particles in DEM and the surrounding clay soil was modelled using FDM. A fast constant strain rate loading was applied to the model which resulted in the governing undrained conditions. The interface between these models was defined specifically, providing proper exchange of data between the models. The numerical simulation results showed that the behaviour of the ordinary and encased stone columns could be accurately captured by the coupled DEM -FDM approach.