Non-linear micro-mechanical behavior of heterogeneous cement-stabilized macadam

Non-linear micro-mechanical behavior of heterogeneous cement-stabilized macadam This study aims to examine the internal meso-mechanical behaviour of cement-stabilised macadams subjected to loading damage following water immersion. By constructing a numerical analysis model, the mechanism of water-immersion weakening of the cement-stabilised macadam material and the internal micro-mechanical behaviour of the specimen were analysed. The numerical simulation results were compared with laboratory test data. The peak stress error was 2.414 %, whereas the maximum error in the strain corresponding to the peak stress of the specimens with different immersion durations was 3.05 %. By comparing the direction distributions of the normal and tangential contact forces in each stage, it was found that the tangential forces at approximately 0º, 90º, 180º, and 270º were much smaller than those at other angles. In the peak stress stage, the normal contact angles between 45–135º and 225–285º were higher than those at other angles. For model M (mean factor) less than 20, the peak stress was significantly affected and increased as a function of M. When M was greater than 20, the peak stress fluctuated within a small range as a function of M. This study provides a new research method for the study of the structural mechanical effect of pavement structures based on cement-stabilised macadam materials damaged following water immersion.


Introduction
Cement-stabilised macadam (CSM) materials have excellent mechanical properties, such as high strength, high stiffness, and good stability. Moreover, they have the ability to substitute aggregates, which can effectively reduce project costs [1,2], thus resulting in their widespread usage as pavement base materials. As a non-homogeneous, discontinuous, artificial, quasi-brittle material, it is an important research direction in solid mechanics. As the bearing layer of the pavement structure, the CSM base layer can transfer the surface layer load evenly to the foundation. Increasing the strength of the base layer can effectively improve the stability of the surface layer [3]. Thus, an in-depth study of the structural mechanical properties of the subgrade is crucial to the service duration and quality of pavements. Therefore, it is very meaningful to study the mechanical properties of CSM materials [4]. Most studies on CSM have focussed on the analysis of the cracking pattern of the subgrade, often based on laboratory tests, to analyse the macro-mechanical properties of the structure. The development of the intrinsic structural model of the pavement has gradually progressed from early linear elastic mechanics to the current stage of elastic-plastic and damage fracture mechanical models [5]. In the numerical simulations of road structures, the main methods include finite, boundary, and discrete element methods. Currently, the first two methods are at a more mature stage, and the model calculations are efficient, thus yielding stable results. The development of macro-mechanics is based on various tests aimed at optimising the mechanical properties of CSM materials. In his study on the tensile and shear damage behaviours of cement-treated substrates subjected to perimeter pressure, Namikawa applied tensile forces to specimens using a drainage triaxial tension test [6]. Zhang performed three-point bending tests on CSM with different cement ratios and curing times and analysed the fracture characteristics of the materials by determining the fracture toughness, fracture energy, and deflection of the CSM [7]. Most numerical analysis models are treated as single homogeneous models to improve computational efficiency. Namikawa developed three finite element models of cementtreated substrates in direct tension, split tension, and bending tests to explain the differences in tensile strength [8]. Karadag created a two-dimensional symmetric finite element model for a flexible pavement that could effectively simulate the interactions between layers [9]. Fracture energy significantly decreased the brittleness index of large particles, as proposed by Tasdemir [10]. Mihashi proposed an uncertainty in the relationship between material strength and aggregate size [11], whereas Chen reported that the strength decreased at increasing particle sizes [12]. Owing to the inability of macro-level tests to determine the mechanism of mechanical interaction between internal mesoscopic particles and to control precisely the internal microcrack distribution, porosity, and other mesoscopic parameters, the macroscopic mechanical behaviour exhibited by the structure cannot be proposed as a rule. Since Roelfstra et al. first analysed the non-homogeneous properties of a numerical model of concrete materials from a mesoscopic perspective, concrete has been studied by dividing it into aggregates, cement mortar, and intersections. Thus, to analyse the mechanical properties of concrete, meso-mechanical studies must be supported by numerical analyses [13]. Regression relationships must be established using macroscopic mechanical behaviours, such as the elastic modulus of materials, number of layers, and number of structures in the subgrade structure, and the linear relationship between macroscopic parameters can be used to change the particle grading structure to improve economic benefits. A detailed analysis of the interaction between the macroscopic parameters of pavement foundation soils, such as the inclusion of mesoscopic mechanical parameters and the mathematical model of parameter influence, will yield a more complete and accurate understanding. In addition, when the grade distribution is adjusted, changes in the internal structural mechanical behaviour provide a more in-depth understanding [14]. The discrete element model creation process considers the effects of porosity, pore distribution, particle size distribution, and aggregate shape on the mechanical properties of the model, thus simulating the development of cracks in asphalt concrete in disk-type compaction and threepoint bending tests. Moreover, the effect of particle size on the fracture process of materials has been observed [15][16][17]. A threedimensional discrete element asphalt mixture model was created to simulate the indirect tensile tests of nonhomogeneous brittle materials, and all achieved good agreement. Consequently, the effects of the internal mesoscopic displacement field distribution and variations in the structure, aggregate gradation, loading speed friction coefficient, and tangential shear stiffness ratio on structural strength were investigated [18,19]. Although the present discrete element model can incorporate the influences of mesoscopic particles between material structures, the CSM, as an artificial quasi-brittle material, must present a random distribution of mesoscopic physical parameters within the material structure. It is evident that the mesoscopic parameters are considered as constant values when performing numerical simulations, which is different from the actual situation. The influence of various factors on the contact diameter between different particle contacts was ignored. Moreover, in the loading phase of the model, the simulated results yield a linear relationship between stress and strain, which does not indicate a gradual increase in the rate of stress increase in the early loading phase owing to pore compaction. However, to enhance pavement service life, only probabilistic predictions of pavement water ingress parameters can be made at the design stage [20]. In addition, there is a serious deficiency in the current analysis of the damage mechanism from a mesoscopic perspective following water immersion of CSM materials. Therefore, this study applied discrete element theory to create a numerical analysis model to simulate uniaxial compression experiments on CSM damaged by water immersion. Subsequently, the radius multiplier Weibull distribution function was added to render a random distribution of the particle contact radius in the model and realise the non-homogeneous characteristics of the material. To realise the non-linear dynamic process of GRAĐEVINAR 75 (2023) 5, 429-438 Non-linear micro-mechanical behavior of heterogeneous cement-stabilized macadam the stress-strain curve during loading, mesoscopic parameter control functions were added to the numerical analysis model. The effect of the discrete element model on the damaged submerged specimens was used to analyse the trend of peak stress and contact area, the distribution of normal and tangential contact force directions, the effect of the homogeneity factor on peak stress, and other mesoscopic mechanical behaviour studies.

Uniaxial compression test of water-damaged
water CSM in laboratory settings

Specimen preparation standards and water immersion damage
According to the standard JTG E51-2009, the standard preparation size for gravel gradation specimens is 150 mm x 150 mm ( Table 1). The specimens were immersed in water at the beginning of the preparation stage, and to discharge internal air bubbles, the interference was increased by 37.5 cm every 2 h. The length of immersion was recorded after 8 h. Each immersion length group was controlled at 3, 7, 15, 20, 25, and 30 days (d).

Uniaxial compression tests
To test accurately the strain change process of the specimen during loading, a strain gauge was pasted at the centre of the specimen (25 mm from the upper and lower bottom surfaces) in the vertical direction and on the flat and smooth surfaces, as shown in Figure 1. To reduce the resistance change caused by the loading process at the connecting contact point between the strain gauge and conductor, and hence produce measurement errors, the strain gauge and conductor were welded, as shown in Figure 2. Subsequently, it was wrapped with insulating tape for insulation.

Creation of aggregate model
According to the grading of the laboratory materials, the aggregate shown in Figure 3 was generated in the model, and the particle diameter was generated in each small range specified in Table  1

Constitutive contact model
A mechanical contact relationship was applied between each aggregate interface for the created model aggregates. The discrete element model selects the parallel contact bonding model [21], as shown in Figure 4, where denotes the bond strength, is the friction angle, C is the normal tensile strength, is the parallel shear stiffness, and is the parallel normal stiffness. Parallel contacts can be classified into bonded and un-bonded states. When the tension of a bond exceeds the bond energy limit, it breaks and degenerates into a linear contact model. Furthermore, when relative displacement occurs between the particles, the contact force of the bonding model is affected by the linear contact force, damping force, and force produced by the parallel bonding contact. The resulting contact moment M c is determined only by using the parallel bond moment M c .

Heterogeneous distribution model of mesoscopic parameters
In the parallel-contact bonding model, the contact radius directly determines the basic parameters of the interparticle tensile and flexural torsional strengths. In traditional discreteelement numerical simulation calculations, the radius multiplier is typically not set (it is set to one by default), as shown in Figure  5. If the contact connects two particles (with radii r (1) and r (2) ), the contact bond radius is considered equal to the small particle radius r (2) . If the contact connects the particle to the wall, the contact radius is considered to be equal to the particle radius r (3) . The parallel bond radius multiplier is the contact radius increase or decrease factor, at which point the contact radius can be expressed using Eq. (3). The model water immersion damage process was based on the assumption of constant particle morphology, while the radius multiplier was discounted. By setting the radius multiplier distribution state, the nonhomogeneous characteristics of the bonding material can be achieved, and the homogeneity of the model can be controlled.
(3) However, as an artificial quasi-brittle material, the distribution of the radius multiplier of the CSM should exhibit a certain randomness, which will be different from the actual situation. Therefore, the Weibull distribution model was introduced in this study. As a continuous distribution function, the probability density function of the Weibull distribution is expressed as, where x represents a random variable, u > 0 represents a scale parameter, and k > 0 represents a shape parameter.

Weibull distribution statistics yield:
Mean value: Variance: The symbol Γ represents the gamma function, Therefore, in the probability distribution of the discrete element model, the created radius multiplier obeys the Weibull distribution, which can be expressed as GRAĐEVINAR 75 (2023) 5, 429-438 Non-linear micro-mechanical behavior of heterogeneous cement-stabilized macadam where λ represents the distribution variable of random distribution of radius multiplier, and m represents the homogenisation factor of model material distribution, which changes with the homogenisation factor. The probability density distribution diagram is shown in Figure 6.

Non-linear control of model loading
When the traditional discrete element model is used to simulate the uniaxial compression process, the stressstrain curve in the additional phase is typically linear, and the stress is proportional to the strain when the peak stress is not reached. Typical stress-strain curves in realistic, uniaxial compression tests on CSM materials are shown in Figure  7. The observations indicate that the strain increases with increasing stress in the pre-loading period and exhibits a non-linear phase of gradually decreasing strain growth rate when the stress increases uniformly, that is, the OA section in the figure. The structure is compressed by the pores, aggregates are gradually occluded and stressed, structure is continuously compacted, and the model stiffness gradually increases with increasing stress. When it reaches a slope equal to the material's elastic modulus, the non-linear phase ends in the early stage. At this point, the stress continues to increase at a uniform loading rate, and the strain begins to increase at a fixed rate in the linear loading phase (section AB in Figure 7). To simulate the non-linear loading process using a numerical analysis model, parameter control functions were added to the model. Consequently, the internal mesoscopic parameters of the model changed at increasing stress values during the pre-loading process until the model entered a linear loading phase.

Model meso-parameter calibration
Currently, the direct measurements of meso-parameters inside a material using laboratory equipment are challenging because the test results obtained through macro-tests yield large calculation errors, and the coupling effect between

Macro-mechanical behaviour data extraction
The stress-strain data obtained based on uniaxial compression tests on the submerged specimen in the laboratory and the data obtained based on numerical simulation tests on the discrete element model are shown in Figure 8. The solid and dotted lines represent test and model simulation data, respectively. The curves in the table indicate that the increases in stress with strain in the model and laboratory test curves were the same. In Figure 9, the blue curve represents the peak stress of the uniaxial compression test for different immersion periods in the laboratory, and the red curve represents the peak stress obtained based on model simulations. Appropriate comparisons revealed that the maximum difference between the two occurred in the non-immersed specimen, with an error of 2.414 %. After fitting, the correlation coefficient was R 2 = 0.997, the fitting effect was good, and the fitting function is, y = 27,79538 -0,97379 · (10) = 27,79538 + 14,80824[exp(-0,06576 · x)-1]

Figure 9. Strain at peak stress
The corresponding strains between the test and numerical simulations when the stress reached the peak stage are shown in Figure 10. The black and red curves indicate the experimental and numerical simulation results, respectively. The maximum error between the two occurred when immersed in water for 7 d, reaching 3.05 %. The peak test strain reached a value of 2.109 x 10 -3 after 30 d of immersion, and the model simulation error was only 0.20 %. The peak strain tended to decrease slightly in the range of 0 to 7 days of immersion, whereas it gradually increased as the immersion time increased. Non-linear micro-mechanical behavior of heterogeneous cement-stabilized macadam

Micro-mechanical behaviour data extraction
During the process of water-immersion weakening of the CSM, the contact area was found to decrease as a function of immersion time by detecting the contact cross-sectional area between the particles in the model. As shown in Figure 11, the following functional relationship was obtained by fitting the contact area to the immersion time, y = 5,33273 -0,9815 · (11) = 5,40428 + 2,20860[exp(-0,04444 · x)-1] The correlation coefficient of the fitted function was R 2 = 0.98963; thus, the fit was good.

Distribution law of tangential and normal contact forces
To study the distribution of the tangential contact force in the internal micro-mechanical behaviour at different stages during the loading process of the CSM materials, the distributions of the tangential and normal contact forces of the discrete element model at the initial stage of loading and the stage when the stress reached the peak were extracted, as shown in Figure 13. The green\ and pink areas represent the normal and tangential contact-force distributions, respectively. According to the contact force distribution diagram at the initial stage of loading, the normal contact force distribution in each direction was relatively uniform and the relative values were small (all below 3 N). In the tangential contact-force distribution, the angle ranges were -15(345)-15°, 75-105°, 165-195° and 225-285°, and the value range was 0.275-0.98647 N, which is significantly smaller than that for the other angles of 1.17-2.29201 N; that is, the horizontal and vertical directions were significantly smaller than other directions. At the peak stage of loading, the distribution law of the tangential contact force was the same, and the range for thehorizontal and vertical directions (11583.24579-17018.27479 N) was significantly smaller than that at other angles (18489.45798-34101.61295 N); the normal contact force range for the angular ranges of 45-135° and 225-285° was 67470.7416-120686.0179 N, which was significantly higher than that for the angular ranges of -45(315)-45° and 135-225°, which was 14973.26942-26597.49541 N. The contact force at the peak stage was significantly greater than that at the initial loading stage.

Comparison of macro-and meso-mechanical behaviours
Based on comparisons of the macro-and meso-mechanical behaviour analysis data, the structure's formula used in the fitting of the meso-contact area and immersion weakening was obtained. The formula is identical to the fitting curve of the stress peak and immersion time, as shown in macro-mechanics. The relationship between the peak stress and contact area is presented in Figure 12. Function fitting revealed a positive linear correlation between the magnitude of the stress peak and the contact area. When the contact area increased, the stress peak increased.

Influence of homogenisation factor
The variation in the stress intensity of the discrete element model with the homogenisation factor in the dry state is shown in Figure 14. A comparison of the curves reveals that the homogenisation factor changes; however, the changing trend of the stress-strain loading curve remains the same. The homogenisation factor had a minimal effect on the strength of the model during the early loading stage. In the later stages of loading, the strength gradually increased with an increase in the homogenisation factor. The change in the stress peak with respect to the homogenisation factor is shown in Figure  15. When m reached 20, the size of the stress peak tended to stabilise.

Conclusions
To determine the water-immersion damage mechanism of the CSM material and the influence of water-immersion damage on the mechanical behaviour of the material subjected to loading, a discrete element model of the CSM at different immersion times was established according to the laboratory material gradation. The non-linear change in macro-stress during loading was controlled by adding a parameter control function, and the Weibull distribution function of the meso-parameter (radius multiplier) was added to realise the heterogeneity of the model. Based on the experimental verification in the laboratory and the analysis of the numerical simulation results, the following conclusions were drawn: Non-linear micro-mechanical behavior of heterogeneous cement-stabilized macadam -By comparing the laboratory uniaxial compression test data and the analysed results of the numerical simulation test, the maximum error of the stress peak was obtained as 2.414 %. Furthermore, the maximum strain error corresponding to the stress peak was 3.05 %, satisfying the error requirements. It was proven that the discrete element model is feasible for simulating the uniaxial compression tests of CSM materials. -The relationships between the peak stress and immersion time and between the contact area and immersion time were extracted and fitted. The two fitting functions had the same structural formulae. Subsequently, the two fitting functions were fitted to the peak stress and contact area, which proved that the proportional linear relationship between them increased as a function of the contact area. Simultaneously, the reduction in the contact area was determined to be the main factor in the water immersion damage and stressbearing capacity reduction of the CSM materials.
-Based on the extraction of the normal and tangential contact force direction distributions in the initial loading stage and the stress peak stage in the numerical simulation process, it was found that the transverse and vertical (-15-15°, 75-105°, 165-195°, 225-285°) shear stresses in each stage of the tangential contact force were lower than those at other angles. The normal contact force between 45-135° and 225-285°at the peak stage was significantly higher than that at the other angles. The contact force at the peak stage was significantly greater than that at the initial loading stage.
-The peak stress intensity of the numerical analysis model increased with the homogenisation factor when the homogenisation factor was less than 20. However, when the homogenisation factor was g.