Numerical Simulation of Fatigue Life of Rubber Concrete on the Mesoscale

Rubber concrete (RC) exhibits high durability due to the rubber admixture. It is widely used in a large number of fatigue-resistant structures. Mesoscale studies are used to study the composition of polymers, but there is no method for fatigue simulation of RC. Therefore, this paper presents a finite element modeling approach to study the fatigue problem of RC on the mesoscale, which includes the random generation of the main components of the RC mesoscale structure. We also model the interfacial transition zone (ITZ) of aggregate mortar and the ITZ of rubber mortar. This paper combines the theory of concrete damage to plastic with the method of zero-thickness cohesive elements in the ITZ, and it is a new numerical approach. The results show that the model can simulate reasonably well the random damage pattern after RC beam load damage. The damage occurred in the middle of the beam span and tended to follow the ITZ. The model can predict the fatigue life of RC under various loads.


Introduction
With rapid economic development, the production of cars has increased, leading to the pollution of many waste tires, which are the primary source of waste rubber [1]. The combination of rubber, an excellent elastic material, and concrete, a brittle material, produces rubber concrete (RC), which has the advantages of low modulus of elasticity, high resistance to deformation, good crack resistance, good flexibility, and good wear resistance [2][3][4]. Liu et al. [2] found that RC improved concrete toughness and fatigue properties. Wang et al. [5] used the sounding technique to study the development of the fatigue damage process in RC at three stress levels-0.6, 0.7, and 0.8. It is a continuous process of the cumulative increase in damage, and it is divided into three processes: crack initiation, stable extension, and destabilization damage. With the rapid development of finite element theory and computer technology, concrete research is no longer limited to experimental studies. The method of finite element numerical simulation has become the primary research tool. Liu et al. [2] studied a mesoscale model of RC and analyzed its compressive properties. However, the model considered factors so simple that the results were unconvincing. Many scholars [6,7] have analyzed RC on microscopic, mesoscopic, and macroscopic scales, but no one has used a finite element model (FEM) to study RC fatigue. For this reason, this paper propose a finite element fatigue model of RC on the mesoscale.

Modeling Methods
The mesoscale model's modeling approach begins by considering the geometry of the model generation, which includes aggregate content, size, and location. Subsequently, the constitutive model of concrete and the ITZ model is considered to make the model feasible.

Mesoscale Model Geometry Generation
In mesoscale studies, concrete is usually considered a three-phase material consisting of aggregate, mortar, and ITZ. The rubber particles in RC replace part of the fine aggregates. Concrete is regarded as a homogeneous material in conventional macroscopic concrete FEM. This visual modeling approach makes it difficult to investigate how the concrete's inhomogeneity affects the macroscopic properties. This assumption of the homogeneity of concrete ignored several vital influences, such as aggregate size, particle size distribution, aggregate shape, and the effect of the ITZ. Zhong et al. [30] investigated the effect of aggregate shape (circular, elliptical, and polygonal) on the results of numerical analysis of the mesoscale model. They compared the stress-strain curves under different conditions with the experimental results. The results showed that the circular aggregate model is optimal for the numerical simulations. In this paper, circular aggregates are used so that meshing is easier and computer solutions are faster. In contrast, irregularly shaped aggregates are very complex to mesh and increase the computational burden.
In this paper, the coarse aggregates in the RC mesoscale model are aggregates of 5 mm or more in diameter and the fine aggregates are included in the mortar. The geometry of the mesoscale model needs to comply with three requirements: firstly, all the particles generated must be within the specified boundaries; secondly, none of the particles can overlap; and thirdly, there must be a gap between each particle, as the aggregates are wrapped in a layer of mortar and have no contact. Fuller's particle size [31] distribution curves are used in this paper. The Fuller curve is widely regarded as the grading curve, which provides an optimum particle size distribution for the working condition of the concrete. The Fuller curve equation is as follows: where P represents the percentage of aggregate passing through sieve hole diameter D 0 , D 0 represents the diameter of the sieve hole, and D max represents the diameter of the largest aggregate. As this paper focuses on the 2D level, it is impossible to deal with the 2D problem directly with the help of the 3D Fuller set matching formula. It was used to obtain the best particle size distribution curve in 2D by applying the Walraven formula [32]. The formula is as follows: where P c represents the percentage of the aggregate area, where size D is smaller than D 0 . P k represents the percentage of the aggregate area of the total area. In this paper, P k is taken to be 0.7. D max represents the diameter of the largest aggregate size, and the maximum diameter is taken as 20 mm. The area of aggregate size distribution in the 550 mm × 150 mm area is listed by Formula (2) in Table 1. The random generation is implemented in Python according to the aggregate area in Table 1, and the generated flowchart is shown in Figure 1.  The random generation is implemented in Python according to the aggregate area in Table 1, and the generated flowchart is shown in Figure 1.

Constitutive Model of Concrete
The finite element software ABAQUS (2021 Version, Dassault systemes, France) and the programming language Python 3.8 are interconnected, and the code generated in Section 2.1 can be imported directly into ABAQUS. Aggregates are developed according to the program, and different property values are assigned to the different components to achieve the actual state of the mesoscale RC aggregates. In this paper, coarse aggregate and rubber are considered homogeneous elastomers, and mortar is modeled numerically as a homogeneous continuum with elasticity. The mortar can be regarded as a lower-strength type of concrete, and its constitutive law uses the concrete-damage-plasticity (CDP) model. The CDP is a continuous, plasticity-based damage model that defines the concrete state by defining two mechanical behaviors: tensile cracking and compression damage. This model assumes that the concrete's uniaxial tensile and compressive response is characterized by plastic damage. The evolution of the yield surface is controlled by two hardening variables, the tensile equivalent plastic strain ε pl t and the compressive equivalent plastic strain ε pl c , which are related to the damage mechanisms under tensile and compressive loading. The uniaxial tensile and compressive stress-strain response is shown in Figure 2.

Constitutive Model of Concrete
The finite element software ABAQUS (2021 Version, Dassault systemes, France) and the programming language Python 3.8 are interconnected, and the code generated in Section 2.1 can be imported directly into ABAQUS. Aggregates are developed according to the program, and different property values are assigned to the different components to achieve the actual state of the mesoscale RC aggregates.
In this paper, coarse aggregate and rubber are considered homogeneous elastomers, and mortar is modeled numerically as a homogeneous continuum with elasticity. The mortar can be regarded as a lower-strength type of concrete, and its constitutive law uses the concrete-damage-plasticity (CDP) model. The CDP is a continuous, plasticity-based damage model that defines the concrete state by defining two mechanical behaviors: tensile cracking and compression damage. This model assumes that the concrete's uniaxial tensile and compressive response is characterized by plastic damage. The evolution of the yield surface is controlled by two hardening variables, the tensile equivalent plastic strain ̃ and the compressive equivalent plastic strain ̃, which are related to the damage mechanisms under tensile and compressive loading. The uniaxial tensile and compressive stress-strain response is shown in Figure 2. Where ̃ and ̃ represent the equivalent plastic strains in tension and compression, and represent the elastic strains corresponding to tension and compression, and represent the two damage variables for elastic stiffness degradation, with the damage variables taking values from 0 to 1, and = 0 means the material is undamaged, = 1 means the material is completely damaged, and 0 represents the initial Young's modulus of the material.
In the case of uniaxial tension, the concrete stress-strain response obeys linear elastic variation up to the time of failure stress 0 , which is used to distinguish between the elastic and plastic phases of concrete. After 0 , the concrete enters the damage phase, and microcracking occurs in the macrostructure. In the case of uniaxial compression, the concrete stress-strain response obeys a linear elastic change to the compressive elastic ultimate stress, and 0 , and 0 used to distinguish the elastic phase from the plastic phase under uniaxial compression. Unlike uniaxial tension, there is a hardening phase to the ultimate compressive stress after where the concrete is softened and microcracked.
When a concrete specimen is unloaded from any point in the strain-softening branch of the stress-strain curve, the elastic stiffness of the material appears to be damaged. The stress-strain relationships for uniaxial tensile and compressive loading are (Equation (3)): Where ε pl t and ε pl c represent the equivalent plastic strains in tension and compression, ε el t and ε el c represent the elastic strains corresponding to tension and compression, d t and d t represent the two damage variables for elastic stiffness degradation, with the damage variables taking values from 0 to 1, and d = 0 means the material is undamaged, d = 1 means the material is completely damaged, and E 0 represents the initial Young's modulus of the material.
In the case of uniaxial tension, the concrete stress-strain response obeys linear elastic variation up to the time of failure stress σ t0 , which is used to distinguish between the elastic and plastic phases of concrete. After σ t0 , the concrete enters the damage phase, and microcracking occurs in the macrostructure. In the case of uniaxial compression, the concrete stress-strain response obeys a linear elastic change to the compressive elastic ultimate stress, and σ c0 , and σ c0 used to distinguish the elastic phase from the plastic phase under uniaxial compression. Unlike uniaxial tension, there is a hardening phase to the ultimate compressive stress σ cu after σ cu where the concrete is softened and microcracked.
When a concrete specimen is unloaded from any point in the strain-softening branch of the stress-strain curve, the elastic stiffness of the material appears to be damaged. The stress-strain relationships for uniaxial tensile and compressive loading are (Equation (3)): This paper deals with the numerical simulation of the three-point bending of concrete, where the general form of damage is tensile damage. After being subjected to cyclic loading, the tensile stiffness after damage needs to be redefined, as shown in Figure 3.
This paper deals with the numerical simulation of the three-point bending of concrete, where the general form of damage is tensile damage. After being subjected to cyclic loading, the tensile stiffness after damage needs to be redefined, as shown in Figure 3. where 0 represents elastic strain and ̃ represents cracking strain.
In the event of damage to the concrete, the cracking strain ̃ is defined by the following equation: ABAQUS automatically converts cracking strain to plastic strain for use: According to the Structural Design Code for Concrete [33], considering the tensile and compressive damage variables of the material, the specific concrete intrinsic model is determined by Young's modulus E and Poisson's ratio λ in the elastic phase and by the non-linear stress-strain equation in the inelastic phase. When the concrete structure is under pressure: Where ε el 0t represents elastic strain and ε ck t represents cracking strain. In the event of damage to the concrete, the cracking strain ε ck t is defined by the following equation: ABAQUS automatically converts cracking strain to plastic strain for use: According to the Structural Design Code for Concrete [33], considering the tensile and compressive damage variables of the material, the specific concrete intrinsic model is determined by Young's modulus E and Poisson's ratio λ in the elastic phase and by the non-linear stress-strain equation in the inelastic phase. When the concrete structure is under pressure: n = E c ε c,r E c ε c,r − f c,r (12) where α c is the parameter value of the falling section of the uniaxial compressive stress−strain curve for concrete, f c,r is the representative value of the uniaxial compressive strength of concrete, ε c,r is the peak compressive strain corresponding to f c,r , and d c is the evolutionary parameter for uniaxial compressive damage to concrete. When the concrete structure is in tension: where α t is the parameter value of the falling section of the uniaxial tension stress−strain curve for concrete, f t,r is the representative value of the uniaxial tension strength of concrete, ε t,r is the peak tension strain corresponding to f t,r , and d t is the evolutionary parameter for uniaxial tension damage to concrete.
In accordance with the Structural Design Code for Concrete [33], specific values are shown in Table 2. In addition to this, the plasticity parameters for CDP are self-contained in ABAQUS, as shown in Table 3. The values in Table 3 have been verified by many academics to be generally consistent. This is a fixed value [10].

CE Model of the ITZ
After the aggregate model has been built, there are two approaches to the ITZ. One is establishing a solid FEM of the ITZ [34]. The advantage of this is that it can reflect the thickness relationship of the interface composed of concrete. However, the ITZ's actual thickness is 10-50 µm, which is difficult to achieve with FEM. Even if the thickness is expanded by a factor of 10 to a range that FEM can calculate, this will result in a dense and small mesh division and a significant increase in computational effort. Secondly, the ITZ is considered a zero-thickness element (ZTE) [35], which retains the relevant mechanical properties of the actual ITZ to achieve the accuracy of the simulation, and all ITZs in concrete can be represented by ZTE. In summary, we selected the ZTE.
The ZTE has three ways of simulating the behavior of the ITZ. Firstly, a layer of the ZTE can be inserted using a shared node, which can be used if the CE is on the same mesh as the surrounding element. Secondly, if the elements of the CE are divided differently from the surrounding mesh or if the CE uses a finer discretization than the adjacent parts, the tie constraint can be used to connect the CE to other parts. Thirdly, in some special cases where the requirements are met, a connected interaction can be added directly to the CE in contact without adding additional elements. Figure 4 shows the three methods of CE processing. ZTE can be inserted using a shared node, which can be used if the CE is on the same mesh as the surrounding element. Secondly, if the elements of the CE are divided differently from the surrounding mesh or if the CE uses a finer discretization than the adjacent parts, the tie constraint can be used to connect the CE to other parts. Thirdly, in some special cases where the requirements are met, a connected interaction can be added directly to the CE in contact without adding additional elements. Figure 4 shows the three methods of CE processing. Based on the random aggregates generated by the simulations in this paper, a cohesive zone will be added to the contact surface of the aggregate and mortar. The first ZTE ( Figure 4a) is chosen to insert a layer of CE using a shared node. The size and location of each aggregate are uncertain, so choosing inserted CE is difficult. A Python program finds the node number of the aggregate place and copies the new node at the node number to create a zero-thickness CE. This fits perfectly with the shared node insertion approach. The ITZ generates CE, as shown in Figure 5. The CE damage is divided into four parts: the linear elastic phase, the damage initiation phase, the damage evolution phase, and complete damage.
The online resilience phase of the damage response of the CE is as follows: Based on the random aggregates generated by the simulations in this paper, a cohesive zone will be added to the contact surface of the aggregate and mortar. The first ZTE ( Figure 4a) is chosen to insert a layer of CE using a shared node. The size and location of each aggregate are uncertain, so choosing inserted CE is difficult. A Python program finds the node number of the aggregate place and copies the new node at the node number to create a zero-thickness CE. This fits perfectly with the shared node insertion approach. The ITZ generates CE, as shown in Figure 5.
ZTE can be inserted using a shared node, which can be used if the CE is on the same mesh as the surrounding element. Secondly, if the elements of the CE are divided differently from the surrounding mesh or if the CE uses a finer discretization than the adjacent parts, the tie constraint can be used to connect the CE to other parts. Thirdly, in some special cases where the requirements are met, a connected interaction can be added directly to the CE in contact without adding additional elements. Figure 4 shows the three methods of CE processing. Based on the random aggregates generated by the simulations in this paper, a cohesive zone will be added to the contact surface of the aggregate and mortar. The first ZTE ( Figure 4a) is chosen to insert a layer of CE using a shared node. The size and location of each aggregate are uncertain, so choosing inserted CE is difficult. A Python program finds the node number of the aggregate place and copies the new node at the node number to create a zero-thickness CE. This fits perfectly with the shared node insertion approach. The ITZ generates CE, as shown in Figure 5. The CE damage is divided into four parts: the linear elastic phase, the damage initiation phase, the damage evolution phase, and complete damage.
The online resilience phase of the damage response of the CE is as follows: The CE damage is divided into four parts: the linear elastic phase, the damage initiation phase, the damage evolution phase, and complete damage.
The online resilience phase of the damage response of the CE is as follows: where t n is the nominal stress in the normal direction, t s is the nominal stress in shear in the first direction, t t is the nominal stress in shear in the second direction, E ij is Young's modulus in each direction, and ε i is the strain in the corresponding direction. The quadratic stress criterion formula is used for damage initiation. Damage initiation occurs when the contact-stress ratio involved reaches 1: Damage evolution by way of traction separation is shown in Figure 6.
where is the nominal stress in the normal direction, is the nominal stress in shear in the first direction, is the nominal stress in shear in the second direction, is Young's modulus in each direction, and is the strain in the corresponding direction. The quadratic stress criterion formula is used for damage initiation. Damage initiation occurs when the contact-stress ratio involved reaches 1: Damage evolution by way of traction separation is shown in Figure 6. where 0 is the separation displacement value at the onset of damage and is the separation displacement at the maximum damage.
The material enters the damage phase judged by the damage value D. When D is 1, the material is completely damaged by the following equation: where D is the damage value and is the additional amount of maximum separation displacement during loading.
In this paper, the fracture energy is calculated using the Benzeggagh-Kenane (BK) criterion: where is the hybrid fracture energy, is the type I fracture energy of the cohesive element, is the type II fracture energy of the cohesive element, is the shear deformation energy, and is the tensile deformation energy. The performance of ITZ is difficult to test on the experimental scale, so the determination of simulation parameters for ITZ is difficult to determine. Usually, the performance of ITZ is approximated by the weak mortar composition, and researchers use the percentage of mortar to study and judge the performance of ITZ. Xiao et al. [36] considered the strength of ITZ to be 80% of the mortar. Kim et al. [37] considered the fracture energy of ITZ to be equivalent to 50% of the mortar. Li et al. [38] considered it to be 80%. It was obvious that different researchers have different opinions on determining the mechanical Where δ 0 m is the separation displacement value at the onset of damage and δ f m is the separation displacement at the maximum damage.
The material enters the damage phase judged by the damage value D. When D is 1, the material is completely damaged by the following equation: where D is the damage value and δ max m is the additional amount of maximum separation displacement during loading.
In this paper, the fracture energy is calculated using the Benzeggagh-Kenane (BK) criterion: where G C is the hybrid fracture energy, G C n is the type I fracture energy of the cohesive element, G C s is the type II fracture energy of the cohesive element, G s is the shear deformation energy, and G T is the tensile deformation energy.
The performance of ITZ is difficult to test on the experimental scale, so the determination of simulation parameters for ITZ is difficult to determine. Usually, the performance of ITZ is approximated by the weak mortar composition, and researchers use the percentage of mortar to study and judge the performance of ITZ. Xiao et al. [36] considered the strength of ITZ to be 80% of the mortar. Kim et al. [37] considered the fracture energy of ITZ to be equivalent to 50% of the mortar. Li et al. [38] considered it to be 80%. It was obvious that different researchers have different opinions on determining the mechanical properties of ITZ. The ultimate purpose is to achieve unity between numerical simulations and experiments, so the performance parameters of the ITZ on numerical simulations are determined by trial and error to determine the optimum values of these relevant parameters. The parameters used in this paper are shown in Table 4.

Experiment
The data for this summary test were obtained from Liu et al. [39]. He investigated the effect of rubber substitution rate and rubber particle size on the fatigue life of rubber concrete. The object of study was an RC beam with dimensions of 150 mm × 150 mm × 550 mm. A static load test was conducted under a three-point bending load, and a fatigue test was conducted under a cyclic load. The fatigue life and related fatigue life curves were obtained for different rubber substitution rates, particle sizes, and stress levels.

Experimental Materials
Material parameters are shown in Tables 5 and 6.

Experimental Test Methods
RC specimens with different rubber replacement rates are first tested by static loading to obtain the corresponding peak loads. The fatigue tests are carried out using models of the same material proportions. The maximum and minimum loads are applied to the RC beams using a load-controlled mode, which is an equal amplitude and uniform load mode.

Building Mesoscale Models
A 150 mm × 550 mm rubber concrete beam element is built according to Table 1, with different dosing of rubber concrete beams as shown in Figure 7, where gray means mortar, red means aggregate, and black means rubber. Polymers 2023, 15, x 11 of 23 The load loading point is in the middle of the upper part, with the bottom left constraint 100 mm from the left boundary and the right constraint 100 mm from the right, as in Figure 8. The load loading point is in the middle of the upper part, with the bottom left constraint 100 mm from the left boundary and the right constraint 100 mm from the right, as in Figure 8.
where the period is T, circle frequency ω = 2π/T, the loading initial time is 0 , and the number of steps parameters The parameters used for RC in this paper are shown in Table 7.

Experimental Versus Simulation
By comparing the results of this study with the three-point bending static load peak load results and fatigue load results from the literature [39], the feasibility of the model is verified.

Peak Load
The peak load tests and simulation results for this model under three-point bending loads at different stress levels are summarized in Table 8. It can be seen that the magnitude of the peak load decreases as the rubber content increases, which is in line with the researchers' judgment on the performance of RC. Comparing the test and simulation for The model is solved using the ABAQUS/Standard. First, apply a displacement constraint of 2 mm at the upper load loading point, stop the calculation when the model does not converge, and obtain the peak load F max for the model. Subsequently, the fatigue life of the model is calculated at different stress levels, still using the same model with cyclic concentrated force constraints applied at the upper load loading points. Load application from minimum load P min to maximum load P max , where P min /P max = 0.1, fatigue load stress levels S = P max /F max . In this paper, S takes the values 0.9, 0.85, 0.8, and 0.75. When S is too high, the fatigue damage results are over in one go. When S is too small, the calculation is too large in the numerical simulation phase. The Fourier series method controls the equivalent mean amplitude load when fatigue loads are applied: where the period is T, circle frequency ω = 2π/T, the loading initial time is A 0 , and the number of steps parameters A 1 , B 1 , A 2 , B 2 , ···, A 0 , A 0 . The parameters used for RC in this paper are shown in Table 7.

Experimental Versus Simulation
By comparing the results of this study with the three-point bending static load peak load results and fatigue load results from the literature [39], the feasibility of the model is verified.

Peak Load
The peak load tests and simulation results for this model under three-point bending loads at different stress levels are summarized in Table 8. It can be seen that the magnitude of the peak load decreases as the rubber content increases, which is in line with the researchers' judgment on the performance of RC. Comparing the test and simulation for peak loads at the same stress levels proved that the simulation and test results agree well. The maximum absolute error is 3.6%. The new numerical model proved reliable for peak loads under three-point bending loads.

Fatigue Life
The results of tests and simulations with different dosing levels of rubber concrete at stress levels S = 0.85 and S = 0.75 are summarized in Table 9. The increase in rubber admixture can improve the fatigue resistance of RC and extend the fatigue life. Due to the large dispersion of the fatigue life results, only the minimum and maximum lives are taken as a reference in the test results. Moreover, the overall life trend improves with increasing rubber doping. As shown in Figure 9, the results obtained from the numerical model of rubber concrete in this paper are all between the maximum and minimum values of the test results and meet the feasibility requirements of the model. This proves the reliability of the new numerical model in fatigue life calculation. There are some differences between the expected life and the experiment, but this is acceptable. Because the experiment phase is a one-off for each test beam, the RC is already destroyed after the experiments with peak load. Although each beam is made to the same size and aggregate content, the mechanical properties are not the same. The different mechanical behaviors of the concrete beam can be observed in [37].

Fatigue Life
The results of tests and simulations with different dosing levels of rubber concrete at stress levels S = 0.85 and S = 0.75 are summarized in Table 9. The increase in rubber admixture can improve the fatigue resistance of RC and extend the fatigue life. Due to the large dispersion of the fatigue life results, only the minimum and maximum lives are taken as a reference in the test results. Moreover, the overall life trend improves with increasing rubber doping. As shown in Figure 9, the results obtained from the numerical model of rubber concrete in this paper are all between the maximum and minimum values of the test results and meet the feasibility requirements of the model. This proves the reliability of the new numerical model in fatigue life calculation. There are some differences between the expected life and the experiment, but this is acceptable. Because the experiment phase is a one-off for each test beam, the RC is already destroyed after the experiments with peak load. Although each beam is made to the same size and aggregate content, the mechanical properties are not the same. The different mechanical behaviors of the concrete beam can be observed in [37].

Analysis of Variables
Finite element static pressure simulations of three-point bending were carried out for RC rubber admixtures of 0, 2.5%, 5%, 7.5%, and 10% to obtain the corresponding peak load and displacement relationships. Based on the stress levels S = 0.85 and S = 0.75 above, add stress levels S = 0.9 and S = 0.8 to the loading method for the RC fatigue simulation to analyze the effect of different rubber doping and stress levels on damage form and fatigue life.

Force-Deflection Curves
The force-deflection curves for RC at different admixtures were obtained from numerical simulations, as shown in Figure 10. The concrete deflection increases as the rubber admixture increases and the peak load tends to decrease significantly. The rising and falling phases of the curve for ordinary concrete are steeper than the gentle curve for RC. The trend becomes more subdued as the amount of rubber added increases. The comparison of the trends of the two curves RC-0 and RC-10 in Figure 10 is exceptionally different. It confirmed the effect of rubber particles on concrete in the mesoscale study. Rubber was able to reduce the extension of concrete damage and increase the toughness of concrete, reducing the brittleness of concrete. This reflects the actual validity of the new numerical model.

Analysis of Variables
Finite element static pressure simulations of three-point bending were carried out for RC rubber admixtures of 0, 2.5%, 5%, 7.5%, and 10% to obtain the corresponding peak load and displacement relationships. Based on the stress levels S = 0.85 and S = 0.75 above, add stress levels S = 0.9 and S = 0.8 to the loading method for the RC fatigue simulation to analyze the effect of different rubber doping and stress levels on damage form and fatigue life.

Force-Deflection Curves
The force-deflection curves for RC at different admixtures were obtained from numerical simulations, as shown in Figure 10. The concrete deflection increases as the rubber admixture increases and the peak load tends to decrease significantly. The rising and falling phases of the curve for ordinary concrete are steeper than the gentle curve for RC. The trend becomes more subdued as the amount of rubber added increases. The comparison of the trends of the two curves RC-0 and RC-10 in Figure 10 is exceptionally different. It confirmed the effect of rubber particles on concrete in the mesoscale study. Rubber was able to reduce the extension of concrete damage and increase the toughness of concrete, reducing the brittleness of concrete. This reflects the actual validity of the new numerical model.

Types of Damage
In this paper, the damage to the RC is shown through stiffness in the form of two factors: one is static compression, and the other is fatigue. The visual form of the damage is represented by the SDEG cloud map output by ABAQUS (SDEG = 0 means no damage to the structure, and SDEG = 1 means complete damage to the structure). Figure 11 shows the damage to an ordinary concrete beam of 150 mm × 550 mm without rubber admixture after a three-point bending static load. As the beam damage occurs in the middle of the span, for ease of observation, the structure is taken in the middle of the beam, as shown in the black box in Figure 11, with a size of 150 mm × 150 mm. The following are screenshots of the damage obtained by this method.

Types of Damage
In this paper, the damage to the RC is shown through stiffness in the form of two factors: one is static compression, and the other is fatigue. The visual form of the damage is represented by the SDEG cloud map output by ABAQUS (SDEG = 0 means no damage to the structure, and SDEG = 1 means complete damage to the structure). Figure 11 shows the damage to an ordinary concrete beam of 150 mm × 550 mm without rubber admixture after a three-point bending static load. As the beam damage occurs in the middle of the span, for ease of observation, the structure is taken in the middle of the beam, as shown in the black box in Figure 11, with a size of 150 mm × 150 mm. The following are screenshots of the damage obtained by this method.
The SDEG damage clouds for 0, 2.5%, 5%, 7.5%, and 10% rubber doping after damage are shown in Figure 12. A form of static pressure damage to rubber concrete was observed in the mesoscale study. The damage was mainly at the mid-span of the beam, with an irregular damage zone extending from the bottom to the top. The damage course follows the edges of the aggregate and rubber and is consistent with existing fracture and damage mechanics theories. The point of damage to the zero rubber-doped concrete is only at the opening of the damage zone, with no damage to the surrounding concrete aggregate, as shown in Figure 12a. Damage points occur not only at the opening of the damage zone but also minor damage to the rubber around the opening, as shown in Figure 12b-e. On the mesoscale, it is observed that the rubber particles take up a small part of the load-bearing capacity under load. Furthermore, with the increase of rubber admixture, the damage point at the bottom of the concrete increases, and the damage zone is influenced by the surrounding rubber particles in the middle of the extension. The 10% and 7.5% rubberdoped concrete leads particularly well, with multiple damage points at the bottom and tiny branches of the damage zone midway through, as shown in Figure 12d,e. Various forms of damage indicate that adding rubber particles to concrete helps to retard concrete damage, which also provides the basis for research into the fatigue resistance of RC. The SDEG damage clouds for 0, 2.5%, 5%, 7.5%, and 10% rubber doping after damage are shown in Figure 12. A form of static pressure damage to rubber concrete was observed in the mesoscale study. The damage was mainly at the mid-span of the beam, with an irregular damage zone extending from the bottom to the top. The damage course follows the edges of the aggregate and rubber and is consistent with existing fracture and damage mechanics theories. The point of damage to the zero rubber-doped concrete is only at the opening of the damage zone, with no damage to the surrounding concrete aggregate, as shown in Figure 12a. Damage points occur not only at the opening of the damage zone but also minor damage to the rubber around the opening, as shown in Figure 12b-e. On the mesoscale, it is observed that the rubber particles take up a small part of the loadbearing capacity under load. Furthermore, with the increase of rubber admixture, the damage point at the bottom of the concrete increases, and the damage zone is influenced by the surrounding rubber particles in the middle of the extension. The 10% and 7.5% rubber-doped concrete leads particularly well, with multiple damage points at the bottom and tiny branches of the damage zone midway through, as shown in Figure 12d,e. Various forms of damage indicate that adding rubber particles to concrete helps to retard concrete damage, which also provides the basis for research into the fatigue resistance of RC.
(a) (b) Figure 11. Damage to ordinary concrete and location of damage sampling. The SDEG damage clouds for 0, 2.5%, 5%, 7.5%, and 10% rubber doping after damage are shown in Figure 12. A form of static pressure damage to rubber concrete was observed in the mesoscale study. The damage was mainly at the mid-span of the beam, with an irregular damage zone extending from the bottom to the top. The damage course follows the edges of the aggregate and rubber and is consistent with existing fracture and damage mechanics theories. The point of damage to the zero rubber-doped concrete is only at the opening of the damage zone, with no damage to the surrounding concrete aggregate, as shown in Figure 12a. Damage points occur not only at the opening of the damage zone but also minor damage to the rubber around the opening, as shown in Figure 12b-e. On the mesoscale, it is observed that the rubber particles take up a small part of the loadbearing capacity under load. Furthermore, with the increase of rubber admixture, the damage point at the bottom of the concrete increases, and the damage zone is influenced by the surrounding rubber particles in the middle of the extension. The 10% and 7.5% rubber-doped concrete leads particularly well, with multiple damage points at the bottom and tiny branches of the damage zone midway through, as shown in Figure 12d,e. Various forms of damage indicate that adding rubber particles to concrete helps to retard concrete damage, which also provides the basis for research into the fatigue resistance of RC. The fatigue simulation of the same RC beam at different stress levels is a unique advantage of the fine-view simulation. The stress levels are guaranteed to be the same peak load each time, something that cannot be achieved experimentally. This paper uses fatigue simulations for four stress levels of 0.9, 0.85, 0.8, and 0.75, with each stress level corresponding to five rubber doping levels. Shown in Figure 13 are four forms of stress level fatigue damage for ordinary concrete. It can be observed that fatigue damage to ordinary The fatigue simulation of the same RC beam at different stress levels is a unique advantage of the fine-view simulation. The stress levels are guaranteed to be the same peak load each time, something that cannot be achieved experimentally. This paper uses fatigue simulations for four stress levels of 0.9, 0.85, 0.8, and 0.75, with each stress level corresponding to five rubber doping levels. Shown in Figure 13 are four forms of stress level fatigue damage for ordinary concrete. It can be observed that fatigue damage to ordinary concrete at different stress levels takes the same form, with the damage zone starting at the same point of failure at the bottom of the concrete. In ordinary concrete, from the start of the damage to the end, only the weakest point within the concrete bears the load, regardless of the force acting. Its fatigue damage is also essentially the same as static pressure damage (Figures 12a and 13), proving that ordinary concrete is relatively homogeneous regarding internal forces when damaged, with the same place bearing the load. The fatigue simulation of the same RC beam at different stress levels is a unique advantage of the fine-view simulation. The stress levels are guaranteed to be the same peak load each time, something that cannot be achieved experimentally. This paper uses fatigue simulations for four stress levels of 0.9, 0.85, 0.8, and 0.75, with each stress level corresponding to five rubber doping levels. Shown in Figure 13 are four forms of stress level fatigue damage for ordinary concrete. It can be observed that fatigue damage to ordinary concrete at different stress levels takes the same form, with the damage zone starting at the same point of failure at the bottom of the concrete. In ordinary concrete, from the start of the damage to the end, only the weakest point within the concrete bears the load, regardless of the force acting. Its fatigue damage is also essentially the same as static pressure damage (Figures 12a and 13), proving that ordinary concrete is relatively homogeneous regarding internal forces when damaged, with the same place bearing the load. show fatigue damage at four stress levels for four doped RC. Unlike ordinary concrete, the fatigue loads do not take the same form of damage at different stress levels when rubber is added. As shown in Figure 14, the damage zone for 2.5% admixture stress levels of 0.85, 0.8, and 0.75 differ, and the damage point is also different at the bottom of the concrete. As shown in Figure 15, the location of the damage zone is different for 5% doping stress levels of 0.9, 0.85, and 0.8, but the location of the bottom damage point is the same for stress levels of 0.85, 0.8, and 0.75. As shown in Figure 16, the 7.5% doping stress level only differs in the damage zone and damage at a stress level of 0.9; the damage zone and damage point are essentially the same at other stress levels. As shown in Figure 17, the orientation of the damage zone and the location of the initial damage point at the bottom stabilize and remain the same when the doping level reaches 10%. The rubber dosing ranges from 0 to 10%, with the damage zone orientation and initial damage point location stabilizing from the beginning, through the disorder of the intermediate dosing, and to final stability. The fatigue damage of RC is different from hydrostatic damage, which is also different from ordinary concrete. The addition of the rubber  Figures 14-17 show fatigue damage at four stress levels for four doped RC. Unlike ordinary concrete, the fatigue loads do not take the same form of damage at different stress levels when rubber is added. As shown in Figure 14, the damage zone for 2.5% admixture stress levels of 0.85, 0.8, and 0.75 differ, and the damage point is also different at the bottom of the concrete. As shown in Figure 15, the location of the damage zone is different for 5% doping stress levels of 0.9, 0.85, and 0.8, but the location of the bottom damage point is the same for stress levels of 0.85, 0.8, and 0.75. As shown in Figure 16, the 7.5% doping stress level only differs in the damage zone and damage at a stress level of 0.9; the damage zone and damage point are essentially the same at other stress levels. As shown in Figure 17, the orientation of the damage zone and the location of the initial damage point at the bottom stabilize and remain the same when the doping level reaches 10%. The rubber dosing ranges from 0 to 10%, with the damage zone orientation and initial damage point location stabilizing from the beginning, through the disorder of the intermediate dosing, and to final stability. The fatigue damage of RC is different from hydrostatic damage, which is also different from ordinary concrete. The addition of the rubber creates a fragile ITZ between the rubber and the mortar, even weaker than the ITZ between the aggregate and the mortar. The model successfully simulated the effect of rubber doping.  Table 10 shows the peak loads and the life of the fatigue loads at the corresponding four stress levels for the five doped rubber concretes. At the same stress level, the fatigue life increases as the rubber content increases, indicating that rubber concrete carries higher cyclic loads than ordinary concrete for a given cyclic load. In the case of rubber concrete, this result is because the rubber particles act as an energy absorber and load cushion in   Table 10 shows the peak loads and the life of the fatigue loads at the corresponding four stress levels for the five doped rubber concretes. At the same stress level, the fatigue life increases as the rubber content increases, indicating that rubber concrete carries higher cyclic loads than ordinary concrete for a given cyclic load. In the case of rubber concrete, this result is because the rubber particles act as an energy absorber and load cushion in  Table 10 shows the peak loads and the life of the fatigue loads at the corresponding four stress levels for the five doped rubber concretes. At the same stress level, the fatigue life increases as the rubber content increases, indicating that rubber concrete carries higher cyclic loads than ordinary concrete for a given cyclic load. In the case of rubber concrete, this result is because the rubber particles act as an energy absorber and load cushion in the concrete. Rubber particles have better elastic properties on the mesoscale level than concrete particles. In the case of concrete suffering from tension and compression, part of the energy is converted into the elastic energy of the rubber particles. The fatigue life of 10% rubber is 7.3, 3.89, 4.45, and 2.77 times greater than that of ordinary concrete at four stress levels. The relationship between rubber doping, stress level, and fatigue life is shown in Figure 18. It is obvious that, within a specific range, an increase in rubber content and a decrease in stress level increase the fatigue life of RC. The increase in fatigue life is a non-linear relationship, as shown in Figure 19. The relationship between rubber doping, stress level, and fatigue life is shown in Figure 18. It is obvious that, within a specific range, an increase in rubber content and a decrease in stress level increase the fatigue life of RC. The increase in fatigue life is a nonlinear relationship, as shown in Figure 19.

Discussion
This model simulation study generates the mesoscale structure of RC through random aggregates, applies the improved properties of CDP to mortar, and combines modeling of the aggregate-mortar ITZ and rubber-mortar ITZ to achieve the mesoscale structure of actual RC. It is a new numerical approach. The model was subjected to a series of three-point bending fatigue loads to analyze the causes of damage forms and fatigue life from a mesoscale.

Causes of Damage Types
The structural form of the model after damage by static pressure and fatigue loading is consistent with reality, and the appearance of concrete damage on the mesoscale can be accurately observed from the mesoscale structure. Subsequent damage develops along the weakness of the ITZ around the aggregate and rubber particles. There are two main reasons for producing an irregular damage band consistent with reality. First, the model adds the damage theory in the CE model, setting a zero-thickness damage zone in the aggre-

Discussion
This model simulation study generates the mesoscale structure of RC through random aggregates, applies the improved properties of CDP to mortar, and combines modeling of the aggregate-mortar ITZ and rubber-mortar ITZ to achieve the mesoscale structure of actual RC. It is a new numerical approach. The model was subjected to a series of three-point bending fatigue loads to analyze the causes of damage forms and fatigue life from a mesoscale.

Causes of Damage Types
The structural form of the model after damage by static pressure and fatigue loading is consistent with reality, and the appearance of concrete damage on the mesoscale can be accurately observed from the mesoscale structure. Subsequent damage develops along the weakness of the ITZ around the aggregate and rubber particles. There are two main reasons for producing an irregular damage band consistent with reality. First, the model adds the damage theory in the CE model, setting a zero-thickness damage zone in the aggregate-mortar and rubber-mortar layers. The material's mechanical properties in the CE are less than those of the mortar. When subjected to forces, the ITZ is more easily damaged than the mortar aggregate. Secondly, the mesoscale model is randomly generated for aggregate size and location, which aligns with the actual aggregate distribution of concrete materials and better reflects the model's realism. The rubber particles are smaller than the coarse aggregate, and it is easier for damage to occur around the rubber than around the coarse aggregate.

Factors Influencing Fatigue Life
The addition of rubber benefits the fatigue properties of rubber concrete. When the model is damaged, there is some damage around the rubber particles not in the damage zone. These share some of the fatigue load, confirming the effect of the rubber particles on the mesoscale level. The mesoscale model can be applied to complex fatigue loads. The model shows fatigue life agreement at stress levels of 0.75 to 0.9 and can simulate the effects of fatigue life due to different doping levels of rubber. Rubber particles have better elastic properties on the mesoscale level than concrete particles. When the concrete is loaded, part of the energy is converted into the elastic energy of the rubber particles. RC life increases with increasing rubber and decreases with increasing stress ratio. As a rule of thumb, the magnitude of the stress ratio is related to the logarithm of the fatigue life [40]. After trying various fitting formulae, the following relationship is assumed: (22) where S is the concrete stress ratio, A, B, and C are constants whose magnitude is related to the concrete admixture, and N is the concrete fatigue life.
According to Formula (22). for curve fitting, as shown in Figure 20, the specific formula and correlation coefficient R results are shown in Table 11, and the fitted results meet the requirements. It is found that the stress level is related to the quadratic function of the logarithm of fatigue life. In addition, the results of this study allow for reasonable extrapolation of the three-point bending fatigue life of rubber concrete at dosing levels between 0.75 and 0.9. This provides a corresponding reference for the test, a model for calculating fatigue life correctly on the mesoscale.
where S is the concrete stress ratio, A, B, and C are constants whose magnitude is related to the concrete admixture, and N is the concrete fatigue life. According to Formula (22). for curve fitting, as shown in Figure 20, the specific formula and correlation coefficient R results are shown in Table 11, and the fitted results meet the requirements. It is found that the stress level is related to the quadratic function of the logarithm of fatigue life. In addition, the results of this study allow for reasonable extrapolation of the three-point bending fatigue life of rubber concrete at dosing levels between 0.75 and 0.9. This provides a corresponding reference for the test, a model for calculating fatigue life correctly on the mesoscale.

Potential Applications and Developments
The mesoscale model proposed in this study can accurately represent the fatigue life of rubber concrete under three-point bending fatigue loading. In addition, analysis of the static pressure load's damage form and the RC's ultimate load reveals that the model is also accurately represented. The study of RC is equally informative for ordinary concrete

Potential Applications and Developments
The mesoscale model proposed in this study can accurately represent the fatigue life of rubber concrete under three-point bending fatigue loading. In addition, analysis of the static pressure load's damage form and the RC's ultimate load reveals that the model is also accurately represented. The study of RC is equally informative for ordinary concrete and other polymer admixture concrete. The difference lies in the polymer's shape, size, location, and properties. Shape, size, and position are solved by Python code, and performance could be solved by setting properties on the polymer. However, it is often difficult to achieve the desired effect, and an interface layer is needed to change the mechanical relationship between the different substances. This study could also be applied to reinforced RC to simulate the location of damage and structural life of specific damaged structures in macrostructures to provide an initial structural performance judgment for actual structures.
The limitations of this research method lie in the 2D structure. When considering the mesoscale design in the 3D structure, the lack of computer performance is challenging to resolve, and the vast number of calculations leads to increased calculation time. Future work could be improved to develop a 3D mesoscale concrete model to calculate fatigue life, achieving the desired accuracy and computational efficiency requirements.

Conclusions
This paper proposed an RAM on the scope of a mesoscale study. The model used plastic damage theory and the insertion of cohesive elements in the ITZ, and is a new numerical model. This paper verifies the model's correctness in peak load and fatigue life. Peak loads were verified for five doping levels of 0, 2.5%, 5%, 7.5%, and 10%, and fatigue life was verified for stress levels of 0.75 and 0.85. After this, the results for stress levels of 0.8 and 0.9 were simulated and analyzed. The peak static pressure load in three-point bending was successfully modeled on a mesoscale as decreasing with increasing rubber doping, and the resulting deflection increased with increasing rubber doping. Static pressure and fatigue forms of damage could be observed in the mesoscale, where the point of damage produced by RC damage is not unique and increases with the amount of rubber admixture. The damage element produced by RC damage shows an order-disorder-order process as the rubber dosing increases. It was observed from the model that when the damage occurred to the RC, the internal rubber took the load. In addition, the model could simulate the three-point bending fatigue life at different stress levels for various rubber doping on a mesoscale. A quadratic function relating stress levels with different rubber doping to fatigue life was fitted, which can predict fatigue life for stress levels between 0.75 and 0.9, providing some reference value for the test. The mesoscale model in this paper satisfied the fatigue life simulation requirements perfectly. The method, with model improvements, could be applied to all RC fatigue structures in the future.

Institutional Review Board Statement: Not applicable.
Data Availability Statement: The data used in the article can be obtained from the author here.