Calibration of Concrete Damaged Plasticity Model parameters for shear walls Calibração da relação constitutiva Concrete Damaged Plasticity para aplicação em pilar-parede

Reinforced concrete structures are relatively complex to analyze, with nonlinear effects like cracking, crushing, steel yielding, aggregate interlock, dowel effect, concrete-rebar interaction and so on. The concrete damaged plasticity CDP model is a consolidated smeared-crack model which accounts for multiaxial behavior with good agreement to experimental results. One particular relevant application which benefits greatly from such feature is the shear wall, as shear stress significantly influences its overall behavior, therefore multiaxial constitutive models and three-dimensional finite elements usage consist in a fitting modeling approach. Reinforced concrete shear walls are structures especially useful for lateral force-resisting systems, as they provide ductility, stiffness and strength. Albeit CDP is widely applied, its parameters are not consensus in the literature, which represents a relevant research gap. The present work considers and compares CDP parameters from relevant literature, in order to calibrate those parameters for the case of reinforced concrete shear walls. To this purpose, four wall experiments related in the bibliography are modeled using solid finite elements for concrete and trusses for rebars using commercial package ABAQUS. All walls are flexure-controlled with aspect ratio greater than 2.0. By varying those parameters and comparing obtained force vs. displacement curves and interesting values attained, like yield lateral force and displacement, stiffness and maximum lateral force, it is settled a set of parameters with acceptable response focusing in the post-peak response based on the lower estimated error of displacement capacity. Those parameters agree reasonably with literature, although it is possible that obtained calibration is restricted to flexure controlled shear walls scope. It is possible that usage of trusses to represent reinforcement does not consider dowel effect, so a suggestion for future studies is to change trusses for elements with transverse stiffness, like beams or solids.


INTRODUCTION
Reinforced concrete modeling has to account for its specific non-linear behavior, such as cracking and crushing. Furthermore, in multiaxial states of stress, material response changes significantly, and it is required an appropriate formulation to account its response appropriately [1,2].
One consolidated modeling approach is the Concrete Damaged Plasticity Model. Derived from the Drucker-Prager criteria, it consists in a multiaxial constitutive model with direct application for structural reinforced concrete. Although very robust and extensively used for analyses, there is no consensus among which parameters should be used in each case by comparison of its applications. This consists in an evident demand for better calibration considering the relatively discrepancy in the parameters suggested in the references.
A case in which a multi-axial state of stress is essential is for shear walls. Reinforced concrete shear walls perform very well as lateral-resisting systems, as they provide displacement capacity, stiffness and lateral strength. Hence, they might assure safety and durability for high-rise or seismic vulnerable buildings [3][4][5].
There are relatively direct models for shear wall analysis which are inexpensive in terms of computing cost and still provide an accurate response while only applying linear elements, e.g. the Wide-Column Model which proposes usage of beam elements that can be considered linear with secant behavior or nonlinear based on a fiber section discretization [6]; the Multiple Vertical Line Element Model can consider non-linear and cyclic material behavior and still accounts for relevant responses, like axial-flexure interaction and neutralaxis migration by using only trusses and a horizontal spring. [4,7].
There are also more robust and complex models. The Finite Element Model FEM opens the possibility for continuum bi-dimensional and three-dimensional elements which account shear displacement and shearflexure interaction inherently to the element formulation, making it very convenient in the assessment of structural walls. Nevertheless, setting up and processing FE models might be very expensive, not only for processing demands but also to obtain the appropriate material parameters for the model as well as in user's time and knowledge to assemble it [4,8].
However, using a more robust formulation is highly beneficial in some cases, like nuclear engineering, special structures, research applications or to avoid experimental methods. Therefore, using the FEM is a good option for reinforced concrete shear wall analysis, allowing for a more accurate response in exchange for more time spent in the modeling process if compared to simplified modeling approaches [8].
The shear displacement is relevant for shear walls due to the fact that they are short elements and therefore more susceptible to shear than linear elements like beams and columns, making it important to account for the shear induced displacement which is negligible in linear elements. Deep beams have also this property and could be well analyzed with FEM.
The constitutive model Concrete Damaged Plasticity CDP is multiaxial, it being a pre-requisite for usage with shell or solid elements, and widely applied to model reinforced and prestressed concrete structures in general. It is available in the software package ABAQUS, in which the entrance parameters are the material's dilation angle, eccentricity, the ratio of biaxial to the uniaxial compressive yield stresses, the ratio of the second stress invariant on the tensile meridian to that on the compressive meridian and the viscosity parameter. It is also required to assume an uniaxial compressive and tensile behaviors and it is recommended to model the damage variables. Among correlate literature, there are significant discrepancies between values found for those quantities.
This paper aims to calibrate the parameters for reinforced concrete shear walls modeling with the Concrete Damaged Plasticity model based on experimental evidences. Although the calibration is focused on reinforced concrete shear walls, it can fit other structural elements with similar properties, like deep beams. This represents a significant contribution considering the wide use of numerical analyzes to estimate diverse structural behaviors [2,3,5,6,7,9,10], so appropriate material simulation is a subject of interest of the research community.

CONSTITUTIVE MODELS
Described in the papers by LUBLINER et al. [11] and extended by LEE and FENVES [12], the Concrete Damaged Plasticity CDP constitutive model is mathematically smooth, continuous and smeared-crack model. It can be used to analyze reinforced concrete behavior considering plasticity, softening and damage, applicable also for cyclic analyses. Its yield function adopts a Drucker-Prager hyperbolic function shown in eq. (1) with eq. (2), (3) and (4) defining its dimensionless parameters.
In which, ̅ : is the hydrostatic pressure; ̅ : is the Mises equivalent stress; ̅ is the stress deviator; ̂ is the maximum principal stress; ̅ and ̅ are the cohesion stresses for compression and tension; 〈 〉 is the Macauley bracket returns zero if x<0, else it returns x.
(σ b /σ c ) is the ratio of the biaxial to the uniaxial compressive strength; K c defines the failure surface in deviatoric plane, which is normal to the hydrostatic axis.
The CDP considers nonassociated potential plastic flow. The potential plastic flow G is a Drucker-Prager hyperbolic function defined as equation (5).
where: σ t0 is the uniaxial tensile stress; is the dilation angle; is the eccentricity of the plastic flow and defines the rate at which the function approach the asymptote.
The yield surface in a deviatoric plane, which is orthogonal to the hydrostatic axis, and also the hyperbolic flow potential is illustrated at Figure 1, exhibiting the influence of the K c parameter, as well as the geometric interpretation of the dilation angle and the eccentricity . The yield surface for concrete subject to biaxial stresses is shown in Figure 2. It also illustrates how the yield function (equation 1) depends on the sign of the maximum principal stress: if there is tension at the stress tensor, thus the maximum stress is positive and the term with is not null. Likewise, the parameter with only appears in triaxial compression [12,14].  [13].
The parameters required for the CDP model in ABAQUS are the dilation angle , eccentricity , the ratio of biaxial to the uniaxial compressive yield stresses (σ b /σ c ), the ratio of the second stress invariant on the tensile meridian to that on the compressive meridian K c and the viscosity parameter μ.
The dilation angle is measured in the p-q plane in high stresses, as shown in Figure 1. It is predominantly assumed as 30º [15,16,17], but there is some divergence between values recommended [8].
The eccentricity is a small positive number defining the rate in which the hyperbolic flux approaches its asymptote. The default value in the literature is 0.1, almost undisputed [18,19].
The ratio of biaxial to the uniaxial compressive yield stresses (σ b /σ c ) was experimentally determined as 1.16 by KUPFER, HOLSDORF and RUSCH [1], extensively assumed as such. One particular work [20] presents this parameter as a function of concrete compressive strength f' c calibrated from experiments reported in the literature.
The ratio of the second stress invariant on the tensile meridian to that on the compressive meridian K c must be greater than 0.5 and limited to 1.0, with the standard value being 0.667 [5,21].
As for the viscosity parameter μ, by assuming a sufficiently small value, it may improve obtaining convergence while preserving accurate results, especially for post-peak softening behavior. It is recommended as null [13].
Summarizing all the needed parameters for CDP modeling, the yield function is determined by (σ b /σ c ) and K c , while the flow potential is defined by and . Additionally, it can be considered a viscosity parameter too. It was found in correlate literature information about values adopted in CDP, and they are shown in Table 1. Furthermore, a stress-strain uniaxial relation is required for compressive and tensile behavior. Herein the CHANG and MANDER [33] envelope is applied. Their work recapitulates numerous reports about concrete stress-strain behavior. Some important constants are calculated from concrete compressive strength f' c , like peak strain, initial elasticity modulus, shape factor and tensile strength.
As for damage, it is calculated by equation (6), and it is considered null up to peak point [18]. Where: d is damage in compression or tension; σ is stress on softening; is the strength.
The CHANG and MANDER [33] envelope establish a stress-strain relation under uniaxial compression or tension and estimates strain at peak, modulus of elasticity, tension strength overall shape all from the concrete strength cylinder under compression f' c . One of its benefits is that it only requires one experimental data to define the tension and compression shape. Figure 3 illustrates the expected behavior.
As for rebars, it is assumed a simple bi-linear stress-strain response considering isotropic hardening, suitable for simple elements like trusses representing rebars.

REFERENCE EXPERIMENTS
For the calibration, it is necessary to compare results from the simulation to those of an experiment. Therefore, this study considers experimental evidence presented by Ghorbani-Renani et al [34] and DAZIO, BEY-ER and BACHMANN [35].
GHORBANI-RENANI et al. [34] describe experimental analysis of four walls, of which two are considered here, walls A1M and A2C. Both have identical geometric and mechanical properties, however wall A1M is subjected to monotonic load and A2C, to cyclic load. Herein, those walls are modeled as one, addressed as A1M/A2C. DAZIO, BEYER and BACHMANN [35] assess mainly the influence of the confining reinforcement, presenting the results of six reinforced concrete walls under a cyclic load protocol. In the present work, it is taken walls WSH2, WSH4 and WSH5 in consideration. They share quite similar geometry and properties.
The cross-sections for walls A1M/A2C, WSH2, WSH4 and WSH5 are presented in Figure 4. The selected walls' key characteristic are presented in Table 2, in terms of concrete compressive strength f' c , reinforcement ratio ρ s , transversal reinforcement ratio ρ sw , height h w , width l w , depth h and axial force ratio ν. All except one of those investigations apply cyclic loads. Walls A1M and A2C show similar (force vs. displacement) curves for monotonic and cyclic behavior up to the peak point of cyclic chart, which may be observed on comparing monotonic and backbone curves.
Extending this assumption to walls WSH2, WSH4 and WSH5 is reasonable, considering they share similar geometry and aspect ratio above 2, defining them as flexure-controlled, and the observed behavior corroborates by displaying good increase in lateral force capacity and ductility. Shear controlled walls instead exhibit continuous loss of lateral load capacity after plasticity.
About results evaluation, it must be taken variables that are both comparable and important wall responses. Thus chosen variables are yield force F y , yield displacement d y , stiffness as ratio of difference of force to difference of displacement between points of 20% and 60% of maximum force K 0.2-0.6 and maximum force F max . Yield point criterion must be equal for all curves, so a yield criterion by Park [38] is considered, applicable to force vs. displacement curves. Quickly reviewing it, it may be applied in three steps; first, a horizontal line which intercepts the curve at maximum lateral force point F max ; then another line passing through origin and 75% F max point reaches the horizontal line; from the interception point, a vertical line is projected to the curve and the point of encounter is considered the yield point with coordinates (d y , F y ).

METHODOLOGY
The Finite Element Method tool used was the ABAQUS commercial package, which offers the CDP formulation and a suitable element library containing the required tools for the simulations.
The steel rebars are modeled as truss elements (T3D2) embedded in the concrete solid elements (C3D4). The modeling approach considering trusses embedded on solids is quite common for reinforced concrete analysis [14,16,17,21,23]. Figure 5 illustrates applied finite elements. The embedment simulates perfect bond between concrete and reinforcement, this is a reasonable approach considering no bond-slip occur. As for the chosen rebar element, unlike beams or solid elements, trusses ignore the transversal behavior that is relevant to consider the dowel effect, an important subject for short elements.
Using defined finite elements and constitutive models, it is possible to simulate wall A1M/A2C behavior using a numerical model. By assembling solids and trusses following the general geometry described by GHORBANI-RENANI et al. [34], a three-dimensional structure is obtained and shown at Figure 6.  For the concrete, it is considered the CDP model, initially assuming standard parameters , , (σ b /σ c ), K c and μ as 38º, 0.1, 1.16, 0.667 and 1e-3. Uniaxial stress-strain relationship is adopted according to CHANG and MANDER [33] proposals for compressive and tensile behavior. Damage is calculated in equation 6.  For the mesh, it was applied a remeshing rule based on the element energy available at the used software and endorsed by appropriate literature [39]. After each remeshing iteration, its results are compared to that of the previous analysis and its mesh accepted as suitable if there is no significant discrepancy, which was considered 2% of the lateral force for the same displacement until softening. From the point when softening occurs, the results are not considered.

RESULTS AND DISCUSSION
Initially, the viscosity parameter is analyzed as it must be small enough for the model to attain convergence. It is tested in compressive behavior an isostatic cube assembly with only uniaxial stresses with concrete which considered as the constitutive model CDP with present work's standard parameters. Based in values for viscosity parameter reported by literature as zero to 0.005, for present test it is considered equals to 1e-5, 1e-4 and 1e-5. The obtained stress-strain for all integration points is presented at Figure 9. As observed, even for a literature recommended value, a relatively discordant stress-strain curve was obtained. To numerically evaluate the results, observed peak stresses are compared to the analytical value, then obtaining the error for each model.
Arbitrarily assuming a tolerable error of 1% in the peak stress and a linear relation between error and viscosity parameter, it can be estimated that the viscosity parameter which leads to 1% error equals 2.15e-5. By carrying out the same assessment in tensile behavior, it is obtained as maximum tolerable viscosity parameter 5.8e-5. Thus, herein is adopted as maximum viscosity parameter 1.0e-5, in order to avoid further interference from excessive viscosity parameter.
Setting default μ as 1.0e-5 for the following simulation, it is still required to evaluate the influence of parameters , , (σ b /σ c ) and K c . Simulating wall A1M/A2C to default values and varying only one of them, the curves for force vs. displacement at the top obtained are exhibited in Figure 10.
By applying similar changes to CDP parameters of wall WSH2, the simulation leads to results as shown in Figure 11.  Initial ascending branch is almost independent of changes in CDP parameters in both simulations, as difference between modeling results are almost unchanged from each other in that part. Then, there is a stage in which stiffness becomes very small, observed in charts as an almost horizontal line, and then CDPs parameters influence are evident as curves diverge. Comparing modeling and experimental results, initial ascending branches agree very well only until subtle loss of stiffness, which is due to rebar yielding.
As Figure 10 indicates, dilation angle parameter value significantly modifies observed behavior; by using greater values for dilation angle, it is observed more displacement capacity and strength; on the other hand, assigning smaller dilation angle leads to lesser ductility and strength.
Regarding the eccentricity parameter, its reduction leads to small decrease in displacement capacity for wall A1M/A2C while wall WSH2 was quite indifferent to it. Considering that even those aggressive alterations did not lead to great changes in response and the literature is almost unanimous, the standard eccentricity is kept as 0.1.
Parameter ratio of biaxial to the uniaxial compressive yield stresses (σ b /σ c ) is also analyzed. If it is increased, it is observed small increase at displacement capacity shows sensible increase and vice-versa. Like the eccentricity parameter, as its alteration leads to minor results change and it is almost consensus in observed literature, (σ b /σ c ) is kept as 1.16.
As for variable K c , its alteration leads to significant response changes: by decreasing K c , observed ductility and lateral force increase significantly for both walls. It is observed great displacement capacity on both walls comparing to any simulation curve presented at Figures 10 and 11, except for simulations considering K c as 0.625 or below for WSH2; while for wall A1M/A1C, neither simulation could properly predict the experimental displacement capacity.
This suggests both parameters and K c should be calibrated together. So by varying those parameters on walls WSH4 and WSH5, obtained relations force vs. displacement at load application height. About the initial branch, once again it is observed great agreement between numerical and experimental responses, pointing good model performance; and once again CDP parameters do not impact initial behavior observed. As for ductility, both simulations required reduction at K c parameter in order to take into account the top displacement capacity suggested by experimental evidence. Wall WSH4 and WSH5 only exhibited satisfying ductility for K c about 0.52 and 0.58.
Another issue is about concrete confining reinforcement influence on simulated wall behavior. This is of paramount importance as its one of the advantages of three-dimensional modeling, because it is possible to assembly a construct with every single rebar detailed. Thus concrete confinement might be implicitly taken into account, as CDP constitutive model considers multiaxial behavior.
This way it is possible to consider material mechanical properties and geometry completely individually. One and two-dimensional finite elements modeling approaches cannot consider confining reinforcement directly, as it is orthogonal to the element line or shell section. Even not being inherent of these, concrete confining effects are very important and there are models for its consideration, like the widely applied Mander, PRIESTLEY and PARK [40] model, which propose changes in constitutive model attributed to specific elements.
To confirm whether confining reinforcement presence in models modify their results, only confining reinforcement is removed from simulation. Also, analyzed parameters are varied to check their impact at the model response without confining reinforcement. Results are presented in Figure 13 for walls A1M/A2C and for wall WSH2. As initial ascending branches do not significantly differ from one another, the origin point is only presented in the first chart for each wall, i.e. Figure 11 for walls A1M/A2C and WSH2, Figure 12 for walls WSH4 and WSH5.
Kc=0.5 Kc=0.5835 Curves from Figure 13 confirm that confining reinforcement elements impact its model results, as walls A1M/A2C shows evident loss of ductility if no boundary stirrups are considered. Furthermore this results shows that parameter lose its impact at the force vs. displacement obtained, which occurs at both analysis shown in Figure 13 and also for wall WSH4, as it has only open-loop stirrups as confining reinforcement, which are kept in the model and leads to results shown at Figure 12. However, decreases in parameter K c still increases observed ductility, thus its change is necessary to account for it in wall WSH4 since it is insensible to any changes.
For wall WSH2, assuming K c around 0.5835 presented a curve with similar displacement capacity as experimental evidence does, as seen in Figure 11. Testing K c equals to 0.57, and varying in concrete constitutive models, obtained force vs. displacement curves are presented in Figure 14 for walls A1M/A2C and WSH2; and Figure 15 for WSH4 and WSH5.  Considering K c as 0.57, results demonstrates overestimation in displacement capacity which leads to unsafe predictions, therefore it is also of interest to analyze K c as 0.58, which leads to the results shown in Figure 16 for walls A1M/A2C and WSH2; and Figure 17 for WSH4 and WSH5.    As graphic examination suggests, assuming K c as 0.58 results to wall WSH5 showing good displacement capacity estimative, while walls A1M/A2C and WSH2 have theirs underestimated and WSH4 ductility is greatly underestimated due to its lack of confining reinforcement.
Finally, for parameter assessment, instead of graphic examination, it is calculated dependent variables for acute adjustment. Values are extracted from all curves presented in Figures 16 and 17, except for wall WSH4. Variables F y , d y , K 0.2-0.6 , and F max as designated dependent variables; it is used the yield point definition criterion recommended by Park [38] for reinforced concrete structures. Comparing their results to that of the mean of the experimental ones, an error can be estimated. Values obtained for all these are presented in Table 3. Wall WSH4 is excluded from this analysis, as it lacks confining reinforcement, which is not expected to occur on everyday engineering praxis and changes significantly wall response.   One way to choose appropriate is summing the errors from each dilation angle value an establishing one relation considering the error as a dilation angle function. Admitting a polynomial regression leads to equation (7): Calculating the smallest value for calculated error, it returns a dilation angle of 46.4º. Reviewing succinctly, the calibration brought up the values for , , (σ b /σ c ) and K c , respectively, 46,4º, 0.1, 1.16 and 0.58. GENIKOMSOU and POLAK [14] report that the dilation angle should reach between 31º and 42º, assuming 40º affirming that the results are almost unchanged for dilation angles from 38º up to 42º. Their paper assesses shear punching on concrete slabs and shows more sensibility to the dilation angle value when compared to present observations. That paper considers Kc as 0.667, as herein is 0.58, and it was observed that for smaller values for Kc makes lesser sensible to the dilation angle. Also, GENIKOMSOU and POLAK [14] apply Kc as 1.0, but do not obtain a brittle curve, while for present models that occur, as seen at Figures 10  and 11.
EARIJ et al. [8] define dilation angle as 40º for their simulations and shows no great change in the results when using 30º, 40º or 50º. EARIJ et al. [8] shows better results using beams elements in detriment to truss elements. Therefore, the present calibration might prove inaccurate while using beam elements as rebars. Also using beams as rebars might not need an aggressive change in standard parameters like present models did.
Although present simulation did not achieve good agreement with experimental behavior by using standard CDP parameters, numerous works do [5,16,17,21,22,23]; most of them varies only the dilation angle, as shown in Table 1.

SUMMARY
The objective of this work was to obtain a calibration for the multiaxial concrete model CDP parameters based on experimental results from reinforced concrete shear walls. Concrete is simulated using solid elements C3D4 and steel rebars as T3D2 trusses. The considered parameters are dilation angle , eccentricity , the ratio of biaxial to the uniaxial compressive yield stresses (σ b /σ c ), the ratio of the second stress invariant on the tensile meridian to that on the compressive meridian K c and the viscosity parameter μ.
The viscosity must be small so it does not impact the reliability of the results while setting it to zero makes numerical convergence more difficult. Simulations show that increasing viscosity values leads to an unwanted growth in uniaxial stress as well. Admitting a reasonable 1% increase in peak-stress for uniaxial stress-strain behavior, the viscosity parameters should not be greater than 2.8e-5 in this study.
As for parameters , , (σ b /σ c ) and K c , comparing curves force vs. displacement from numerical modeling and experimental evidences of reinforce concrete shear walls, some important remarks are: 1. The initial ascending branches of experimental and numerical analyses agree very well and observed behavior is insensible to changes in any of these CDP parameters. Therefore, the modeling approach was able to accurately simulate cracked stiffness; 2. It is in the plateau of the curves that CDP parameters have great impact in overall behavior. Using parameters recommended by the literature leads to a response with loss of bearing strength at earlier displacements when compared to experimental data. 3. It was possible to change plateau behavior of numerical analyses by varying CDP parameters.
Still, only parameters and K c impacted it; the influence of parameters and (σ b /σ c ) were negligible. By increasing or decreasing K c obtained force vs. displacement plateau presents more displacement capacity at the same level of strength. 4. Finally, comparing numerical and experimental evidences, one way to simulate RC shear walls and similar structures using CDP material model and account for ductility, parameters , , (σ b /σ c ), K c and μ are recommended equal to 46,4º, 0.1, 1.16, 0.58 and 1e-5. With or without adjustments, FEM with CDP proved a powerful nonlinear analysis representing a good correlation to experimental results, especially for the initial ascending branch.
It is possible that considering the dowel effect would improve observed plateau response, as this work utilizes trusses for rebar modeling and therefore did not account for it as trusses lack transverse stiffness. That could be corrected by changing the element type used to model rebars. This would be possible using linear beam or solid elements, which could be costly for processing demands; still, further research is then required and is underway.