Calibration of a New Concrete Damage Plasticity Theoretical Model Based on Experimental Parameters

The introduction of concrete damage plasticity material models has significantly improved the accuracy with which the concrete structural elements can be predicted in terms of their structural response. Research into this method's accuracy in analyzing complex concrete forms has been limited. A damage model combined with a plasticity model, based on continuum damage mechanics, is recommended for effectively predicting and simulating concrete behaviour. The damage parameters, such as compressive and tensile damages, can be defined to simulate concrete behavior in a damaged-plasticity model accurately. This research aims to propose an analytical model for assessing concrete compressive damage based on stiffness deterioration. The proposed method can determine the damage variables at the start of the loading process, and this variable continues to increase as the load progresses until complete failure. The results obtained using this method were assessed through previous studies, whereas three case studies for concrete specimens and reinforced concrete structural elements (columns and gable beams) were considered. Additionally, finite element models were also developed and verified. The results revealed good agreement in each case. Furthermore, the results show that the proposed method outperforms other methods in terms of damage prediction, particularly when damage is calculated using the stress ratio.


Introduction
Due to its durability and availability of constituents, concrete is one of the most utilized materials in construction [1][2]. Many models attempted to describe the concrete behavior under different loading types during the long usage life, such as elastic and plastic models. A perfect elastic or plastic model cannot capture the actual behaviour of concrete since they neglect the fact of concrete brittleness that appears as concrete cracking or crushing. Therefore, it was necessary to find new models that take into account the actual behavior of concrete. One modern model is the Concrete Damage Plasticity (CDP) model, in which concrete damage (cracking or crushing) is combined with the traditional plastic softening model.
Currently, the CDP model is one of the most popular constitutive models to simulate the behavior of concrete. This model is a modification of the Drucker-Prager strength hypothesis and was theoretically described by Lubliner et al. in 1989 and developed by Lee and Fenves in 1998 [3,4]. Inelastic deformation, stiffness and strength degradation, as well as hardening and softening of materials can all be modeled using continuum models based on the theory of plasticity and continuum damage mechanics. Damage-plasticity models are particularly useful for cohesive-frictional materials such as concrete or other cementitious materials, as they can accurately predict the material's behavior in both tension and compression. Plasticity and damage model combinations are usually based on isotropic hardening combined with isotropic or anisotropic damage variables. Anisotropic damage assumes a different strength and stiffness in different directions in concrete material. This type of coupling is complex and is not straightforward for structural analysis. Isotropic damage assumes that the strength and stiffness are degraded equally in different directions [5]. In this study, isotropic damage is considered. Figure 1 depicts the combination of plasticity and damage models-the loading branches are represented by thick lines, while the unloading/reloading branches are represented by thin lines. The nonlinear behavior of concrete can be displayed through either a plasticity or damage model. Nevertheless, none of these models alone is capable of representing amply real concrete behaviour. Plastic models (refer to Figure 1-a) might represent the observed deformation but do not capture the stiffness deterioration. On the other hand, as shown in Figure 1-b, the damage models can describe the stiffness degradation but are not suitable for capturing the deformation observed in experiments. Thus, the coupling between plasticity and damage models is widely accepted as essential in describing the nonlinear response of concrete, as shown in Figure 1-c [6].  Recently, scholars performed a numerical analysis for the seismic behaviour of retrofitted RC beam-column joints. Whereas, they utilized the damage plasticity model for concrete and the bilinear plasticity model for steel reinforcement [7]. Others studied the effect of viscosity parameter on numerical simulation of fire damaged concrete columns. The concrete damage plasticity model was adopted for the calculation of constitutive concrete material used in hybrid concrete columns [8], reinforced concrete columns, and the viscosity coefficient [9]. Fire-resistant concrete slabs with nonlinear layered shell elements were also studied [10]. The authors suggested a model for simulation of deflections and stress states and the concrete damage-plasticity which was utilized at elevated temperature to effectively model concrete response. The damaged plasticity models were adopted in the assessment of the structural response of gravity beams were done by Nabi and Rahim (2021) [11]. That was done numerically using common material and damage models. The model's numerical structural response was compared to the experimental response. The comparison revealed that numerical models accurately predict the structural response of gravity beams. The behavior of column base-plates with different configurations by finite element simulation was also investigated [12]. The numerical model was calibrated with an experimental test on six different base plate configurations. Concrete Damage Plasticity (CDP) was used in Abaqus to simulate concrete. Menh et al. used a concrete damage plasticity model to predict and simulate the behavior of high-strength concrete types C60, C80, and C110 under static and dynamic loading conditions [13]. The model was developed using data from previous studies to improve and develop the original CDP model in the ABAQUS software.
The results indicate that the CDP model developed in this study agrees reasonably well with experimental results for both static and dynamic loading conditions, with a high degree of reliability. Figure 1(c) shows that the unloaded part seems to be weakened because the initial stiffness appears to be degraded or damaged. The degradation of the initial stiffness is identified by damage variables (d), which can take values from zero for undamaged material to one for wholly damaged [14]. Damage evolution laws play an essential role, especially when these laws are imposed. Many researchers have suggested different approaches to define the damage parameters. Most of their rules were derived based on splitting damage into tensile and compressive parts, and each one is calculated separately. In some researches, the total damage was determined with some combination rules such as [5,15]. According to the general formula for computing damage in tension and compression [16]. In most of the available damage calculating methods, the variables require experimental calibration [17]. Other methodologies used empirical approaches with some hypotheses on stresses invariants or iterative techniques to fit experimental results [18][19]. In ABAQUS software, the stress-strain relations are governed by a single degradation variable, compressive damage variable (dc), or tensile damage variable (dt), as clarified in the next section. If the damage variables are not specified, the model behaves as a plasticity model [20].
The aim of the present study can be summarized in the following points, where the study provides valuable information about:  The differences between each of the damage, plastic, and CDP models.
 Propose a new method to compute the compressive damage parameters and verify the numerical results predicted by this method with those experimentally and numerically obtained.
 Influence of FE mesh size on numerical results in CDP model for different methods of damage variables calculating.

Damage Parameters
Concrete undergoes a series of damage stages, ranging from micro-cracking to complete failure as it is subjected to more and more stress. In order to accurately predict the residual capacity of damaged structures to withstand additional loads, this process must be mathematically simulated. Over the course of history, numerous damage models have been put forth. Generally speaking, none of these can accurately predict the residual strength of damaged members. The dc, dt are the key parameters which were proved by other scholars as they represent the degradation in stiffness of concrete. The stress-strain relationship used in the CDP model adopted by Abaqus can written in the following formulation [21]: where, is the nominal (Cauchy) stress tensor, d is the Stiffness degradation variable (It is also defined as a damage parameter), Eo and E are Elastic stiffness of undamaged and damaged material, respectively. and are total and plastic strain tensors. On the other hand, the effective stress ̅ is defined as: The term (1-d) indicates the ratio between the effective load carrying area and the total area, i.e., the overall area minus the damaged area to the overall cross-sectional area. When no damage occurred, i.e., d = 0, the effective stress ( ̅) equals the nominal stress σ [22]. The effective stressed cross-sectional area of the model is obtained using the method proposed by Al-Rub and Kim [5]. Generally, the damage parameters (dc and dt) represent the degradation of elastic stiffness, as illustrated in Figure 2 under both compressive and tensile loading. The unloading response is weakened if the concrete specimen is unloaded from any point of the softening branch of the stress-strain curve, i.e., the elastic stiffness of concrete seems to be damaged. Two damage variables characterize this process. These variables range from zero for an undamaged material to one for total damage, representing a total loss of strength [21].  The damage parameters play an essential role to describe the actual behavior of the material. Some researchers define compressive damage as the ratio between the inelastic strain to the total strain ( ⁄ ). Similarly, the tensile damage variable is the ratio between cracking strain to total strain ( ⁄ ) [20]. At the same time, Kmiecik (2011) obtained damage values by calculating the ratio of the stress for the descending part of the stress-strain curve to the ultimate strength of concrete ( ⁄ ) [4]. In comparison, Arya used the two above definitions to determine the damage parameters [23]. Also, some researchers have assumed the ratio of plastic strain to inelastic strain to find the damage parameters [24,25]. Besides, some researchers did not specify any damage parameters with input data in the Abaqus program, which leads to the concrete model behaving like a plasticity model, not as a damage plasticity model [26].
Evaluation of damage parameters according to the ratio of inelastic or cracking strain to total strain is not accepted programmatically in Abaqus program because this process leads to the existence of non-plastic strains, which is apparent in the plastic strain equation (Equation 6), i.e., ( ) equal to zero by replacing dc by ( / ) Ratio. However, it is theoretically adequate for simulation of concrete as a damage model, as shown in Figure 3 case (a), i.e., the total strain was recovered immediately after the load removal. On the other hand, the main shortcoming of the second definition concerns the stress ratios ( / ) from which the damage parameters are determined after the maximum stress is reached. However, cracking is started with a stress level less than the ultimate stress. This response shown in Figure 3 case (b), which indicates perfect plastic behaviour up to the maximum strength of concrete; after that, the damage grows to complete failure. Finally, case (c) plots as a horizontal line to represent the modelling of concrete with a pure plasticity model, where the plastic strain equals the inelastic strain, which is increasing with zero compressive damage value. In the proposed method, the compressive damage increases with the increasing of the plastic strain, gradually, up to complete damage. This response obtains when the damage finds by the stiffness ratio as shown in Figure 3 case (d). This method is described in the next section.

Proposed Method for Calculating dc
The concrete response under uniaxial compression and cyclic loading must be known for obtaining the compressive damage parameters. The flowchart shown in Figure 4 illustrates the research methodology and in the following subsections, a complete description is presented.

Concrete Stress-Strain Relation for Uniaxial Compression
In the CDP model, the definition of concrete stress-strain behavior under compression needs to know stresses, and inelastic strain corresponds to stress values and compressive damage (dc) in tabular format. In general. The inelastic strain (or crushing strain) ( ), which occurs due to the opening of micro-cracks then merging to macro cracks, is determines by subtracting the elastic strain ( ) from the total strain ( ). Finally, the inelastic strain is converted to plastic strain ( ) according to damage parameter dc, as follows: According to Equation 1 and Equation 2, compressive damage parameter can be written as follows: Then, the calculation of dc needs to estimate the value of secant stiffness (E) at any point in the compression stressstrain curve, i.e., stiffness of damaged material. The value of (E), at any point, could be found by defining the cyclic response of concrete. In this study, two constitutive relationships were used to describe the compression behaviour of concrete. Those relations are presented in the next section.

Concrete Stress-Strain Relation under Cyclic Loading
Generally, the strength and deformation of concrete exposed to cyclic loading differ significantly from those subject to monotonic loading. Figure 5(a) shows the stress versus strain curves of concrete subjected to cyclic loading. Figure  5(b) shows the response of the concrete under one loading cycle and unloading cycle, where (oa) and (ab) are the loading and unloading curves, respectively. When stresses dropped to zero after point (a), most of the overall strain corresponding could be recovered during unloading instantaneously ( ′ ). Besides, a small portion of strain referred to as elastic hysteresis can also be recovered after some time ( ′′ ). The residual strain is called the unrecovered strain ( ′ ) [27].  There are numerous constitutive models for describing the behaviour of concrete under cyclic loadings [28][29][30]. This study used the constitutive relationship proposed by Sima et al. [30] in this model. The secant stiffness could be shown in the following expression: where ε' is the strain at peak stress. Using Equation 9 requires Excel tabular format calculations. The accuracy of the proposed approach shows in the following sections.

Verification of the Proposed Method
This section describes the modelling of different structural specimens with different mechanical properties to check the accuracy of the proposed equation. The FE results obtained using Equation 9 are compared with those recorded by experimental and numerical work. Briefly, these specimens were listed below:  Standard cylinder with dimension 150×300 mm exposed to compression.
The dimensions and details of RC specimens are shown in Figure 6 and Table 1. Also, the characteristics of the FE models are displayed in Table 2.   It should be noted that the tensile behaviour of concrete was modelled by using the constitutive relationship suggested by Belarbi & Hsu (1994) [33] depicted in Figure 6. This model expresses mathematically as: σ t = E cm ε t for ε t ≤ ε cr (10) σ t = f cr (ε cr ε t ⁄ ) o.4 for ε t > ε cr (11) where; is the average tensile strain at which the concrete start cracking and fcr is the tensile strength of concrete at . Figure 7 shows a linear elastic relation up to ultimate tensile strength. When the concrete reaches the cracking state, a phenomenon called tension stiffening will occur. This phenomenon should be considered when the tension zone is reinforced with steel or other reinforcement material. This state results due to the bond between concrete and steel bars between cracks and the aggregate interlocking. The tension stiffening phenomenon is responsible for giving the plasticity nature to the behavior of concrete under tension. Accordingly, the stress in reinforced tensioned zones does not descend rapidly but gradually [4]. Consequently, the concrete tensile damage (dt) calculates as Equation 12:

Case Study 1 -Concrete Cylinders
The problem of uniaxial compression of concrete cylindrical specimens with 150 mm diameter and 300 mm height was modelled in the Abaqus program. The cylinder is discretized with a uniform mesh of three-dimensional 8-node hexahedron solid elements (C3D8R). Figure 8 displays four schemes of the cylinder discretized with a coarse mesh (M50, 78 elements), a medium-mesh (M35, 198 elements), fine mesh (M25, 576 elements), and very-fine mesh (M15, 2360 elements), respectively. This analysis aims to establish the allegedly low sensitivity to mesh size on the behaviour of concrete cylinder specimens under uniaxial compression. In which, load-deformation control was applied with an axial velocity of -80 mm/sec at the top surface of the cylinder. The cylinder's top and bottom surface displacements could move freely except for the axial displacement at the top surface. Besides, all rotations about principal axes were allowed. The analysis was carried out using the explicit dynamic method with a one-second step time. Plain concrete was assumed to have a compressive strength of 30 MPa for the analysis.  [34] has been adopted in Abaqus to represent the compressive behavior of concrete. This relation was taken as a reference model. The other three relationships were based on the compressive damage calculation method. Obviously, by implementing the plasticity model, the compressive and tensile damage will not be specified. In Figure 9, the numerical stress-strain curves were plotted together with the  t  t f cr  cr reference model obtained with Euro-code. This figure shows the accepted accuracy of the proposed method compared with other outcomes which indicates a superior ability of this model to capture the deterioration in stiffness. Additionally, the effect of mesh size is limited, however, all predicted data per the proposed approach are conservative with certain margins.

Figure 9. Compressive behavior of concrete with different element size
The mesh sensitivity on stress-strain curves is also considered as presented in Figure 10. In comparing coarse and medium meshes with the reference model, the maximum differences in ultimate strength were 2.5 and 12% when the compressive damage is calculated by the proposed method and method of stress ratio and as show in Figure 10-aand Figure 10-b, respectively. Nonetheless, this percentage increased to 13.5% when the plasticity method modulates the concrete, as shown in Figure 10-c.

Case Study 2 -RC Columns
Two RC columns tested by [31] were modelled via the Abaqus program. These columns were tested under biaxial loaded and pinned end conditions. Each column has a slenderness ratio equal to 30, the steel reinforcement details and dimensions as shown in Figure 6. Numerically, both column ends are modelled as pinned-pinned conditions using two reference points (RP) fixed at both ends of columns. The coordinate of these points was specified based on the point of load application, i.e., the eccentricity of load (e). The load was applied with an axial velocity of -80 mm/sec at the top end with biaxial eccentricity equal to 30 mm and 50 mm for C1 and C23, respectively. Additionally, the columns were discretized with a uniform mesh C3D8R, as shown in Figure 11-a. In the case of C1 column, the cracking occurs at the upper half, as displayed in Figure 11-b. In the C23 column, the cracks appeared at the tensile side at the mid-height region, as shown in Figure 11-c. Accordingly, the failure mode in C23 column is a more frequent flexural failure than C1 column. The compressive damage in concrete elements in C23 column shows in Figure 11-d. The FE load-deflection curves were compared with those obtained from the test results, as displayed in Figures 12-a and 12

Case Study 3 -RC Gable Beams
Case three shows a study based on the experimental test of simply supported gable roof beams done by [32]. The original paper provides the influence of openings on the flexural behaviour of RC beams. Six RC beams were tested, one solid (GB) and the other 5 with openings of different shapes. The test was executed under monotonic static load at mid-span until complete failure, i.e., concentrated load. In this study, two RC beams were selected to examine the accuracy of the proposed method with the test results. The experimental arrangement of selected beams and meshed models are displayed in Figures 13-a and 13 those numerically predicted, as shown in Figures 13-c and 13-d. Besides, the experimental load-mid-deflection curves were compared with those gained by FE analysis with plasticity and damage-plasticity model as shown in Figure 14.
Generally, it can be seen that the load-deflection curves predicted by the proposed method have good agreement with experimental results as compared with the plasticity model. When compared to a plasticity model, it can be seen that load-deflection curves predicted by the proposed method are more accurate. Solid and hollow beam ultimate loading predicted by the proposed model differed slightly from the results of the plasticity model by 3% and 1%, indicating the model's ability to accurately predict data. When compared to the experimental results, both solid and hollow beams were slightly overestimated in strength by 4 and 7%, respectively, in the proposed approach.

Conclusions
It is evident from the cases examined in this study that the proposed method is effective in determining compressive damage parameters for a variety of reasons, including:  The damage begins as soon as the model is loaded, and it doesn't stop until it's completely destroyed.
 Comparing the results obtained by the plasticity model and those obtained by stress ratio, when the damage parameters were calculated using the proposed method, the ultimate strength was only marginally influenced by the element size that was used in the modeling process in comparison to the results obtained by either model. When the compressive damage is calculated using the proposed method or the stress ratio method, the maximum differences in ultimate strength between coarse and medium meshes were 2.5 percent and 12 percent, respectively, when the compressive damage is calculated using the proposed method or the stress ratio method. However, when the concrete is modulated using the plasticity method, this percentage rises to 13.5 percent, indicating a significant increase.
 The stress ratio's evaluation of compressive damage may be flawed due to some factors, such as the fact that the damage begins after the maximum load has been reached. Cracks are also formed simultaneously as the ultimate stress in the material is reached but before the ultimate stress is reached in the material. On the other hand, the coupling between the plasticity and damage models, on the other hand, is not observed in the ascending portion of the loading curve. The damage begins as soon as the model is loaded, and it doesn't stop until it's completely destroyed.

Data Availability Statement
The data presented in this study are available in the article.