Abstract
A nonlocal damage-plasticity coupled area-weighted model was proposed to investigate the fracture process of the ductile fracture materials with softening behavior, such as asphalt concrete materials. The model can overcome the mesh sensitivity problem in the analysis. A subroutine of ABAQUS software was developed to apply this nonlocal damage-plasticity coupled area-weighted model into finite element analysis. Then, the subroutine was adopted in finite element models to simulate several loading conditions. The results indicate that the nonlocal damage-plasticity coupled area-weighted model and the subroutine describe the softening behavior of asphalt concrete materials better than traditional softening models. And the robustness and stability of the proposed model were also justified. The proposed model provides a concrete and accurate alternate to predict the damage condition of asphalt concrete.
1. Introduction
The local deformation band occurs when the structures and/or materials are in unsteady states, which is the initiation of the ultimate collapse. These phenomena are usually observed in engineering, such as the shear band of the metal and local strain zone of asphalt concrete [1–3]. It is well known that the strength of a structure is determined by the strength of the weakened region rather than the mean strength. When the local deformation zone occurs in the weakened region, microcracks will form and propagate in microcosmic level. The essence of local deformation band reveals the complex nonlinearity of the broken behavior of materials. Tremendous attention has been paid to the broken behavior of various materials. The fracture mechanics is usually assumed an effective tool to simulate cracking propagation. However, the fracture mechanics can be merely applied on the condition where macrocracks exist. The damage mechanics has been widely used to solve the cracking propagation problems for decades [4]. Nevertheless, traditional damage models are usually not nonlocal models and are facing the mesh sensitivity phenomenon, which leads to irrational results when adopted in finite element methods to simulate damage behaviors of material and structure.
The mesh sensitivity phenomenon in traditional damage models is believed to be due to the fact that they can only use elliptic equations rather than hyperbolic equations as their control equations. When the transformation from elliptic equations to hyperbolic equations is encountered, the extra supplementary equations are needed to obtain the closure equations. How to find appropriate supplementary equations is a worldwide famous puzzle [5]. The supplementary equations need to not only cover the statistic characteristics of materials in microlevel, but also have a succinct expression. Great efforts have been put into finding supplementary equations [6]. However, there are no conclusive solutions for this problem yet. With the development of numerical science, some nonlocal damage models have been proposed to deal with softening behaviors of materials and widely utilized to describe the damage behaviors of materials and structures without the mesh sensitivity phenomenon.
Pijaudier-Cabot and Benallal [7] proposed a widely used nonlocal damage-plasticity coupled area-weighted model. They defined a variable called nonlocal strain, the mean value of all equivalent strains in continuum medium points over whole local deformation zone: where = nonlocal strain; = equivalent strain in a continuum medium point; is the characterization volume.
For homogeneous materials, . And can be defined as is assumed to be a Gauss distribution function, and can be regarded as a regularized factor: Since the parameter represents the influence range of the nonlocal region, it was named as the characteristic length.
The nonlocal damage-plasticity coupled area-weighted model [7] has been proved to have good accuracy and preciseness on simulating softening behaviors and ductile fracture behaviors of materials.
Saanouni et al. [8, 9] assumed the damage factor as the nonlocal area-weighted variable and applied it to analyze the creep fracture of the metal.
The two kinds of nonlocal area-weighted models mentioned above [7–9] are based on the traditional plastic theory, and the traditional yield criteria, flow rules, and consistency condition can be easily modified and applied with them. Meanwhile, the mesh sensitive phenomenon is not a problem in such two models. However, due to the quasi brittle fracture characteristic and softening behavior of concrete materials before the ultimate failure, the two nonlocal area-weighted damage models are not suitable for it.
Aifantis and his coworkers [10–12] proposed a gradient-dependent plastic nonlocal model. The loading conditions in plastic theory were modified, and a Laplacian term was introduced for the equivalent strain and a factor corresponding to the Laplacian term. The modified loading condition is listed as follows: where is the equivalent stress; is the equivalent strain; is the plastic strength, which is the function of the equivalent strain when the material is plastic hardening; is factor of the Laplacian term; and is the Laplacian term of , which can induce the effect of the dislocation, when dealing with the metal material.
de Borst et al. [13] modified the Aifantis model and proposed the de Borst-Pamin-Geers gradient-enhanced nonlocal constitutive model. Peerlings and Geers are well recognized for their contribution in the gradient-dependent plastic nonlocal model. They further developed the nonlocal constitutive model from mathematical perspectives around 2000 and proposed a rigorous nonlocal model named Peerlings-Geers-de Borst-Brekelmans gradient-enhanced nonlocal constitutive model [14].
However, a huge defect exists in the gradient-enhanced nonlocal constitutive models [12–14] mentioned above; that is, the models were not based on the second law of thermodynamics, and the Clausius-Duhem inequality is not satisfied. Svedberg and Runesson [15] derived a more rigorous gradient-enhanced nonlocal constitutive model based on the second law of thermodynamics. It should be noted that the Svedberg and Runesson model [15] contains an internal variable from the thermodynamics frame, which is developed from the Valanis theory [16, 17]. Therefore, this model is complicated comparing to the Aifantis nonlocal model [10–12]. Nedjar [18] proposed a gradient dependence nonlocal constitutive model based on the thermodynamics equations including a nonlocal damage variable . When there is no damage, and when the material loses its loading capacity, . Liebe et al. [19] proposed a plastic gradient dependence coupled with damage nonlocal constitutive model as well. Both the two models above [18, 19] assume there is energy dissipation in local region during the unsteady state of structures and materials.
It was reported that nonlocal models might be applicable for the analysis of cement concrete material. However, asphalt concrete material is not similar to the cement concrete in forming material strength in microlevel. The aggregates in asphalt concrete form the framework to bear the loading, and the asphalt mastic provides cohesion between aggregates and mastic. The asphalt concrete material always manifests the viscoelastic characteristic which is due to the rolling and compression of molecular chain in asphalt mastic. The strength of mastic is mainly dominated by the Van der Waals force in microlevel, which is much smaller than the electrostatic force. Therefore, the asphalt mastic is usually confronted with permanent deformation, even when the external loading is small. Meanwhile, the ultimate fracture of asphalt concrete is mainly determined by the permanent deformation. Therefore, the damage behavior and essence of asphalt concrete material are more involved compared with those of the cement concrete material. In continuum mechanics, the asphalt concrete material manifests the ductile fracture behavior like metal materials. Moreover, it manifests the softening behavior after strength peak like cement concrete. Based on the above phenomena, we need to develop a nonlocal model to simulate the damage behavior of asphalt concrete. Although the Pijaudier-Cabot-Bazant area-weighted nonlocal model is suitable for simulating the nonlocal ductile fracture behavior of asphalt concrete, a damage factor for simulating the softening behavior of asphalt concrete material after the strength peak is still necessary. In this paper, a modified Pijaudier-Cabot-Bazant area-weighted nonlocal model coupled with damage theory was proposed to fit the above requirements.
2. Issues for Traditional Damage-Plasticity Coupled Models
There is a very severe mesh sensitivity phenomenon leading to distortion of numerical results when using traditional local models to simulate the softening behavior of material and structure. When the material or structure has manifested the softening behavior, the deformation will happen inside one element. This phenomenon cannot be erased by decreasing element size, which conflicts with the essence of the finite element method. Figure 1 gives a mesh sensitive phenomenon happening in an X-type beam under uniaxial tension loading, when the traditional damage-plasticity coupled model is applied in the finite element model. It can be seen that, by merely changing the mesh density, the damage field for the model with sparse grids (Figure 1(a)) is quite different from the model with dense grids (Figure 1(b)). The damage field for the model with sparse grids (Figure 1(a)) focuses on two elements in the middle with a peak of the damage value of 0.993, while the damage field for the model with dense grids (Figure 1(b)) covers more finite elements with a peak of damage value of 0.785. That is to say, when the traditional damage-plasticity coupled models are utilized to simulate the softening behavior of materials, the numerical solution is not converging or stable.
(a) Sparse grids
(b) Dense grids
The reaction forces on the two models in Figures 1(a) and 1(b) are plotted in Figure 2. In the elastic phase, the reaction force curves are identical with each other. However, when softening behavior occurs, it can be obviously observed that significant difference exists between the reaction forces from the two models. It must be admitted that traditional plasticity-damage coupled models cannot provide an accurate and converged solution for the softening behavior of materials.
3. Nonlocal Damage-Plasticity Coupled Area-Weighted Model
3.1. Theoretical Framework
Based on the nonlocal area-weighted model, the constitutive equation of local variables stress tensor and strain tensor at the material point is listed as follows: where is the four-order stiffness tensor.
Assume there is an area covering point ; the relation between the nonlocal strain tensor and the local strain tensor at point can be expressed as where is a material point adjacent to point in area and is the regularized weight function and can be written as where is the Gauss type weighting function, is the vector radius, and .
In this model, is defined as where is the size of the nonlocal zone. The permanent deformation occurs in asphalt concrete under loading. Therefore, the plastic property should be defined. Based on the small deformation assumption, the strain is decomposed as elastic strain and plastic strain, as shown in where is the whole strain tensor, is the elastic strain tensor, and is the plastic strain tensor.
According to the flow rule, the plastic strain tensor is written as follows: where is the plastic multiplier, is the loading condition, and the rate of (10) is as follows: The superscript denotes the rate form of the variable, and the deviatoric stress tensor is defined as , and the Von-Mises equivalent effective stress is as follows: Assume the initial yield stress is , and according to the plastic flow theory, the plastic loading condition is given in the power series hardening mode [20]: where is the power series hardening parameter, is the dimensionless elastic modulus parameter, and is equivalent effective plastic strain and can be expressed as follows: According to the consistency condition, , , and , the rate form of the plastic multiplier is derived as Since the plastic deformation leads to no increase on stress, the softening behavior is not considered in traditional plastic models. However, asphalt concrete has softening behavior. A damage factor coupled with the plastic strain has been defined to explain the softening behavior of asphalt concrete. Based on the principle of equivalent strain from Lemaitre, the damage constitutive equation in effective stress space is expressed as follows: The superscript “” denotes the variable in the effective stress space, is the initial elastic four-order tensor without damage, and and are elastic strain tensor and plastic strain tensor, respectively. The initialization and evolution damage variable is determined by the damage evolution function, which is proposed by Saanouni et al. [21] and developed by Shen et al. [22] and Cao [23]. The rate form of the damage evolution function is listed as follows: where , , and are the material parameters.
3.2. Development on Subroutine of the Nonlocal Damage-Plasticity Coupled Area-Weighted Model
Based on the software platform ABAQUS, the nonlocal damage-plasticity coupled area-weighted model proposed in this paper has been developed via the user defined material subroutine (UMAT) of ABAQUS using FORTRAN language. The flowchart of developing the UMAT is listed in Figure 3.
4. Numerical Tests to Validate the Proposed Model
4.1. Simulation on the Uniaxial Tension Test of X-Type Beam
The X-type beam with uniaxial tension was simulated with the proposed model. The geometric parameters of the X-type beam are shown in Figure 4. The thickness of the beam is 3. Basic mechanical properties of this beam are summarized in Table 1. Three different mesh densities were applied on the beam and the element sizes were 900, 1200, and 1800.
In order to monitor the softening scale of the material, a user defined variable (SDV) was defined as the damage factor via UMAT. When there is no damage in material, and when the material totally loses its capacity, .
Figures 5 and 6 provide the damage contours and tension strain contours of models with different mesh density. It can be found from Figures 5(a), 5(b), and 5(c) that the damage distributes in a fixed region regardless of the mesh density. The same phenomenon can be observed in the tension strain distribution contours. It is worth mentioning that the mesh density should remain in some degree to obtain rational results when using this proposed model for the fact that there is some disordered deformation in some elements when the element number is 900 (Figure 5(a)).
(a) Number of elements = 900
(b) Number of elements = 1200
(c) Number of elements = 1800
(a) Number of elements = 900
(b) Number of elements = 1200
(c) Number of elements = 1800
The reaction forces for models with the three different mesh densities are shown in Figure 7. They are almost overlapped with each other in the whole loading procedure except the model with 900 elements after some time. The dissimilarity is mainly induced by the huge element deformation.
(a)
(b)
4.2. Simulation on the Uniaxial Tension Test of Asymptotic Beam
When uniaxial loading is imposed on the asymptotic beam, there will exist a local deformation zone. The geometric parameters of the asymptotic beam are 150 in height, 20 in width, and 5 in width (without dimension), and the central section area is about 90% of the boundary one.
Since the area defect of the central section is not as severe as the rest of the regions, the stress intensity phenomenon is not comparably strong. The constraints were employed in the whole boundary node to form a fixed support, so that the stress intensity would not appear in the boundary region. The finite element models were divided into 15, 20, and 30 elements in length, respectively, and 3 elements in thickness. The damage distributions of the asymptotic beams with different mesh densities are shown in Figure 8. No mesh dependence phenomenon can be observed among the beams. The reaction force curves on beams with different mesh densities were almost identical (Figure 9), indicating that no mesh sensitivity problem exists. Additionally, no disordered deformation was observed. Therefore, the nonlocal model proposed was proven to be robust and stable.
(a) Number of elements = 900
(b) Number of elements = 1200
(c) Number of elements = 1800
4.3. Simulation on Three-Point Bending Beam
Bending is the most common condition happening in structure, which induces the ultimate collapse. The following part is to validate the proposed model in the bending condition. A three-point bending beam was selected with 150 in height, 20 in width, and 5 in width (without dimension). Three different mesh densities were generated on the finite element models, including 1200, 1600, and 1800 elements, respectively. The basic mechanical properties of the beam are shown in Table 1. The damage distributions of the models with different mesh densities are very similar, as shown in Figure 10, indicating that the proposed model is validated in bending conditions of materials.
(a) Number of elements = 1200
(b) Number of elements = 1600
(c) Number of elements = 1800
5. Conclusion
In order to simulate the damage process of asphalt concrete material, especially the softening behavior, a modified Pijaudier-Cabot-Bazant area-weighted nonlocal model coupled with damage theory was developed and the proposed nonlocal model was integrated into finite element method using the ABAQUS/Standard and ABAQUS/UMAT. The conclusion can be summarized as follows.
When utilizing the proposed nonlocal model to simulate the ductile fracture behavior and softening behavior of materials, the local deformation region can be manifested without the mesh dependence problem, and closed-form solution can also be obtained.
A series of loading condition numerical tests were conducted in finite element models with the proposed nonlocal model, including the uniaxial tension numerical tests and three-point bending numerical test. It indicates that, compared with traditional damage-plasticity coupled theory, the proposed nonlocal model successfully avoids the mesh sensitivity problem and is more robust and stable when simulating the ductile fracture and softening behavior of the material and structure.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.