Microscopic mechanism of subgrade vibration compaction based on discrete element method

1.   Introduction

Compaction is an important means to ensure the performance and longevity of the subgrade and pavement structures. To enhance quality control of compaction, many scholars have conducted research on the mechanisms of vibratory compaction. Traditional methods mostly employ finite element analysis to study the macroscopic mechanical response of the compaction system. However, the essence of compaction lies in the rearrangement and redistribution of material particles under external forces, which is often difficult to capture using macroscopic mechanical methods, thus limiting the understanding of the compaction mechanism. Consequently, some research studies have started to introduce micro-mechanical methods based on discrete elements, aiming to analyze the behavior of internal microstructures and reveal the compaction mechanism of particle systems at a deeper level.

The discrete element method can reflect the compositional characteristics of non-uniform materials and is particularly adept at simulating the micro-behavior of heterogeneous materials during large deformations. It has already found some application in the field of road engineering. Salazar et al. ^[1] employed a rolling friction model to establish a three-dimensional discrete element model, simulating the direct shear test of coarse sand. They improved computational efficiency using the radius expansion method and demonstrated that the stress paths of sand particles can be appropriately reproduced. Ma et al. ^[2] used PFC2D/3D software to simulate the effects and depth of compaction reinforcement in crushed stone soil foundations. By analyzing the dynamic hysteresis effect, contact model parameters, and the influence of the shape of the compaction hammer on the computation results, they concluded that pore changes induced by particle rotation must be considered when analyzing such foundations. Michael et al. ^[3] optimized the algorithm for discrete element simulation of crushed stone materials, generating samples with realistic grain shapes, sizes, and relative densities. They successfully validated the program through indoor heap soil tests, compaction tests, consolidation tests, triaxial tests, and multistage shear tests. Ma et al. ^[4,5] analyzed the structural characteristics of reinforced retaining walls using discrete elements, and by analyzing the particle displacement, stress distribution, force chain distribution and potential damage surface of stepped reinforced soil retaining wall (RSRW), it was concluded that the potential damage surface of stepped RSRW can be determined by observing the changes in the force chain structure, and the numerical simulation results of the potential damage surface match well with the theoretical calculation results. The aforementioned research data suggests that employing the discrete element method to simulate the mechanical behavior of soils is feasible.

Vibratory compaction is a complex mechanical process, and there is currently limited experimental evidence for its theoretical validation ^[6,7,8] due to operational issues. Therefore, using discrete element software to analyze the particle motion response and changes in the mechanical response between particles during the vibratory compaction process can effectively optimize existing theories. Wang ^[9] verified the variations in porosity and mechanical response during the vibratory compaction process by establishing a "steel wheel-soil" discrete element model. They also analyzed the variations in microscopic parameters between particles from a dynamic perspective. Wang ^[10] used porosity as a quantitative evaluation index for compaction quality and applied the three-dimensional discrete element analysis method to study the variations in compaction quality with particle shape parameters and vibration force parameters. Wang ^[11] employed a three-dimensional discrete element method for modeling and concluded that the dynamics of the system during vibratory compaction primarily depend on factors such as vibration frequency and amplitude, particle friction, and initial conditions.

At present, the research on the compaction mechanism of subgrade soil is only limited to the change rule between the compaction parameters of subgrade soil and the various types of responses of subgrade soil. The traditional macro-mechanical methods and finite element simulation methods are unable to better interpret the compaction mechanism of subgrade soil from the material point of view. In addition, both methods analyze the subgrade soil material as a homogeneous material, thus ignoring the material properties of the subgrade soil. In contrast, the use of discrete elements to study the compaction mechanism of the subgrade soil from a fine viewpoint not only reflects the inhomogeneous material properties of the subgrade soil, but also analyzes the response mechanism of each index of the subgrade soil in the compaction process from a microscopic viewpoint.

2.   Discrete element model of vibratory compaction system

2.1.   Generation of particles

A vibratory compaction model for sandy soil was established using the PFC2D5.0 software. For sand materials, the bonding between particles can be ignored and the load-bearing capacity is mainly provided by the interlocking and sliding friction between particles. However, the basic unit of traditional discrete element methods is a sphere model, which lacks angularity and cannot reflect the interlocking characteristics between particles. Therefore, this section attempts to explore methods for generating irregularly shaped particle models.

The mainstream method for generating irregularly shaped particles is through the use of image processing techniques. By setting a grayscale threshold, the image of the mixture specimen is processed into a binary image, distinguishing the aggregate from other materials and subsequently generating a two-dimensional digital specimen in PFC software ^[12], as shown in Figure 1. Figure 1(a) shows the image of the asphalt mixture after gray scale binarization, where white represents aggregate and black represents asphalt. In Figure 1(b), to ensure that the sandy material has more distinctive features in the discrete elements, the image is processed so that black is aggregate and white is asphalt. It should be noted that in PFC-based simulations, using image-reconstructed particles is a more accurate method as it can address the limitations of the "ball generate" command in PFC software in terms of computational capabilities and porosity. However, currently, the models generated using image processing techniques are mostly of small size, primarily for indoor mechanical testing requirements. However, vibratory compaction models typically have larger dimensions, and if this method is used, the computational requirements may exceed the capabilities of the software itself. Therefore, such methods are not applicable to vibratory compaction problems.

Figure 1.  The application of image processing method in Particle Flow Code (PFC): (a) grayscale image of the specimen; (b) Discrete particle models.

DownLoad: Full-Size Img PowerPoint

This study used the disk packing method. Using the "Geometry" command, geometric shapes can be created or imported as the external contour of the particles. With the help of the "clump" module, the program can automatically fill disk particles inside the contour and fit them to the external contour. The filling principle is to find a disk among N possible disks that maximizes the additional coverage area between particles. The fitting accuracy of the contour boundary can be controlled by setting the "bubblepack distance" keyword between 0 and 180, which represents the tangent angle at the intersection of adjacent disks. A larger value results in a smoother boundary after filling. Additionally, the "ratio" keyword can control the minimum and maximum radius ratio of the filling disks.

We used AutoCAD to draw the two-dimensional graphics of aggregates, and then imported them to generate DXF files to create a geometry template. The outline of the templates is shown in Figure 2. Based on these templates, the particles were filled. After multiple tests, the parameter for bubblepack distance was determined to be 150, and the minimum to maximum radius ratio was 0.5. As a result, particles filled by six discs were obtained (Figure 2(c), (d)). Table 1 provides the area parameters for the two types of particles.

Figure 2.  Contour of the particle template: (a) stone 1; (b) stone 2; Particle shapes created by the disc filling method: (c) stone 1; (d) stone 2.

DownLoad: Full-Size Img PowerPoint

Table 1.  Particle area parameters.

Particle size/mm	Particle area/mm2
	Stone1	Stone2
1	1.4282	1.6685

---

 | Show Table

DownLoad: CSV

After determining the particle shape, the corresponding particle group model can be generated based on the gradation. In this study, the soil from an in-situ compaction test was subjected to sieve analysis, and the gradation obtained is shown in Table 2. In this table, particles larger than 0.5 mm account for less than 50% of the total mass, indicating that the soil is classified as medium sand. By calculation, the coefficient of curvature (${C}_{c} = 1.914$) is within the empirical required range of 1–3, and the coefficient of uniformity (${C}_{u} = 36$) is greater than the empirical required value of 5. Therefore, this type of soil belongs to a well-graded soil.

Table 2.  Gradation of sandy soil.

Particle size/mm	10.0	5.0	2.0	1.00	0.5	0.25	0.075
Passing rate/%	100	72.5	42.7	36.5	21.2	18.4	7.8

---

 | Show Table

DownLoad: CSV

The unit sieve percentage is a percentage obtained by dividing the mass mi retained in each sieve by the total mass M. When the density of the selected material is constant, the unit sieve percentage is actually the percentage of the volume Vi in a certain particle size interval to the total volume V. Therefore, the volume Vi of particles in a given particle size interval can be calculated under the condition that the volume V of the desired generated aggregate and the unit sieve percentage p in each particle size interval are known. Since the PFC software often uses the number of generated particles when generating particles, it is assumed that the particle radius distribution in each particle size interval follows a normal distribution. The number of particles in the interval is obtained by dividing the volume Vi of the particles in the interval by the mean particle volume $\stackrel{-}{{V}_{i}}$. The specific calculation method is shown in Eqs (1)–(4).

As can be seen from the above calculation formulas, when the gradation is known, the number of particles can be obtained as long as the average volume in each particle size interval is obtained. Under 2D conditions, the above volumes are represented by areas. Based on Table 1, the shape of the particles was proportionally scaled to obtain the particle area of different particle sizes, and the number of particles in different gradation intervals was calculated accordingly.

Based on the above method, a sandy soil layer of 4 m ^* 0.5 m was established with an initial porosity of 0.20. In order to reduce computational complexity, the particles with a diameter larger than 2 mm were generated using the two templates mentioned above in a 1:1 ratio, while the particles smaller than 2 mm were generated using a spherical model, but the calculation method for the number of particles remained the same as before. Additionally, the particles smaller than 0.25 mm were omitted ^[13,14] because these particles mainly serve to fill voids and do not significantly affect the mechanical properties of the structure.

2.2.   Establishment of contact model

For the contact between the drum and the particles, as well as between particles, a linear contact model with rolling resistance was used in this study. This model is based on the linear contact model (Figure 3) with the addition of stiffness to resist rolling and is suitable for systems where relative rotation between particles is more pronounced. The model has six main parameters: The stiffness and damping coefficients in the tangential and normal directions (${k}_{s}$, ${k}_{n}$, ${\beta }_{s}$, and ${\beta }_{n}$), the friction coefficient (μ) and the rolling resistance coefficient (${\mu }_{r}$) ^[15].

Figure 3.  Schematic of linear contact model.

DownLoad: Full-Size Img PowerPoint

The contact parameters between particles were calibrated through triaxial compression tests. The simulated stress-strain curve was fitted to the stress-strain curve obtained from laboratory tests to determine the parameters. In PFC2D, a particle group with dimensions of 50 mm ^* 100 mm was generated to meet the requirements of the triaxial compression test size and the maximum particle size not exceeding 1/10 of the sample size specified in the standard. Since the walls in PFC cannot directly apply forces, the confining pressure was applied using a servo mechanism. By controlling the velocity of the wall with servo control, the particle group was subjected to a confining pressure of 1MPa to achieve consolidation and stability. After stabilization, the top and bottom walls were closed with servo control and compressed toward each other. The compression speed of the top and bottom walls was related to the height of the sample, with lower heights corresponding to slower speeds to ensure that the sample was not rapidly destroyed due to excessive compression speed. Additionally, during compression, contact forces and displacements of the two walls were recorded, and stress and strain were calculated accordingly to obtain the corresponding test curve.

Figure 4 shows the three-axis compression test and simulation models. The stress-strain curve shown in Figure 5 was simulated using the contact parameters in Table 3. Comparing the simulation results with the experimental results, it can be seen that the simulated curve is in good agreement with the experimental values. Therefore, this set of contact parameters can be used for subsequent simulation analyses.

Figure 4.  Triaxial compression test and PFC2D simulation image.

DownLoad: Full-Size Img PowerPoint

Figure 5.  Comparison of simulation and experimental results of triaxial compression test.

DownLoad: Full-Size Img PowerPoint

Table 3.  Contact parameters between particles.

Contact parameters	${k}_{s}$	${k}_{n}$	${\beta }_{s}$	${\beta }_{n}$	μ	${\mu }_{r}$
Values	1.0 * 108	0.5 * 108	0.7	0.7	0.4	0.2

---

 | Show Table

DownLoad: CSV

3.   Simulation tests and analysis of vibratory compaction

3.1.   Arrangement of measuring cycles

A simulation model of vibratory compaction with a width and height of 4 m ^* 0.5 m was constructed using PFC2D. The working parameters of the drum are specified in Table 4, and the contact parameters between the drum and particles are shown in Table 5.

Table 4.  Working parameters of the drum.

Diameter/m	Mass/kg	Eccentric force/kN	Vibratory frequency/Hz
1.6	10,000	150	29

---

 | Show Table

DownLoad: CSV

Table 5.  Contact parameters between the drum and particles.

Contact parameters	${k}_{s}$	${k}_{n}$	${\beta }_{s}$	${\beta }_{n}$	μ	${\mu }_{r}$
Values	1011	1011	0.7	0.7	0.4	0.2

---

 | Show Table

DownLoad: CSV

Before simulating the vibratory compaction, static rolling was first simulated to achieve a certain degree of soil compaction, and then the vibrating wheel was centered for single-point vibratory compaction. It should be noted that in practice, the compactor is rolling and a certain point or area of the compacted layer will not undergo such continuous compaction. However, in the process of mechanism research, changes in various aspects of the soil at different times need to be explored, so the single-point loading was adopted to target a limited area for greater accuracy.

To observe the mechanical behavior of different parts of the soil under the action of vibration, 25 measurement circles were arranged in the center (x = 2 m) or on both sides below the drum, with each circle having a radius of 1 cm. The changes in horizontal stress (stressxx), vertical stress (stressyy), and shear stress (stressxy) are monitored through the measurement circles. The specific positions of the measurement circles and the corresponding monitoring indicators are numbered as shown in Table 6. The vibration model is shown in Figure 6.

Table 6.  Measurement circle information.

Numbers of the measuring cycles	Positions		Numbers of the monitoring indicators
	x(m)	y(m)	stressxx hist ID	stressyy hist ID	stressxy hist ID
1	1.00	0.45	5	30	55
2	1.00	0.35	6	31	56
3	1.00	0.25	7	32	57
4	1.00	0.15	8	33	58
5	1.00	0.05	9	34	59
6	1.50	0.45	10	35	60
7	1.50	0.35	11	36	61
8	1.50	0.25	12	37	62
9	1.50	0.15	13	38	63
10	1.50	0.05	14	39	64
11	2.00	0.45	15	40	65
12	2.00	0.35	16	41	66
13	2.00	0.25	17	42	67
14	2.00	0.15	18	43	68
15	2.00	0.05	19	44	69
16	2.50	0.45	20	45	70
17	2.50	0.35	21	46	71
18	2.50	0.25	22	47	72
19	2.50	0.15	23	48	73
20	2.50	0.05	24	49	74
21	3.00	0.45	25	50	75
22	3.00	0.35	26	51	76
23	3.00	0.25	27	52	77
24	3.00	0.15	28	53	78
25	3.00	0.05	29	54	79

---

 | Show Table

DownLoad: CSV

Figure 6.  Vibration simulation model.

DownLoad: Full-Size Img PowerPoint

3.2.   Stress distribution below the center of the drum

3.2.1.   Vertical stress

Figures 7 and 8 show the vertical stress change curves at depths of 0.05, 0.25, and 0.45 m below the center of the vibrating wheel. It can be observed that:

Figure 7.  Vertical stress below the center of the vibrating wheel.

DownLoad: Full-Size Img PowerPoint

Figure 8.  Stress wave propagation phenomenon.

DownLoad: Full-Size Img PowerPoint

1) In Figure 8, there is a time lag between the stress state of the upper soil and that of the lower soil, which is a proof of the stress wave theory. In addition, fluctuations in the waveform can be observed at the peaks and valleys, which is a phenomenon of dispersion and indicates that the fluctuations are the result of the combined action of harmonics at different frequencies. This pattern is present in all subsequent curve diagrams.

2) Within 0–1 s, the stress on the top layer is much greater than that on the middle and bottom layers. This is mainly due to two factors: Smaller contact area and great energy dissipation. In the initial stage of vibratory compaction, the contact area between the vibrating wheel and the soil surface is small, resulting in high initial stress. The energy dissipation of stress waves is large due to insufficient soil density in the initial state, but the force transmitted to the lower soil is mainly elastic waves with smaller amplitudes. The different amplitudes in the curve in Figure 8 can prove this point.

3) After 1 s of vibration, the stress on the upper soil decreases because the soil undergoes a certain degree of compression after vibration, causing the vibrating wheel to sink and increasing the contact area between the wheel and soil, leading to a decrease in stress. In addition, the vibration force on the middle and lower soils keeps increasing because the soil density continues to increase and the energy dissipation of the stress wave attenuation gradually decreases, resulting in an increase in force on the lower structure.

3.2.2.   Horizontal stress

The changes and distribution characteristics of horizontal stress are more complex. Figure 9 shows that:

Figure 9.  Horizontal stress below the center of the vibrating wheel.

DownLoad: Full-Size Img PowerPoint

1) The horizontal stress on the surface layer is compressive stress, and the magnitude is relatively large. However, with the continuous vibration, the compressive stress gradually decreases, which is because of the same reason for the change in vertical stress on the surface layer. The sharp changes in the curve values in the figure are due to particle rearrangement, which should be excluded when analyzing the regularity of the continuous change of stress.

2) The stress of the measurement circle at the bottom shows a gradually increasing trend throughout the process, indicating that as the vibratory compaction proceeds, the upper soil gradually becomes denser, and both the elastic stiffness and plastic limit increase continuously. The attenuation of plastic waves in the vibration wave decreases, resulting in greater vibration force.

3) After 1 s of vibration, the distribution trend of horizontal stress is different from the traditional mechanical understanding. For general pavement structure stress, it is usually believed that the maximum horizontal stress is at the bottom of the layer and under tension. Because the structural layer is regarded as a plate, and when subjected to vertical load, it will undergo relative bending and tension, resulting in maximum tensile stress at the bottom of the layer. However, under the vibration of sandy soil structure layers, their stress characteristics are different. The areas at the top and bottom are subject to greater compressive stress horizontally, while those in the middle have lower compressive stress horizontally. This trend is similar to the distribution of horizontal force in the vertical direction in splitting tests, and such a distribution trend only occurs in the later stage of vibration. This indicates that in the early stage, the entire soil is relatively loose, and the transfer of horizontal stress in the vertical direction still shows a decreasing trend. However, in the later stage, the overall integrity of the pavement structure increases. Under vibration, vertical compression produces a tendency of lateral expansion. This tendency to expand decreases the compressive stress in the transverse direction, resulting in such distribution for horizontal stress.

3.2.3.   Shear stress

From the shear stress distribution curve in Figure 10, it can be observed that:

Figure 10.  Shear stress below the center of the vibrating wheel.

DownLoad: Full-Size Img PowerPoint

1) The magnitude of shear stress is relatively small compared to vertical and horizontal stresses. The overall shear stress shows a decreasing trend within the depth range. The shear stress at the top is relatively large, and the direction of force changes sharply, resulting in significant shear deformation. The shear stress at the bottom is relatively small, and the direction of force has basically not changed, indicating that the shear deformation of the particles at the bottom is small.

2) In the later stage of vibration, the amplitude of shear stress on the top and middle soils decreases, indicating a reduction in shear deformation behavior and relatively stable particle groups.

3.3.   Stress distribution of the soil surface

3.3.1.   Vertical stress

As shown in Figure 11, the distribution trend of vertical stress of the soil surface is relatively simple. Similar to before, the reduction in the peak value of stress is due to the increase in the contact area of the vibrating wheel. The vertical stress is larger in the middle and smaller on both sides, indicating that the transfer of vertical stress along the horizontal direction is weak and has little effect on the soil far away horizontally. However, since the surface layer is a free interface without tensile force constraint, under the alternating action of vertical stress, particles are prone to jumping. Therefore, the degree of compaction is usually not high.

Figure 11.  Vertical stress of the soil surface.

DownLoad: Full-Size Img PowerPoint

3.3.2.   Horizontal stress

As shown in Figure 12, compared with vertical stress, horizontal stress of the soil surface has a greater transfer distance and less attenuation. This indicates that when the vibrating wheel is in action, there is a strong pushing effect on the surface soil in front and behind the vibrating wheel, causing significant deformation and uplift of the surface soil. This is also a significant feature under vibration. Other changes in the magnitude of forces follow the same pattern as before and will not be repeated.

Figure 12.  Horizontal stress of the soil surface.

DownLoad: Full-Size Img PowerPoint

3.3.3.   Shear stress

As shown in Figure 13, the horizontal shear stress of the soil surface is relatively small compared to vertical and horizontal stresses. Its transfer distance is similar to that of horizontal stress. Within 1 second after the start of vibration, the shear deformation of the soil mainly occurs around the vibrating wheel, and the shear force on the soil far away is small. After 1 second, the shear force on the soil far away becomes significant, accompanied by significant uplift of the soil.

Figure 13.  Shear stress of the soil surface.

DownLoad: Full-Size Img PowerPoint

4.   Mechanism of vibratory compaction

4.1.   Vibratory compaction theory

Regarding the mechanism of vibration compaction, four mainstream theories have been proposed in previous studies ^[16,17]:

1) Resonance theory. Since the vibration compaction output is the eccentric force with a certain vibration frequency, when the vibration frequency is close to the natural frequency of the compacted material, resonance occurs and particles are rearranged, resulting in higher compaction density.

2) Repeated loading theory. It mainly explains that under periodic compression movement, materials are subjected to repeated loading caused by vibration to achieve compaction.

3) Alternating shear strain theory. This theory believes that the alternating shear stress generated by vibration causes soil to undergo shear deformation, resulting in particle rearrangement and increased compaction density.

4) Reduction of internal friction theory. This theory believes that the shock wave generated by the vibration load can continuously propagate and diffuse along the depth direction of compacted material, making the static friction between particles into dynamic friction, causing movement and rearrangement. At this time, the internal friction between particles under motion state decreases sharply, making compaction easier.

4.2.   Validation and analysis of vibratory compaction theory

According to the simulation results of vibration compaction, the four vibratory compaction theories are verified, and the following conclusions can be obtained:

First, the resonance theory has certain limitations. The main idea of this theory is that the particle rearrangement is caused by resonance. However, from the stress-strain curve of the bottom layer shown in Figure 6, it can be observed that the stress at the bottom shows a relatively smooth wave curve only in the later stage. This indicates that there are many influencing factors in the early stage, while the overall structural integrity improves in the later stage, so the wave frequency becomes relatively stable. Therefore, even if resonance occurs, the resonance of the entire structure layer should occur in the later stage of compaction. In fact, during the compaction process, the parameters of the soil change with the increase of compaction degree, which will inevitably change the natural frequency of the soil. Therefore, the occurrence of resonance is a temporary state and is not applicable to the entire compaction process.

The repeated loading theory and the alternating shear strain theory can explain the law of vibration compaction. The particle rearrangement is due to sliding and rolling between particles, which is related to shear force. Therefore, cyclic shear strain is the cause of material compaction. In addition, there is a similar rule in the above curves: After stress fluctuates at a certain level for a period of time, it changes. This can explain that under repeated loading, the original stable structure is destroyed, and particles rearrange to form a new stable structure. This is also one of the reasons why vibration compaction is better than static compaction. Vibration loading can destroy the original stable structure, but under static loading due to the self-locking effect, increasing load does not have a significant effect on enhancing compaction.

The internal friction reduction theory is not clearly reflected in this model. According to the above results, during compaction, the vertical shear stress gradually decreases with time, but at the same time, the settlement of the vibration roller also gradually decreases. Therefore, it is difficult to directly verify the promoting effect of internal friction reduction on compaction.

5.   Conclusions

To further elucidate the compaction mechanism of the subgrade material during compaction, this paper established a discrete element particle model in line with the actual state of the subgrade soil, and analyzed the force characteristics of the area below the center of the vibrating wheel and the surface area of the soil body under vibration. Finally, the compaction mechanism of the subgrade soil under the corresponding state was derived through the distribution characteristics of various stresses of the subgrade soil. The conclusions are as follows:

1) The distribution of vertical stress below the drum center shows a decreasing trend in depth range. During vibration, the stress amplitude of the lower structure continues to increase because the plastic wave part of the vibration wave has less loss in the upper structure after gradual compaction. The plastic limit dynamic strain of the soil is elevated, so it can withstand higher loads. However, the stress of the upper structure decreases due to the increase in the contact area of the wheel.

2) The distribution of horizontal stress below the drum center is similar to the stress distribution of splitting tests, which is different from the distribution of horizontal force under traditional static loads. This indicates that non-cohesive soil does not have plate characteristics during compaction.

3) The soil at the surface is subjected to significant pushing effect. The horizontal stress transmission distance is greater than that of vertical stress. This feature, coupled with shear action, causes a certain degree of uplift deformation of the soil around the drum. This phenomenon increases the contact area between the vibration wheel and soil, and the subsequent vibration force decreases accordingly.

4) Regarding the theory of compaction process, the resonance theory has certain limitations and cannot reflect the characteristics of the entire compaction process. The repeated loading theory and alternating shear strain theory can explain the law of vibration compaction. The internal friction reduction theory is not clearly reflected in this model. Further studies are needed to investigate the compaction mechanism of different types of subgrade soils under actual compaction conditions.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

This work is supported by The Fundamental Research Funds for the Central Universities: 2242023K30042, The Fundamental Research Funds for the Central Universities: 2242023K40017 and the Natural Science Foundation of Shanxi Province, China: 2022JQ-743.

Conflict of interest

The authors declare there is no conflict of interest.