Discrete element modeling and mesoscopic parameter calibration of bottom voiding
The calculation model for the bottom voiding of the concrete pavement is illustrated in Fig. 2. The dimensions of the computational model were referenced from the numerical model sizes presented in literature. This selection was made to ensure computational efficiency while simultaneously avoiding model size effects during the simulation. This choice guarantees the accuracy and generalizability of the results. The concrete panel has a length of 5 m, a width of 4 m, and a thickness of 0.26 m. The base layer is 0.5 m thick, and the foundation is expanded to dimensions of 7 × 6 × 6 m. The voiding area is positioned at the bottom of the concrete panel, with dimensions of 1 × 1 × 0.01 m. The primary calculation parameters for each structural layer are provided in Table 1.
The vehicle load is represented by a single wheel pressure with a 20 × 20 cm contact area, and the wheel spacing is 30 cm. Through calculations and comparison of different loading positions, it was observed that when the pavement slab edge is grouted, the wheel load symmetrically applied at the bottom plate on one side of the slab is most detrimental to the stress distribution of the pavement. The loading arrangement is shown in Fig. 3.
In coupled simulations, when the discrete element method (DEM) is used to model macroscopic engineering, the number of particles required is large and the calculations are time-consuming. The continuous-discrete coupling method can help alleviate this problem by reserving a void region in the continuous element core section, where DEM modeling is applied. Communication between continuous and discrete elements takes place through an interaction interface, thus ensuring both computational efficiency and accuracy.
The principle of FEM-DEM coupling is illustrated in Fig. 4. This approach uses walls/shells as interaction interfaces. Velocity data from the continuous element nodes are passed to the wall/shell unit nodes, where they are converted into unbalanced forces, which are then distributed to the DEM particles via shape functions. Likewise, the unbalanced forces in the DEM are converted into velocity and sent back to the continuous element nodes through the wall/shell unit nodes, continuing the cycle until dynamic equilibrium is reached. Once coupling is activated, the calculation must run under large deformation conditions, and to ensure the balance of forces and moments at the coupling interface, full computation mode should be activated.
The contact bonding constitutive model in the particle flow method can effectively simulate the generation and development of cracks in the specimen. In the program, the cohesion between particles is determined by setting the tangential and normal bonding strengths and the friction coefficient. If the stress in any direction of the particle model exceeds its corresponding bonding strength during loading, the bonding between particles will fail, thus generating cracks. Cracks can only form during the loading process in the particle flow program if a cohesive contact model is applied between the particle bodies. Thus, the microscopic parameters set when assigning the contact model will significantly affect the number and position of the cracks that form in the model. The cohesion between particles can be considered as a cylindrical surface in the normal direction within the model plane, as shown in Fig. 5. The crack representation in this model is as follows:
Assume that the two particles generating the crack are A and B, the thickness of the crack cylindrical surface can be expressed as:
In the formula: d represents the distance between the two particles; represents the radius of particle A; represents the radius of particle B;
The central position of the cylindrical surface can be expressed as:
In the formula: represents the normal direction from to .
The radius of the cylindrical surface can be represented as:
As described above, in the particle flow program, cracks in the model can be represented by parameters such as thickness, radius, normal direction, and center point location. The thickness can be represented by the gap between the two particles, the radius can be represented by the length of the mid-plane of the connecting cylindrical surface between the particles, the crack normal direction aligns with the direction of the line connecting the centers of the two particles, and the center point is the intersection of the centerline between the particles and the positions of the two particle centers. Therefore, in numerical simulations, by recording the type of damage when the contact bond breaks, it is possible to overcome the limitations of laboratory and field in-situ tests, which are often unable to observe and differentiate between tensile damage and shear damage.
In the discrete element method, the macroscopic mechanical properties of the material are governed by the micro-parameters of the particles and the contact constitutive model. In this study, uniaxial compression tests are conducted to obtain the physical and mechanical property indicators of the discrete element model. The uniaxial compression test was conducted using a servo-hydraulic testing machine with a maximum loading capacity of 3000 kN. The machine is equipped with a high-precision load cell and a linear variable differential transformer (LVDT) for measuring axial deformation, with a resolution of 0.001 mm. Standard concrete test blocks were selected to conduct uniaxial compressive mechanical tests. The standard cubic test blocks with the size of 100 × 100 × 100 mm were used as the key specimens for the uniaxial compressive tests. During the test, the axial load and corresponding axial deformation were continuously recorded until the specimen failed. For each of the three groups, the compressive strength, elastic modulus, and peak strain were calculated based on the recorded data. To ensure the reliability and reproducibility of the results, the final test results (including compressive strength, elastic modulus, and peak strain) were determined as the average value of the three groups of test results. This averaging method helps to eliminate random errors caused by minor variations in specimen preparation and testing processes, providing more representative and robust data for subsequent analysis. The test model of the standard concrete test block and the discrete element calculation model are shown in Fig. 6, which provides a simulation basis for the subsequent in-depth numerical analysis and study of mechanical properties. The particle size and micro-parameters of the uniaxial specimen are consistent with the foundation parameters. Rigid wall boundaries are used as the upper and lower loading plates for the specimen, with a loading rate of 0.01 mm/s, and the normal and tangential contact stiffness are identical to the specimen's stiffness. After multiple trial calculations, the mesoscopic parameters are as shown in Table 2, and the stress-strain curves of different models are as shown in Fig. 7.
Comparison and validation between DEM and laboratory tests in the research on the concrete surface layer, uniaxial compression tests were conducted both through DEM simulations and laboratory experiments, yielding the corresponding stress-strain curves as shown in the Fig. 7a. A comparative analysis of the curve profiles reveals a high degree of consistency between the DEM simulation curve and the laboratory test curve. During the stress development phase, both curves exhibit a similar pattern: stress increases progressively with strain, reaches a peak value, and then declines. Focusing specifically on the peak strength, both the laboratory tests and DEM simulations yield a value of approximately 30 MPa, demonstrating excellent numerical agreement. This consistency strongly validates the rationality and effectiveness of the mesoscopic parameters employed in the DEM simulations. In DEM modeling, mesoscopic parameters are fundamental for constructing the model and accurately representing the material's mechanical behavior. The close match between the simulation results and the actual laboratory test data, particularly in such a critical mechanical peak strength and the overall curve trend, indicates that these parameters can precisely capture the mechanical response characteristics of the concrete surface layer under uniaxial compression. Consequently, they provide a reliable data foundation and parameter basis for subsequent DEM-based simulations under more complex conditions, enabling in-depth investigations into the evolution of the concrete surface layer's mechanical properties under various conditions.
For the base layer, stress - strain curves from DEM simulation and lab uniaxial compression test are shown in the Fig. 7b. As strain rises from 0, both curves DEM: solid line; dashed line show stress increasing, reflecting internal force mobilization under compression. Both DEM and lab tests give a peak strength of 5 MPa, with peak strain at 0.32%. This high consistency validates the rationality of mesoscopic parameters in DEM.
The FEM-DEM coupled model for the bottom void of the concrete pavement is shown in Fig. 8. The subgrade consists of particles with a radius of 3 cm, and the concrete slab consists of particles with a radius of 2.5 cm. The model is composed of a large number of particles, totaling 125,900, with 59,587 particles in the concrete surface layer and 66,314 particles in the subgrade. This model is developed using the DEM and requires validation before use to ensure its reliability for further research. Choosing the appropriate particle radius ensures the accuracy of the simulation results. A larger particle radius reduces computation time but may ignore some microscopic mechanical behaviors, while a smaller radius provides more accurate simulation of micro-mechanical behaviors at the cost of increased computational cost. The PFC is used to simulate the behavior of granular materials, reflecting the mechanical properties of the material through particle interactions. The subgrade is modeled with a finite element mesh, which helps improve computational efficiency. The parameters are shown in Table 3. The FEM-DEM coupling method combines the benefits of the FLAC and the PFC. The FEM is advantageous for continuous media and large deformation problems, while the DEM is excellent for simulating granular materials and crack propagation. This coupled method enables a more comprehensive simulation of complex geotechnical engineering problems.
For the discrete element model, the boundaries are defined as free boundaries. This means that when the model is subjected to external loads, the particles within the model are allowed to move freely in the spatial domain without any constraints, which accurately simulates the actual mechanical behavior of the material under unconstrained conditions. In terms of the mesh boundaries, specific constraints are applied to ensure the stability of the computational domain during the simulation. The surrounding boundaries of the mesh are constrained in the x and y directions, preventing any displacement in these two directions, while the z direction remains unconstrained. The bottom boundary of the mesh is fully constrained in all x, y, and z directions to fix it in space, providing a stable foundation for the entire simulation system.
When constructing the coupling model, the core research area is set as a discrete domain, and DEM is used to depict mesoscopic cracking and block movement; the area with gentle deformation and uniform mechanical response is set as a continuous domain, and FEM is used to efficiently calculate stress and strain, achieving "refinement of the core and high efficiency of the periphery". To ensure the reliability of the finite element (FEM) results for the foundation layer (a key part of the FEM-DEM coupling model) and balance computational accuracy and efficiency, a mesh sensitivity analysis was conducted as a critical preliminary step before formal simulation.
The mesh sensitivity analysis in this study focuses on the foundation layer (FEM domain), with the settlement of the most central mesh after the upper particles exert pressure on the mesh under self-weight as the reference index. We have designed four groups of schemes with different mesh sizes to compare the computational accuracy (difference between relative error and theoretical value) and computational efficiency (computational time). The specific parameters of each scheme are detailed in Fig. 9; Table 4.
Based on the above comprehensive evaluation, with the results of Group 4 as the reference standard, Group 3 demonstrates a relatively small relative error (0.66%) while significantly reducing computational resource consumption (with 66.7% fewer finite elements and 57.1% less computational time compared to Group 4). This balance between accuracy and efficiency is particularly critical for large-scale FEM-DEM coupling simulations, as it avoids excessive computational burdens caused by over-pursuit of precision while ensuring that the simulation results meet engineering analysis requirements. Therefore, considering both the reliability of research data and the practicality of computational implementation, this study ultimately selects the mesh size of Group 3 as the optimal unit size for the foundation layer.
The establishment of the void model employs an efficient approach. Initially, a wall is created in the void region. The role of this wall is to maintain the model's stability at the initial stages of model setup. Once the model reaches equilibrium and appropriate microscopic parameters are assigned, the wall is deleted to create the void region. The benefit of this approach is that after the model achieves equilibrium, the wall is removed to create the void region. At this stage, the model is stable, allowing for further calculations to be carried out based on this stability, thus reducing the number of iterations and computational load.
This paper uses a rigid wall as the loading plate for the subgrade model. Based on the research results of Zhang et al., after several trial calculations, a loading rate of 0.005 mm/s was determined. This speed guarantees that the model undergoes quasi-static loading during the loading process, thus optimizing computational efficiency.
Figure 10 shows displacement contour plot and force chain distribution of the DEM model after self-weight equilibrium. From the figure, it can be observed that the settlement displacement shows a roughly horizontal distribution. This implies that in this soil model, the settlement amounts at various horizontal layers are fairly uniform. This phenomenon is consistent with the property of uniform vertical compression of the soil under gravity. In actual engineering, this kind of horizontally distributed settlement displacement typically occurs in homogeneous soil layers. Under the influence of gravity, the soil layers compress uniformly downward, resulting in similar settlement amounts across each layer. The force chain shows a dense distribution in the lower section and a more sparse distribution in the upper section. This distribution pattern aligns with the actual soil's initial stress distribution under the influence of gravity. Due to gravity, the lower layers of soil must bear the weight of the upper soil. As a result, the force chain between particles in the lower soil is more compact. The upper soil, however, experiences less pressure, resulting in a relatively sparse force chain.
The force chain shows a dense distribution in the lower section and a sparse distribution in the upper section. This distribution pattern aligns with the initial stress distribution in the soil under the influence of gravity. In the soil, due to gravity, the lower layers must support the weight of the upper layers. As a result, the force chain between particles in the lower soil is more compact. In contrast, the upper soil experiences less pressure, and the force chain is relatively sparse.