Investigation of the dynamic response of subgrade vibration compaction based on the finite element method

1.   Introduction

Compaction is essential during road construction, affecting the pavement's service quality and durability. Quality control of the compaction can improve the strength and stiffness of the pavement material. Hence, the compacted material's overall stability and impermeability can be significantly improved. This would avoid excessive plastic deformation or strength loss of the road base pavement material under cyclic wheel loads and external environmental conditions such as temperature, scouring, freezing, and thawing ^[1]. Field construction data shows that the proper compaction of road materials may ensure the strength and stability of the pavement structure ^[2].

Accurate, fast, and comprehensive compaction measurement is crucial to ensure compaction quality ^[3,4,5]. However, the traditional methods mainly rely on field sampling tests (such as the cutting ring method, sand replacement method and core-drilling inspection), which are "point control" and "delay control" ^[6,7,8]. Achieving timely and comprehensive compaction quality control is challenging because the results are only available after compaction. Many scholars have studied intelligent compaction technology to overcome the shortcomings of traditional methods. This technology requires real-time feedback on the degree of compaction of the compacted material through the use of a continuous detection device to adjust the working parameters of the roller according to the feedback information.

Intelligent compaction technology is based on the monitoring and use of dynamic responses. Meehan et al. ^[9] calculated the compaction index CMV based on the measured acceleration signal and analyzed the correlation between this index and the commonly used in-situ test results. White et al. ^[10] transformed the dynamic responses of the roller into the compaction meter value and machine drive power, and carried out regression analysis between these roller-integrated compaction values and the in-situ measurements. The results showed statistical significance relationships. Sivagnanasuntharam et al. ^[11] summarized the correlation coefficients between different roller-integrated compaction values and spot test measurements (obtained from various intelligent compaction demonstrations for asphalt pavements), but the results were not consistent. Shi et al. ^[12] believed that the existing compaction indexes are not suitable for evaluating the compaction quality of coarse-grained geomaterials, so they proposed a new index called ECP, which needs to monitor the responses of roller travel, vibration and transmission system. Vennapusa et al. ^[13] studied the method of reflecting compaction uniformity with roller-integrated values and analyzed its feasibility. To sum up, previous studies mainly focused on transforming dynamic responses into compaction indexes and studying their application effect, while at the mechanism level, the interaction of the "vibrating wheel-soil" system ^[14] and the resulting dynamic responses have not been well studied ^[15,16]. Therefore, this paper investigates the wheel-soil interaction though the use of the finite element method. The dynamic stress-strain response is analyzed to reflect the response characteristics and mechanism of the soil under compaction conditions.

2.   Development of finite element model

The mechanism of the internal interaction of the system "vibrating wheel-soil" is extremely complicated, while accurate prediction is required in engineering practices. Therefore, it is difficult to meet the practical needs by using a simple dynamic model analysis method. In this paper, based on summary of previous research, a finite element method is used to establish a three-dimensional model of the rolling condition of the vibrating wheel to explore the interaction between the compacted material and the vibrating wheel. The stress state, plastic deformation of the soil body, and other characteristics are analyzed. The finite element model development is presented in the next subsections.

2.1.   Parameters of the wheel model

The roller model is usually simplified to a single vibratory wheel in the simulation for vibratory compaction. Due to a rubber damper between the steel wheel and the frame results in relatively small frame vibrations, its inertial effects can be ignored, and only its gravitational static load assigned to the steel wheel is considered. However, for the vibrating wheel, it is necessary to model it explicitly to reflect its vibratory characteristics. In this paper, the rigid-body assumption is adopted to describe the vibrating wheel, and the relevant working parameters are designed with reference to the XS223JS single-steel wheel roller of XCMG, as shown in Table 1.

In Table 1, some parameters are rounded. In addition, the centrifugal force and vibratory frequency correspond to the operating mode of low frequency and high amplitude to reflect the working condition of high compaction efficiency in the early stage.

In this model, the horizontal component of the excitation force is not considered, which is consistent with the working condition of the vertical vibratory roller. In addition, it should be noted that in terms of the realization of the steel wheel drive, the single steel wheel vibratory rollers on the market can be mainly divided into mechanically driven and fully hydraulic driven types. The former is not equipped with a travel drive system on the steel wheel that needs to generate rolling force through the friction between the steel wheel and the soil and the traction of the frame to form a combined torque. The latter is also equipped with a drive system on the front wheel, which can provide rolling force coupling for the steel wheel by itself and no longer depends on the frictional resistance of the soil body. In contrast, the latter has significantly better performance than the former, because the horizontal force of the steel wheel on the soil produces a horizontal pushing effect on the front soil, which causes the arching phenomenon of the front soil. This phenomenon is more obvious in the finite element, because the roadbed model in the finite element analysis is a continuous slab structure, and its ability to transfer shear stress is much stronger than that of the discontinuous loose soil in the actual project, so the bulge deformation is more obvious.

Moreover, the internal structure of the uplifted soil in the actual project has mostly been destroyed, and the strength is not high, so it can be easily flattened again. But, in the finite element method, it is difficult to reflect this change in the intrinsic model of the material. Therefore, during the analysis, the stiffness and strength of the uplifted soil units due to shear stress do not actually change when compared with those before uplifting, so the compaction work will be more than that of the actual project. To weaken the impact of this effect and reduce the gap with the actual project, the steel wheel is simulated in this model with a matching rolling angular velocity, along with the translational velocity, to simulate the travel drive system on the wheel.

2.2.   Establishment of finite element model

2.2.1.   Constitutive model and parameters of subgrade soil

The soil model established in this paper is divided into two layers: upper and lower. The upper layer is the soil to be compacted, and the linear Drucker-Prager model (D-P model) is used. The lower layer is compacted soil, and a linear elastic model is used. The linear D-P model in ABAQUS is an extension of the classical D-P model, whose projection on the meridian plane remains linear. Still, the projection on the π-plane is no longer circular, so it can reflect the properties of geotechnical materials with different tensile and compressive strengths (Figure 1). The model has simple parameters and is highly adaptable to geotechnical materials, so it is widely used in the simulation analysis for geotechnical materials ^[17,18,19].

The expression for the yield criterion in the meridional plane for the linear D-P model is as follows ^[20]:

where F is the yield surface function, and t is the deviatoric stress parameter:

where p, q, r are the three stress invariants; d and β denote the cohesion and friction angle of the D-P model, respectively; K is the ratio of the triaxial tensile yield stress to the triaxial compressive yield stress of the material, which controls the shape of the yield surface. When K = 1, the yield trajectory of the linear D-P model in the π-plane is the Mises circle.

The flow law of the linear D-P model is related to the dilation angle $\psi$. The model assumes that the plastic flow direction of the material in the $p-t$ stress plane is at an angle $\psi$ to the t axis, so its plastic potential surface function is expressed as follows.

where G is the plastic potential surface function. In general, $\psi \ne \beta$ is taken to correspond to the non-associative flow law. In addition, the linear D-P model adopts the hardening law of isotropic reinforcement. But, in this model, the hardening properties of the material are ignored.

In summary, the linear D-P model has six parameters to be determined, which are the soil elastic modulus E, Poisson's ratio μ, density ρ, yield stress ${\sigma }_{c}^{0}$, friction angle β and shape parameter K (the dilation effect is ignored, so the dilation angle is set to zero). By comparison, the linear elastic model is simpler and contains only the first three parameters. To ensure the accuracy of the finite element analysis results, it is necessary to take reasonable values for the above parameters. In this paper, the values are selected with reference to the research data of Li ^[21]. The elastic moduli of the upper and lower soil layers are 20.4 and 50.0 MPa, Poisson's ratios are 0.20 and 0.30, and the soil densities are 1750 and 2050 kg/m³, respectively. The yield stress and friction angle in the linear D-P model can be generally determined by triaxial tests. However, in the absence of experimental data, the parameters can also be transformed by the cohesive force and friction angle of the Mohr-Coulomb model. Under the condition of triaxial compression, the parameter matching relationship between the two model is as follows:

where $c$ and $\varphi$ are the cohesion and friction angle in the Mohr-Coulomb model.

In the linear D-P model, $K\ge 0.778$ is required to ensure the yield surface convexity. According to the requirements of Eq (6), the internal friction angle $\varphi \le {22}^{\circ }$. However, many materials' actual internal friction angles are greater than that value. Currently, the calculation is usually performed by using $K = 0.778$. However, if the actual friction angle is much larger than ${22}^{\circ }$, the parameters of the D-P model obtained through this method do not fit well with the parameters of the Mohr-Coulomb model. Therefore, this model only applies to materials with a small friction angle. In addition, the soil's cohesion and friction angle will be attenuated under the effect of the vibratory load. Sun ^[22] pointed out in his study that due to the effect of vibration, the friction angle of soil generally decreases by ${1}^{\circ }~{3}^{\circ }$, and the cohesion can be attenuated by more than 50%. The attenuation is more serious when the soil is soft. Considering the above factors, the final cohesive force of 18.0 kPa and the internal friction angle of ${22.3}^{\circ }$ is plugged into Eqs (4) and (5) for parameter conversion to obtain the linear D-P model yield stress and friction angle. The final values of the soil parameters are organized as shown in Table 2.

2.2.2.   Mesh size and artificial boundary

The mesh size of finite elements significantly impacts the computational accuracy, efficiency, and convergence, so it must be carefully considered when selecting the mesh size. For dynamics problems, the discretization of the continuum produces two adverse effects: "low-pass" and "dispersion." The low-pass effect refers to a cutoff frequency in the discrete structure, above which the vibrational components cannot propagate. The dispersion effect refers to the fact that the spectrum of the discrete system is narrower than that of the continuous structure, resulting in a weaker concentration of vibrational energy in time and space ^[23].

Both effects are related to the size of the grid. The denser the element division and the closer to the continuum, the weaker the adverse effect of discretization. According to the research results of Zong ^[23], and the grid size should satisfy the following requirements:

where $l$ is the maximum grid size; ${\lambda }_{min}$ is the minimum wavelength selected, ${\lambda }_{min} = {v}_{s}/{f}_{max}$; ${v}_{s}$ is the shear wave velocity (${v}_{s} = \sqrt{G/\rho }$); ${f}_{max}$ is the maximum intercepted frequency.

Using the basic parameters of the soil listed in Table 2, and then, according to the conversion relationship between the moduli, the shear modulus of the upper soil $G = 8.5$ MPa can be calculated, from which the shear wave velocity can be obtained as 69.1 m/s. The vibratory frequency applied in this model is 30 Hz, and the highest fluctuation frequency of 120 Hz is employed according to the relevant literature, so the minimum wavelength ${\lambda }_{min} = 0.576$ m. Substituting the wavelength into Eq (7), we can get $l\le 0.07$ m. On this basis, the final grid size of the soil is taken as $l = 0.05$ m. The grid size of the steel wheel is equivalent to the soil, which is designed to avoid the large and unreasonable bulge deformation of the soil elements.

Artificial boundary processing is mainly purposed to solve the reflection problem of vibratory waves. The road foundation is a semi-infinite spatial body, but only a finite part can be intercepted for calculation in the finite element analysis. As a result, the outgoing waves in the calculation area cannot propagate to infinity through the artificially intercepted boundary. Instead, they are reflected to the calculation area, which affects its accuracy. Therefore, when performing foundation dynamics analysis, the intercepted artificial boundary must be specially processed to fit the boundary conditions of the infinite domain in actual project.

This model adopts the idea of a viscoelastic artificial boundary ^[24,25], which adds damping at the boundary nodes to absorb the energy of vibratory waves and adds springs to ensure the model's elastic recovery and avoid the model's overall drift phenomenon (Figure 2) ^[26]. In addition, to facilitate the implementation of finite elements, the method of replacing the spring-damper elements at the nodes with equivalent solid cells ^{[27,28,29,30,31]} is adopted. This method extends another layer of hexahedral cells at the model boundary normal to the edge interface. These cells' stiffness and damping matrixes are used to equate the spring-damping elements distributed at the boundary. The moduli and damping coefficients of the equivalent cells are calculated as follows:

where G is the material's shear modulus, $\stackrel{~}{G}$ and $\stackrel{~}{E}$ are the shear modulus and elastic modulus of equivalent boundary cells, $\stackrel{~}{\mu }$ and $\stackrel{~}{\eta }$ are the Poisson's ratio and damping coefficient of equivalent elements, h is the equivalent cell thickness, R is the distance from the wave source to the artificial boundary, ${\alpha }_{T}$ and ${\alpha }_{N}$ are the tangential and normal viscoelastic artificial boundary correction factors, ${c}_{s}$ and ${c}_{p}$ are the wave velocities of S and P waves, respectively.

where ρ is the material's density, E and μ are the material's elastic modulus and Poisson's ratio, respectively.

The finite element model of the soil established has a size of $10\times 6\times \left(0.5+0.3\right)$ m, where the thickness of the upper layer is 0.3 m and that of the lower layer is 0.5 m. There are five kinds, and nine artificial boundaries need to be processed (the relative lateral boundaries belong to the same kind). The R-values of the boundaries were determined based on the principle of the shortest perpendicular distance. According to the basic parameters of the soil provided in Table 2, the elastic moduli and damping coefficients of the artificial boundary cells were calculated. The results are shown in Tables 3 and 4.

Based on the above description, the modeling of the vibrating wheel and the soil can be completed. However, the contact properties need to be set to establish the connection between them. In this model, the normal contact property is defined as a "hard contact", meaning that the contact force only exists when the contact gap is zero. The contact properties in the tangential direction were determined using the Coulomb friction model ^[32] with an empirical friction coefficient of 0.3. The contact calculation method uses the penalty function method. This method allows a small penetration of the master surface nodes to the slave surface. And, it induces a penalty function proportional to the penetration depth as the interface contact force. The final completed finite element model is shown in Figure 3.

3.   Analysis of the dynamic response of vibratory compaction

3.1.   Vertical stress distribution characteristics of the soil under the wheel

The vertical stress of the soil under the wheel is the main reason for the compacted process, and its value can reflect the compaction effect on the soil to a certain extent. Generally, in the early stage of compaction, the greater the vertical stress and the deeper the transmission, the better the compaction effect. In the following, the stress distribution of the soil under the wheel is analyzed according to the finite element simulation data.

3.1.1.   Distribution of vertical stresses of the soil along the transverse direction of the wheel

Figure 4 gives a cloud diagram of the vertical stress distribution of the soil under the wheel along the transverse direction of the vibrating wheel. It can be found that the stress distribution of the soil under the wheel is not uniform along the transverse direction of the vibrating wheel, but presents a "saddle shape" with large sides and a small middle. This phenomenon is more obvious near the surface of the soil body. As the stress transfer deepens, the depressing part in the middle of the stress contour gradually filled up and finally leveled with the two sides, which was because the stiffness of the vibrating wheel is much larger than that of the soil. No bending deformation will be produced during compaction, so the soil under the wheel must produce a uniform settlement to ensure contact with the steel wheel. To make the settlement uniform, the soil on both sides must bear more pressure, thus producing a "saddle-shaped" stress contour.

Accordingly, a simple speculation about the compaction effect on the soil under the wheel can be obtained. In the surface layer of the soil, due to the saddle-shaped stress distribution, the compression density of the soil on both sides was higher than that in the middle. In the deeper layer of the soil body, on the contrary, the vertical stress of the soil body on both sides was lower than the middle, so the compression deformation was smaller. This may allow a uniform cumulative settlement in the surface layer of the soil. At the same time, a transitional intermediate layer tends to form between the two, where the vertical stress and compression deformation are uniformly distributed, and the compaction effect is optimal, which is consistent with the actual engineering experience.

In addition, the soil sample simulated in this model is cohesive soil. If the compacted material is sandy soil, the shape of the stress contour may be opposite to that of the cohesive soil, which means that a parabolic shape with a high center and low sides. The sandy soil has no cohesive force, so the sand particles at the edge of the steel wheel can easily extrude to the sides, thus transferring part of the stress they should bear to the middle. This could also lead to poor compaction on both sides of the sandy soil, which is treated in actual engineering practice by overlapping the rolling width.

For cohesive soils, the overlapping of the rolling width is carried out in practice in order to ensure the compaction effect on the soil at the edge of the vibrating wheel. However, the value of the overlap width is determined by engineering experience and lacks theoretical guidance. This section analyzes the reasonable value of the overlap width from the perspective of soil stress distribution along the transverse side of the drum.

Figure 5 provides the distribution of peak vertical stresses at different depths along the transverse direction of the steel wheel. At 0.80-1.0 m from the center of the wheel, the vertical stresses at different depths are significantly reduced, and the curve's inflection point keeps approaching the center of the wheel as the depth deepens. For example, at the soil depth of 5 cm, when the distance from the center of the wheel reaches 0.90 m, the discounting of the curve begins to increase significantly. The vertical stress of the soil at this point is 96.1% of the vertical stress at the center. However, at 0.95 m, the vertical stress of the soil decays to 85.5% of the center value, and at 1.0 m, the vertical stress is only 51.7% of the center value. Similarly, at the soil depth of 15 cm, the curve's inflection point appears at 0.85 m from the center of the wheel, and the vertical stress amplitude at this point accounts for 95.4% of the central value. By advancing 0.05 m towards the edge, the stress amplitude decays to 84.1% of the center value. At the edge of the drum, the stress amplitude is only 45.8% of the center value. Finally, at the bottom of the soil with a depth of 25 cm, the inflection point of the curve shifts to the left at 0.80 m. The stress amplitude at this point accounts for 94.2% of the center value, and the vertical stresses decay to 84.3, 71.6, 57.9 and even 43.9% of the center value in 0.5 m steps from this point to the edge. In this paper, the stress curve at a depth of 25 cm is selected as the reference base (to be conservative) when determining the overlap width of the rolling. The overlap width is about 20 cm, i.e., about 1/5 of the half-width of the vibrating wheel. This value can be used as a reference for determining the overlap width in actual projects.

3.1.2.   Distribution of vertical stresses of the soil along the longitudinal direction of the compaction path

The cloud diagram of the vertical stress distribution of the soil along the longitudinal direction of the compaction path is shown in Figure 6. It can be found that the center of the vertical stress distribution of the soil under the wheel has an offset from the mass center of the drum. In this model, the difference between the two was about 4 cm. This phenomenon was because the contact surface between the soil and the vibrating wheel is not horizontal. Due to the sinking of the wheel and the uplifting of the soil in front of it, the actual contact surface between the two formed an inclination in front, and the force on the soil located on this inclination is greater than that on the soil directly below the drum.

The above findings can illustrate two issues. First, the compressive deformation of the soil under the wheel is not only the result of the vertical stress but also the effect the horizontal shear stress. However, this point has not received much attention in the previous research. Researchers often ignore the effect of horizontal forces when conducting the analysis of mathematical models, or simply attribute the horizontal forces to the friction between the wheel and the soil ^[33,34]. However, the actual stress state of the soil is much more complex than the assumptions researchers make, which is one of the reasons why the existing theoretical models may not be applicable to engineering practices. The influence of this effect needs to be considered in the subsequent research work.

Second, these findings indicate that the horizontal action force is unavoidable in the vibratory compaction process. The traditional compaction theory considers that the horizontal force between the wheel and the soil is formed by friction. And, the horizontal force had a pushing effect on the front soil, which causes the front soil to bulge and affects the compaction quality. To weaken this effect, a hydraulic steel wheel roller was used to replace the traditional mechanical roller. The steel wheel was equipped with driving equipment, which can weaken the effect of frictional resistance, and even change the direction of frictional resistance. Still, the horizontal force of the soil in front of the vibrating wheel is unavoidable, because the horizontal force, in addition to the frictional resistance, includes the tangential component of the contact force, which exists together with the contact force and cannot be eliminated. Due to the small inclination of the soil body in front of the wheel, the tangential component of the contact force is not large. However, since the soil is a friction-type material, the value of its shear strength is related to the vertical pressure. The vertical pressure of the soil in front of the vibrating wheel was very small, and its shear strength was not high enough to resist the tangential component. Therefore, there was still a sizable bulge in front of the wheel. This phenomenon is detrimental to the improvement of compaction degree and compaction uniformity, so, in practice, pre-compaction is required before vibratory compaction to improve the initial strength of the material and reduce the sinking of the wheel.

To further investigate the influencing factors of the soil uplift effect in front of the wheel, the relationship between the soil uplift height and the forward speed of the wheel is summarized as shown in Figure 7 (the uplift height is calculated as the average value of the maximum uplift height of the six nodes in the center area of the subgrade surface). With the increased driving speed, the uplift effect of the soil in front of the steel wheel obviously weakens. This is because the increase in driving speed will shorten the time that the soil is subjected to the action of the vibrating wheel, then the accumulation time of soil uplift deformation will naturally be shortened. However, although increasing the driving speed contributes to overcoming the bulging effect of the soil in front, it may not benefit the overall compaction effect, because the settlement of the soil under the wheel will also be reduced, resulting in a lower compaction increment in a single pass. In practice, how to adjust the working parameters of the roller, so that the uplift of the soil in front of the steel wheel and the settlement of the soil under the wheel can achieve the most optimal balance, still needs further research.

3.2.   Contact analysis of "vibrating wheel-soil"

The contact analysis of "vibrating wheel-soil" can be used to obtain some parameters in the contact action, which are essential for the theoretical analysis of intelligent compaction and the calculation of mathematical models. There needs to be more experimental data on the contact parameters of vibratory compaction, but these parameters can be extracted more easily through the finite element method. Three units were chosen on the central axis of the compaction path of the roller, and the contact stress data of the nodes of the three elements were extracted. The soil surface nodes' contact stress time history curve was obtained as shown in Figure 8.

According to Figure 8, it can be found that the maximum value of contact stress at this time is about 460 kPa. With the roller's vibration, the contact pressure's value also fluctuates. When the roller was near the node (the time when the forced vibration is the most intense), the fluctuation range of the contact stress was about 100 kPa, and the stress value at the trough was about 350 kPa. In addition, observing the distribution of the contact stress time history curve, it can be found that the curve is not a symmetric structure with the horizontal coordinate of the peak as the axis. The time duration on the left side of the peak is slightly shorter than the time duration on the right side. Moreover, when observing the three waveforms in the middle of the time course curve, the peak value of the left waveform is generally higher than that of the right waveform. This is caused by the tilt of the contact surface in front of the wheel. The tilt of the contact surface causes the center point of the contact stress to be shifted to the front of the wheel so that the nodal stress reaches its peak earlier than the time when the mass center of the vibrating wheel reaches directly above the mode, which leads to the asymmetry between the left and right sides of the curve. The higher peak of the left waveform is also due to the fact that the inclined surface is subjected to a more significant contact stress. At the moment of the left waveform peak, the node is located on the micro plane with a higher inclination, while at the moment of the right waveform peak, the inclination of the micro plane has dropped a lot, so the contact stress is less.

In addition, the grounding width of the vibrating wheel can be estimated according to Figure 8 to make up for the shortage of experimental data. It can be found that the time interval between the node starting to bear the contact pressure and the force disappearing was about 0.27 s, which meant the contact time between the node and the wheel. The forward speed of the drum in this model was 0.8 m/s so that the grounding width could be calculated as 22 cm. However, it should be noted that although the concept is clear, the total grounding width may not be very useful in the mechanical analysis because the contact pressure is not uniformly distributed within the total width. This is demonstrated in Figure 8, where the shape of the contact stress curve can reflect the stress variation on the time scale and be taken as the distribution of stress on the spatial scale to some extent. Therefore, according to this figure, it is easy to find that the contact stress decays rapidly in the lateral direction, which means that only some of the contact stresses effectively force the soil to produce strong forced vibrations and yield within the grounded width. In general, only the contact stresses in the middle three waveforms significantly affect the soil. In the other contact areas, the compressive deformation of the soil is not dominant. Instead, the uplift deformation due to shear stresses is more pronounced. Therefore, it is recommended to define the width range corresponding to the middle three waveforms as the effective contact width, which characterizes the compaction capacity of the vibratory roller along the compaction path direction. In the subsequent "vibratory wheel-soil" contact analysis, it is suggested to consider the effective vertical stresses in this width range rather than the average contact stress in the total grounding width.

The calculation of the effective width can be performed by assuming effective contact stresses. Since only the contact stress of the middle three waveforms during vibratory compaction has a significant effect on the vertical deformation of the soil, it is reasonable to assume the valley value of the middle three waveforms, i.e., 350 kPa, as the minimum effective contact stress for this model. According to the data in Figure 8, the time occupied by the interval with the stress greater than 350 kPa can be obtained as about 0.09 s. The value was converted into the distance traveled by the wheel in space, which was calculated as 7 cm, i.e., the effective contact width is about 7 cm, accounting for about a third of the total grounding width.

4.   Conclusions

By establishing a three-dimensional finite element model of vibratory compaction, the elastic-plastic characteristics of the compacted soil were considered, and the soil under the wheel's vertical stress and contact force characteristics were analyzed. The main conclusions are as follows:

1) When the compacted material is cohesive soil, along the lateral direction of the vibrating wheel, the vertical stress contour of the surface layer of the soil is "saddle-shaped". With the diffusion of stress, the shape of the stress contour would switch to a parabolic type, which is caused by the stiffness of the vibrating wheel being greater than that of the soil. The soil's vertical stress attenuated at the wheel's edge. Inflection points can be found by observing the vertical stress distribution curves at different depths. And, the deeper the soil layer, the closer is the location of the inflection point to the center of the drum. To ensure the overall compaction quality, the overlap width of the vibrating wheel in actual project can be determined according to the position of the inflection point of the stress distribution curve of the bottom soil, which was about 1/5 of the half-width of the wheel.

2) Along the compaction path's longitudinal direction, the soil's vertical stress distribution shows a typical stress bubble shape. The center point of the stress bubble is located in front of the mass center of the vibrating wheel, which is due to the uplift of the soil ahead of the wheel and the inclination of the contact surface between the wheel and the soil. This phenomenon makes the contact pressure between the wheel and soil produce a horizontal component, thus intensifying the shear stress, which is unfavorable to the vibratory compaction. Increasing the speed of the vibrating wheel could overcome this uplift effect, but the balance with the compaction effect of a single pass needs to be considered.

3) The total grounding width of the vibrating wheel can be estimated from the time course curves of the contact stress at the element nodes. However, the contact forces are not uniformly distributed within this width. The effective contact width of the vibrating wheel can be defined based on the loading time through significant compressive deformation of the element nodes, which was about a third of the total contact width. Within this region, the soil deformation is dominated by vertical compressive deformation. In contrast, the soil deformation in the other region ranges is dominated by plastic uplift caused by shear stress. Using the effective contact width when describing the wheel-soil contact load is more reasonable.

Acknowledgments

This paper was supported by the Department of Transportation of Shandong Province (No. 2018B51).

Conflict of interest

The authors declare that there is no conflict of interest.