Abstract
In this paper, the elastoplastic damage constitutive model of concrete is established based on the damage energy release rate to analyze the failure of concrete structures. Based on the user-defined material program VUMAT interface provided by the ABAQUS, the explicit algorithm is used to secondary develop the concrete model. The model is verified from mechanical behavior, mesh-dependence, and stress state under reciprocating loading. The results show that the model can reasonably simulate the mechanical state of concrete under compression and tension and can effectively solve the mesh-dependence problem in finite element problems. The simulation results were still consistent when the mesh size was 2.5 mm. The low cycle reciprocating simulation test of the concrete column also proves that the model has better hysteretic and bidirectional coupling performance. It is verified that the VUMAT material subroutine can meet the analysis requirements of reinforced concrete structures and can be used to simulate the dynamic response of structures with strong nonlinear stages.
1. Introduction
Earthquake is one of the great natural disasters in the world, many parts of the world suffer from severe earthquake disaster loss. How to reduce the loss caused by the earthquake and ensure safety of people’s life and property is the current need. With the continuous development of urban construction, there are more and more complex high-rise building structures, which put forward new requirements for engineering technology which requires advanced technical experiments, analysis, and simulation results to ensure stability and safety of the structure [1].
The investigation results of previous earthquake disasters show that the collapse of structures is the leading cause of casualties and property losses [2, 3]. Structures usually enter the elastoplastic stage under large earthquakes, and the analysis method of linear elasticity of structures can no longer reflect structures’ accurate and effective stress. With the significant progress of the constitutive relation theory, numerical analysis algorithms, and computer software and hardware, nonlinear whole-process analysis methods for engineering structures have gradually become popular. The nonlinear analysis method can reflect the component’s internal force and deformation state at each moment under earthquake and other disasters, find the position in the stress and deformation, and judge the yield mechanism, weak link, and possible failure type of the structure. Reinforced concrete (RC) structure is still the main structural form of building structure engineering at present and for a long time, so conducting nonlinear dynamic elastoplastic analysis of the RC structure under earthquake action is essential.
Concrete is made of coarse and fine aggregate by gelation. Due to the inhomogeneity of its internal structure and cohesion-friction characteristics, it shows complex nonlinear elastoplastic characteristics under a multiaxial stress state. Concrete is the most widely used engineering material in civil engineering, used in building structures, bridges, dams, and other complex projects. So, the establishment and realization of the constitutive relation model of concrete material is the critical link in structural analysis and simulation, but it also has higher requirements for theory and technology. The constitutive relation refers to the functional relationship between stress and strain rate or a material’s stress tensor and strain tensor. Because the constitutive relation describes the basic mechanical properties of materials, it forms the basis for studying the deformation and movement of components and structures under external action and plays a vital role in the mechanical behavior of the structures. For the problem of structural damage, the constitutive relation should include the description of material damage and failure from the material level. The purpose of the study of concrete constitutive relation is to analyze the mechanical behavior of engineering structures under external and self-action to provide a scientific basis for the rational design of concrete structures.
With the wide application of finite element software such as ABAQUS, ANSYS, and SAP2000 in structural analysis [4, 5], many scholars began to use finite element software for dynamic elastoplastic analysis of structures. Xinzheng team of Tsinghua University developed a subprogram of concrete constitutive model based on the implicit algorithm for the applicability of the concrete constitutive model of the fiber beam element in finite element software [6]. Many scholars rely on the secondary development of software computing platforms to improve the material properties in the software so that the calculation can save time and get relatively high precision calculation results [7–10]. The Bézier method and Galerkin method are widely used in numerical simulation algorithm applications which have the advantages of high accuracy and wide applicability. Kabir et al. developed a nonlinear vibration and postbuckling configuration of Euler–Bernoulli composite beams augmented with graphene nanoplatelets (GnP) based on the Bézier method [11, 12]. Belardi et al. based on the Galerkin method for bending analysis of linear orthogonal anisotropic composite plate, the improved model was proved to have good consistency by comparison [13, 14]. Chen et al. [15] proposed an anisotropic critical state soil model to simulate soil’s plastic behavior, which significantly impacts the swelling response of anisotropy caused by plastic deviational strain. In recent years, elastoplastic dynamic time history analysis has been successfully applied in engineering, and experts have widely recognized its analysis results. The finite element software has a rich material library and powerful preprocessing and postprocessing function, especially the explicit analysis module can avoid the nonconvergence problem of the implicit algorithm when dealing with highly nonlinear problems. At the same time, the dynamic elastoplastic analysis of large complex structures can be carried out using this technique, saving a lot of calculation storage space and improving the speed of analysis and simulation by using parallel computing technology.
The standard methods for solving engineering structures are Newmark-β method and central difference method, in which the Newmark-β method is an implicit algorithm and the central difference method is an explicit algorithm. The implicit algorithm must introduce the Jacobian matrix K in the solution process. The implicit algorithm is robust and ensures the calculation accuracy, but it needs to compute the inverse matrix of K at each iteration, which is prone to numerical convergence problems such as stiffness mutation and negative stiffness. The explicit algorithm does not need to invert the matrix K in the solution process, which can avoid the nonconvergence problem encountered in the strong nonlinear stage. However, the disadvantage of the explicit algorithm is that it is generally a conditionally stable algorithm, and the required step size of the solution is several orders of magnitude smaller than the implicit algorithm. Most concrete constitutive models adopt implicit algorithms, which are not ideal for simulating large deformation and collapse of structures. However, the concrete damage plastic (CDP) model of ABAQUS, as an empirical constitutive model, cannot sufficiently describe the damage evolution process of concrete and cannot be applied to the beam unit.
This paper establishes a double scalar elastoplastic damage constitutive relationship of concrete by introducing the energy equivalent strain based on the damage energy release rate. Based on the VUMAT subroutine interface of the ABAQUS, the solution process of the concrete constitutive model is divided into three calculation steps: elastic prediction, shaping correction, and damage correction by the operator separation algorithm. The concrete uniaxial compression and tension models were established to verify the accuracy of the constitutive model, respectively, and the mesh-dependence analysis was carried out on the concrete constitutive model. The low cycle reciprocating test model of reinforced concrete specimens is established to verify whether the model can reasonably and accurately grasp the elastoplastic stress characteristics of reinforced concrete beam-column and its effectiveness in nonlinear response analysis of structures under earthquake.
2. The Establishment of Concrete’s Elastoplastic Damage Constitutive Model
The nonlinear behavior of concrete under stress includes the influence of microcrack and plastic flow, which results in the following significant characteristics. An unstable region after peak stress is accompanied by apparent stiffness degradation and strength softening. The hysteretic characteristics during loading and unloading are as follows: when the deformation exceeds the fixed threshold value, the concrete has unrecoverable deformation after unloading entirely. The strength and ductility of the material increase obviously when there are confinements (such as biaxial compression stress state). Unilateral effect: tensile strength and compressive strength are obviously different [16]. The damage has the characteristic of directionality, and the antidirection loading will close the damage crack, which leads to the recovery of concrete stiffness.
The stress-strain theory diagram of concrete under uniaxial compression and tension is shown in Figure 1. In the tensile stage, the stress and strain of concrete basically obey the linear growth relationship. When the stress reaches 50%–70% of the tensile strength, the deformation growth accelerates, and the stress-strain relationship deviates from the linear and begins to show nonlinear characteristics. Subsequently, cracks appear on the surface of the concrete tensile specimen, and finally, the damage is caused by the stress concentration at the crack, which gradually accumulates and eventually leads to the failure of the member. Similar to tension, the uniaxial compression failure of concrete also produces cracks and defects inside the specimen, which develop due to stress and lead to the gradual accumulation of damage.
(a)
(b)
Using the theory of plastic mechanics to obtain the mechanical behavior of concrete can reflect the stress state of reciprocating loading to a certain extent. However, it is difficult to reflect the complex characteristics of concrete, such as stiffness degradation, strength softening, and material anisotropy. After long-term research, it is found that the damage theory is a reasonable choice for the constitutive relation modeling of concrete pseudo-brittle materials [17]. Currently, most of the constitutive relations of concrete are based on empirical constitutive models [18–20] to simulate the stress and strain of concrete under tension and compression. The research shows that the constitutive model of the damage evolution law described by continuum mechanics and irreversible thermodynamics can better express the constitutive relation of concrete based on the damage energy release rateWith the development of damage mechanics, the second-order microdamage tensor for anisotropic materials has been proposed, and the nonuniform growth of microdamage tensor components has been proved. Various coupled thermal-viscoplastic damage theories are derived, which provide a reference for the theoretical model of elastoplastic damage mechanics of concrete [21–23]. Therefore, this paper adopts concrete damage theory, energy variable theory, thermodynamic theorem, and energy dissipation principle to deduce the concrete damage constitutive model, introduces the concept of energy equivalent strain, and establishes damage evolution law and damage criterion in effective stress space.
2.1. Thermodynamic Fundamentals and Damage Constitutive Relations
Damage is an irreversible energy dissipation process, so the damage mechanics of damage and its mechanical effects must involve irreversible thermodynamics. According to the first law of thermodynamics, the energy change rate in the system concerning time is equal to the sum of the power done and the change rate of the heat energy transmitted to the system by the external force. Its differential form can be expressed as follows: where is density, is the differential of internal energy density with respect to time, is the heat generated per unit time by unit mass of heat source, is the effective stress tensor, is the effective strain rate, and is the heat outflow through unit area in unit time.
The first law of thermodynamics establishes the conversion and conservation relationship between mechanical energy and heat energy. The second law of thermodynamics explains the automatic realization of this conversion and the nature and development direction of the conversion process. The second law of thermodynamics applies to irreversible thermodynamic systems, stating that heat cannot be transferred from a cold body to a hot body without causing other changes. The second law of thermodynamics is usually expressed by Clausius–Duhem inequality. Considering that in the isothermal thermodynamic process, the temperature is constant and the influence of internal heat source and heat flow is not taken into account, the first law of thermodynamics and the second law of thermodynamics can be simplified as
According to formula (3), for any allowable isothermal pure mechanical process, part of the work done by the external force is converted into the Helmholtz free energy of the object, and the other part is used for energy dissipation in the deformation process, and its magnitude is nonnegative.
For isothermal pure mechanical processes, another commonly used free energy is Gibbs free energy , which satisfies the following Legendre transformation with Helmholtz .
Substituting the above equation into the formulas for the first and second laws of thermodynamics yields the following equation:
The material damage process is an irreversible thermodynamic process. In this process, heat dissipation is not considered in weak, that is, only isothermal pure mechanical process is considered, then the Helmholtz free energy of the material is the state function of strain tensor and damage variable , which is as follows:
Take the derivative of this with respect to time.
Substitute the above formula into the formula for the fundamental principles of thermodynamics.
Due to the arbitrariness of the strain rate, the necessary and sufficient condition to satisfy the above inequality is
Letwhere is the damage energy release rate.
From the formula (9), it is easy to derive the following:
Formula (12) is the stress-strain constitutive relation with damage variables. In the classical deterministic damage mechanics theory, the damage development equation is given according to the principle of maximum dissipation energy.where is the derivative of damage variable to time; is the rate of damage evolution, and is a kind of damage potential function or dissipation function.
By combining (12) and (13) with specific damage criteria, a complete material damage evolution process can be solved.
2.2. Double Scalar Damage Constitutive Relation
In the isothermal adiabatic state, it can be assumed that the elasticity of the material is decoupled from the plastic Helmholtz free energy. Therefore, the Helmholtz free energy of damaged materials can be decomposed into elastic and plastic parts. The damage of materials is only related to the elastic deformation, and the damage further affects the plastic deformation of the materials. Accordingly, the total strain tensor in the nonlinear development process of material stress should be expressed as the sum of elastic strain tensor and plastic strain tensor, namely
In order to describe the irreversible energy dissipation process of damage evolution, tensile damage variables and shear damage variables are introduced to describe the Helmholtz free energy degradation of the damaged material under positive effective forces and negative effective forces action, that is, the elastoplastic Helmholtz free energy of the damaged material is expressed asWhere and are the elastic and plastic Helmholtz free energy of the material, respectively.
Based on the damage energy release rate, the damage evolution laws of tensile damage variables and shear damage variables are established. They can be obtained according to formula (11).
Note that in the above model, only the elastic part of the Helmholtz free energy is considered, and the damage energy release rate thus obtained is only related to the initial elastic Helmholtz free energy of the material, i.e., .
Because the damage energy release rate controls the damage evolution, the equivalent transformation does not affect the nature of the damage development. In practical engineering, concrete can be regarded as an isotropic material in theory, and only the first invariants and the second invariants need to be considered without considering the full triaxiality. After derivation and simplification, the tension damage energy release rate formula and the shear damage energy release rate can be expressed aswhere is the shear damage parameter, and are the first invariant of the effective stress component and the second invariant of deviatoric stress .
With the development of damage mechanics, many scholars established orthogonal flow theory based on damage surface by analogy with classical plastic flow theory. Since the damage energy release rate reflects the change rate of damage energy consumption on damage development, a natural choice is to express the damage surface with the function of damage energy release rate. The damage function is established and the determined surface is called the damage surface. Then the occurrence criteria of tensile and shear damage can be expressed aswhere is any increasing scalar function of a scalar; determines the size of the damage surface at the current moment and is expressed aswhere is the damage energy release rate threshold of tensile damage and shear damage.
When the damage energy release rate exceeds the initial damage energy release rate threshold , the material will enter nonlinear state. When the damage energy release rate exceeds the historical maximum damage energy release rate , further damage will occur.
According to the orthogonal flow rule, the evolution rule of damage variables can be obtained as
After derivation, the expression of damage variable can be obtained as
3. Numerical Realization of Damage Constitutive Relation
In this paper, through the VUMAT subroutine interface of ABAQUS, the user-defined material subroutine of the elastoplastic damage constitutive model is compiled by the Fortran language. Due to the need for an iterative algorithm to calculate the solution of stress , incremental strain loading is used to complete the numerical realization of the constitutive relation. In fact, due to the intrinsic incremental nature of the evolutionary equation, the incremental variable loading method is used to solve the constitutive equation.
The effective stress space and the damage space are decoupled from each other, so the updating algorithm of the effective stress is similar to the stress updating algorithm in general elastoplastic mechanical constitutive calculation, and the operator separation algorithm is adopted. By combining the operator decomposition system with the implicit solution method, the basic solution process of the elastoplastic damage constitutive model can be obtained. Firstly, the maximum possible effective stress increment and Cauchy stress increment can be obtained by the elastic prediction. Then, the plastic strain was obtained by the Newton iteration method, and the effective stress was reflected on the yield surface. Finally, the Cauchy stress is modified by calculating two damage variables respectively. The flow chart of realizing the overall constitutive relation by the numerical algorithm is shown in Figure 2.
Considering the development of multidimensional damage, the effective stress is decomposed into positive and negative components during damage correction. Decomposing the stress tensor into the sum of positive and negative stresses:where the positive and negative stress tensors are respectively defined aswhere and are principal stress and corresponding direction vector of stress tensor , respectively; is the slope function, and is expressed as
In the iterative process of damage correction, if both and are satisfied, no damage loading occurs. Otherwise, damage loading occurs and damage variables are updated.
Finally, the final stress is calculated
4. Numerical Simulation Model Validation
In order to verify the accuracy of the VUMAT subroutine in simulating the concrete model, tests were carried out on mechanical behavior, mesh-dependence, and stress state under reciprocating loading. The results were compared with ABAQUS’s concrete damage plastic (CDP) material model.
4.1. Mechanical Behavior
The uniaxial compressive and uniaxial tensile test models of concrete were established to verify the mechanical behavior of the models. The analysis step uses dynamic, explicit procedure in ABAQUS/EXPLICIT module. For the uniaxial compression test, the component size is 150 ∗ 150 ∗ 300 mm. The uniaxial tensile test adopts dog bone type specimen, and the cross section size is 100 ∗ 100 mm. C3D8R mesh type were used in both models, and the mesh size was 10 mm globally. The loading mode is displacement loading, and the amplitude type is set to smooth. Three concrete strengths of C25, C30, and C35 were used for simulation. The comparison between simulation results and the CDP model is shown in Figure 3.
It can be seen from the figure that the VUMAT model can describe the stress of concrete under uniaxial compression and tension well. As the concrete damage energy release rate is used as the damage evolution standard in the model, it can be seen that the model is gentler in the descending section. The VUMAT model is consistent with the CDP model in simulating the uniaxial tensile state of concrete. The failure threshold of the material element is set in VUMAT, so when the damage evolution of the element reaches a certain degree, the element fails, and the fracture failure occurs in the simulation process.
4.2. Mesh-Dependence
The efficiency of finite element calculation is closely related to the mesh. Based on ensuring the accuracy, the larger the element size, the higher the computational efficiency. The concrete constitutive model based on the damage theory is applied to the integral points of elements. When the damage reaches the threshold, the deletion of elements is the deletion of integral points of elements. In practical engineering, the component’s macroscopic minimum visible crack size is 0.1 mm. Therefore, to better simulate the damage development and crack propagation, the element size of the structure along the crack propagation path will be pretty small. At the same time, it shows that the dynamic algorithm has a stability limit, and the smallest unit in the model determines its stability step size. The smaller the unit, the smaller the stability step size, and the higher the calculation cost.
This paper uses the concrete uniaxial compression test to verify the mesh-dependence of the VUMAT concrete elastoplastic damage subroutine. The mesh is divided into 15 mm, 10 mm, 5 mm, and 2.5 mm specifications to verify the mesh-dependence. Figure 4 shows the test results. When 2.5 mm mesh scale is used, the model contains 727,936 elements, significantly reducing computational efficiency. As seen from the enlarged area in the figure, the denser the grid, the faster the damage of the model develops, but the overall trend is still stable. It is proved that the numerical results of the model in the uniaxial compression test are consistent, which indicates that the VUMAT program can better solve the dependence of the grid size and has universality in engineering structure simulation.
4.3. Low Cycle Reciprocating Test Simulation of Reinforced Concrete
In order to thoroughly verify the validity of the model, the low cycle reciprocating test model of the reinforced concrete column was established. The model size and concrete column reinforcement are shown in Figure 5. The loading rules adopted in the low-cycle reciprocating test are shown in Figure 6. Concrete material parameters were simulated using the C30 standard of the VUMAT model and the C30 standard of the CDP model. To prevent stress concentration during loading, a reference point was established at the top center of the specimen, and coupling constraints were adopted with the loading surface. The embedded region is used to embed rebars into concrete. Since the explicit dynamic method is used for quasi-static analysis, it is necessary to reduce the influence of inertia force on the whole component. Therefore, by reducing the loading rate and the calculation time step, the acceleration can be approached to a small size to ignore the influence of inertial force.
Figure 7 shows the calculated hysteresis curves of the reinforced concrete column. Under the action of low cyclic reciprocating load, cracks, and other damage phenomena occur in concrete, and the damage will gradually accumulate, which will lead to the gradual reduction of the stiffness of concrete, that is, the phenomenon of stiffness degradation, which is best illustrated in the figure. Compared with the hysteresis curve obtained by the CDP model, the hysteresis curve of the concrete VUMAT subroutine is fuller and can better reflect the bidirectional bending coupling performance of reinforced concrete columns. However, compared with the actual test situation, the simulated hysteretic curve is relatively full and does not show the reinforcement slip well, especially in the late loading stage. Moreover, the reinforcement constitutive model adopted is the linear strengthening model, which has some errors with the constitutive relation of reinforcement in engineering practice, and the energy dissipation capacity is more potent. The simulation results show that the curve conforms to the concrete constitutive model, indicating that the compiled constitutive subroutine of the concrete VUMAT material is correct and can be used for the dynamic elastoplastic analysis of structures.
5. Discussion
In this paper, the elastoplastic damage constitutive model of concrete is established, and the failure state of concrete is analyzed by the secondary development. It is verified that the model can better express the stress state of concrete and effectively solve the problem of mesh-dependence in finite element calculation. Through a low cycle reciprocating simulation test, the model can better show the hysteretic behavior of concrete columns, which verifies the correctness of the model. Through this study, further research can be carried out in the following aspects:(1)As a kind of composite material, the inherent defects of concrete have the property of random distribution. The initial damage distribution and the subsequent damage development are inevitably characterized by randomness. Therefore, the random field of damage evolution of composite concrete can be established in the subsequent research, and the random damage constitutive model of concrete can be developed based on the deterministic model.(2)The concrete subroutine under the explicit algorithm can be used to analyze and calculate the structural dynamic response using the central difference method, which can avoid the nonconvergence problem in the strong nonlinear stage compared with the implicit algorithm. In the case of multiple strong nonlinear structure calculations, such as structural collapse simulation, the concrete constitutive model under the explicit algorithm is more advantageous.
6. Conclusion
In this paper, based on the VUAMT interface of the ABAQUS platform, the concrete constitutive model is developed by the explicit algorithm, and simulation tests verify the model. Given the subroutine’s influence on reinforced concrete’s hysteretic behavior, the reinforced concrete column’s low cycle reciprocating simulation test is established. The following conclusions are drawn:(1)Based on previous studies, the elastoplastic damage constitutive model of concrete is summarized, and the concrete constitutive equation based on the damage energy release rate is established, which can better describe the damage failure state of concrete.(2)The explicit algorithm is used for secondary development of the model, and the relevant calculation model is established for verification. Compared with the test results of the CDP model, this model can better describe the mechanical state of concrete under compression and tension.(3)Because of the mesh-dependence problem of the finite element model, several different mesh specifications were used for model verification. The results show that the model can improve the mesh-dependence in the finite element simulation.(4)Compared with the CDP model, the hysteretic curve of the VUMAT subroutine is fuller, which can better describe the hysteretic performance of the reinforced concrete model and be applied to the dynamic elastoplastic analysis of structures.
Data Availability
All data included in this study are available upon request from the corresponding author.
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.
Acknowledgments
The writers gratefully acknowledge the financial support of the National Natural Science Foundation of China (no. 52168072 and no. 51808467) and “Ten Thousand Talents” High-level Talent Support Project of Yunnan Province (2020). The corresponding author is Dewen Liu and the address is [email protected].
Supplementary Materials
Part of the VUMAT program code is provided in the supplementary file which mainly describes the calculation of damage variable. (Supplementary Materials)