# Viscoelastic-Thermodynamic Model and Crack Propagation in Asphalt Concrete Pavement

Authors: Ruijia Yuan; Liqin Zhao
Journal: Security and Communication Networks
 (2022)

---

## Abstract

Asphalt concrete pavement has a wide range of applications in the construction of high-grade highways. Asphalt concrete pavement has elastic properties, and its viscoelasticity and viscoplasticity are obviously affected by temperature. To determine whether a structure can continue to be used safely, it is most important to determine whether microscopic or macroscopic cracks present in the structure continue to propagate and cause structural failure. The purpose of this study was to study the effects of different ambient temperatures and mobile vehicles on the stress around the crack tip and the propagation path of different types of asphalt pavements. Based on the theory of fracture mechanics and viscoelasticity, this work studies the problem of crack propagation from a macroscopic perspective. It proposes to use the finite element method to analyze the crack problem in the pavement structure. With the help of ABAQUS software, a curved beam model with cracks was established for calculation, and moving loads and ambient temperature fields were added for analysis. The experimental results in this study show that with the decrease in the ambient temperature, the stress at the crack tip increases continuously. As the analysis time increases, the stress also increases. When the crack is smaller, the stress value at the tip is larger. Whena = 15 mm, the stress value at the crack tip is the largest. However, when a = 20 mm and a = 25 mm, the stress at the crack tip is significantly reduced. It shows that the crack is easier to crack when the crack value is smaller. After the crack has grown to a certain length, the growth rate will slow down.


---

## Body

## 1. Introduction

Asphalt pavement is favored for its excellent driving performance and has become the first choice for road construction in various countries. However, asphalt pavement also has a serious problem of early failure. That is, before reaching the design life, its good driving performance will be lost due to reflection cracks, temperature cracks, rutting, peeling, oil flooding, water damage, and other reasons. Among them, cracking and rutting are the most common and serious. The existence of these forms of damage has caused great pressure on pavement maintenance work.The high modulus asphalt concrete (HMAC) base has very good rutting and fatigue resistance, but they increase the risk of low-temperature cracking compared with conventional asphalt concrete (AC). During fracture testing, especially at lower test temperatures, it has become very common to quantify low-temperature cracks. Because sealants are exposed to high pressure for long periods of time, their viscoelastic properties at low temperatures play a crucial role in preventing bond failure.The main task of pavement design is to ensure that unacceptable damage does not occur during its life. This is a common goal of different design approaches. Experimentation is an effective way to screen and evaluate different design methods and is widely used by researchers. However, due to the limitations of the test itself (time-consuming, labor-intensive, and little data obtained), it is impossible for researchers to conduct similar experiments on a large scale. The methods of mechanical analysis and finite element simulation are relatively simple and effective analysis tools. It plays an important role in guiding the experimental design, epitaxial test results, and analyzing the mechanism of reflection cracks, the expansion process, and the selection of anticountermeasures. Therefore, it is very necessary to choose a suitable analysis method to analyze the cracking process of reflection cracks in the asphalt surface.Accurate quantitative information on crack length and width is important for assessing crack severity and making accurate pavement maintenance decisions. However, due to noise and fluctuations at the crack edges of the pavement image, it is difficult to precisely determine the edge correspondence of cracks necessary to measure the crack width at a particular pixel because different automated or manually performed measurements may yield different results on the same pixel.With the development of industry and technology, the previous methods can no longer meet the current environment. The innovation of this study is to use the finite element method to analyze the characteristic constitutive model and the theoretical basis of fracture mechanics based on the thermal viscoelastic-thermodynamic behavior of asphalt pavement materials. It was studied with the help of the finite element commercial simulation software ABAQUS, and a model with reflection cracks was used for numerical simulation.


## 2. Related Work

With the development of science and technology, many scholars have carried out related research on it. Among them, Bodin et al. proposed a nonlocal damage model to predict the fatigue crack behavior of pavement [1]. Gichul proposed a method to determine the dynamic modulus and phase angle of hot mix asphalt (HMA). This method is a simpler, nondestructive impact resonance (IR) test. It tested 35 HMA blends with different aggregate and binder properties using IR and traditional dynamic modulus testing. The comparison of the dynamic moduli and phase angles obtained from the IR test with those obtained from the conventional dynamic modulus test shows the reliability of the proposed procedure. The main advantage of these combined test methods is that they are simpler and require less expensive test equipment and setup [2]. Rys et al. presented a comparison of these two road materials in terms of low-temperature cracking, using a statistical approach based on an ordered logistic regression model. The analysis is based on the results of a field survey. The survey was carried out on 80 selected road sections that are normally used in Poland. The strength of the low-temperature crack is an analytical parameter. The results show that the pavement with high modulus asphalt foundation is 2.45 times more likely to be in the cracked pavement group than the pavement with the conventional asphalt concrete foundation [3]. Vega-Zamanillo et al. presented the performance of different types of asphalt mixtures containing recycled asphalt pavement (RAP) and industrial by-products (such as cupola slag and green sand). In this study, conventional asphalt was used as the binding material, while limestone was used as the natural aggregate. The Marshall method was used to design asphalt mixtures and void percentages. Both water sensitivity and wheel tracking tests were performed to evaluate the performance of the mixture. Although the only natural aggregate is limestone, this greatly hinders the performance of these mixtures. Laboratory results show that they can be successfully used in the construction of low-traffic intensity roads as bonding layers [4].However, most of the research methods of the above scholars are based on theoretical aspects. Practical applications are rare, and previous analytical methods of testing have often used a single number (fracture energy or fracture toughness) to quantify crack resistance. These tests cannot capture both the initiation and propagation of cracks. Resistance curves or R curves are widely used in many fields such as metals, polymers, and composites. The R curve treats crack resistance as a function of crack growth, including growth and propagation. Therefore, Yang and Braham tested three asphalt concrete mixtures, including hot mix, hot mix with recycled asphalt pavement, and hot mix with RA at two temperatures, three aging levels, and two humidity conditions through semicircular bending fracture. It turns out that the R curve can be a useful tool for quantifying the crack resistance of asphalt concrete in crack initiation and propagation [5]. Li et al. selected two types of hot cast sealants and two types of silicone sealants for extension tests with reference to standard adhesion tests. The results show that hot poured sealants may be more suitable for use in areas with long-term low temperature in winter [6]. Instead of measuring crack width on a pixel-by-pixel basis, Weng et al. provided an accurate and reliable segmentation-based method to measure crack width. Extensive experiments on synthetic and real pavement images show that the method can accurately and reliably quantify various cracks with an average accuracy of 93.7% in crack width [7]. These studies have deepened the analysis of crack propagation of asphalt concrete pavement to a certain extent, and combined with the actual calculation method, the general law of crack propagation can be obtained, which plays a role in road safety prevention. However, these calculation methods are complicated or the predictions are not accurate enough. The viscoelastic-thermodynamic model and crack propagation analysis method proposed in this study can be more accurate and efficient.


## 3. Theoretical Basis of the Viscoelastic-Thermodynamic Model

### 3.1. Fracture Mechanics Theory

Fracture mechanics is a new discipline developed over the past 30 years. It is used to study the equilibrium, propagation, instability, and strength of cracks under various environmental conditions (such as load effects, temperature changes, and humidity changes) and to predict structural stability [8, 9]. Its main task is to determine the stress intensity factor (linear elastic fracture mechanics) or the integral and the maximum displacement of the crack surface opening (elastoplastic fracture mechanics) and then determine the stress-strain field at the crack tip. If the required surface energy per unit area of the new surface is Γ, then once the crack covers the unit area, the required surface energy is 2Γ due to the formation of the top and bottom surfaces. Next, some materials undergo a certain amount of plastic deformation before breaking. This requires a certain amount of plastic deformation work [10]. If the plastic deformation work consumed when the crack extends to a unit area is UP, then the total energy required for the crack to extend to a unit area is R:(1)R=2Γ+UP.Assuming that the testing machine and the specimen together form a system, the energy that the system can provide when the crack propagates per unit area. Then, the conditions for crack propagation should be as follows:(2)G≥R.Letu1 and u2 be the potential energy before and after the crack propagation system, and then, the change in potential energy is as follows:(3)−Δu=u2−u1.Assuming that the crack in a crack body expands the areaΔA, the energy consumed by the crack expansion is as follows:(4)GI=ΔuΔA=−∂u∂A.From the above formula, the critical condition forI-type crack propagation is as follows:(5)GI=2Γ+UP.Among them, the subscriptI of GI represents a I-type crack. Under critical conditions, if 2Γ+UP is represented by GIC, the fracture criterion can be written uniformly as follows:(6)GI=GIC.Among them,GIC represents fracture toughness.If an infinite plate contains a central through crack, as shown in Figure1, the plate thickness is B, and the crack length is 2a, then the crack propagation area is A = 2Ba. There is a uniform tensile stress σ perpendicular to the crack acting at infinity. Substituting this formula into formula (4) has the following:(7)GI=−∂u∂A=−12B∂u∂a.Figure 1
Load distribution.From the knowledge of elasticity theory, it can be known that(8)u=πσ2a2EB.Therefore,(9)GI=πσ2aE.If plastic deformation is not considered, there are as follows:(10)GI=πσ2aE=2Γ.Therefore, when the plate is subjected to forceσ, the maximum allowable half-crack length is as follows:(11)ac=2EΓπσ2.When the plate has a central crack of length2a, the maximum tensile stress it can withstand is as follows:(12)ac=2EΓπσ.Figure2 shows a plane crack. The coordinate origin O is selected at the crack tip. r and θ are polar coordinates, and x and y are Cartesian coordinates.Figure 2
Schematic diagram of the crack tip.Then, the asymptotic stress field and displacement field of the crack tip in polar coordinates and Cartesian coordinates are as follows:(13)σrr=KΙ2πr54cosθ2−14cos3θ2+K∐2πr−54sinθ2+34sin3θ2.(14)σθθ=KΙ2πr34cosθ2+14cos3θ2+K∐2πr−34sinθ2−34sin3θ2.(15)σrθ=KΙ2πr14sinθ2+14sin3θ2+K∐2πr14cosθ2+34cos3θ2.(16)σx=KΙ2πrcosθ21−sinθ2sin3θ2−K∐2πrsinθ22+cosθ2cos3θ2.(17)σy=KΙ2πrcosθ21+sinθ2sin3θ2+K∐2πrsinθ2cosθ2cos3θ2.(18)τxy=KΙ2πrcosθ2sinθ2cos3θ2+K∐2πrcosθ21−sinθ2sin3θ2.(19)u=KΙ4Gr2π2χ−1cosθ2−cos3θ2+K∐4Gr2π2χ+3sinθ2+sin3θ2.(20)v=KΙ4Gr2π2χ−1sinθ2−sin3θ2−K∐4Gr2π2χ−3cosθ2+cos3θ2.In the formula,G—shear modulus, G=E/21+μ, χ=3−μ/4+μ (plane stress), and χ=3−4μ (plane strain).


### 3.2. Viscoelasticity Theory

Viscoelasticity means that asphalt has both solid properties, such as elasticity, strength, and thus stability, and liquid properties, such as flow properties that vary with time, temperature, and load size and speed [11, 12]. It is manifested as mechanical relaxation, that is, the change in mechanical properties over time [13, 14].In ABAQUS, the time-incremental finite element model is used for finite element calculation. Viscoelastic three-dimensional constitutive relations include differential constitutive equations and integral constitutive equations. In ABAQUS, the viscoelastic effects of the stress spherical tensor and the stress deviatoric tensor are considered, respectively. The viscoelastic integral formula can be expressed as follows:(21)σt=∫0t2Gτ−τ′e˙dt′+I∫0tKτ−τ′ϕ·dt′.In the formula,e and ϕ are the deviatoric and bulk strain tensors, respectively, K is the bulk modulus, and G is the shear modulus. Both K and G are functions of relaxation time τ, and t′ represents a different sign. I is the unit tensor.The relationship between relaxation timeτ and actual time t can be expressed by the following formula:(22)τ=∫0tdt′Aθθt′,dτdt=1Aθθt.In the formula,θ represents the temperature, and Aθ represents the displacement factor. So, if Aθ=1, then τ=t. Asphalt mixtures have thermorheologically simple material properties over a wide temperature range, satisfying the Williams– Landel–Ferry (WLF) formula. ABAQUS can convert any temperature to a reference temperature using the WLF formula.(23)−lgAθ=hT=C1θ−θ0C2+θ−θ0,where θ0 is the reference temperature, and C1 and C2 are the standard constants at the reference temperature.The bulk relaxation modulusKt and shear relaxation modulus Gt are expressed in the basic form of the Prony series as follows:(24)Gt=G0α∞G+∑i=1nGαiGexp−tτiG.(25)Kt=K0α∞K+∑i=1nKαiKexp−tτiK.In the formula, shear modulusGi, G∞: G∞=G0α∞G and Gi=G0αiG, bulk modulus Ki, K∞: K∞=K0α∞K and Ki=K0αiK, relative modulus αiG, αiK: αiG=Gi/G0 and αiK=Ki/K0, and G0 and K0 are the transient moduli of the viscoelastic material, respectively, and are defined as follows:(26)G0=Gt=0=G∞+∑i=1nGGi.(27)K0=Kt=0=K∞+∑i=1nKKi,where K∞ and G∞ are the long-term bulk modulus and long-term shear modulus. Delay times τiK and τiG are independent of each other and can be different quantities. τi=τiK=τiG is assumed in ABAQUS. On the other hand, the number of terms nK of the Prony series of bulk modulus (the shear modulus G(t) and bulk modulus K() of a viscoelastic material are expressed as time-variable Prony series. It is a generalized Maxwell model that can be used to describe relaxation effects through Prony series), and the number of terms nG of the Prony series of shear modulus may be different. In many practical cases, if the material can be deformed in other ways, there will never be a change in volume, so nK=0 can be assumed. The following focuses on the shear corresponding behavior and can be derived similarly for the iterative derivation of bulk modulus.The deviatoric stress integral formula is as follows:(28)S=∫0t2G∞+∑i=1nGGiexpτ′−ττie·dt′=∫0τ2G∞+∑i=1nGGiexpτ′−ττidedτ′dτ′.The above formula can be written as follows:(29)S=2G0e−∑i=1nαiei.In the formulaG0=G∞+∑i=1nGi, G0 are called transient shear moduli; αi=Gi/Go, αi are the quantities related to the ith Prony series, and(30)ei=∫0τ1−expτ′−ττidedτ′dτ′,where ei is the viscous strain (creep) of the ith Prony series. In the finite element analysis, this formula is rewritten as a time-incremental finite element equation. It is assumed that the strain varies linearly with time within a time increment step, that is, de/dτ′=Δe/Δτ.The above formula is decomposed into two parts:(31)ein+1=∫0τn1−expτ′−τn+1τidedτ′dτ′+∫τnτn+11−expτ′−τn+1τidedτ′dτ′.Now, the following is observed:(32)1−expτ′−τn+1τi=1−exp−Δττi+exp−Δττi1−expτ′−τn+1τi,(31) is brought to get(33)ein+1=1−expΔττi∫0τndedτ′dτ′+expΔττi∫0τn1−expτ′−τnτidedτ′dτ′+ΔeΔτ∫τnτn+11−expτ′−τn+1τidτ′.The first and last integrals in the formula can be evaluated quickly. From the following formula, the second integral represents the viscoelastic strain during the onset of the increment, so the change in viscous strain is as follows:(34)Δei=1−exp−Δττienexp−Δττi−1ein+Δτ−τi1−exp−ΔττiΔeΔτ=τiΔτΔττi+exp−Δττi−1Δe+1−exp−Δττien−ein.IfΔτ/τi⟶0 (i.e., Δτ/τi<10−7), the above formula can be simplified to(35)Δei=Δττi12Δe+en−ein.Therefore, the new viscous strain can be calculated by the above two formulas in one-time increment step. The new deviatoric stress increment is then obtained by (2)–(29).The tangent modulus can be directly obtained from these differential deviatoric stress increments, namely,(36)ΔS=2G0Δe−∑i=1nGαiein+1−ein.The deviatoric stress incrementΔS is derived from the deviatoric strain increment Δe to obtain the shear relaxation modulus. The equation is linear. The value of the modulus depends on the relaxation time.(37)GT=∂ΔS∂Δe=G01−∑i=1nαiτiΔτΔττi+exp−Δττi,Δττi>10−7G01−∑i=1n12αiΔττiΔττi<10−7.Similarly, the bulk modulus can be derived as follows:(38)KT=∂Δσ∂Δε=K01−∑i=1nαiτiΔτΔττi+exp−Δττi−1,Δττi>10−7K01−∑i=1n12αiΔττiΔττi<10−7.Finally, it is also necessary to establish a relationship between the relaxation time increment△τ and the real-time increment △t. Aθ is the nonlinear equation for temperature, and the direct approximation of Aθ was found to have large errors. However, hθ is a smooth and smooth function of temperature change, a nearly linear temperature increment function. Assuming further that temperature is a linear function of time within a time increment step, the following relationship is obtained:(39)hθ=−lnAθθt=a+bt,or(40)Aθ−1θt=expa+bt.Among them,(41)a=1△ttn+1hθn−tnhθn+1,b=1△thθn+1−hθn.Then, the relationship between△τ and △t is as follows:(42)△τ=∫tntn+1expa+btdt=1bexpa+btn+1−expa+btn.It can also be written as follows:(43)△τ=Aθ−1θn+1−Aθ−1θnhθn+1−hθn△t.Elastic modulusETt,t, shear modulus GTt,t, bulk modulus KTt,t, and Poisson’s ratio μTt,t are all related to temperature and time changes. Since the temperature change in μTt,t is small with time, it is assumed that it is a constant value, which is similar to the elastic constitutive, and has the following equation:(44)KTt,t=GTt,t21−μTt,t31−2μTt,t.Relaxation test—generalized Maxwell to ABAQUS-Prony series. First, the txt file of the stress relaxation test results is imported into origin for processing and fit into a generalized Maxwell model, namely, the relaxation master curve familyETt,t, and then, GTt,t and KTt,t are obtained.(45)GTt,t=ETt,t21+μ,KTt,t=ETt,t31−2μ.


### 3.3. Moving Loads

From the actual situation and previous studies, it can be analyzed that the pressure and position of the vehicles on the road are constantly changing due to the different sizes and types of vehicles [15, 16]. Therefore, it is necessary to simplify the vehicle into four forms: wave load, impact load, horizontal movement vertical load, and random load according to the usual method [17, 18]. The random load describes the driving trajectory more objectively and is close to reality. However, when the load moves to describe the force of the road surface, it is difficult to accurately describe the random load in a mathematical way. The moving load form can reflect the motion nature of the vehicle load. It mainly uses the equivalent uniform load. That is, the maximum pressure and amplitude of the load on the road surface are unchanged during the movement process. At present, when using the finite element method to solve the response of the pavement under the action of the dynamic load, the form of the moving load is usually used. When the vehicle is driving on a relatively smooth road, the vibration impact effect of the vehicle on the road caused by the unevenness of the road can be ignored. When the wheels are in contact with the road surface, the wheels exerted on the road surface are complex loads that are perpendicular to the road surface and parallel to the road surface [19, 20]. However, for the load parallel to the road surface, because of its rapid attenuation, the damage to the road surface is small and can be ignored. For simplicity, only moving loads in the direction perpendicular to the road surface are considered. When a vehicle passes a certain part of the road surface, deflection occurs here under the load perpendicular to the road surface. If the load is too large, the short-term deflection of the pavement will be too large, resulting in fracture [21–23].The distribution of the moving load is shown in Figure3. The load F moves at a constant speed c, and its expression is as follows:(46)qx,y,t=Qda−d120else≤x−ct≤da+d12∥dw−d22≤y≤dw+d22.Figure 3
Load distribution.In the formula,c is the moving speed of the load, d1 is the length of the load, d2 is the width of the load, da is the distance between the center of the area before and after the load, and dw is the distance between the center of the left and right areas of the load.When analyzing the vehicle load on the road surface, in addition to analyzing the form of load application, the way the wheel contacts the road surface should also be analyzed. Vehicle loads transmit forces through the interaction between the tires and the road surface. The size of the vehicle-mounted pressure on the road surface is directly affected by the stiffness of the tire itself, the form of the tire acting on the road surface, and the size of the vehicle pressure on the tire. When the tires act on the road surface, the shape of the mutual contact surface can be regarded as an ellipse, but for heavy-duty vehicles, the shape of the contact surface can be regarded as a rectangle, and the calculation of the elliptical contact shape is troublesome. To facilitate the simulation calculation, the contact shape of the vehicle wheel load involved in the simulation analysis in this study adopts an equivalent rectangle. That is, the vertical uniform rectangular pressure distribution load is used to simulate the load applied by the tire to the road surface [24, 25].For one of the single wheel loading, according to the requirements of the “Code for Design of Highway Cement Concrete Pavement” (it has played an important role in guiding the design of cement concrete pavement in China’s highways and ensuring the quality of the pavement) [26, 27], the contact surface of the tire when it acts on the road surface can be regarded as a combination of a rectangle of 0.6 L × 0.4 L and two semicircles with a diameter of 0.6 L. L is the length of vehicle load action along the direction of vehicle travel. The contact surface between the single wheel load and the road surface can be equivalent to a rectangle with a length of 0.87 L and a width of 0.6 L. The equivalent conversion of the contact shape of the single wheel load is shown in Figure 4.Figure 4
Equivalent area of single wheel load contact.The resulting equivalent contact area is as follows:(47)A=0.4L×0.6L+π×0.3L2,F=A×P.In the formula,A is the contact area between each tire and the road surface, F is the load applied by each tire to the road surface, and P is the ground contact pressure of the tire.


## 3.1. Fracture Mechanics Theory

Fracture mechanics is a new discipline developed over the past 30 years. It is used to study the equilibrium, propagation, instability, and strength of cracks under various environmental conditions (such as load effects, temperature changes, and humidity changes) and to predict structural stability [8, 9]. Its main task is to determine the stress intensity factor (linear elastic fracture mechanics) or the integral and the maximum displacement of the crack surface opening (elastoplastic fracture mechanics) and then determine the stress-strain field at the crack tip. If the required surface energy per unit area of the new surface is Γ, then once the crack covers the unit area, the required surface energy is 2Γ due to the formation of the top and bottom surfaces. Next, some materials undergo a certain amount of plastic deformation before breaking. This requires a certain amount of plastic deformation work [10]. If the plastic deformation work consumed when the crack extends to a unit area is UP, then the total energy required for the crack to extend to a unit area is R:(1)R=2Γ+UP.Assuming that the testing machine and the specimen together form a system, the energy that the system can provide when the crack propagates per unit area. Then, the conditions for crack propagation should be as follows:(2)G≥R.Letu1 and u2 be the potential energy before and after the crack propagation system, and then, the change in potential energy is as follows:(3)−Δu=u2−u1.Assuming that the crack in a crack body expands the areaΔA, the energy consumed by the crack expansion is as follows:(4)GI=ΔuΔA=−∂u∂A.From the above formula, the critical condition forI-type crack propagation is as follows:(5)GI=2Γ+UP.Among them, the subscriptI of GI represents a I-type crack. Under critical conditions, if 2Γ+UP is represented by GIC, the fracture criterion can be written uniformly as follows:(6)GI=GIC.Among them,GIC represents fracture toughness.If an infinite plate contains a central through crack, as shown in Figure1, the plate thickness is B, and the crack length is 2a, then the crack propagation area is A = 2Ba. There is a uniform tensile stress σ perpendicular to the crack acting at infinity. Substituting this formula into formula (4) has the following:(7)GI=−∂u∂A=−12B∂u∂a.Figure 1
Load distribution.From the knowledge of elasticity theory, it can be known that(8)u=πσ2a2EB.Therefore,(9)GI=πσ2aE.If plastic deformation is not considered, there are as follows:(10)GI=πσ2aE=2Γ.Therefore, when the plate is subjected to forceσ, the maximum allowable half-crack length is as follows:(11)ac=2EΓπσ2.When the plate has a central crack of length2a, the maximum tensile stress it can withstand is as follows:(12)ac=2EΓπσ.Figure2 shows a plane crack. The coordinate origin O is selected at the crack tip. r and θ are polar coordinates, and x and y are Cartesian coordinates.Figure 2
Schematic diagram of the crack tip.Then, the asymptotic stress field and displacement field of the crack tip in polar coordinates and Cartesian coordinates are as follows:(13)σrr=KΙ2πr54cosθ2−14cos3θ2+K∐2πr−54sinθ2+34sin3θ2.(14)σθθ=KΙ2πr34cosθ2+14cos3θ2+K∐2πr−34sinθ2−34sin3θ2.(15)σrθ=KΙ2πr14sinθ2+14sin3θ2+K∐2πr14cosθ2+34cos3θ2.(16)σx=KΙ2πrcosθ21−sinθ2sin3θ2−K∐2πrsinθ22+cosθ2cos3θ2.(17)σy=KΙ2πrcosθ21+sinθ2sin3θ2+K∐2πrsinθ2cosθ2cos3θ2.(18)τxy=KΙ2πrcosθ2sinθ2cos3θ2+K∐2πrcosθ21−sinθ2sin3θ2.(19)u=KΙ4Gr2π2χ−1cosθ2−cos3θ2+K∐4Gr2π2χ+3sinθ2+sin3θ2.(20)v=KΙ4Gr2π2χ−1sinθ2−sin3θ2−K∐4Gr2π2χ−3cosθ2+cos3θ2.In the formula,G—shear modulus, G=E/21+μ, χ=3−μ/4+μ (plane stress), and χ=3−4μ (plane strain).


## 3.2. Viscoelasticity Theory

Viscoelasticity means that asphalt has both solid properties, such as elasticity, strength, and thus stability, and liquid properties, such as flow properties that vary with time, temperature, and load size and speed [11, 12]. It is manifested as mechanical relaxation, that is, the change in mechanical properties over time [13, 14].In ABAQUS, the time-incremental finite element model is used for finite element calculation. Viscoelastic three-dimensional constitutive relations include differential constitutive equations and integral constitutive equations. In ABAQUS, the viscoelastic effects of the stress spherical tensor and the stress deviatoric tensor are considered, respectively. The viscoelastic integral formula can be expressed as follows:(21)σt=∫0t2Gτ−τ′e˙dt′+I∫0tKτ−τ′ϕ·dt′.In the formula,e and ϕ are the deviatoric and bulk strain tensors, respectively, K is the bulk modulus, and G is the shear modulus. Both K and G are functions of relaxation time τ, and t′ represents a different sign. I is the unit tensor.The relationship between relaxation timeτ and actual time t can be expressed by the following formula:(22)τ=∫0tdt′Aθθt′,dτdt=1Aθθt.In the formula,θ represents the temperature, and Aθ represents the displacement factor. So, if Aθ=1, then τ=t. Asphalt mixtures have thermorheologically simple material properties over a wide temperature range, satisfying the Williams– Landel–Ferry (WLF) formula. ABAQUS can convert any temperature to a reference temperature using the WLF formula.(23)−lgAθ=hT=C1θ−θ0C2+θ−θ0,where θ0 is the reference temperature, and C1 and C2 are the standard constants at the reference temperature.The bulk relaxation modulusKt and shear relaxation modulus Gt are expressed in the basic form of the Prony series as follows:(24)Gt=G0α∞G+∑i=1nGαiGexp−tτiG.(25)Kt=K0α∞K+∑i=1nKαiKexp−tτiK.In the formula, shear modulusGi, G∞: G∞=G0α∞G and Gi=G0αiG, bulk modulus Ki, K∞: K∞=K0α∞K and Ki=K0αiK, relative modulus αiG, αiK: αiG=Gi/G0 and αiK=Ki/K0, and G0 and K0 are the transient moduli of the viscoelastic material, respectively, and are defined as follows:(26)G0=Gt=0=G∞+∑i=1nGGi.(27)K0=Kt=0=K∞+∑i=1nKKi,where K∞ and G∞ are the long-term bulk modulus and long-term shear modulus. Delay times τiK and τiG are independent of each other and can be different quantities. τi=τiK=τiG is assumed in ABAQUS. On the other hand, the number of terms nK of the Prony series of bulk modulus (the shear modulus G(t) and bulk modulus K() of a viscoelastic material are expressed as time-variable Prony series. It is a generalized Maxwell model that can be used to describe relaxation effects through Prony series), and the number of terms nG of the Prony series of shear modulus may be different. In many practical cases, if the material can be deformed in other ways, there will never be a change in volume, so nK=0 can be assumed. The following focuses on the shear corresponding behavior and can be derived similarly for the iterative derivation of bulk modulus.The deviatoric stress integral formula is as follows:(28)S=∫0t2G∞+∑i=1nGGiexpτ′−ττie·dt′=∫0τ2G∞+∑i=1nGGiexpτ′−ττidedτ′dτ′.The above formula can be written as follows:(29)S=2G0e−∑i=1nαiei.In the formulaG0=G∞+∑i=1nGi, G0 are called transient shear moduli; αi=Gi/Go, αi are the quantities related to the ith Prony series, and(30)ei=∫0τ1−expτ′−ττidedτ′dτ′,where ei is the viscous strain (creep) of the ith Prony series. In the finite element analysis, this formula is rewritten as a time-incremental finite element equation. It is assumed that the strain varies linearly with time within a time increment step, that is, de/dτ′=Δe/Δτ.The above formula is decomposed into two parts:(31)ein+1=∫0τn1−expτ′−τn+1τidedτ′dτ′+∫τnτn+11−expτ′−τn+1τidedτ′dτ′.Now, the following is observed:(32)1−expτ′−τn+1τi=1−exp−Δττi+exp−Δττi1−expτ′−τn+1τi,(31) is brought to get(33)ein+1=1−expΔττi∫0τndedτ′dτ′+expΔττi∫0τn1−expτ′−τnτidedτ′dτ′+ΔeΔτ∫τnτn+11−expτ′−τn+1τidτ′.The first and last integrals in the formula can be evaluated quickly. From the following formula, the second integral represents the viscoelastic strain during the onset of the increment, so the change in viscous strain is as follows:(34)Δei=1−exp−Δττienexp−Δττi−1ein+Δτ−τi1−exp−ΔττiΔeΔτ=τiΔτΔττi+exp−Δττi−1Δe+1−exp−Δττien−ein.IfΔτ/τi⟶0 (i.e., Δτ/τi<10−7), the above formula can be simplified to(35)Δei=Δττi12Δe+en−ein.Therefore, the new viscous strain can be calculated by the above two formulas in one-time increment step. The new deviatoric stress increment is then obtained by (2)–(29).The tangent modulus can be directly obtained from these differential deviatoric stress increments, namely,(36)ΔS=2G0Δe−∑i=1nGαiein+1−ein.The deviatoric stress incrementΔS is derived from the deviatoric strain increment Δe to obtain the shear relaxation modulus. The equation is linear. The value of the modulus depends on the relaxation time.(37)GT=∂ΔS∂Δe=G01−∑i=1nαiτiΔτΔττi+exp−Δττi,Δττi>10−7G01−∑i=1n12αiΔττiΔττi<10−7.Similarly, the bulk modulus can be derived as follows:(38)KT=∂Δσ∂Δε=K01−∑i=1nαiτiΔτΔττi+exp−Δττi−1,Δττi>10−7K01−∑i=1n12αiΔττiΔττi<10−7.Finally, it is also necessary to establish a relationship between the relaxation time increment△τ and the real-time increment △t. Aθ is the nonlinear equation for temperature, and the direct approximation of Aθ was found to have large errors. However, hθ is a smooth and smooth function of temperature change, a nearly linear temperature increment function. Assuming further that temperature is a linear function of time within a time increment step, the following relationship is obtained:(39)hθ=−lnAθθt=a+bt,or(40)Aθ−1θt=expa+bt.Among them,(41)a=1△ttn+1hθn−tnhθn+1,b=1△thθn+1−hθn.Then, the relationship between△τ and △t is as follows:(42)△τ=∫tntn+1expa+btdt=1bexpa+btn+1−expa+btn.It can also be written as follows:(43)△τ=Aθ−1θn+1−Aθ−1θnhθn+1−hθn△t.Elastic modulusETt,t, shear modulus GTt,t, bulk modulus KTt,t, and Poisson’s ratio μTt,t are all related to temperature and time changes. Since the temperature change in μTt,t is small with time, it is assumed that it is a constant value, which is similar to the elastic constitutive, and has the following equation:(44)KTt,t=GTt,t21−μTt,t31−2μTt,t.Relaxation test—generalized Maxwell to ABAQUS-Prony series. First, the txt file of the stress relaxation test results is imported into origin for processing and fit into a generalized Maxwell model, namely, the relaxation master curve familyETt,t, and then, GTt,t and KTt,t are obtained.(45)GTt,t=ETt,t21+μ,KTt,t=ETt,t31−2μ.


## 3.3. Moving Loads

From the actual situation and previous studies, it can be analyzed that the pressure and position of the vehicles on the road are constantly changing due to the different sizes and types of vehicles [15, 16]. Therefore, it is necessary to simplify the vehicle into four forms: wave load, impact load, horizontal movement vertical load, and random load according to the usual method [17, 18]. The random load describes the driving trajectory more objectively and is close to reality. However, when the load moves to describe the force of the road surface, it is difficult to accurately describe the random load in a mathematical way. The moving load form can reflect the motion nature of the vehicle load. It mainly uses the equivalent uniform load. That is, the maximum pressure and amplitude of the load on the road surface are unchanged during the movement process. At present, when using the finite element method to solve the response of the pavement under the action of the dynamic load, the form of the moving load is usually used. When the vehicle is driving on a relatively smooth road, the vibration impact effect of the vehicle on the road caused by the unevenness of the road can be ignored. When the wheels are in contact with the road surface, the wheels exerted on the road surface are complex loads that are perpendicular to the road surface and parallel to the road surface [19, 20]. However, for the load parallel to the road surface, because of its rapid attenuation, the damage to the road surface is small and can be ignored. For simplicity, only moving loads in the direction perpendicular to the road surface are considered. When a vehicle passes a certain part of the road surface, deflection occurs here under the load perpendicular to the road surface. If the load is too large, the short-term deflection of the pavement will be too large, resulting in fracture [21–23].The distribution of the moving load is shown in Figure3. The load F moves at a constant speed c, and its expression is as follows:(46)qx,y,t=Qda−d120else≤x−ct≤da+d12∥dw−d22≤y≤dw+d22.Figure 3
Load distribution.In the formula,c is the moving speed of the load, d1 is the length of the load, d2 is the width of the load, da is the distance between the center of the area before and after the load, and dw is the distance between the center of the left and right areas of the load.When analyzing the vehicle load on the road surface, in addition to analyzing the form of load application, the way the wheel contacts the road surface should also be analyzed. Vehicle loads transmit forces through the interaction between the tires and the road surface. The size of the vehicle-mounted pressure on the road surface is directly affected by the stiffness of the tire itself, the form of the tire acting on the road surface, and the size of the vehicle pressure on the tire. When the tires act on the road surface, the shape of the mutual contact surface can be regarded as an ellipse, but for heavy-duty vehicles, the shape of the contact surface can be regarded as a rectangle, and the calculation of the elliptical contact shape is troublesome. To facilitate the simulation calculation, the contact shape of the vehicle wheel load involved in the simulation analysis in this study adopts an equivalent rectangle. That is, the vertical uniform rectangular pressure distribution load is used to simulate the load applied by the tire to the road surface [24, 25].For one of the single wheel loading, according to the requirements of the “Code for Design of Highway Cement Concrete Pavement” (it has played an important role in guiding the design of cement concrete pavement in China’s highways and ensuring the quality of the pavement) [26, 27], the contact surface of the tire when it acts on the road surface can be regarded as a combination of a rectangle of 0.6 L × 0.4 L and two semicircles with a diameter of 0.6 L. L is the length of vehicle load action along the direction of vehicle travel. The contact surface between the single wheel load and the road surface can be equivalent to a rectangle with a length of 0.87 L and a width of 0.6 L. The equivalent conversion of the contact shape of the single wheel load is shown in Figure 4.Figure 4
Equivalent area of single wheel load contact.The resulting equivalent contact area is as follows:(47)A=0.4L×0.6L+π×0.3L2,F=A×P.In the formula,A is the contact area between each tire and the road surface, F is the load applied by each tire to the road surface, and P is the ground contact pressure of the tire.


## 4. Simulation Calculation of Crack Propagation Based on Finite Element

### 4.1. Determination of Parameters and Methods

#### 4.1.1. Determination of the Action Mode and Size of the Moving Load

In the analysis of the model built in this study, the realization of moving loads in ABAQUS (ABAQUS is a set of powerful finite element software for engineering simulation. It solves problems ranging from relatively simple linear analysis to many complex nonlinear problems. ABAQUS includes a rich library of elements that can simulate arbitrary geometric shapes) is a very important step. The moving loads distributed as shown in Figure3 are applied using the custom load subroutine method. The DLOAD subroutine (the DLOAD subroutine is a subroutine in ABAQUS for defining complex loads that vary with time and space position or other variables, written in Fortran 77) is called in the implicit analysis. It compiles a program that conforms to the size and distribution of the actual load applied in this study on the basis of the existing fixed commands. The distributed loads used in this study are four rectangular surface loads distributed on the four corners of the rectangle, the distance between the front and rear load centers is 1 m, and the distance between the left and right load centers is 0.4 m. The area where the moving load passes on the model surface forms two moving load bands, as shown in Figure 5.Figure 5
Moving load band.It minimizes the size of the model while meeting the analysis requirements. So when writing subroutines, time (1)-L is used to ensure the coordination of time and loading area. time(1) represents the moving distance of the preceding load. L represents the distance between the front and rear load centers. In this way, the front load is applied to the model surface first. With the increase in time(1), the latter load also starts to be applied to the model surface after the former load moves by L. When the front load is out of the model, the latter load will continue to increase with time until it leaves the model. By this method, the established model can be reduced, thus saving computing resources and time.In this study, the single axle load is 28 kN, and the ground pressure of the tire is assumed to be 0.35 MPa. The front and rear axles all adopt a single wheel and single axle structure, assuming that the axle loads assigned to the front and rear axles are the same. The pressure surface load was selected as the load, and two moving tire load zones with a width of 10 cm were divided on the surface of the established asphalt pavement model. It applies moving loads to the model in the same manner as the moving load subroutine described earlier.


#### 4.1.2. Model Establishment and Finite Element Mesh Division

The first step in finite element simulation analysis is to build a model that is close to reality. This study uses thexyz coordinate system. The model size is 600 mm ∗ 100 mm ∗ 100 mm, and the asphalt concrete bending beam models with initial cracks of 5 mm, 10 mm, 15 mm, 20 mm, and 25 mm, respectively, are shown in Figure 6. This model can better simulate the failure process of pavement with initial cracks and has been widely used in many simulation analyses.Figure 6
Asphalt concrete pavement beam model with initial reflection cracks (a = 15 mm).A crucial part of finite element analysis in ABAQUS is the division of elements. The rationality of the unit division has a great influence on the simulation results. According to the analysis principle of finite element method mentioned above, the precision of numerical calculation result of finite element often depends on the division quality of element mesh. However, the smaller the element size is selected, the denser the grid, the larger the memory required for model analysis, and the higher the requirements for the computer. An increase in the number of cells taken also increases the calculated intake error of the response. Therefore, in the model calculation, the division technique, the size of the element, and the number of elements contained in the model should take into account the performance of the computer. Reasonable mesh division can improve the efficiency of model analysis and calculation. The quality of meshing will affect the accuracy of the results, and there are several control factors.The size of the pavement model determines the size of the mesh, the element shape is hexahedron, and the division technique is structural division. The configuration of the vehicle load, such as the wheel spacing, and the mesh is refined at the wheel. The processing of details also has a certain impact on the accuracy of the results. If the mesh is too small, the computation will be too large, and if the mesh is too coarse, the result will be inaccurate. This model is an eight-node linear hexahedron element, with reduced integration, hourglass control, and a total of 8400 elements, as shown in Figure7.Figure 7
Finite element meshing of an asphalt concrete beam with initial cracks (a = 15 mm).


#### 4.1.3. Determination of Boundary Conditions, Loads, and Material Properties

Boundary conditions also have a significant impact on the finite element analysis of pavement structural responses. According to different analysis needs, this study adopts the symmetric, antisymmetric, completely fixed boundary conditions and vehicle moving load of the initial analysis step. In the later analysis step, the displacement boundary load is set, and a set of reference points is established, so that the force of the structure can reach equilibrium. Elastic and viscoelastic materials that vary with temperature have different mechanical properties at different temperatures. It defines the time-domain viscoelastic material. Under the time list, the input method of viscoelastic material parameters is selected according to your needs. It enters the Prony series, selects Prony, and adds the corresponding shear relaxation modulus, bulk relaxation modulus, and relaxation time at different temperatures. At the same time, the corresponding temperature elastic modulus, Poisson’s ratio, and temperature values are added to the elastic material parameters.


### 4.2. Simulation of Crack Propagation in Different Environments

#### 4.2.1. Stress Analysis of Crack Tip Corresponding to Different Viscoelastic Material Parameters

A viscoelastic asphalt concrete beam with an initial reflection crack is subjected to a moving load at a certain speed. It sets the ambient temperature to 20°C and changes the material parameters of asphalt concrete, which not only changes the proportion of asphalt and aggregate, etc., thereby increasing or decreasing the viscosity. The viscoelasticity of material I-materialV gradually decreases, and the reduction range is 5%. The stress change around the crack tip of the asphalt concrete beam model is shown in Figure 8.Figure 8
Stress cloud diagram around the crack of viscoelastic pavement under moving load (a = 25 mm, V = 80 dm/s, T = 20°C).It can be observed from Figure8 that material I has strong viscoelasticity, and the stress value around the crack tip is relatively high. The stress value during the whole moving load application process is slightly higher than the other four materials. With the change in material, the maximum stress value is slightly lower than that of material I, and the difference in stress value during vehicle operation is relatively large. Therefore, when the asphalt material ratio is within a certain range, the stress value around the crack remains stable. When the viscoelasticity of the asphalt mixture is high, the stress value around the reflection crack is also high, and the crack has the risk of further expansion. As the viscoelasticity of the asphalt mixture decreases, even the stress values around the crack tip decrease. However, in the case of hot weather, it is not conducive to the stability of the actual road surface. These are all factors that should be considered in practical engineering.


#### 4.2.2. Stress Analysis around the Crack Tip under Moving Loads with Different Speeds

In the process of driving on the road, different road sections have different speed limits. Different types of vehicles also have their own different travel speeds. In the simulation analysis of road cracks, the speed of the vehicle is a factor that cannot be ignored. The application of the dynamic load is much closer to the actual load than the static load, which also makes the analysis results more realistic and accurate. Therefore, it is also necessary to explore the effect of different vehicle speeds on the cracks of asphalt pavement. Figure9 shows the stress cloud diagram of the beam model under the action of moving loads with different velocities and the comparison of the results.Figure 9
Stress cloud diagram around the crack of the pavement under moving load (a = 25 mm, V = 60 dm/s, T = 20°C).As shown in Figure10, moving loads at different speeds were applied and the effect of speed on the stress around the crack tip was observed. From the numerical results, it can be seen that under the environmental condition of 20°C, the stress values generated by the loads in different speed ranges are similar. That is, the change in velocity in one cycle has little effect on the stress at the crack tip.Figure 10
Stresses around cracks in viscoelastic pavement at different driving speeds (a = 25 mm, T = 20°C).


#### 4.2.3. Stress Analysis of Crack Tip When Different Initial Temperature Fields Are Applied

In the research process of asphalt concrete pavement, because the asphalt concrete material is greatly affected by temperature, the crack state, propagation phenomenon, and the stress around the tip of the asphalt pavement with initial crack will be different at different ambient temperatures. Figure11 shows the stress nephogram around the crack of viscoelastic asphalt pavement at 20°C. In recent years, there have been many studies on the effect of temperature on asphalt mixtures. Many research results also have great practical application value. In this study, ambient temperature fields with different temperature values are added, and displacement loads equivalent to the actual ones are applied to explore the mechanical properties of asphalt concrete pavements with initial reflection cracks of different sizes at different ambient temperatures. Tables 1–5 provide the simulation results of analysis and comparison at different temperatures with 5 mm, 10 mm, 15 mm, 20 mm, and 25 mm initial reflection cracks, respectively.Figure 11
Stress cloud diagram around cracks of viscoelastic asphalt pavement at 20°C (a = 15 mm).Table 1
Stress around cracks of viscoelastic asphalt pavement under different temperature conditions (a = 5 mm).
20°C30°C40°C50°C60°CTime (s)Misesmax (kPa)Time (s)Misesmax (kPa)Time (s)Misesmax (kPa)Time (s)Misesmax (kPa)Time (s)Misesmax (kPa)00000000000.16660.15000.14200.14000.13850.213320.210000.28400.28000.27700.319980.315000.312600.312000.311550.426640.420000.416800.416000.415400.533300.525000.521000.520000.519250.639960.630000.625200.624000.62310Table 2
Stresses around cracks of viscoelastic asphalt pavement under different temperature conditions (a = 10 mm).
20°C30°C40°C50°C60°CTime (s)Misesmax (kPa)Time (s)Misesmax (kPa)Time (s)Misesmax (kPa)Time (s)Misesmax (kPa)Time (s)Misesmax (kPa)00000000000.17660.15800.14800.14400.14030.214320.210800.29000.28400.27880.320980.315800.313200.312400.311730.427640.420800.417400.416400.415580.534300.525800.521600.520400.519430.640960.630800.625800.624400.62328Table 3
Stresses around cracks of viscoelastic asphalt pavement under different temperature conditions (a = 15 mm).
20°C30°C40°C50°C60°CTime (s)Misesmax (kPa)Time (s)Misesmax (kPa)Time (s)Misesmax (kPa)Time (s)Misesmax (kPa)Time (s)Misesmax (kPa)00000000000.112660.110800.19800.19400.19030.219320.215800.214000.213400.212880.325980.320800.318200.317400.316730.432640.425800.422400.421400.420580.539300.530800.526600.525400.524430.645960.635800.630800.629400.62828Table 4
Stresses around cracks of viscoelastic asphalt pavement under different temperature conditions (a = 20 mm).
20°C30°C40°C50°C60°CTime (s)Misesmax (kPa)Time (s)Misesmax (kPa)Time (s)Misesmax (kPa)Time (s)Misesmax (kPa)Time (s)Misesmax (kPa)00000000000.117160.115300.114300.113900.113530.223820.220300.218500.217900.217380.330480.325300.322700.321900.321230.437140.430300.426900.425900.425080.543800.535300.531100.529900.528930.650460.640300.635300.633900.63278Table 5
Stresses around cracks of viscoelastic pavement at different temperatures (a = 25 mm).
20°C30°C40°C50°C60°CTime (s)Misesmax (kPa)Time (s)Misesmax (kPa)Time (s)Misesmax (kPa)Time (s)Misesmax (kPa)Time (s)Misesmax (kPa)00000000000.127110.120300.118400.117700.116530.233770.225300.222600.221700.220380.340430.330300.326800.325700.324230.447090.435300.431000.429700.428080.553750.540300.535200.533700.531930.660410.645300.639400.637700.63578It can be seen from Tables1–5 that under different ambient temperature conditions, when other conditions are the same, the stress around the crack tip of different sizes is different. As the ambient temperature decreases, the stress at the crack tip increases, and as the analysis time increases, the stress also increases.It can be seen from Figure12 that the impact on the asphalt pavement with initial crack elastic and viscoelastic materials is more obvious under the low-temperature environment. That is to say, the damage mostly occurs when the ambient temperature is low. For example, in winter, the numerical analysis results are consistent with reality. It shows that the numerical simulation process of viscoelastic materials can more effectively simulate the crack propagation process and mechanical behavior of asphalt mixture trabeculae under different ambient temperatures.Figure 12
Stresses around cracks in viscoelastic asphalt pavement under different temperature conditions (v = 80 dm/s, T = 20°C).Figure13 is an analysis diagram of the tip stress results of numerical simulation of initial cracks of different sizes at the same temperature. It can be seen from the figure that when the crack size is different, the stress at its tip is also different. When the crack is smaller, the stress value at the tip is larger. When a = 15 mm, the stress value at the crack tip is the largest. However, when a = 20 mm and a = 25 mm, the stress at the crack tip is significantly reduced. It indicates that the crack may be easier to crack when the crack value is smaller. After the crack has grown to a certain length, the growth rate will slow down.Figure 13
Surrounding stress of viscoelastic pavement with initial cracks of different lengths (T = 20°C).


## 4.1. Determination of Parameters and Methods

### 4.1.1. Determination of the Action Mode and Size of the Moving Load

In the analysis of the model built in this study, the realization of moving loads in ABAQUS (ABAQUS is a set of powerful finite element software for engineering simulation. It solves problems ranging from relatively simple linear analysis to many complex nonlinear problems. ABAQUS includes a rich library of elements that can simulate arbitrary geometric shapes) is a very important step. The moving loads distributed as shown in Figure3 are applied using the custom load subroutine method. The DLOAD subroutine (the DLOAD subroutine is a subroutine in ABAQUS for defining complex loads that vary with time and space position or other variables, written in Fortran 77) is called in the implicit analysis. It compiles a program that conforms to the size and distribution of the actual load applied in this study on the basis of the existing fixed commands. The distributed loads used in this study are four rectangular surface loads distributed on the four corners of the rectangle, the distance between the front and rear load centers is 1 m, and the distance between the left and right load centers is 0.4 m. The area where the moving load passes on the model surface forms two moving load bands, as shown in Figure 5.Figure 5
Moving load band.It minimizes the size of the model while meeting the analysis requirements. So when writing subroutines, time (1)-L is used to ensure the coordination of time and loading area. time(1) represents the moving distance of the preceding load. L represents the distance between the front and rear load centers. In this way, the front load is applied to the model surface first. With the increase in time(1), the latter load also starts to be applied to the model surface after the former load moves by L. When the front load is out of the model, the latter load will continue to increase with time until it leaves the model. By this method, the established model can be reduced, thus saving computing resources and time.In this study, the single axle load is 28 kN, and the ground pressure of the tire is assumed to be 0.35 MPa. The front and rear axles all adopt a single wheel and single axle structure, assuming that the axle loads assigned to the front and rear axles are the same. The pressure surface load was selected as the load, and two moving tire load zones with a width of 10 cm were divided on the surface of the established asphalt pavement model. It applies moving loads to the model in the same manner as the moving load subroutine described earlier.


### 4.1.2. Model Establishment and Finite Element Mesh Division

The first step in finite element simulation analysis is to build a model that is close to reality. This study uses thexyz coordinate system. The model size is 600 mm ∗ 100 mm ∗ 100 mm, and the asphalt concrete bending beam models with initial cracks of 5 mm, 10 mm, 15 mm, 20 mm, and 25 mm, respectively, are shown in Figure 6. This model can better simulate the failure process of pavement with initial cracks and has been widely used in many simulation analyses.Figure 6
Asphalt concrete pavement beam model with initial reflection cracks (a = 15 mm).A crucial part of finite element analysis in ABAQUS is the division of elements. The rationality of the unit division has a great influence on the simulation results. According to the analysis principle of finite element method mentioned above, the precision of numerical calculation result of finite element often depends on the division quality of element mesh. However, the smaller the element size is selected, the denser the grid, the larger the memory required for model analysis, and the higher the requirements for the computer. An increase in the number of cells taken also increases the calculated intake error of the response. Therefore, in the model calculation, the division technique, the size of the element, and the number of elements contained in the model should take into account the performance of the computer. Reasonable mesh division can improve the efficiency of model analysis and calculation. The quality of meshing will affect the accuracy of the results, and there are several control factors.The size of the pavement model determines the size of the mesh, the element shape is hexahedron, and the division technique is structural division. The configuration of the vehicle load, such as the wheel spacing, and the mesh is refined at the wheel. The processing of details also has a certain impact on the accuracy of the results. If the mesh is too small, the computation will be too large, and if the mesh is too coarse, the result will be inaccurate. This model is an eight-node linear hexahedron element, with reduced integration, hourglass control, and a total of 8400 elements, as shown in Figure7.Figure 7
Finite element meshing of an asphalt concrete beam with initial cracks (a = 15 mm).


### 4.1.3. Determination of Boundary Conditions, Loads, and Material Properties

Boundary conditions also have a significant impact on the finite element analysis of pavement structural responses. According to different analysis needs, this study adopts the symmetric, antisymmetric, completely fixed boundary conditions and vehicle moving load of the initial analysis step. In the later analysis step, the displacement boundary load is set, and a set of reference points is established, so that the force of the structure can reach equilibrium. Elastic and viscoelastic materials that vary with temperature have different mechanical properties at different temperatures. It defines the time-domain viscoelastic material. Under the time list, the input method of viscoelastic material parameters is selected according to your needs. It enters the Prony series, selects Prony, and adds the corresponding shear relaxation modulus, bulk relaxation modulus, and relaxation time at different temperatures. At the same time, the corresponding temperature elastic modulus, Poisson’s ratio, and temperature values are added to the elastic material parameters.


## 4.1.1. Determination of the Action Mode and Size of the Moving Load

In the analysis of the model built in this study, the realization of moving loads in ABAQUS (ABAQUS is a set of powerful finite element software for engineering simulation. It solves problems ranging from relatively simple linear analysis to many complex nonlinear problems. ABAQUS includes a rich library of elements that can simulate arbitrary geometric shapes) is a very important step. The moving loads distributed as shown in Figure3 are applied using the custom load subroutine method. The DLOAD subroutine (the DLOAD subroutine is a subroutine in ABAQUS for defining complex loads that vary with time and space position or other variables, written in Fortran 77) is called in the implicit analysis. It compiles a program that conforms to the size and distribution of the actual load applied in this study on the basis of the existing fixed commands. The distributed loads used in this study are four rectangular surface loads distributed on the four corners of the rectangle, the distance between the front and rear load centers is 1 m, and the distance between the left and right load centers is 0.4 m. The area where the moving load passes on the model surface forms two moving load bands, as shown in Figure 5.Figure 5
Moving load band.It minimizes the size of the model while meeting the analysis requirements. So when writing subroutines, time (1)-L is used to ensure the coordination of time and loading area. time(1) represents the moving distance of the preceding load. L represents the distance between the front and rear load centers. In this way, the front load is applied to the model surface first. With the increase in time(1), the latter load also starts to be applied to the model surface after the former load moves by L. When the front load is out of the model, the latter load will continue to increase with time until it leaves the model. By this method, the established model can be reduced, thus saving computing resources and time.In this study, the single axle load is 28 kN, and the ground pressure of the tire is assumed to be 0.35 MPa. The front and rear axles all adopt a single wheel and single axle structure, assuming that the axle loads assigned to the front and rear axles are the same. The pressure surface load was selected as the load, and two moving tire load zones with a width of 10 cm were divided on the surface of the established asphalt pavement model. It applies moving loads to the model in the same manner as the moving load subroutine described earlier.


## 4.1.2. Model Establishment and Finite Element Mesh Division

The first step in finite element simulation analysis is to build a model that is close to reality. This study uses thexyz coordinate system. The model size is 600 mm ∗ 100 mm ∗ 100 mm, and the asphalt concrete bending beam models with initial cracks of 5 mm, 10 mm, 15 mm, 20 mm, and 25 mm, respectively, are shown in Figure 6. This model can better simulate the failure process of pavement with initial cracks and has been widely used in many simulation analyses.Figure 6
Asphalt concrete pavement beam model with initial reflection cracks (a = 15 mm).A crucial part of finite element analysis in ABAQUS is the division of elements. The rationality of the unit division has a great influence on the simulation results. According to the analysis principle of finite element method mentioned above, the precision of numerical calculation result of finite element often depends on the division quality of element mesh. However, the smaller the element size is selected, the denser the grid, the larger the memory required for model analysis, and the higher the requirements for the computer. An increase in the number of cells taken also increases the calculated intake error of the response. Therefore, in the model calculation, the division technique, the size of the element, and the number of elements contained in the model should take into account the performance of the computer. Reasonable mesh division can improve the efficiency of model analysis and calculation. The quality of meshing will affect the accuracy of the results, and there are several control factors.The size of the pavement model determines the size of the mesh, the element shape is hexahedron, and the division technique is structural division. The configuration of the vehicle load, such as the wheel spacing, and the mesh is refined at the wheel. The processing of details also has a certain impact on the accuracy of the results. If the mesh is too small, the computation will be too large, and if the mesh is too coarse, the result will be inaccurate. This model is an eight-node linear hexahedron element, with reduced integration, hourglass control, and a total of 8400 elements, as shown in Figure7.Figure 7
Finite element meshing of an asphalt concrete beam with initial cracks (a = 15 mm).


## 4.1.3. Determination of Boundary Conditions, Loads, and Material Properties

Boundary conditions also have a significant impact on the finite element analysis of pavement structural responses. According to different analysis needs, this study adopts the symmetric, antisymmetric, completely fixed boundary conditions and vehicle moving load of the initial analysis step. In the later analysis step, the displacement boundary load is set, and a set of reference points is established, so that the force of the structure can reach equilibrium. Elastic and viscoelastic materials that vary with temperature have different mechanical properties at different temperatures. It defines the time-domain viscoelastic material. Under the time list, the input method of viscoelastic material parameters is selected according to your needs. It enters the Prony series, selects Prony, and adds the corresponding shear relaxation modulus, bulk relaxation modulus, and relaxation time at different temperatures. At the same time, the corresponding temperature elastic modulus, Poisson’s ratio, and temperature values are added to the elastic material parameters.


## 4.2. Simulation of Crack Propagation in Different Environments

### 4.2.1. Stress Analysis of Crack Tip Corresponding to Different Viscoelastic Material Parameters

A viscoelastic asphalt concrete beam with an initial reflection crack is subjected to a moving load at a certain speed. It sets the ambient temperature to 20°C and changes the material parameters of asphalt concrete, which not only changes the proportion of asphalt and aggregate, etc., thereby increasing or decreasing the viscosity. The viscoelasticity of material I-materialV gradually decreases, and the reduction range is 5%. The stress change around the crack tip of the asphalt concrete beam model is shown in Figure 8.Figure 8
Stress cloud diagram around the crack of viscoelastic pavement under moving load (a = 25 mm, V = 80 dm/s, T = 20°C).It can be observed from Figure8 that material I has strong viscoelasticity, and the stress value around the crack tip is relatively high. The stress value during the whole moving load application process is slightly higher than the other four materials. With the change in material, the maximum stress value is slightly lower than that of material I, and the difference in stress value during vehicle operation is relatively large. Therefore, when the asphalt material ratio is within a certain range, the stress value around the crack remains stable. When the viscoelasticity of the asphalt mixture is high, the stress value around the reflection crack is also high, and the crack has the risk of further expansion. As the viscoelasticity of the asphalt mixture decreases, even the stress values around the crack tip decrease. However, in the case of hot weather, it is not conducive to the stability of the actual road surface. These are all factors that should be considered in practical engineering.


### 4.2.2. Stress Analysis around the Crack Tip under Moving Loads with Different Speeds

In the process of driving on the road, different road sections have different speed limits. Different types of vehicles also have their own different travel speeds. In the simulation analysis of road cracks, the speed of the vehicle is a factor that cannot be ignored. The application of the dynamic load is much closer to the actual load than the static load, which also makes the analysis results more realistic and accurate. Therefore, it is also necessary to explore the effect of different vehicle speeds on the cracks of asphalt pavement. Figure9 shows the stress cloud diagram of the beam model under the action of moving loads with different velocities and the comparison of the results.Figure 9
Stress cloud diagram around the crack of the pavement under moving load (a = 25 mm, V = 60 dm/s, T = 20°C).As shown in Figure10, moving loads at different speeds were applied and the effect of speed on the stress around the crack tip was observed. From the numerical results, it can be seen that under the environmental condition of 20°C, the stress values generated by the loads in different speed ranges are similar. That is, the change in velocity in one cycle has little effect on the stress at the crack tip.Figure 10
Stresses around cracks in viscoelastic pavement at different driving speeds (a = 25 mm, T = 20°C).


### 4.2.3. Stress Analysis of Crack Tip When Different Initial Temperature Fields Are Applied

In the research process of asphalt concrete pavement, because the asphalt concrete material is greatly affected by temperature, the crack state, propagation phenomenon, and the stress around the tip of the asphalt pavement with initial crack will be different at different ambient temperatures. Figure11 shows the stress nephogram around the crack of viscoelastic asphalt pavement at 20°C. In recent years, there have been many studies on the effect of temperature on asphalt mixtures. Many research results also have great practical application value. In this study, ambient temperature fields with different temperature values are added, and displacement loads equivalent to the actual ones are applied to explore the mechanical properties of asphalt concrete pavements with initial reflection cracks of different sizes at different ambient temperatures. Tables 1–5 provide the simulation results of analysis and comparison at different temperatures with 5 mm, 10 mm, 15 mm, 20 mm, and 25 mm initial reflection cracks, respectively.Figure 11
Stress cloud diagram around cracks of viscoelastic asphalt pavement at 20°C (a = 15 mm).Table 1
Stress around cracks of viscoelastic asphalt pavement under different temperature conditions (a = 5 mm).
20°C30°C40°C50°C60°CTime (s)Misesmax (kPa)Time (s)Misesmax (kPa)Time (s)Misesmax (kPa)Time (s)Misesmax (kPa)Time (s)Misesmax (kPa)00000000000.16660.15000.14200.14000.13850.213320.210000.28400.28000.27700.319980.315000.312600.312000.311550.426640.420000.416800.416000.415400.533300.525000.521000.520000.519250.639960.630000.625200.624000.62310Table 2
Stresses around cracks of viscoelastic asphalt pavement under different temperature conditions (a = 10 mm).
20°C30°C40°C50°C60°CTime (s)Misesmax (kPa)Time (s)Misesmax (kPa)Time (s)Misesmax (kPa)Time (s)Misesmax (kPa)Time (s)Misesmax (kPa)00000000000.17660.15800.14800.14400.14030.214320.210800.29000.28400.27880.320980.315800.313200.312400.311730.427640.420800.417400.416400.415580.534300.525800.521600.520400.519430.640960.630800.625800.624400.62328Table 3
Stresses around cracks of viscoelastic asphalt pavement under different temperature conditions (a = 15 mm).
20°C30°C40°C50°C60°CTime (s)Misesmax (kPa)Time (s)Misesmax (kPa)Time (s)Misesmax (kPa)Time (s)Misesmax (kPa)Time (s)Misesmax (kPa)00000000000.112660.110800.19800.19400.19030.219320.215800.214000.213400.212880.325980.320800.318200.317400.316730.432640.425800.422400.421400.420580.539300.530800.526600.525400.524430.645960.635800.630800.629400.62828Table 4
Stresses around cracks of viscoelastic asphalt pavement under different temperature conditions (a = 20 mm).
20°C30°C40°C50°C60°CTime (s)Misesmax (kPa)Time (s)Misesmax (kPa)Time (s)Misesmax (kPa)Time (s)Misesmax (kPa)Time (s)Misesmax (kPa)00000000000.117160.115300.114300.113900.113530.223820.220300.218500.217900.217380.330480.325300.322700.321900.321230.437140.430300.426900.425900.425080.543800.535300.531100.529900.528930.650460.640300.635300.633900.63278Table 5
Stresses around cracks of viscoelastic pavement at different temperatures (a = 25 mm).
20°C30°C40°C50°C60°CTime (s)Misesmax (kPa)Time (s)Misesmax (kPa)Time (s)Misesmax (kPa)Time (s)Misesmax (kPa)Time (s)Misesmax (kPa)00000000000.127110.120300.118400.117700.116530.233770.225300.222600.221700.220380.340430.330300.326800.325700.324230.447090.435300.431000.429700.428080.553750.540300.535200.533700.531930.660410.645300.639400.637700.63578It can be seen from Tables1–5 that under different ambient temperature conditions, when other conditions are the same, the stress around the crack tip of different sizes is different. As the ambient temperature decreases, the stress at the crack tip increases, and as the analysis time increases, the stress also increases.It can be seen from Figure12 that the impact on the asphalt pavement with initial crack elastic and viscoelastic materials is more obvious under the low-temperature environment. That is to say, the damage mostly occurs when the ambient temperature is low. For example, in winter, the numerical analysis results are consistent with reality. It shows that the numerical simulation process of viscoelastic materials can more effectively simulate the crack propagation process and mechanical behavior of asphalt mixture trabeculae under different ambient temperatures.Figure 12
Stresses around cracks in viscoelastic asphalt pavement under different temperature conditions (v = 80 dm/s, T = 20°C).Figure13 is an analysis diagram of the tip stress results of numerical simulation of initial cracks of different sizes at the same temperature. It can be seen from the figure that when the crack size is different, the stress at its tip is also different. When the crack is smaller, the stress value at the tip is larger. When a = 15 mm, the stress value at the crack tip is the largest. However, when a = 20 mm and a = 25 mm, the stress at the crack tip is significantly reduced. It indicates that the crack may be easier to crack when the crack value is smaller. After the crack has grown to a certain length, the growth rate will slow down.Figure 13
Surrounding stress of viscoelastic pavement with initial cracks of different lengths (T = 20°C).


## 4.2.1. Stress Analysis of Crack Tip Corresponding to Different Viscoelastic Material Parameters

A viscoelastic asphalt concrete beam with an initial reflection crack is subjected to a moving load at a certain speed. It sets the ambient temperature to 20°C and changes the material parameters of asphalt concrete, which not only changes the proportion of asphalt and aggregate, etc., thereby increasing or decreasing the viscosity. The viscoelasticity of material I-materialV gradually decreases, and the reduction range is 5%. The stress change around the crack tip of the asphalt concrete beam model is shown in Figure 8.Figure 8
Stress cloud diagram around the crack of viscoelastic pavement under moving load (a = 25 mm, V = 80 dm/s, T = 20°C).It can be observed from Figure8 that material I has strong viscoelasticity, and the stress value around the crack tip is relatively high. The stress value during the whole moving load application process is slightly higher than the other four materials. With the change in material, the maximum stress value is slightly lower than that of material I, and the difference in stress value during vehicle operation is relatively large. Therefore, when the asphalt material ratio is within a certain range, the stress value around the crack remains stable. When the viscoelasticity of the asphalt mixture is high, the stress value around the reflection crack is also high, and the crack has the risk of further expansion. As the viscoelasticity of the asphalt mixture decreases, even the stress values around the crack tip decrease. However, in the case of hot weather, it is not conducive to the stability of the actual road surface. These are all factors that should be considered in practical engineering.


## 4.2.2. Stress Analysis around the Crack Tip under Moving Loads with Different Speeds

In the process of driving on the road, different road sections have different speed limits. Different types of vehicles also have their own different travel speeds. In the simulation analysis of road cracks, the speed of the vehicle is a factor that cannot be ignored. The application of the dynamic load is much closer to the actual load than the static load, which also makes the analysis results more realistic and accurate. Therefore, it is also necessary to explore the effect of different vehicle speeds on the cracks of asphalt pavement. Figure9 shows the stress cloud diagram of the beam model under the action of moving loads with different velocities and the comparison of the results.Figure 9
Stress cloud diagram around the crack of the pavement under moving load (a = 25 mm, V = 60 dm/s, T = 20°C).As shown in Figure10, moving loads at different speeds were applied and the effect of speed on the stress around the crack tip was observed. From the numerical results, it can be seen that under the environmental condition of 20°C, the stress values generated by the loads in different speed ranges are similar. That is, the change in velocity in one cycle has little effect on the stress at the crack tip.Figure 10
Stresses around cracks in viscoelastic pavement at different driving speeds (a = 25 mm, T = 20°C).


## 4.2.3. Stress Analysis of Crack Tip When Different Initial Temperature Fields Are Applied

In the research process of asphalt concrete pavement, because the asphalt concrete material is greatly affected by temperature, the crack state, propagation phenomenon, and the stress around the tip of the asphalt pavement with initial crack will be different at different ambient temperatures. Figure11 shows the stress nephogram around the crack of viscoelastic asphalt pavement at 20°C. In recent years, there have been many studies on the effect of temperature on asphalt mixtures. Many research results also have great practical application value. In this study, ambient temperature fields with different temperature values are added, and displacement loads equivalent to the actual ones are applied to explore the mechanical properties of asphalt concrete pavements with initial reflection cracks of different sizes at different ambient temperatures. Tables 1–5 provide the simulation results of analysis and comparison at different temperatures with 5 mm, 10 mm, 15 mm, 20 mm, and 25 mm initial reflection cracks, respectively.Figure 11
Stress cloud diagram around cracks of viscoelastic asphalt pavement at 20°C (a = 15 mm).Table 1
Stress around cracks of viscoelastic asphalt pavement under different temperature conditions (a = 5 mm).
20°C30°C40°C50°C60°CTime (s)Misesmax (kPa)Time (s)Misesmax (kPa)Time (s)Misesmax (kPa)Time (s)Misesmax (kPa)Time (s)Misesmax (kPa)00000000000.16660.15000.14200.14000.13850.213320.210000.28400.28000.27700.319980.315000.312600.312000.311550.426640.420000.416800.416000.415400.533300.525000.521000.520000.519250.639960.630000.625200.624000.62310Table 2
Stresses around cracks of viscoelastic asphalt pavement under different temperature conditions (a = 10 mm).
20°C30°C40°C50°C60°CTime (s)Misesmax (kPa)Time (s)Misesmax (kPa)Time (s)Misesmax (kPa)Time (s)Misesmax (kPa)Time (s)Misesmax (kPa)00000000000.17660.15800.14800.14400.14030.214320.210800.29000.28400.27880.320980.315800.313200.312400.311730.427640.420800.417400.416400.415580.534300.525800.521600.520400.519430.640960.630800.625800.624400.62328Table 3
Stresses around cracks of viscoelastic asphalt pavement under different temperature conditions (a = 15 mm).
20°C30°C40°C50°C60°CTime (s)Misesmax (kPa)Time (s)Misesmax (kPa)Time (s)Misesmax (kPa)Time (s)Misesmax (kPa)Time (s)Misesmax (kPa)00000000000.112660.110800.19800.19400.19030.219320.215800.214000.213400.212880.325980.320800.318200.317400.316730.432640.425800.422400.421400.420580.539300.530800.526600.525400.524430.645960.635800.630800.629400.62828Table 4
Stresses around cracks of viscoelastic asphalt pavement under different temperature conditions (a = 20 mm).
20°C30°C40°C50°C60°CTime (s)Misesmax (kPa)Time (s)Misesmax (kPa)Time (s)Misesmax (kPa)Time (s)Misesmax (kPa)Time (s)Misesmax (kPa)00000000000.117160.115300.114300.113900.113530.223820.220300.218500.217900.217380.330480.325300.322700.321900.321230.437140.430300.426900.425900.425080.543800.535300.531100.529900.528930.650460.640300.635300.633900.63278Table 5
Stresses around cracks of viscoelastic pavement at different temperatures (a = 25 mm).
20°C30°C40°C50°C60°CTime (s)Misesmax (kPa)Time (s)Misesmax (kPa)Time (s)Misesmax (kPa)Time (s)Misesmax (kPa)Time (s)Misesmax (kPa)00000000000.127110.120300.118400.117700.116530.233770.225300.222600.221700.220380.340430.330300.326800.325700.324230.447090.435300.431000.429700.428080.553750.540300.535200.533700.531930.660410.645300.639400.637700.63578It can be seen from Tables1–5 that under different ambient temperature conditions, when other conditions are the same, the stress around the crack tip of different sizes is different. As the ambient temperature decreases, the stress at the crack tip increases, and as the analysis time increases, the stress also increases.It can be seen from Figure12 that the impact on the asphalt pavement with initial crack elastic and viscoelastic materials is more obvious under the low-temperature environment. That is to say, the damage mostly occurs when the ambient temperature is low. For example, in winter, the numerical analysis results are consistent with reality. It shows that the numerical simulation process of viscoelastic materials can more effectively simulate the crack propagation process and mechanical behavior of asphalt mixture trabeculae under different ambient temperatures.Figure 12
Stresses around cracks in viscoelastic asphalt pavement under different temperature conditions (v = 80 dm/s, T = 20°C).Figure13 is an analysis diagram of the tip stress results of numerical simulation of initial cracks of different sizes at the same temperature. It can be seen from the figure that when the crack size is different, the stress at its tip is also different. When the crack is smaller, the stress value at the tip is larger. When a = 15 mm, the stress value at the crack tip is the largest. However, when a = 20 mm and a = 25 mm, the stress at the crack tip is significantly reduced. It indicates that the crack may be easier to crack when the crack value is smaller. After the crack has grown to a certain length, the growth rate will slow down.Figure 13
Surrounding stress of viscoelastic pavement with initial cracks of different lengths (T = 20°C).


## 5. Conclusions

Based on the study and analysis of fracture mechanics and viscoelasticity theory, this study simulates and analyzes the asphalt concrete pavement beam model with initial reflection cracks and the coexistence model of multiple cracks and verifies it with the actual life. The following conclusions can be obtained:(1)
The viscoelastic asphalt concrete beam with the initial reflection crack is subjected to the actual moving load. The load moves at a constant speed. When the ambient temperature remains unchanged, the stress around the crack does not change much with the change in speed; that is, the crack propagation is little affected by the speed.(2)
A viscoelastic asphalt concrete beam with an initial reflection crack is subjected to a moving load at a certain speed. It keeps the ambient temperature unchanged, changes the material parameters of asphalt concrete, and changes the proportion of asphalt and aggregate. Therefore, when the asphalt material ratio is within a certain range, the stress value around the crack remains stable. When the viscoelasticity of the asphalt mixture is high, the stress value around the reflection crack is also high, and the crack has the risk of further expansion. As the viscoelasticity of the asphalt mixture decreases, even if the stress value around the crack tip decreases, it is not conducive to the stability of the actual road surface in the case of hot weather. These are all factors that should be considered in practical engineering.(3)
A viscoelastic asphalt concrete beam with an initial reflection crack is subjected to a moving load at a certain speed. It changes the surrounding ambient temperature field and should be the time-temperature equivalence of the viscoelastic material. Different temperature values make their material properties different. Therefore, the stress values around the reflection crack at different temperatures vary greatly. With the decrease in the ambient temperature, the stress value is higher; that is, the crack is more likely to expand in the low-temperature environment.



---

*Source: 6603130-2022-06-24.xml*