Experimental and Analytical Modeling of Ground Displacement Induced by Dynamic Compaction in Granular Soils

Abstract

Dynamic compaction (DC) is a ground treatment method that achieves soil densification effects using impact forces. The ground displacement of a crater induced by a hammer is often used for the determination of densification, but less attention has been paid to internal displacement in the ground. To establish an overall understanding of the displacements caused by DC, a laboratory experiment was conducted with sand. The experiment included four energy levels by changing the falling height of the hammer. Meanwhile, a calculation model based on stochastic media theory was proposed to calculate the displacement in the soil. The relationship between the geometric characteristics of the crater and the internal displacement of the soil was established in the model based on the experimental results. The ranges of the relevant parameters were determined, and the feasibility of the calculation model was verified. The model showed good consistency with the experimental data. By selecting the critical settlement, the model could be used to estimate the specific densification scope, including the reinforcement depth and radius. This method can provide a reference for the calculation and optimization of DC.

1. Introduction

Dynamic compaction (DC) is an engineering method that is widely used in ground treatments in many great engineering projects, such as dams and airport construction. Compared with other treatment methods, DC can economically and efficiently improve the bearing capacity and resistance of a ground settlement. This technology can be simply described as dropping a tamper (10~40 t) from a height of 10~20 m to impact the ground. Generally, several blows are required until the ground achieves a specific degree of improvement [1,2]. The DC method can be applied to a wide range of soil types, including collapsible loess, sand, and gravel soil, as well as unsaturated soil [3,4,5,6,7,8]. However, DC was shown to be more suitable for granular soils, such as soil–rock mixtures (SRM) and sands [9,10].

It is challenging to calculate the effects of DC because of its dynamic characteristics. Therefore, recent research on DC has primarily focused on tests and simulations [2,11,12,13]. For engineering design or control aspects, a simple calculation method is still needed. Menard and Broise [9] were the first to propose an equation to calculate the densification depth $D=\sqrt{WH}$ containing two parameters: the weight of the hammer W and the falling height of the hammer H. Mayne et al. [14] added a coefficient n to this equation to adapt it to complex working conditions, and many subsequent studies were conducted to refine the value of this coefficient [2,15,16,17,18]. However, the shortcomings of this equation are obvious in practical applications. For example, the calculation results for the depth produce a large range that cannot be used directly, and the values of the coefficient n vary from 0.2 to 1, leading to large uncertainties in the calculation results of the densification range. In addition, the soil properties are not included, and the horizontal densification radius cannot be determined using this method. From the perspective of engineering, numerical simulation methods are complex for the design and optimization of site construction and difficult for engineers to apply.

In the study of DC, understanding the state change of the soil, including the surface and inner soil of the ground, is essential. The surface changes in the ground include the crater and ground deformation [18,19]. Ground deformation generally occurs when the DC tamping numbers exceed the limit tamping numbers; therefore, DC is often terminated before ground deformation is obvious during actual construction. A crater is formed by the impact of the tamper during DC. In general, crater settlement is used as an evaluation standard for the densification effect [20,21]. For instance, in China’s Technical Code for the Ground Treatment of Buildings (JGJ 79-2012), a settlement following two tampings of < 5 cm is used as an indication to stop tamping. When the settlement change is small, this indicates that a critical value inside the soil has reached the maximum effect under the current energy level; accordingly, it can reflect the change inside the soil to a certain extent. The depth of a crater is often changed with the tamping numbers and energies, and thus, the geometric characteristics of a crater can also reflect the construction information. The inner changes of the ground can therefore reflect the DC densification effect. Multiple studies examined the inner change through aspects of the stress, soil characteristics (standard penetration test (SPT), pressuremeter test (PMT), and density logging), and deformation using field tests, laboratory tests, and numerical simulations [2,22,23,24]. Compared with other parameters, deformation can be observed more easily in soil without using complex equipment. In addition, the occurrence of deformation indicates that the soil state has changed and can therefore be used to characterize the densification of the soil.

In this study, a laboratory experiment was conducted to analyze the displacement induced by DC. A calculation model was proposed in this study based on the displacement of the surface and inner soil. In this model, a relationship between the crater parameters and the internal displacement parameters of the soil was established so that the densification range, including the depth and the horizontal radius, could be determined via this method simply by using the parameters of the crater. This method can provide a reference for the calculation and optimization of DC.

2. Laboratory Experiment Investigating DC

Based on the experimental design in Chen’s work [25], the experimental soil sample was sandy soil, and the main properties of the soil are presented in

Table 1. The Gaussian fitting results for three energy levels.

Type	Dry Density	Moisture Content	Interior Friction Angle	Cohesive Force
Sand	1.5 g/cm3	10%	18°	0

. The gradation curves of the tested sands are presented in Figure 1, and the type of sand was classified as medium sand based on the Standard for Geotechnical Testing method (GB/T 50123-2019). To ensure the uniformity of the soil sample density during sample preparation, 3 cm-thick soil layers were loaded in a 2 m × 2 m reservoir (Figure 2a) and were separated by flexible colored paper, as shown in Figure 2b,c. The paper used in the test was very narrow and soft. During the experiment, a small amount of water was sprayed on the paper, making it easier to damage when it received an external force, and thus, the influence on the experiment was very small. After the layout of the model was completed, the surface settlement was measured during DC. The tamper had a weight of 5 kg and a radius of 10 cm (Figure 2d). The tests were conducted in four groups with falling heights of 10 cm, 20 cm, 40 cm, and 60 cm, which were referred to as T1, T2, T3, and T4, respectively. After each test, each section was excavated, and the settlement characteristics were drawn according to the shape of each layer of flexible colored paper (Figure 2e,f). The position of each layer was measured by using a millimeter ruler and image recognition technology [26,27], and the error was within 1 mm.

3. Results and Analysis

3.1. Relationship between the Tamping Settlement and Tamping Numbers

Figure 3 shows the relationship curves between the tamping settlement and the tamping numbers (i.e., the numbers of tamping blows) for four different tamping energies. Figure 3 shows that the tamping settlement gradually increased with increasing tamping numbers. However, the single tamping settlement gradually decreased with increasing tamping numbers. The settlements of the soil surface during the last four tamping blows were basically unchanged, indicating that the soil had reached a relative equilibrium state given the energy level. As the energy level increased, the tamping amount increased and the number of blows for the tamping convergence amount changed. For example, the tamping amount of the T1 energy level was basically unchanged at the beginning of the 7th blow, while the T4 energy level still has a large settlement value at the 22nd blow. Therefore, the depth of the crater was a manifestation of the energy level, which obviously reflected the different effects of the energy level changes.

3.2. Changes in the Crater with Different Tamping Numbers

Figure 4 shows the change in the crater given different numbers of tamping blows. For the same energy level, as the tamping numbers increased, the depth of the crater increased continuously, and the slope of the crater became steeper with increasing tamping numbers. For different energy levels, when the number of tamping blows was the same, the shape of the crater differed with respect to the slope and the depth of the crater. From the geometry of the crater, the deformation mainly developed downward, and the depth increased with increasing tamping numbers, as shown in Figure 3. The ground heave was not obvious compared with the whole crater. It indicated that most of the soil in the crater was pushed downward. Therefore, from this point of view, different tamping energy levels resulted in different crater shapes and different numbers of tamping blows led to different crater parameters.

3.3. Inner Soil Displacement Induced by DC

Due to the DC treatment, the soil particles were forced to move downward by the impact of the tamper. As shown in Figure 4, the most obvious displacement of each layer occurred along the central axis of the crater. The displacement decreased with the increase in depth and distance from the central axis. The displacements caused by different energy levels were also different. Take the displacement of the 36 cm layer, i.e., the red box in Figure 5a–c, as an example: the maximum displacements of three energy levels were about 0.8, 2, and 3.2 cm, respectively. Therefore, the displacement induced by DC could also reflect the energy levels in terms of the geometry of the crater.

3.4. Analysis of the Experimental Results

In Figure 5, the displacement of each layer was at a maximum in the middle and gradually decreased with the increase in the horizontal distance, which was similar to a Gaussian distribution. In Figure 6, the displacement of the T1 level was fitted using the Gaussian distribution equation, which is simply expressed as

$$y-y_0=a{e}^{-b{\left(x-x_c\right)}^2}$$

(1)

where y₀ and y are the heights before and after the DC of each layer, respectively; x is the horizontal distance of each point; x_c is the horizontal distance of the central axis; and a and b are constant parameters.

The fitting results in Figure 6 were commendably consistent with the experimental data. The fitting results of the three energy levels are presented in

Table 2. The Gaussian fitting results for three energy levels.

Energy Level	Height (cm)	a	b	Coefficient of Determination
T1 (H = 10 cm)	41.6	−1.67	0.025	0.983
	38	−0.81	0.019	0.955
	35.4	−0.64	0.01	0.825
	32	−0.38	0.006	0.93
T2 (H = 20 cm)	39	−2.73	0.017	0.978
	36	−2.13	0.013	0.992
	32.6	−1.768	0.0125	0.952
	30	−1.678	0.01	0.894
	26	−0.51	0.0085	0.81
T3 (H = 40 cm)	36.4	−3.5	0.013	0.967
	33.4	−3.03	0.011	0.971
	29.8	−1.23	0.007	0.868
	26.3	−0.324	0.005	0.78
	22.6	−0.529	0.003	0.92

, which shows a, b, and the fitting coefficient of each layer. The fitting coefficients varied from 0.78 to 0.992, which indicated that the displacement of each layer caused by the DC conformed to the Gaussian distribution.

4. Modeling Analysis of Displacement Induced by DC

Through the analysis of the experimental results, the displacement of each layer caused by the DC conformed to the characteristics of a Gaussian distribution. This phenomenon is similar to the basic assumption of stochastic media theory, which is often used to calculate the subsidence induced by excavation, such as mining and tunneling. As shown in Figure 7a, the excavation induced the movement of the upper medium, which spread gradually as the distance increased. If the probability of movement is 1 for the particle on the excavation, then the movement probability of the two particles immediately above is 1/2 for each. By analogy, the farther the distance, the smaller the probability and the wider the scope. The movement of soil induced by DC can also be expressed using this perspective. The movement of the particles in the crater will push the soil under the crater downward, as shown in Figure 7b, which is similar to Figure 7a. The definition of stochastic media in the theory refers to granular soils, such as sands and soil–rock mixtures. The material in this test was sand, and thus, the displacement induced by DC was analyzed using this theory.

4.1. Computational Model of DC

The displacement caused by the downward movement of a unit cell in the crater can be divided into two parts: vertical and horizontal displacements. The vertical displacement is more convenient to detect than the horizontal displacement. Therefore, in this study, only the change in the vertical displacement was considered and the horizontal displacement was neglected in the subsequent calculation process.

In the Cartesian coordinate system, a unit line segment in the crater is denoted as dx. The downward force from dx will induce the displacement of a layer at a certain depth. The attribution of displacement is similar to the Gaussian attribution results of the experiment, and the displacement is symmetric about the central axis of the crater. Therefore, it is assumed that the deformation probability density function $A\left(x_A,y_A,z_A\right)$ in Figure 8 induced by dx is f(x²) and that the curve of the function is symmetric about the z-axis.

Therefore, the deformation probability of point A induced by dx can be expressed as

$$P\left(dx\right)=f\left(x^2\right)dx,$$

(2)

For an infinitesimal microplane ds = dxdy in the crater, the deformation probability of point A can be expressed as

$$P\left(ds\right)=P\left(dx\right)P\left(dy\right)=f\left(x^2\right)dxf\left(y^2\right)dy=f\left(x^2\right)f\left(y^2\right)ds,$$

(3)

Rotating the coordinate system along the z-axis from (x, y, z) to (x₁, y₁, z), as shown in Figure 8, Equation (3) can be transformed into

$$P\left(ds\right)=P\left(d{x}_1\right)P\left(d{y}_1\right)=f\left({x_1}^2\right)dxf\left({y_1}^2\right)dy=f\left({x_1}^2\right)f\left({y_1}^2\right)ds,$$

(4)

Assuming that the axis x₁ passes through point A in the new coordinate system, then

$${x_1}^2=x^2+y^2,y_1=0,$$

(5)

Therefore,

$$f\left({x_1}^2\right)f\left({y_1}^2\right)ds=f\left(x^2+y^2\right)f\left(0\right)ds=cf\left(x^2+y^2\right)ds,$$

(6)

$$f\left(x^2\right)f\left(y^2\right)=f\left({x_1}^2\right)f\left({y_1}^2\right)=cf\left(x^2+y^2\right),$$

(7)

where c is a constant.

Equation (8) can be obtained by differentiating both sides of Equation (7):

$$\begin{array}{c}f\left(y^2\right)\frac{\partial f\left(x^2\right)}{\partial x^2}=c\frac{\partial f\left(x^2+y^2\right)}{\partial \left(x^2+y^2\right)}\frac{\partial \left(x^2+y^2\right)}{\partial x^2}=c\frac{\partial f\left(x^2+y^2\right)}{\partial \left(x^2+y^2\right)}\\ f\left(x^2\right)\frac{\partial f\left(y^2\right)}{\partial y^2}=c\frac{\partial f\left(x^2+y^2\right)}{\partial \left(x^2+y^2\right)}\frac{\partial \left(x^2+y^2\right)}{\partial y^2}=c\frac{\partial f\left(x^2+y^2\right)}{\partial \left(x^2+y^2\right)}\end{array}\},$$

(8)

Then,

$$\begin{gathered} f\left(y^2\right)\frac{\partial f\left(x^2\right)}{\partial x^2}=f\left(x^2\right)\frac{\partial f\left(y^2\right)}{\partial y^2}. \\ \mathrm{or} \\ \frac{1}{f\left(x^2\right)}\frac{\partial f\left(x^2\right)}{\partial x^2}=\frac{1}{f\left(y^2\right)}\frac{\partial f\left(y^2\right)}{\partial y^2}, \end{gathered}$$

(9)

Equation (9) shows that the calculations on both sides of the equation are independent of x and y.

Let

$$\frac{1}{f\left(x^2\right)}\frac{\partial f\left(x^2\right)}{\partial x^2}={k_1}^2,$$

(10)

The solution to this equation can be obtained through integration:

$$f\left(x^2\right)=k_2{e}^{-{k_1}^2{x}^2},$$

(11)

where k₁ and k₂ are constants.

The deformation probability of point A induced by ds is thus

$$P\left(ds\right)=f\left(x^2\right)f\left(y^2\right)dxdy={k_2}^2{e}^{-{k_1}^2\left(x^2+y^2\right)}dxdy,$$

(12)

The soil mass will be compressed during DC. Assuming that the amount of ground surface heave is negligible during DC, the total deformation should be equal to the compressed volume V_e of the unit cell in the crater after compression. The initial volume V₀ of the unit is assumed to be 1, and the compression coefficient of the soil mass is η₁. Therefore,

$$V_e={\displaystyle {\int }_{+\infty }^{+\infty }{\displaystyle {\int }_{+\infty }^{+\infty }f\left(x^2\right)f\left(y^2\right)dxdy}}={\eta }_1{V}_0={\eta }_1,$$

(13)

By substituting Equation (12) into Equation (13), the following relation can be obtained (for details, please see Appendix A):

$${k_2}^2={\eta }_1\frac{{k_1}^2}{\pi },$$

(14)

Equation (11) can be transformed into

$$f\left(x^2\right)=k_2{e}^{-{k_1}^2{x}^2}=\sqrt{{\eta }_1\frac{{k_1}^2}{\pi }}{e}^{-{k_1}^2{x}^2},$$

(15)

The distribution density W_e of the deformation caused by the unit subsidence in the crater can thus be expressed as

$$W_e=f\left(x^2\right)f\left(y^2\right)={\eta }_1\frac{{k_1}^2}{\pi }{e}^{-{k_1}^2\left(x^2+y^2\right)},$$

(16)

It is assumed that the ground surface is perfectly flat without heaving adjacent to the crater’s outer edge, that the tamper is smooth and round, and that the crater formed by the impact is circular. The DC process can therefore be simplified to an axisymmetric problem. Taking the z-axis as the axis of symmetry, the displacement distribution density W_A of A(x_A, y_A, z_A) in Figure 8 caused by any unit dx_edy_edz_e in the crater can be expressed as

$$W_{\mathrm{A}}={\eta }_1\frac{{k_1}^2}{\pi }\exp \left[-{k_1}^2\left[{\left(x-x_e\right)}^2+{\left(y-y_e\right)}^2\right]\right]d{x}_{e}d{y}_{e}d{z}_e,$$

(17)

The displacement $W\left(x_{\mathrm{A}},y_{\mathrm{A}},z_{\mathrm{A}}\right)$ of point A caused by the whole crater can thus be obtained via integration:

$$W\left(x,y,z\right)={\displaystyle {\iiint }_{\Omega }{\eta }_1\frac{{k_1}^2}{\pi }\exp \left[-{k_1}^2\left[{\left(x-x_e\right)}^2+{\left(y-y_e\right)}^2\right]\right]d{x}_{e}d{y}_{e}d{z}_e},$$

(18)

where Ω represents the area in the crater.

In Equation (18), the compression coefficient η₁ is related to soil properties, which is set to a constant for homogeneous soil. The parameter k₁ can control the distribution of W(x,z), which varies with the change of depth; therefore, k₁ is a parameter that is related to depth z. The distribution after DC is shown in Figure 5 and the shape of the distribution is controlled by a and b in Equation (1). Based on the results in Table 2, the relationship between b and 1/z is linear, as shown in Figure 9, and can be expressed as Equation (19). The coefficient of determination for the fitting line in Figure 9 was 0.843, and thus it can be considered that a linear relationship existed between b and 1/z.

$$b=10.1/z$$

(19)

From the expressions of Equations (1) and (18), b and k₁ have similar effects. Therefore, k₁ is assumed to be linear with 1/z, which can be expressed as

$$k_1={\eta }_2/z$$

(20)

The solution of Equation (18) can be calculated in a rectangular coordinate system. The coordinate system of the crater is converted from a rectangular coordinate system to a cylindrical coordinate system, and the unit representation in the crater is changed from $d{x}_{e}d{y}_{e}d{z}_e$ to $\rho d\theta d\rho d\zeta$. As shown in Figure 10, the geometric parameters of the crater are the crater depth h and the dummy variable of depth value ζ in the range of (0, h). The cross-section of the crater is a trapezoid, and the radii at different depths can be expressed as $\rho \left(\zeta \right)$:

$$\rho \left(\zeta \right)={\rho }_1-\left({\rho }_1-{\rho }_2\right)\zeta /h$$

(21)

where ${\rho }_1$ and ${\rho }_2$ are the radii of the top and bottom of the crater, respectively.

To simplify the calculation, the vertical plane of the x-axis was chosen to be analyzed so that the y coordinates in this plane are 0. Combined with Equation (20), the displacement W(x,z) generated at any depth can be calculated using Equation (22):

$$\begin{gathered} W\left(x,z\right)={\eta }_1{{\displaystyle \int }}_0^h{{\displaystyle \int }}_0^{\rho \left(\zeta \right)}{{\displaystyle \int }}_0^{2\pi }\frac{{k_1}^2}{\pi }\exp \left[-{k_1}^2\left[{\left(x-\rho \cos \theta \right)}^2+{\left(0-\rho \sin \theta \right)}^2\right]\right]\rho d\theta d\rho d\zeta \\ ={\eta }_1{{\displaystyle \int }}_0^h{{\displaystyle \int }}_0^{\rho \left(\zeta \right)}{{\displaystyle \int }}_0^{2\pi }\frac{{{\eta }_2}^2}{\pi {\left(z-\zeta \right)}^2}\exp \left[-\frac{{{\eta }_2}^2\left[{\left(x-\rho \cos \theta \right)}^2+{\left(\rho \sin \theta \right)}^2\right]}{{\left(z-\zeta \right)}^2}\right]\rho d\theta d\rho d\zeta \end{gathered}$$

(22)

The deformation value at a certain depth can thus be calculated using Equation (22), where the unknown quantities are η₁ and η₂.

4.2. Determination of the Parameters

In the model, the parameter η₁ is used to quantify the volume compression of the soil in the crater before and after DC. The density of sand is 1.5 g/cm³ before DC and 1.71 g/cm³ after DC. Thus, η₁ can be calculated using ${\eta }_1=V_{final}/V_{initial}=1.5/1.71=0.88$, where V_initial and V_final are the volumes before and after DC, respectively.

The parameter η₂ is related to the shape of the displacement distribution. With the increase in depth, the peak value decreases and the range expands, and thus, this parameter is mainly used to describe the diffusion effect of the displacement. η₂ can be obtained using the data from the experiment, which can be determined using the least-squares method:

$$R\left({\eta }_2\right)={\displaystyle \sum {\left(W_i-W_i^{*}\right)}^2}$$

(23)

where W_i and $W_i^{*}$ are the calculated and observed displacements, respectively. A smaller R(η₂) indicates that the parameter η₂ is closer to the target value.

The relationship between η₂ and R(η₂) are calculated using Equation (23), which is presented in Figure 11. The results show that the η₂ values of three energy levels were 1.85, 1.7, and 1.9, respectively. The comparison between the tested and calculated results are presented in Figure 12; it shows that the calculated data were very close to the test data, and thus, the feasibility of this method was demonstrated.

5. Discussion

To verify the validity of the calculation results, a T4 model test level was used based on the test results. In the verification process, for convenience, the average value of η₂ was taken to be 1.8. Figure 13 was obtained by introducing two parameters into the calculation of Equation (22). The results showed that the calculated data were very close to the test data, demonstrating the feasibility of the method.

Using this method, the displacement of each position in the soil could be calculated simply from the shape parameters of the crater. In Figure 14, the relationship between the density ratio and the displacement is presented, showing that the density of the sand was improved, that is, increased, with increasing displacement. When the displacement was less than 1.5 cm or greater than 3 cm, the increase in density was not obvious. When it was between 1.5 cm and 3 cm, the density ratio and displacement basically showed a linear relationship (the red box in Figure 14). The displacement reflected the change in the internal state of the soil, that is, the densification effect of the soil. This phenomenon can be used to determine the boundary of the densification area. In other words, the densification range inside the soil can be determined from the corresponding soil displacement. Figure 13 shows a contour map of the T4 level displacement distribution. Three areas—areas I, II, and III—are presented in Figure 13, with boundary displacements of 2 cm, 1.5 cm, and 1 cm, respectively. The distributions of these three areas were bulb-shaped, which was consistent with the rule of the densification scope. If the dry density of the ground needs to be improved to 1.6 g/cm³ after DC, then the corresponding deformation is approximately 2 cm, and the densification area can be determined from area I in Figure 13. The densification depth and radius were approximately 15 cm and 12 cm, respectively, for the T4 energy level. Figure 13 shows that this method can not only calculate the densification depth but also determine the densification width or the entire densification area. Therefore, this method can also be applied to the optimization of DC and the determination of the densification range.

6. Conclusions

In this study, a laboratory experiment was conducted to investigate the displacement of sand induced by DC; the experiment was conducted in four groups with falling heights of 10 cm, 20 cm, 40 cm, and 60 cm, which were referred to as T1, T2, T3, and T4, respectively. A calculation model was established to analyze the surface and inner displacement of the soil. The feasibility of the model for calculating the internal soil settlement was verified using the results of a laboratory experiment. The following conclusions were obtained:

(1)

The shape of the crater was related to the energy level and the tamping number, and the slope and depth of the crater could directly reflect the change in tamping parameters; when the depth of the crater was large, the corresponding internal displacement was also large. Therefore, the geometric parameter of the crater can be seen as the link between the tamping parameters and internal displacement.

(2)

The displacement of soil induced by DC was similar to a Gaussian distribution; therefore, a model based on stochastic media theory was proposed to mathematically describe the displacement. The model contains only two parameters: η₁ for the compression effect of soil and η₂ for the diffusion effect in the ground. The deformation of any point in the soil affected by DC can be determined using these two parameters combined with the parameters of the crater.

(3)

Compared with traditional methods, such as the Menard equation and numerical simulation, this model has the following advantages when calculating the DC densification range. First, the shape of the crater is considered. This reflects the energy level and soil properties of DC, and the shape parameters can be easily obtained through in situ measurements. Second, the calculation model can determine the densification depth and width of DC at the same time, rather than simply calculating the densification depth. Third, the calculation model is relatively simple and can be directly calculated using MATLAB or Python.