“Structural engineering is the art of molding materials we don't wholly understand, into shapes we can't fully analyze, so as to withstand forces we can't really assess, in such a way that the community at large has no reason to suspect the extent of our ignorance”

James E. Amrhein

## Resumen

Primera parte. En esta parte se desarrolla un nuevo modelo basado en la mecánica del medio continuo para simular el comportamiento monotónico y cíclico de estructuras de hormigón armado. El modelo desarrollado combina una nueva metodología para el cálculo de las variables del daño en el Modelo de Daño Plástico del Hormigón “CPDM”, y un nuevo enfoque para integrar el CPDM con un modelo de interface de 3-D desarrollado en otra investigación. También se presenta un nuevo esquema para implementar la interfaz del modelo en un modelo FEM continuo de regiones con armaduras que se cruzan en varias direcciones. La precisión, la fiabilidad y la insensibilidad a la malla del modelo propuesto se verifican simulando varias pruebas incrementales y cíclicas; los resultados obtenidos se comparan con experimentales, lográndose un ajuste satisfactorio.

## Summary

Under severe seismic excitation, structural behavior of buildings and other constructions is highly complex. It involves, among other issues, soil-structure interaction, large strains and displacements, damage, plasticity, and near-collapse behavior. Moreover, in reinforced concrete structures, there are several coupled degradation and failure modes: cracking, crushing and spalling of concrete, yielding and pull-out of tensioned reinforcement, yielding and buckling of compressed reinforcement. Furthermore, another circumstance makes the situation more alarming: given the increasing awareness and concern on the huge worldwide seismic risk, earthquake engineering has experienced in last years substantial advances. New design and analysis strategies have been proposed, leading to relevant developments. These developments rely on extensive testing and numerical simulation mainly based on oversimplified models referred in this work as structural component-based models, as a result of their moderate computational cost. Therefore, there is a strong need of verifying the reliability of the new developments by comparison with analyses performed using more advanced simulation tools and with experiments. This work is organized in two parts. First part presents an accurate model, while the second part deals with a more simplified model, although highly computational efficient.

First part. This research clarifies the aforementioned issues by developing a new continuum mechanics-based model for simulating the monotonic and cyclic behavior of reinforced concrete structures. The developed model combines a new methodology for calculating the damage variables in Concrete Plastic Damage Models “CPDM”, and a new approach to integrate CPDM with a 3-D interface bond-slip model developed by other researchers. A new scheme to implement the interface model in a continuum FEM model of regions with crossing reinforcement bars is also presented in this research. Mesh-insensitivity, accuracy and reliability of the proposed model are verified by simulating several monotonic and cyclic tests; the obtained results are compared with experimental ones, satisfactory agreement has been accomplished.

Second part. The developed model in the First Part is compared with oversimplified structural component-based models that are commonly used in earthquake engineering; results have shown the superiority of the proposed model to predict the actual behavior of highly damaged RC elements and frames, capturing strength reduction, stiffness degradation and pinching phenomena. However, some of the structural component-based models have shown an acceptable performance considering the law computationally cost in comparison with the advanced continuum mechanics-based model. After this conclusion, this part presents a numerical study on the relation among the non-simulated deterioration modes of the elements in non-ductile RC frames and their final capacity. A structural component-based model has been developed for simulating the nonlinear dynamic behavior of non-ductile reinforced structures, accounting for flexure, shear and axial deterioration modes. The developed model is numerically efficient, thus being suitable for day use in earthquake engineering. The capacity of the developed model is verified by simulating the nonlinear dynamic behavior of an existing non-ductile building and the prototype building. Obtained results show that the developed model, despite its moderate computational cost, detects and reproduces accurately the nonlinear dynamic behavior of non-ductile RC structures, as well, capturing the deterioration modes that are blind to the simplified models. Comparison with results from more simplified models highlights the importance of hidden failure modes in the behavior of each element and in the overall collapse mechanisms. The relation between the non-simulated failure modes and the so-called "Structural Resurrection" is addressed.

## List of symbols

### Roman letters. Lower case

ac / at / bc / bt: dimensionless coefficients in equations

b: ${\textstyle {\epsilon }_{c}^{pl}/{\epsilon }_{c}^{ch}}$ ratio

c / c1 / c2: cohesion / coefficients in the uniaxial tensile behavior of concrete

d / dc / dt: damage variable / compression damage variable / tension damage variable

b´´/d´´: are the width and depth of the confined core respectively

f / fcm / ftm / fc0 / ft0 / fck: stress strength / concrete compressive stress strength / concrete tensile stress strength / limit stress of linear compressive branch / limit stress of linear tensile branch / characteristic value of concrete compressive strength

f / fb0 / fc0 / fcm / ftm / fc0 / ft0 / fck/f'c: stress strength / biaxial compressive yield strength / uniaxial compressive yield strength / concrete compressive stress strength / concrete tensile stress strength / limit stress of linear compressive branch / limit stress of linear tensile branch / characteristic value of concrete compressive strength/concrete compressive strength.

fy / fu: longitudinal reinforcing bars stresses yield / ultimate

fyh / fu: transverse reinforcing bars yield stress

gc / gt: compressive / tensile energies per unit volume dissipated by damage along entire deterioration process

${\textstyle {h}_{c}}$ / ${\textstyle {h}_{t}}$: weighting factors accounting for stiffness recovery

leq: mesh size (finite element characteristic length)

r*: stress state; for uniaxial stress r*(σ11) = 1 for tension and r*(σ11) = 0 for compression

sc / st: coefficients accounting for stress state and stiffness recovery effects

p: hydrostatic pressure stress

q: Von Mises-equivalent effective stress

r*: stress state; for uniaxial stress r*(σ11) = 1 for tension and r*(σ11) = 0 for compression

sc / st: coefficients accounting for stress state and stiffness recovery effects

w / wc: crack opening / crack opening at fracture

### Roman letters. Upper case

Db: reinforcement bar diameter

E / E0 / Eci: modulus of deformation / undamaged modulus of deformation / tangent modulus of deformation of concrete for zero stress

Es / Esh: steel modulus of elasticity / slope of hardening branch

EIflex : column effective flexural rigidity

G / Gch / GF: flow potential / crushing energy per unit area / fracture energy per unit area

H: Mohr-Coulomb yield surface function

I1: first invariant of stress tensor

J2 / J3: second / third invariants of deviatoric stress tensor

Kc: ratio of second stress invariants on tensile and compressive meridians

Kdeg: shear degradation stiffness (Elwood 2005 mode [1])

Cslip: Stiffness of linear slip spring

L: column height

L: mesh size for the spring bond-slip model.

MSP: moment at spalling of concrete

Fbs : is the equitant bond force

S1 / S2 / S3: relative normal displacement / longitudinalslip / transverse slip

Speak / ${\textstyle {S}_{max}^{+}\,/\,{S}_{max}^{-}\,/{\,S}_{cum}}$: peak slip / maximum positive slip / maximum negative slip / cumulated slip

SR:clear spacing between bar ribs

S: is the center to center spacing of the hoops

As´´: is the cross-sectional area of the hoop bar

Ag/Aeff: are the gross section area/ the effective shear section area.

Rc : is the pinching parameter of [124] model.

M: fourth-order tensor for anisotropic damage model

I identity tensor

### Greek letters

ε / εc / εt / εel / εpl / εcm / εtm: strain / compression strain / tensile strain / elastic strain / plastic strain / strain at compressive strength / strain at tensile strength

${\textstyle {\epsilon }_{c}^{pl}}$ / ${\textstyle {\epsilon }_{t}^{pl}}$ / ${\textstyle \,{\epsilon }_{c}^{el}}$ / ${\textstyle \,{\epsilon }_{t}^{el}}$ / ${\textstyle \,{\epsilon }_{c}^{ch}}$ / ${\textstyle \,{\epsilon }_{t}^{ck}}$ / ${\textstyle {\epsilon }_{0c}^{el}}$ / ${\textstyle \,{\epsilon }_{0t}^{el}}$: strains at Figure 27; subindexes “c”, “t”, “0c” and “0t” and refer to compression, tension, undamaged compression and undamaged tension, respectively; superindexes “pl”, “el”, “ch” and “ck” and refer to plastic, elastic, crushing and cracking, respectively

εs / εysh u: axial steel strain / yielding steel strain / hardening steel strain / ultimate steel strain

${\textstyle {\epsilon }_{c}^{ch}}$ / ${\textstyle \,{\epsilon }_{t}^{ck}}$: crushing compressive strain / cracking tensile strain

${\textstyle {\epsilon }_{c}^{pl}}$ / ${\textstyle \,{\epsilon }_{t}^{pl}}$: plastic compressive / tensile strain

εsh / εu: steel strain that corresponds to onset of hardening / ultimate strain

ε50u 50c: are the strain corresponds to 50% of concrete maximum compressive strength for unconfined and confined concrete respectively

ρb,s / ρf,s: reduction factor in bar bearing resistance accounting for steel yielding in tension / reduction factor in bar friction resistance accounting for steel yielding in tension

ρb,c / ρf,c: reduction factor in bar bearing resistance accounting for slip history / reduction factor in bar friction resistance accounts for slip history

ρn: reduction factor accounting for opening of splitting cracks

ρsh: is the transverse reinforcement ratio

ρ´´ : is the volumetric ratio of confining hoops to volume of concrete core measured to the outside of the perimeter hoops.

ϵ: eccentricity of the plastic potential surface

ϕ: friction angle

θ: Lode similarity angle

θ: is the critical crack angle (Elwood 2005 mode [1])

θcap: is the plastic rotation at onset of strength loss (Ibarra et al model [2])

θpc: is the post-capping stiffness (Ibarra et al model [2])

λ: is the energy dissipation parameter Ibarra et al model [2])

σ / σ11 / σt0 / σc(1) / σc(2) / σc(3): stress / first principal uniaxial stress / uniaxial tensile stress at failure / concrete compressive stress at first / second / third segment

σ' /α' : the deviatoric parts of the stress and backstress tensors σ /α respectively

σ0 : the effective stresses tensor

ε/ εe/ εpl: are tensors of (total strain, elastic strain, plastic strain)

${\textstyle {\overline {{\sigma }_{c}}}}$ / ${\textstyle {\overline {{\sigma }_{t}}}}$ : effective compressive / tensile cohesion stress

σ1:normal bond stress

σ / σ11: normal stress / first principal uniaxial stress

τ2 / τmax / τb/ τf:longitudinalshear bond stress / peak bond strength / due to bar bearing resistance/ due to friction resistance

τ3:transverseshear bond stress

ξ: distance from origin of stress space to stress plan

ѱ: dilatancy angle

${\textstyle {C}_{m\,}and\,{\gamma }_{m\,}}$: Steel material parameters

Φy: Yield curvature

Δyflexslipshear: are the yield drift displacement/ displacement due to flexure/ displacement due to slip/ displacement due to shear.

Δy/L, Δy/L: Are the drift limit angle at flexure-shear and flexure-shear-axial deterioration (Elwood 2005 mode [1])

## General introduction to earthquake engineering

This section presents a brief introduction and a concise historical review of earthquake-resistant design of structures and the nonlinear analysis of RC structures.

The first seismic analysis methods appear on the year 1923 in Japan (after the earthquake in Kanto [3] and can be included within the package of so-called Earthquake Analysis Methods Based on Resistance. These procedures were intended to provide buildings with lateral (horizontal) resistance; it was believed that if the structure of the building had enough lateral resistance it should be capable to survive the design earthquake. This resistance is guaranteed by designing the structure to be able to withstand horizontal forces applied at each floor level and in each direction of the building (usually two orthogonal directions). Figure 1 illustrates this concept.

 Figure 1. Lateral forces that are equivalent to a seismic input.

In Figure 1, V is the sum of the forces acting at each floor level; in other words, the horizontal interaction force between the ground and the building. V is also known as base shear. Obviously, the value of V quantifies the severity of the earthquake effect on the building.

In the firstly developed earthquake-resistant design methods, horizontal forces represented in Figure 1 were obtained by multiplying the weight of each floor by a constant coefficient. This ratio between the horizontal and vertical forces was called seismic coefficient and in the first 1923 Japanese Seismic Code [3] it was estimated as 0.1. This value gradually increased as it was experienced that structures designed with this resistance value failed when an earthquake stronger than expected occurred. This ratio took to the values of 0.10, 0.15 and 0.20 until, thanks to the development of computers and by having more and more seismic experiences, it was concluded that structures that had been designed with a certain lateral resistance, did not reach collapse but could suffer damage in the case of a larger earthquake. After that, resistance was not the primary goal and everybody started paying more attention to the ductility; it can be roughly defined as the ability of a given structure to resist after the onset of damage. The ductility of a given building can be estimated from observed damages. The regulations began to introduce the concept of ductility by quantifying it with a response reduction factor, which reduces the equivalent lateral forces shown in Figure 1; it was mentioned in the 1957 American design code [4]. Thus, this approach has been incorporated to the current worldwide regulations. In summary, most of the earthquake-resistant regulations require to provide buildings with a certain level of lateral resistance. This resistance is obtained by dividing the resistance that a given building should have to remain in the elastic range under the design input by the aforementioned response reduction factor. This factor should obviously be equal to or greater than the unity. This coefficient is represented by different symbols in each standard; in the European standard [5] it is named q, in the United States [6]it is known as R. It is remarkable that, in fact, this ratio does not take into account only the ductile behavior of the structure but also includes the over-resistance of the building due to the conservative considerations that are regularly considered (safety factors, among others) and the increase of the material resistance under dynamic inputs (“strain rate effect”).

In any case, it should be kept in mind that in these methods the effect of the earthquake on the structure is characterized by means of equivalent static forces Figure 1; they are determined as those that generate a lateral displacement equal to the maximum one that would occur along the duration of the earthquake. However, another possible strategy is to represent the seismic action by a much more direct way: as input accelerograms. In this case, dynamic analysis must be performed to determine the time-history responses; then, the maximum values will be selected, they would represent the design demands. This formulation is often referred to as earthquake-resistant design based on dynamic calculations. This strategy seems appropriate and has apparently shown to be quite capable of simulating the actual seismic behavior of structures with great accuracy and reliability; however, there are some drawbacks that hinder the use of such formulations: (1) the information about the earthquakes that may occur for a particular structure during its lifetime is limited, which severely impairs the accuracy of the study, (2) for economic reasons, structures are designed to behave nonlinearly during the design earthquake (the most severe earthquake expected with a reasonable probability) and, hence, nonlinear dynamic analyses are a must. Dynamic analyses in the nonlinear regime are much more complex than the, already complex, dynamic linear calculations.

Earthquake-resistant design based on spectra:

In general terms, these methods are based on estimating the equivalent static forces (which characterize the effect of the seismic action) in terms of the fundamental period of the structure. This is done by using response spectra; they are plots whose ordinates are certain response magnitudes and whose abscissas are the natural periods of SDOF systems that represent the structure. Up to date, three types of spectra have been basically proposed: absolute acceleration, relative displacement, and energy spectra. In the absolute acceleration spectra, the ordinates are the ratio between the maximum absolute acceleration in the top of the building and the maximum input acceleration in the base of the building. In the relative displacement spectra, the ordinates are the ratio between the maximum relative displacement between the top and the base of the building and the maximum input relative displacement. In the energy spectra, the ordinates are the input energy introduced by the seismic input in the building. It is noteworthy that each of these three spectra considers a meaningful response magnitude: the relative displacement is an indicator of the apparent structural damage level (i.e. not cumulative), the absolute acceleration is related the human perception of the motion and the damage to the facilities (and, more generally, to all the non-structural elements) and the energy reports on the accumulated structural damage.

Linear spectra plot the ratio between the maximum values of the response of an elastic single-degree-of-freedom system and of the input acceleration. Figure 2 shows an elastic model of a single-degree-of freedom system undergoing a horizontal ground motion zg.

 Figure 2. Elastic single-degree-of-freedom systems.

In Figure 2, m, c and k are the mass, damping and stiffness coefficients, respectively, y is the relative displacement between the mass and the base (degree-of-freedom) and zg is the displacement of the ground. Noticeably that the model in Figure 2 can incorporate the relative displacement between the foundation and soil, this is termed as Soil Structure Interaction SSI . Yet this formulation is commonly applied to horizontal motion, can be also considered for vertical vibrations.

The equation of motion of the system described in Figure 2, is given by

 ${\displaystyle m{\ddot {y}}+c{\overset {\cdot }{y}}+ky=-m{\ddot {z}}_{g}}$
(1)

By dividing both sides by m, (1) becomes

 ${\displaystyle {\ddot {y}}+2z{w}_{0}{\overset {\cdot }{y}}+{w}_{0}^{2}y=-{\ddot {z}}_{g}}$
(2)

In this relationship, ${\textstyle {w}_{0}}$ is the undamped natural frequency of the system and ζ is the critical damping factor. These coefficients are given by

 ${\displaystyle {w}_{0}={\sqrt {\frac {k}{m}}}}$ ${\displaystyle z={\frac {c}{2m{w}_{0}}}}$
(3)

The damped natural frequency ${\textstyle {w}_{d}}$ is related to ${\textstyle {w}_{0}}$ and to ζ by

 ${\displaystyle {w}_{d}={w}_{0}{\sqrt {1-{z}^{2}}}}$
(4)

It is remarkable that, unless the damping ζ takes extremely high values, ${\textstyle {w}_{0}}$ and ${\textstyle {w}_{d}}$ are nearly coincident.

The acceleration, velocity and displacement spectra are obtained, for each input zg(t), as the maximum values of the absolute acceleration ${\textstyle \,{\ddot {x}}\,}$ (where ${\textstyle {\ddot {x}}=}$${\displaystyle {\ddot {y}}+{\ddot {z}}_{g}}$), relative velocity ${\textstyle {\overset {\cdot }{y}}\,}$ and relative displacement y. They depend on the natural period T (T = 2 π/ ω0) and on the damping factor ζ. These quantities are obtained by the following linear relationships [7]:

 ${\displaystyle y\left(t\right)=-{\frac {1}{m{w}_{d}}}\int _{0}^{t}m{\ddot {z}}_{g}\left(t\right)\mathrm {sin} \,\left({w}_{d}\left(t-t\right)\right){e}^{-z{w}_{0}\left(t-t\right)}dt}$
(5)
 ${\displaystyle {\overset {\cdot }{y}}\left(t\right)=\int _{0}^{t}{\ddot {z}}_{g}\left(t\right)\mathrm {cos} \,\left({w}_{d}\left(t-t\right)\right){e}^{-z{w}_{0}\left(t-t\right)}dt-}$${\displaystyle {\frac {z{w}_{0}}{{w}_{d}}}\int _{0}^{t}{\ddot {z}}_{g}\left(t\right)\mathrm {sin} \,\left({w}_{d}\left(t-t\right)\right){e}^{-z{w}_{0}\left(t-t\right)}dt}$
(6)
 ${\displaystyle {\ddot {x}}\left(t\right)=\mathrm {\left({\frac {2{z}^{2}{w}_{0}^{2}}{{w}_{d}}}-{\frac {{w}_{0}^{2}}{{w}_{d}}}\right)\int _{0}^{t}{\ddot {z}}_{g}cos} \,\left({w}_{d}\left(t-t\right)\right){e}^{-z{w}_{0}\left(t-t\right)}dt-}$${\displaystyle 2z{w}_{0}\int _{0}^{t}{\ddot {z}}_{g}\left(t\right)\mathrm {cos} \,\left({w}_{d}\left(t-t\right)\right){e}^{-z{w}_{0}\left(t-t\right)}dt}$
(7)

Figure 3.a.b.c show, relative displacement, relative velocity and absolute acceleration spectra, respectively. Such spectra correspond to the accelerogram registered in the ICA2 station (E-W component) during the Pisco earthquake, 15 august 2007.

 1) Relative displacement spectra 2) Relative velocity spectra 3) Absolute acceleration spectra Figure 3. Displacement, velocity and acceleration spectra for Pisco earthquake 2007- ICA2-EW

Figure 3.a.b.c show that the spectral ordinates decrease with the increasing damping ratio; this shows that damping has a beneficial effect, since it contributes to reduce relevant response magnitudes (relative displacement, relative velocity and absolute acceleration). Moreover, the spectrum corresponding to zero damping exhibits sharper peaks than the spectra for non-zero damping; it means damping contributes to smoothen the spectra, e.g. making it less sensitive to small period changes.

It has been demonstrated [8] that for small values ​​of damping and not too long periods (under 10 seconds), the velocity spectra are obtained by multiplying the acceleration spectra by T / 2π and that the displacement spectra are obtained in the same way from the velocity ones:

 ${\textstyle S_{v}=S_{a}(T/2\pi )\qquad S_{d}=S_{v}(T/2\pi )=S_{a}(T/2\pi )^{2}}$ (8)

Absolute acceleration response spectra

As discussed in the previous subsection, the absolute acceleration response spectra are curves that represent, in ordinates, the ratio between the maximum values of the absolute acceleration of the SDOF system that represents the dynamic behavior of the structure in a given vibration mode and the ground acceleration. The design spectra are smoothed envelopes obtained from a number of individual records. Figure 4 shows the design spectrum of the Spanish regulation [NCSE-02 2002].

 Figure 4. Design acceleration spectrum [NCSE-02 2002].

The spectrum shown in Figure 4 consists of three branches: a linearly increasing one, a constant one and a hyperbolically decreasing one. Periods TA and TB depend on the characteristics of the soil, being higher as it has less stiffness; in some codes, the spectral ordinate (e.g. the height of spectrum) also grows as the flexibility of the soil does. The interpretation of each of these branches in terms of the effect of the earthquake on the structure is quite clear: (1) short-period structures are very rigid and tend to behave as the surrounding soil, but its motion is amplified as its rigidity decreases, (2) in the medium period range, the ground motion reaches its highest amplification inside the building and, (3) in the long periods range, structures are flexible enough so that its stiffness is not capable of overcoming the high inertia forces. This interpretation helps us to understand the influence of the soil stiffness in TA and TB: for stiff soil the range of building periods whose motion is highly amplified (in between TA and TB) is narrow, while this range widens and encompasses higher rise buildings as the soil becomes less stiff. This spectrum is commonly presented in dimensionless form (the ordinates Sa are dimensionless).

Performance-based earthquake-resistant design

The objective of the current seismic design codes is to prepare the structure to resist the design seismic input only under ultimate limit state; in other words, the structure is intended to resist the design earthquake with an acceptable level of serious damage but without collapse (in other words, avoiding at all costs the loss of human lives). Remarkably, that approach does not include any requirement about the behavior under seismic actions with lower or higher level of severity; this contrasts with the usual strategy against other type of actions (gravity, for example) where two types of limit states (ultimate and service) are considered. This approach is broadly valid and has been used for decades but was in shortage especially after the Northridge earthquake in 1994 and Kobe in 1995; after these highly severe earthquakes it was found that some structures, even those relatively new and that had been designed with the latest seismic standards, did not collapse (and in them there were no human casualties), but the damage to buildings (both structural and non-structural) was very serious. In the Kobe earthquake, some hospitals had been so intensely reinforced that effectively its structure did not collapse but absolute accelerations in the building were so high that it damaged the installations and were unusable at the time of greatest need (a few hours after the earthquake). After these events, the earthquake engineering was directed not only to prevent loss of human lives but also to quantify, reduce and prevent the damage. Depending on the damage, we are able to accept when an earthquake occurs, different solutions can be proposed. This strategy is commonly known as “Performance Based Design”; it is mainly described in the references [9], [10], [11], [12]. These documents present different seismic design methodologies oriented to control and to quantify the level of structural damage due to seismic action and to design structures that do not exceed the corresponding level.

Based on structural and non-structural damage, the following four levels of performance (“Performance States”) are defined according to SEAOC 1995 [9]:

• Fully Operational. Uninterrupted service. Negligible structural and non-structural damage.
• Operational. Most of the activities can be resumed immediately. The structure is safe and can be inhabited. The essential activities are maintained while the non-essential ones are interrupted. Repairs are necessary to resume the non-essential activities. Slight damage.
• Life Safe. Moderate damage, the structure remains safe. Some elements or components of the building may be protected to avoid damage. The risk of loss of life is low. The building may need to be evacuated after the earthquake. The repair is possible, but can be economically unfeasible.
• Near Collapse. Severe damage, but without risk of collapse. Possible fall of non-structural elements.

More recently, another similar classification is considered according to ATC-40 [10] and FEMA 356 [12];

• Immediate Occupancy. Achieve essentially elastic behavior by limiting structural damage (e.g., yielding of steel, significant cracking of concrete, and non-structural damage.). Important services are not uninterrupted. The period of lack of functionality (“down time”) is about 14 hours.
• Damage Control. Slight structural damage. Achievable occupants’ safety. The essential activities are repairable. Moderate overall damage. The period of lack of functionality (“down time”) is about 2 or 3 weeks.
• Life Safety. Limit damage of structural and non-structural components but no collapse. No risk from falling non-structural elements. The evacuation of the occupants can be done without risk. Possibility of irreparable building.
• Collapse Prevention. Severe structural damage, with risk of collapse. Likely fall of non-structural elements. The evacuation of the occupants may involve risk. Building likely irreparable.

These four levels are often represented by their initials: IO, DC, LS and CP. The three levels IO, LS and CP are the most commonly used for seismic design; Figure 5 presents a graphical and easily understandable way, the practical significance of these levels and their relationship with the percentage of damage. The case “operational” in this case refers to a building without any damage.

 Figure 5. Damage levels [13]

Regarding the seismic action, four levels of severity a defined as specified in Table 1.

Table 1. Severity levels of the seismic inputs
 Design Earthquake Return Periods (years) Probability of Occurrence Frequent 43 50% in 30 years Occasional 72 50% in 50 years Rare 475 10% in 50 years Very rare 970 10% in 100 years

Table 1 shows that the severity of the earthquakes is quantified in terms of their return period; it is understood as the average of the elapsed time among earthquakes with the same magnitude or, almost equivalently, as the inverse of the probability of occurrence in one year. In some cases, seismic actions more severe than those contained in Table 1 are considered; the so-called MCE (“Maximum Considered Earthquake”) corresponds to a return period of about 2475 years. The relationship between the return period T and the probability pn of being exceeded n years is given by the expression ${\textstyle T=}$${\displaystyle -n/ln\left(1-{p}_{n}\right)}$ ; it is often used to indicate the severity of an earthquake by the probability p50 to be exceeded in 50 years, for example, in the case of MCE

is ${\textstyle {p}_{50}=1-{e}^{-{\frac {50}{2475}}}=0.02}$ and in the case of an earthquake “Rare” is ${\textstyle {p}_{50}=1-}$${\displaystyle {e}^{-{\frac {50}{475}}}=0.10}$.

Table 2 shows the demand levels for each of the four performance levels previously described [9] , when the earthquakes that have the probability of occurrence specified in Table 1 occur.

Table 2 shows three levels of protection (expressed by the three represented diagonals): less intense for systems of moderate importance (“Basic Facilities”), more intense for major facilities (“Essential / Hazardous Facilities”) and even more intense for crucial facilities (“Safety critical Facilities”). For example, in “Essential / Hazardous Facilities” (diagonal terms) it is required that for an earthquake of return period of 75 years the building remains fully operational, for an earthquake of return period of 475 years the building keeps operating in its major functions and for a return period of 970 years the building is able to preserve the lives of its occupants.

Table 2. Severity levels of the seismic inputs and damage levels

Nonlinear static analysis (Pushover)

The method of earthquake-resistant design based on nonlinear static analyses consists basically of comparing the capacity of the structure, characterized by a capacity curve representing its behavior under pushing incremental forces, with the effect of the design earthquake, characterized by a demanding spectrum. The intersection between both curves is termed as target drift or target displacement or the most common “the performance point”; in other words, that point indicates the effect produced by the earthquake on the structure ATC-40 1996 [10]. The capacity curve is usually expressed by representing on the vertical axis the base shear force V and on the horizontal axis the displacement of the top floor. The analysis that generates this curve is static, monotonic, and obviously is nonlinear, being commonly known as push-over. Figure 6 shows an example of a capacity curve obtained from a push-over analysis.

 Figure 6. Capacity curve obtained from push-over analyses ATC-40 1996 [10]

In the push-over analyses the base shear force is distributed along the floors according to certain patterns; the most commonly used are the first modal shape, uniform or linear (“triangular”) distributions. The push-over analyses are made incrementally, in other words, the lateral forces are increased progressively. For small values ​​of V, the behavior of the structure is linear and as V increases, the structure is becoming gradually more damaged; the stiffness of the structure decreases and its capacity curve becomes more flat. The smallest slope of the capacity curve with the increasing displacement illustrates clearly the elongation of the natural period of the structure.

Some researchers [14], [15] have proposed techniques to modify the distribution of the lateral forces among the floors to take into account the variation of the modal properties and the contribution of higher modes " Modal Push-Over Analysis".

The demand is characterized by the design spectrum for the considered level of seismic action Table 1; to be able to intersect it with the capacity curve, it is represented as the absolute acceleration spectrum Sa (vertical axis) vs. the relative displacement spectrum Sd (horizontal axis). This type of representation is commonly known as “Acceleration-Displacement Response Spectra” (ADRS).

The methods mostly used to obtain the target displacements are: Capacity Spectrum Method ATC-40 1996 [10], Displacement Coefficient Method FEMA-356 2000 [12], Equivalent linearization method FEMA-440 2005 [16], Modified displacement coefficient method FEMA-440 2005 [16], Modified Capacity Spectrum ATC-40 1996 [10].

The analysis push-over characterizes the nonlinear dynamic behavior of the structure by means of increasing static forces. The main drawback of this strategy is that the response of the structure to a given input is not incremental but cyclical and the push-over analyses cannot take into account the accumulated plastic strain, in other words, the cumulated damage. Therefore, we cannot establish a clear relationship between the maximum displacement of the structure and the energy stored during plastic deformation cycles. When the structure enters the inelastic range, deterioration occurs by the accumulation of plastic incursions; that can produce the complete breakdown of structural elements for deformations smaller than those that could be resisted under monotonic forces. This type of failure is called low cycle fatigue or plastic fatigue. Another drawback of the earthquake design strategy based on displacements is that the hysteretic behavior is interpreted as an equivalent viscous damping (ζeq); this introduces a relevant error, especially for significant levels of damping. Moreover, such identification is not based on any physical principle that justifies, in inelastic systems, the existence of a direct relationship between the energy corresponding to the maximum displacement and the equivalent viscous damping.

Dynamic analysis:

This procedure evaluates the effect of earthquakes on buildings based on determining the dynamic response (commonly known as “time history”) to the expected accelerograms. The most relevant response quantities are; the maximum relative displacements in between consecutive floors (inter-story drifts), and the maximum absolute accelerations. The maximum relative displacements report about the experienced level of structural damage, the maximum absolute accelerations are directly correlated with the non-structural damage and the human comfort conditions. Since the dynamic calculations take into account the performance of buildings under seismic inputs in a more direct way than in the methodologies based on response spectra, in general the dynamic analyses are able to provide more accurate results. In particular, the comparison between the nonlinear static methods (push-over) and the nonlinear dynamic methods is clearly favorable to them because, besides being more accurate in general, they have two important advantages: (i) by considering the cyclical behavior they are able to reproduce the accumulated plastic damage and (ii) the consideration of the effect of damping (both the present in the undamaged structure and the generated for increasing damage) is more direct.

In zones of medium or high seismicity, buildings are often designed by accepting a given level of structural damage under the design earthquake. Accordingly, in these cases the dynamic analyses should be nonlinear, in other words, must be able to reproduce the behavior of the structure when it been damaged and therefore has experienced significant reductions in its strength and rigidity. Moreover, second-order analyses may be necessary because of the significant relative horizontal displacements, this being another source of complexity and increased computational cost. Although the nonlinear dynamic analyses are increasingly used in the earthquake-resistant design of important structures, this procedure is rarely used in the design of ordinary structures, this is due to the high computational cost involved and to the effort required to properly interpret the large amount of generated information.

Incremental Dynamic Analysis IDA

With the main purpose of alleviating the problem derived from the fact that the push-over analysis cannot take into account the accumulated plastic strain, the so-called Incremental Dynamic Analysis IDA has been proposed [17]. It is remarkable that the incremental dynamic analyses require making several nonlinear dynamical calculations, which are expensive in computational time; on the other hand, it may be necessary to perform second-order analyses. However, the incremental dynamic analyses, especially when applied to several earthquakes, constitute powerful formulations, which may provide greater and more useful information than the rest of approaches that have been described in this section.

The results of these procedures are usually represented by the so-called IDA curves. These representations consist of capacity curves similar to the result of the push-over analyses; on the horizontal axis, an index related to the magnitude of the response is usually represented and the vertical axis usually contains an index related to the severity of excitation. Figure 7 shows the results of this kind; Figure 7.(a) corresponds to a single record and Figure 7.(b) corresponds to multiple (30) records. In both representations the severity of the seismic action is quantified by the ordinate of acceleration response spectrum for the first mode Sa(T1, 0.05) and the magnitude of the response is quantified by the maximum value (along the duration of the earthquake) of the relative displacement between floors (inter-story drift). Figure 7.(a) shows both increases and decreases of the damage on the upper floors with increasing severity of excitation, this effect is obviously due to the “protection” provided by the lower floors. None of the other methods described in this section are able to predict this phenomenon so clearly. Figure 7.(b) shows the remarkable variability in the response of a determined structure to records that have, in first approximation, a comparable level of severity.

 1) Single Register 2) Thirty Registers Figure 7. Examples of IDA curves [17]

Usually the damage thresholds IO, LS and CP are related with certain values ​​of the index that quantify the magnitude of the excitation (ordinate in Figure 7); in this way performance-based analyses can be made from incremental dynamic calculations.

Deterioration modes in RC structures

While buildings are usually designed for seismic resistance using elastic analysis, most will experience significant inelastic deformations under large earthquakes. Deterioration modes for RC frame components have been identified in [18] based on a review of experimental tests, published information and observation from past earthquakes. They are classified into six modes depending on the type of structural element and the physical behavior associated with deterioration shown in Figure 8. For each mode, currently available nonlinear components models are rated from low to high confidence “1 to 5” for their ability to simulate the deterioration.

 Figure 8. Deterioration modes for reinforced concrete frame components [18]

Progressive deterioration can cause structural collapse, this can be generated in two forms; Local Collapse; identified as part of structure has been collapsed in result of ground shaking causes elements deterioration modes form a collapse mechanism, the localized part in RC frame structures is one floor that is usually the so called “soft story”, vertical loads still can be transmited by the non-damaged component. Second form is the Global Collapse; the entire structure collapsed and cannot transmit any vertical neither lateral loading. In regards to structure collapse, two possible mechanisms “scenarios” are defined; Sidesway Collapse induced by dynamic instability as result of unlimited story drifts amplified by a combination of P-Δ effects and deterioration in strength of the structural components. Vertical Collapse; this kind of collapse occurs as a result of the loss of one or several elements the vertical load carrying capacity causes a local collapse, and eventually entire vertical collapse due to a developed progressive collapse mechanism. These possible collapse mechanisms are principally dependent on the developed deterioration modes at the components level. Contribution of the predefined deterioration modes on the collapse mechanism "scenario" is defined in Figure 9 according to [18].

   Figure 9. Deterioration modes for reinforced concrete frame components [18]

Nonlinear Analysis

Modern performance-based design methods require ways to determine the realistic behavior of structures under strong earthquakes. Enabled by advancements in computing technologies and available test data, nonlinear analyses provide the means for calculating structural response beyond the elastic range, including strength and stiffness deterioration associated with inelastic material behavior and large displacements. As such, nonlinear analysis plays an important role in the seismic assessment of new and existing buildings. In contrast to linear elastic analysis and design methods that are well established, nonlinear inelastic analysis techniques and their application to design are still evolving and may require engineers to develop new skills. Nonlinear analyses require thinking about inelastic behavior and limit states that depend on deformations as well as forces, different degradation modes among their interaction at the structural components must to be captured.

Depending on the structural configuration, the results of nonlinear analyses can be extremely sensitive to the assumed input parameters and of implemented numerical models, these models must confirm the locations of expected or observed inelastic deformations and trace the structural behavior up to the onset of collapse. This requires sophisticated models that are validated against physical tests to capture the highly nonlinear response approaching collapse.

Capturing the aforementioned deterioration models and collapse mechanisms requires that the numerical model has the capacity to represent; the behavior of concrete and steel materials at high level of monotonic and cyclic loading, as well the interaction between both materials during the loading process characterized by the bond-slip behavior. Numerical models can be differentiated by the way that plasticity is distributed through the member cross sections and along its length. Figure 10 illustrates several idealized model types for simulating the nonlinear response of beam-column elements, same concepts can be applied for other element (e.g. braces and shear walls).

 Figure 10. Beam-Column nonlinear models [19]

Concentrated Plasticity Models are the simplest formulation, representation of such models are shown in Figure 10-a,b. these models are based on concentrating the inelastic deformations at a specific location where the damage is expected to be located (e.g. ends of beams and columns). Numerically this is represented by attaching zero-length or spring element at location of the inelastic deformation, the behavior of this nonlinear element is governed by the so-called “backbone” or capacity curve that can be represented by moment-rotation or force-deformation relation for ductile and brittle behavior respectively. For cyclic loading, a hysteretic behavior is also incorporated in the behavior of the nonlinear element. Both the capacity curve and the hysteretic properties in these models are usually calibrated after experimental results and observation of past earthquakes. This formulation can be used to reflect the effect of some degradation phenomena such as bond-slip or bar buckling.

Distributed Plasticity Models are more complex to represent the distribution of the inelastic deformation. The distribution varies as well between the simple and the more complex. The Finite Length Hinge Model Figure 10-c is an efficient distributed plasticity formulation with designated hinge zones at the member ends. Cross sections in the inelastic hinge zones are characterized through either nonlinear moment-curvature relationships or explicit fiber-section integrations that enforce the assumption that plane sections remain plane. The inelastic hinge length may be fixed or variable, as determined from the moment-curvature characteristics of the section together with the concurrent moment gradient and axial force. Integration of deformations along the hinge length captures the spread of yielding more realistically than the concentrated hinges, while the finite hinge length facilitates calculation of hinge rotations.

Another more complex distributed plasticity formulation is the Fiber Section Model Figure 10-d. This model distributes plasticity by numerical integrations through the member cross sections and along the member length. The section is discretized by several fibers “portions”. Uniaxial material models are defined to capture the nonlinear hysteretic axial stress-strain characteristics in the cross sections, each fiber is assigned to a material (e.g. steel, confined and unconfined concrete). The plane-sections-remain-plane assumption is enforced, where uniaxial material “fibers” are numerically integrated over the cross section to obtain stress resultants (axial force and moments) and incremental moment-curvature and axial force-strain relations. The cross section parameters are then integrated numerically at discrete sections along the member length, using displacement or force interpolation functions. The most complex models are the Continuum Models Figure 10-e. These models discretize the continuum along the member length and through the cross sections into small (micro) finite elements with nonlinear hysteretic constitutive properties. This fundamental level of modeling offers the most versatility since the inelastic behavior can be captured in all the location and for different material, moreover the relation between concrete and steel “bond-slip” can be treated also by a constitutive model. The big challenge for Continuum formulation is providing integrated efficient model for RC simulation.

From the previous literature; two main modeling formulations can be defined to simulate the complex inelastic behavior of RC structures, the definition is based on the method each model presents the degradation at the components and material levels. The most commonly used formulation in earthquake engineering and performance based design is the Structural Component-Based Formulation represented by the models in Figure 10-a,b,c,d, and the most complex and general formulation is the Continuum-Based Formulation represented by Figure 10-e. In this dissertation, the simple and advanced formulations are verified and used to develop an appropriate tool for simulating the nonlinear dynamic behavior of RC frame structures.

### Objectives

#### Main objective

To establish the requirements of developing a computationally efficient model for nonlinear cyclic behavior for RC structures. This study will consider advanced models derived from continuum mechanics-based formulations, and the simplified structural component-based models. This research will focus on RC frames behavior.

#### Particular objectives

This global objective will be pursued through the following particular objectives:

• Review the numerical models commonly used for describing the nonlinear dynamic behavior of RC frames under earthquake excitation. This study encompasses the most spread models based on structural component-based approaches (e.g. Lumped Plastic, Distributed Plastic).
• Review the more sophisticated numerical models derived from continuum mechanics formulations. This study encompasses the models for concrete, reinforcing bars and the models to describe the relation between both materials under monotonic and cyclic loading.
• Development of 3D continuum mechanics-based model for simulating the nonlinear cyclic behavior of RC structures, the model uses concrete plastic damage model and phenomenological interface bond-slip model.
• Verifying the capacity of the developed 3D continuum mechanics-based model, as well the simplified components-based models with experimental results of RC structures subjected to monotonic and cyclic loading.
• Development of an advanced and day use structural component-based model to simulate the nonlinear dynamic behavior of non-ductile RC structures.
• Performing Incremental Dynamic Analysis (IDA) study using the developed structural component-based model. The study investigates the importance of hidden failure modes as well as the relation between these modes and the so called “Structural Resurrection Phenomena”

### Outline

This dissertation is organized in two parts, description of each part as well the new contribution is presented as follow:

Part I focuses on the "numerical models of RC structure that are based on continuum mechanics-based formulations". Two major contributions are developed in this part. First contribution is a new methodology for calculating the damage variables "damage evolution" in Concrete Plastic Damage Models (CPDM), implementation of this methodology for CPDM of ABAQUS finite element analysis program is described as well. The second contribution is a new modeling scheme for 3D implementation of an interface bond-slip model developed by [20] and implemented in ABAQUS. Combining the methodology for calculating damage evolution in CPDM, as well the modeling scheme of the bond-slip model, results in an integrated continuum-based model for monotonic and cyclic simulation of RC structure. Proposed model is implemented in ABAQUS and verified with experimental results. Part I is organized in seven chapters as follow:

Chapter 1 provides an introduction to Part I, motivation and new contribution are also described. Chapter ‎2 describes in general the modeling of concrete material, models derived from plasticity, damage and coupled plastic-damage are described in this part, as well as the used concrete plastic damage model. Chapter 3 presents the proposed methodology for calculating the damage evolution in concrete plastic damage models. Chapter 4 describes the used model of steel material, while Chapter ‎5 describes the bond-slip phenomena and the implemented interface model in this work, the new modeling scheme is described in the subsection ‎5.4. Chapter ‎6 presents numerical simulation of real experiments conducted on reinforced concrete structures, experimental results are compared with ones obtained by the proposed models, finally Chapter ‎7 addresses a general concluding remarks about Part I.

Part II focuses on numerical models that are based on simplified formulations referred as "structural component-based approaches" models. In this part the oversimplified models and commonly used in earthquake engineering are verified with experimental results and compared with the advanced model derived in Part I. The main contribution in this part is developing an advanced structural component-based model for nonlinear dynamic simulation of non-ductile RC framed structures. A numerical study is conducted using the developed model on the relation among the non-simulated deterioration modes of the elements of non-ductile RC frames and their final capacity. Part II is organized in six chapters as follow:

Chapter 1 provides an introduction to Part II, motivation and new contribution are also described. Chapter ‎2 describes the structural component-based models commonly used in the simulation of reinforced concrete structures. Chapter ‎3 verifies the capacity of the models presented in chapter ‎2 to capture the nonlinear behavior of the simulated experiments in Part I, obtained results are also compared with the ones from the continuum mechanics-based model developed in Part I. Chapter ‎4 presents the developed model to simulate the nonlinear dynamic behavior of non-ductile reinforced concrete structures. Chapter ‎5 verifies the capacity of the presented model in chapter ‎4 by performing numerical simulation of a two non-ductile structures, obtained results are also compared with experimental ones. Conclusions of this part are addressed in Chapter ‎6.

In the final conclusions chapter, the achievements of the present study are exposed along with concluding remarks and future works derived from this dissertation.

## 1 Introduction and motivation

Under severe seismic excitation, structural behavior of buildings and other constructions is highly complex. It involves, among other issues, soil-structure interaction, large strains and displacements, damage, plasticity, and near-collapse behavior. Moreover, in reinforced concrete structures, there are several coupled degradation and failure modes: cracking, crushing and spalling of concrete, yielding and pull-out of tensioned reinforcement, yielding and buckling of compressed reinforcement. Numerical model must be able to predict the global behavior of the elements among the interaction of the forces acting on each element. In other words, the model must capture the interaction between axial, shear, torsion and flexural moment, as well the subsequent degradation modes represented in strength reduction, stiffness degradation and pinching. Models based on simplified theories referred in this work as structural component-based models are commonly used in earthquake engineering, as a result of their moderate computational cost. Such models might serve to predict the seismic behavior of RC structures in some particular cases. However, in many cases the oversimplified models are not capable to capture adequately the entire range of deformation till the structure reaches its collapse state.

Therefore, in earthquake engineering, advanced numerical simulations based on continuum mechanics are strongly necessary. These models are commonly used in small size simulation due to their expensive computationally cost. The motivation in this part is to developed computationally efficient model based on advanced formulations. The developed model can be used for simulating small size structures or to investigate particular joints or some parts in reinforced concrete structures. Another application of the proposed model is to improve the performance of existing simplified models, as well to calibrate new models that can be used afterwards to simulate large scale models.

Using continuum mechanics-based formulations to describe the behavior of a mixed material such concrete implies that each material is assigned to a set of elements, the relation between materials must be also described by incorporating new elements or by imposing some bond conditions. In the developed model, concrete material is simulated by solid finite elements, steel bars are simulated by truss finite element and bond between steel and concrete is simulated by either; an interface elements or by assuming perfect bond conditions for cyclic and monotonic loading respectively. This part of the work describes shortly the bases of concrete constitutive modeling, as well the used damage-plastic model. Regarding steel bars, the basic plasticity model is described with its implementation. The complex relation between concrete and steel bars "bond-slip" is revised under monotonic and cyclic loading, modeling these phenomena in a 3D continuum model is explained.

The main contributions in this part are: Developing a new methodology for calculating the damage variables "damage evolution" in Concrete Plastic Damage Models "CPDM" [21] [22] [23] [24] [25] [26] [27] [28] [29] [30], this methodology is described in section ‎3. As well as a new scheme to implement interface bond-slip model in 3D FEM simulation described in subsection ‎5.4. As a final contribution, a new approach to integrate CPDM and bond-slip model in simulating the cyclic and monotonic behavior of RC structures is presented in subsection ‎6.1

## 2 Concrete modeling

### 2.1 Introduction

Quasi-brittle materials, as concrete, exhibit nonlinear stress-strain response mainly because of micro-cracking. Cracks are oriented as the stress field and generate the failure modes. In tension, failure is localized in a narrow band; stress-strain behavior is characterized by sudden softening accompanied with reduction in the unloading stiffness. In compression, failure begins usually in the outside and is more complex, involving volumetric expansion, strain localization, crushing, inclined slipping and spalling; stress-strain behavior involves ductile hardening followed by softening and reduction in the unloading stiffness. In mixed stress states, failure depends usually on the ratio between the principal stresses. Nonlinear concrete response can be represented using plasticity or damage theory. However, none of these formulations alone is able to describe adequately this phenomenon. Plastic models [31] [32] might represent realistically the observed deformation in high confined concrete but do not capture the stiffness degradation observed in experiments [21]. Damage-based models [33] [34] [35] [36] are based on gradual reduction of the elastic stiffness; they can describe the stiffness degradation in tension and low confined compression, but are not suitable to capture the irreversible deformations observed in experiments and the inelastic volumetric expansion in compression. In addition, fracture propagation can be represented by embedded crack models, where standard FEM interpolations are enriched with strain or displacement discontinuities [37] [38] [39]. These models can be used for high strain localization problems (fracture). The mechanical behavior among the constitutive modeling of concrete are described in the next sections.

### 2.2 Mechanical behavior

The fundamental characteristics of concrete behavior are established through experimental testing of plain concrete specimens subjected to specific, relatively simple load histories. Continuum mechanics provides a framework for developing an analytical model that describe these fundamental characteristics. Experimental data provide additional information for refinement and calibration of the analytical model. Some important mechanical features of concrete are summarized in this section. This furnishes a background for the review and further study on the constitutive modelling of concrete in the following subsections.

#### 2.2.1 Uniaxial behavior

The behavior of concrete is highly nonlinear in both uniaxial tension and compression. Under uniaxial compression stress state, five different deformational stages can be observed as shown in Figure 11.a. For axial stresses up to about 30% of the maximum compressive stress, the uniaxial compressive behavior of concrete can be considered linear with existing micro-cracks in the material remaining nearly unchanged (Zone A). The second stage is between 30 % and 50% of the maximum compressive stress, results in some reduction in the initial stiffness due to the significant increase in crack imitation (Zone B). noticeably that crack does not continue to grow under constant load. Third stage is between 50 % and 75% of maximum compressive stress, results also in more reduction in the initial stiffness due to further initiation and growth of crack (Zone C). At this stage crack continue to grow under constant load. For stress state more than 75% of the maximum compressive strength, results in increasing compressive strain under constant loading (Zone D). After reaching the maximum compressive stress, softening behavior is observed with crack localization that varies according to the boundary conditions (Zone E). The softening under compression is more complicated than tension one and no general model is yet accepted such the crack model in tension. This is because of many factors such as the size and shape of the specimen, concrete strength and boundary conditions.

 a) Stress-strain relation [40] b) normalized stiffness reduction evolution versus normalized plastic strain under cyclic loading [41] Figura 11. Uniaxial concrete response under monotonic and cyclic compression loading

In uniaxial tension stress state, observed deformation process is different from that in compression. The low tensile strength of concrete is primarily due to the low tensile strength of the aggregate-mortar interface, which has a significantly lower strength than the mortar. This interface is known to be the weakest link in this composite material, with cracks usually occurring at the interface. Figure 12 shows typical concrete stress stain response under cyclic loading [42]. The response is essentially linear before reaching the maximum tensile strength capacity, corresponds to stable microcracks. After reaching the maximum strength, softening is observed due to the development of localized continuous crack system mainly perpendicular to the stress direction. Unloading and reloading after reaching the maximum capacity showed also degradation in the initial stiffness.

#### 2.2.2 Multiaxial behavior

Since plain concrete in RC elements is usually subjected to multi-dimensional stresses, developed constitutive models are not only based on the uniaxial behavior of concrete. A number of researchers have studied the biaxial and triaxial behavior of concrete experimentally and analytically.

Investigations by [43] and [44] mainly focused on biaxial loading and developed experimentally a failure surface as shown in Figure 13. The 2D failure surfaces were extended by data from [45] by studying the effect of low levels of confinement pressure in the third direction. Results obtained by [45] showed that relatively small confinement pressure can significantly increase the concrete strength in the direction perpendicular to the confinement pressure up to 2.5 times the uniaxial compressive strength Figure 13.

 1) Biaxial compressive failure surface for concrete. Data from [43,44,45] and presented by [46] 2) Concrete tensile and compressive meridians. Data from [47-50] and presented by [51] Figure 13. Experimental failure surfaces for concrete

Tow envelopes must be considered to describe the behavior of concrete; elastic limit surface defining the elastic region and the failure surface characterizing the maximum-strength envelope as shown in Figure 14.a. It's been wildly accepted that concrete behavior can be considered as isotropic, therefore both surfaces can be described in terms of stress invariants; I1 the first invariant of the stress tensor, J2 and J3 the second and third invariants of deviatoric stress tensor respectively. Haigh-Westergaard space can be used to define the failure surface in principal stress space Figure 14.a.b, in which the position of a stress point is determined by three coordinates ρ , ξ and θ, where ρ is the octahedral radius, ξ is the distance from the origin of stress space to the stress plan, θ is the Lode similarity angle.

 ${\displaystyle {I}_{1}\,={\sigma }_{11}+{\sigma }_{22}+{\sigma }_{33}}$
(9)
 ${\displaystyle {J}_{2}\,={\frac {1}{6}}\left[{\left({\sigma }_{11}-{\sigma }_{22}\right)}^{2}+\right.}$${\displaystyle \left.{\left({\sigma }_{22}-{\sigma }_{33}\right)}^{2}+{\left({\sigma }_{33}-{\sigma }_{11}\right)}^{2}\right]+}$${\displaystyle {{\sigma }_{12}}^{2}+{{\sigma }_{13}}^{2}+{{\sigma }_{23}}^{2}}$
(10)
 ${\displaystyle {J}_{3}\,={\frac {1}{27}}\left(2{I}_{1}^{3}-9{I}_{1}{I}_{2}+27{I}_{31}\right)}$
(11)
 ${\displaystyle \xi \,={\frac {{I}_{1}}{\sqrt {3}}}\qquad \qquad }$ ${\displaystyle \rho \,={\sqrt {2{J}_{2}}}}$
(12)

 a) In principal stress space [52] b) In deviatoric plans [53] c) In the meridian sections [53] Figure 14. Concrete failure surfaces in different views

From experimental results in Figure 13, it can been seen that the failure surface has an open shape. The shapes of the failure surface in the deviatoric and meridians plans is shown respectively in Figure 14.b.c, it can be seen that the deviatoric plans are different in size and shape according to the level of the hydrostatic pressure. ρs, ρt and ρc correspond respectively to shear, tension and compression meridians and are determined according to Lode angle.

### 2.3 Concrete constitutive models

In principle, it is desired that the above-mentioned macroscopic features of the

material behavior be reflected in any constitutive models dedicated to concrete

modelling. However, it is quite difficult to incorporate all of these aspects of material

behavior in a constitutive model. Those experimentally observed features are all of

macroscopic nature, which can only be characterized through some material and

structural quantities and cannot always represent what truly happens at the microscopic

level. This is the disadvantage of the macroscopic approach to constitutive modelling. However, developing an analytical model for concrete has been the aim of many researchers in the last thirty years. Early models were based on elastic theory, more recent use plasticity, damage and fracture mechanics. Majority of the models capture particular aspects concrete response under different type of loading, however, most of the models have acceptable accuracy and efficiency.

#### 2.3.1 Models based on plasticity

Plasticity in material can be defined as the accumulation of uncovered deformation upon loading beyond the yield limit. Experimental results showed that concrete exhibits plastic strain when it's subjected to compression and tension, therefore, constitutive model must incorporate plastic theory.

Developing a plasticity-based constitutive model requires defining; rule for decomposition of the total stain, elastic material constitutive relationship, yield/failure surfaces law and flow rule.

The total stain is usually assumed to be the sum of elastic and accumulated plastic strain

 ${\displaystyle {\overset {\cdot }{\boldsymbol {\varepsilon }}}\,={\overset {\cdot }{{\boldsymbol {\varepsilon }}^{e}}}+}$${\displaystyle {\overset {\cdot }{{\boldsymbol {\varepsilon }}^{pl}}}}$
(13)

Elastic constitutive relationship follows Hook's law

 ${\displaystyle {\boldsymbol {\quad \quad \,\,}}{\overset {\cdot }{\boldsymbol {\sigma }}}=\,{\boldsymbol {D}}_{0}^{el}\,:{\overset {\cdot }{{\boldsymbol {\varepsilon }}^{e}}}=}$${\displaystyle {\boldsymbol {D}}_{0}^{el}\,:\left({\overset {\cdot }{\boldsymbol {\varepsilon }}}-{\overset {\cdot }{{\boldsymbol {\varepsilon }}^{pl}}}\right)}$
(14)

${\textstyle {\boldsymbol {D}}_{0}^{el}}$ is the elastic stiffness tensor, ${\textstyle {\boldsymbol {\sigma }}}$ and ${\textstyle {\boldsymbol {\epsilon }}}$ are the stress and strain tensors.

For concrete, the available material data facilitated definition of the yield surface in stress space and it is most appropriate to consider a yield surface that evolves as a function of the load history. A hardening rule defines the evolution of the yield surface. The flow rules define the evolution of a set of internal variables that uniquely define the material state. In particular, a flow rule defines the orientation of plastic strain which may be associated, defined as normal to the yield surface, or non-associate. Proposed plastic models for concrete vary in the definition of the yield surface, the hardening rules and the flow rules.

Regarding the distinction of the yield surface and the failure surface Figure 14.a, we can see that these two surfaces coincide in plasticity theory. In other words, a single loading surface acts as a yield-failure surface in plasticity theory. This combined surface is often a scaled down version of the failure envelope of the material. Numerous forms of yield surfaces have been proposed and can be classified based on either the number of model parameters [53] or on the shape of the surface in principal stress space.

The first yield criteria to characterize the behavior of plain concrete were Mohr-Columb and Druker-Prager [42] shown in Figure 15, these criteria were developed to describe the response of material such as sand, rock and concrete for which hydrostatic pressure affects the material yield and failure strengths.

The Mohr-Coulomb criterion is defined as follows:

 ${\displaystyle H\left({I}_{1},\,{J}_{2},\theta ,\phi ,c\right)=\,{\frac {1}{3}}{I}_{1}\,sin\,\phi +}$${\displaystyle {\sqrt {\,{J}_{2}}}sin\,\left(\theta +{\frac {\pi }{3}}\right)+\,{\frac {\sqrt {\,{J}_{2}}}{\sqrt {\,3}}}cos\left(\theta +\right.}$${\displaystyle \left.{\frac {\pi }{3}}\right)sin\phi -c\,cos\phi =0}$
(15)
 ${\displaystyle \theta =\,{\frac {1}{3}}\mathrm {acos} \,\left({\frac {3{\sqrt {3}}}{2}}{\frac {\,{J}_{3}}{\,{{J}_{2}}^{\frac {3}{2}}}}\right)}$
(16)

In the previous, ϕ is the friction angle, and c is the material cohesion. Equations (15) and (16) represent a straight line of variable slope in the meridian plane and an irregular hexagon in the ${\textstyle \pi }$ plane as shown in Figure 15.

 Figure 15. Mohr-Coulomb and Drucker-Prager Failure Criteria [42]

The Drucker-Prager criterion represents moderately well the response of plain concrete

subjected to multi-axial compression and provides a smooth yield surface Figure 15. This criterion is defined as follows:

 ${\displaystyle H\left({I}_{1},\,{J}_{2},\alpha ,y\right)={\sqrt {\,{J}_{2}}}+\alpha {I}_{1}+y=}$${\displaystyle 0}$
(17)

In the previous criterion, α and y are material constants. The values of theses constants significantly change the yield surface as shown in Figure 16.

 Figure 16. Drucker-Prager Failure Criterion Compared with Experimental Data as presented by [51]

As indicated in Figure 16, Imran and Pantazopoulou 1996 [48] propose α =0.3 for characterizing the response of concrete subjected to triaxial compression. The response of concrete subjected to biaxial compressive loading Kupfer et al. 1969 [43] and Yin et al. 1989 [44] is characterized well by α =0.1. All of the presented yield criteria are calibrated to predict the observed uniaxial compressive strength.

Comparison of the Drucker-Prager criterion with experimental data shows that while

the criterion may be used to represent the response of concrete subjected to multi-axial

compression, the model over-estimates the capacity of concrete subjected to compression-tension or tension-tension type loading. Variation in concrete response under various load regimes has been addressed by a number of researcher through the use of multi-surface plasticity models. Murray et al. 1997 [54] model proposes a three surface model to characterize the response of plain concrete subjected to biaxial loading, this approach was extended to concrete loaded in three-dimensions by Chen and Chen 1975 [31] and Lubliner et al. 1989 [30].

The main shortcoming of Mohr-Columb and Druker-Prager surfaces is that they assume a linear relationship between J2 and I1 (between ρ and ξ in the meridian plan), although this relationship has been experimentally shown to be nonlinear Figure 14.c. More recent failure criteria with nonlinear relation between J2 and I1 with incorporating Lode angel θ in the formulation have incorporated the [53], [55], [51], [56]. The typical deviatoric and meridian sections of those failure surfaces are shown in Figure 17

 Figure 17. Deviatoric and meridian sections of two typical failure surfaces

In plasticity theory, the definition of a yield surface, the shape of which is usually similar to that of the failure surface (i.e. the yield surface by Grassl et al. 2002 [56]), is required. However, as pointed out by [53] yield surfaces as scaled down versions of failure surfaces at maximum loading are inadequate for concrete modelling. The open shape of such yield surfaces does not reflect the true behavior of concrete under hydrostatic loading. A solution for this is the use of an additional “cap surface” for the behavior of the model under hydrostatic compressive pressure.

The volumetric expansion of concrete under compression makes the application of the associated flow rule for concrete inappropriate. A non-associated flow rule, which is defined by the plastic potential other than the yield function should be used instead. This feature is included in the model developed by [53], [22], [55]. Deterioration in the initial stiffness beyond unloading and reloading has a significant impact on the response of RC structures. Such aspect should be captured by the analytical model that aims to represent the behavior under cyclic loading. Conventional plasticity models are not able to represent this degradation as plastic theory was invented for metallic materials.

#### 2.3.2 Models based on damage mechanic

The damage "stiffness degradation" of a continuous solid is an alteration of the elastic properties during load application due to a decrease of the effective strength area. This effective area loss is normally caused by the increase of voids and/or micro fractures [57]. The continuous damage theory was first introduced by Kachanov in 1958 [58], further contribution was given by Rabotnov in 1963 [59] by introducing the concept of effective stresses. However, the basic development of continuum damage mechanics only began in the 1970s and then in the 1980s with a more rigorous basis, based on thermodynamics and micromechanics. Since then there have been numerous continuum damage mechanics models proposed for the constitutive modelling of materials in general and concrete in particular.

The quantities of continuum mechanics are defined at a mathematical point. However, from the physical point of view, and accounting for the heterogeneity of the material in reality, these quantities should be considered in two scale [60] to have been averaged over a certain volume called a “Representative Volume Element” whose size depends on each material [61]. To define the material damage at a mathematical point M in Figure 18, let us consider a Representative Volume Element (RVE) oriented by a plane defined by its normal ${\textstyle {\overset {\rightarrow }{n}}}$ and its abscissa x along the direction ${\textstyle {\overset {\rightarrow }{n}}}$.

 Figure 18. Definition of damage variable [61]

The damage value D(M, ${\textstyle \,{\overset {\rightarrow }{n}}}$, x) at point M in the direction ${\textstyle {\overset {\rightarrow }{n}}}$ and at abscissa x is defined as:

 ${\displaystyle D\left(M,{\overset {\rightarrow }{n}},x\right)={\frac {\delta {S}_{{D}_{x}}}{\delta S}}}$
(18)

δS is the area of the intersection considered plan and the RVE, and ${\textstyle \delta {S}_{{D}_{x}}}$ is the effective area of intersections of all microcracks and microcavities in δS as illustrated in Figure 18. The anisotropic aspect of damage can be understood by assuming that the failure of the RVE depends on the direction on ${\textstyle {\overset {\rightarrow }{n}}}$. If microcracks and microcavities are uniformly distributed in the RVE, it's adequate to assume isotropic damage and damage variable D(M, ${\textstyle \,{\overset {\rightarrow }{n}}}$, x) doesn't depend on the direction, thus and scalar damage will control the deterioration.

To understand better the damage concept, let's consider the case of uniaxial tension force F with scalar damage variable. After reaching the maximum tensile strength and after occur of damage, the cross sectional area is reduced and becomes the effective cross sectional area SSD. Where S is the original cross sectional area and SD is the total area of microcracks. The stress is no longer σ = F/ S but replaced by the effective stress σ0 = F / (SSD) = σ (1− D) ≥ σ. This concept can be extended to multiaxial stress state.

Physically, the degradation process of the virgin materials is due to the presence and growth of small fractures and micro voids. This growth process can be simulated within the context of the mechanics of continuous media, taking into consideration the theory of internal state variables, introducing an internal variable of damage represented by a scalar, vector or tensor. This internal variable of damage characterizes the level of damage in the material and transforms the stress real tensor into the so called "effective stress tensor" as follow:

 ${\displaystyle {\boldsymbol {\sigma }}_{\mathit {\boldsymbol {0}}}={\boldsymbol {M}}^{-1}\,:{\boldsymbol {\sigma }}}$
(19)

Mis a fourth-order tensor for anisotropic damage model. It's well approved that damage in concrete can be accepted to be isotropic. However, damage models using scalar damage variables are still preferred because of their simplicity in the formulation.

For isotropic damage mode, the material degradation is alike in all direction and the tensor M is reduced to M = (1-d) I, where I is an identity tensor. The equation (19) can be written for scalar-damage model:

 ${\displaystyle {\boldsymbol {\sigma }}_{\mathit {\boldsymbol {0}}}={\frac {\boldsymbol {\sigma }}{\left(1-d\right)}}}$
(20)

where d is an internal scalar variable range between 0 for no damage state and 1 for full distortion, σ0 is effective stress tensor measured in the non-damaged space. The effective stress concept was formulated for the first time in connection with the equivalent deformation hypothesis by [62] as follow " The deformation associated to a damaged state subjected to a stress σ is equivalent to the deformation associated to the non-damage state subjected to an effective stress σ0). A graphical representation of the effective and real space is illustrated in Figure 19.

 Figure 19. Graphic representation of the effective stress hypothesis presented by [57]

Continuum damage mechanics approach has been proved by many authors to be appropriate for constitutive models of concrete; different damage models have been proposed by researchers, among others [63], [64], [34], [65], [66], [67] [36].

Material deterioration in damage-based model is controlled by damage variable. The evolution law of damage, which plays a very important role in any damage-based model, is different for many continuum damage mechanics models. However, it is possible to group almost all existing approaches into three categories: one when the damage evolution law is imposed (e.g. [68,65,67]). Another approach is based on obtaining the damage evolution from a dissipation potential [69,61]. Third approach is based on using implicitly defined damage evolution law (e.g. [23,24,25], ).

#### 2.3.3 Coupling between damage and plasticity

It is been widely accepted that coupling between damage and plasticity models is essential to capture the nonlinear behavior of concrete [26]. Coupled damage and plasticity models for concrete differ mainly in the coupling method and the damage evolution law.

In the implicit coupling, damage and plasticity are implicitly embedded in the yield and damage criteria [23,24,25], with the material strength being a decreasing function with respect to the damage variable. This implicit coupling characterizes the strength reduction due to the material deterioration and is equivalent to introducing effective instead of nominal stress into the yield function. Therefore, the concept of effective stress is still applicable in this case. This way of introducing coupling enables the constitutive modelling to use separate yield and failure criteria. The corresponding internal variables (damage variable and plastic strains for the coupled model) of the model do not explicitly depend on each other. Nevertheless, the parameter identification becomes more difficult as the model responses in this case are governed by all the internal variables used in the coupled model.

Other researchers describe coupling using a single function, in which only one loading function is specified and used to control the dissipation process. This function can be a damage loading function [68] or a yield function [22,61], the concept of effective stresses is applied by defining the yield function in the effective stresses space. In the first case with a damage loading function governing the dissipation process, an evolution law for the plastic strain is required [68]. For the use of a yield function, damage measure can be based on some criteria or by postulating damage variables law [22,61]. Figure 20 shows a representation of coupling between plasticity and damage [30].

 (a) Plasticity Model (b) Damage Model (c) Plastic Damage Model Figure 20. Representation of coupling between damage and plasticity

In coupled models, plasticity for concrete can be described with isotropic hardening; however, damage in many cases is not isotropic but has preferential directions [70]. Some plastic isotropic damage models have been proposed, e.g. [21-23,25-29]; these models have shown good performance in capturing concrete behavior in tests on full-scale structures [70] [26]. Anisotropy can be added to the damage to capture the anisotropy feature of concrete both for compression and tension. Although anisotropic damage models are complex and coupling with plasticity in the application to practical engineering is not straightforward, researchers have investigated this issue and proposed plastic anisotropic damage models, among others [64,71-80].

Advantages and limitations of some coupled plastic damage models for concrete are presented briefly as follow:

Yazdani and Schreyer 1990 [81]: In this approach, the pressure-dependent damage surface is used as a failure surface, and enhanced with a Von Mises yield surface. In principal stress space, this yield surface covers the damage surface in the positive

hydrostatic axis, and lies almost entirely inside the damage surface in the opposite

direction Figure 21. Attention was also paid to the multi-axial behaviour of the

model, with the failure surface being calibrated to account for the strength enhancement and ductility under increasing lateral confinement Figure 21.

 Figure 21. Coupled damage-plasticity model by Yazdani and Schreyer 1990 [81]

Compared to other coupled damage-plasticity models, this model is able to capture the irrecoverable strains using inelastic damage mechanism to account for the misfit of crack surfaces. Nevertheless, this capability is restricted to compressive mode of cracking only [81]. As a consequence, with a pure damage mechanism activated in tension, the model exhibits an inability to capture the observed permanent deformation in tensile loading.

Lee and Fenves 1998 [22]: This model has been developed for the numerical analysis of concrete dam, it was originally proposed for monotonic, cyclic and dynamic behavior by Lubliner et al.1989 [30]. In this model, a combined yield-failure surface is defined and evolves with damage variables, for which the evolution laws are postulated. The yield surface in two-dimensional principal stress space is shown in Figure 22.

 Figure 22. Yield surface in plane stress space Lee /Lubliner [22,30]

Yield criterion is based on a modified version of Drucker-Prager function. Stiffness deterioration is described by a scalar damage that is calculated as a function of the uniaxial damage variables in tension and compression. The opening/closing of microcracks also features in the model. This model has no damage loading function, and the evolution of the yield surface is governed by the damage variables.

Although the whole model is defined from several “pieces”, its thermodynamic

consistency can be readily ensured. However, the model showed very good performance in simulating concrete behavior under uniaxial and biaxial stress state, but is not recommend for high multiaxial stress state. This model is adopted in this research. The constitutive behavior, as well the implementation is presented in the next subsection.

Faria et al. 1998 [68]:

Similarly to the model of Lee and Fenves 1998 [22], this model possesses only one loading surface and was formulated on the basis of an effective stress tensor. However, in contrast with that by [22], a composite damage loading surface, instead of a yield surface, is used to govern the constitutive behavior of the model.

 a) Initial 2D elastic domain b) Uniaxial cyclic behavior Figure 23. Coupled damage-plasticity model by Faria et al. 1998 [68]

The composite damage loading surface is formed by two separated damage surface, these function are used to define the strength envelope shown in Figure 23a. The evolution of damage variables in tension and compression are defined by law for each that contains three material parameters to be defined experimentally. This model showed good performance for biaxial and multiaxial loading, but fails to capture the permanent deformation in tension as shown in Figure 23b.

### 2.4 Concrete Plastic Damage Model

Since concrete behavior is brittle, but, under stress reversal, tensile cracks might close, then broken parts being reassembled. Therefore, concrete behavior can be better described with models that combine damage and plasticity. These models are particularly well suited for reproducing failure modes that are based on tensile cracking and compression crushing.

After the previous literature, in this work the coupled model of Lee and Fenves 1998 [22] is adopted. This model is described with a multiaxial model that considers parallel combination of scalar (isotropic) damaged elasticity and non-associated multi-hardening plasticity. This model is termed as “Concrete Plastic Damage Model” (CPDM) along this research. Was proposed for monotonic, cyclic and dynamic behavior by [30] and was further developed by [22]. This model shows good performance in primarily uniaxial and biaxial stress states, but should not be used in case of significant triaxial compressive stresses. Since the purpose of this research is to simulate the monotonic and cyclic behavior of RC structures where no serious triaxial compressive stresses are expected, hence, this model can be considered appropriate. The fundamental bases of this model are described in this subsection, necessary parameters to use the predefined model in the FEM software ABAQUS are also described. This subsection furnishes the base for the next subsection where a new methodology to define damage evolution is presented.

This models assumes that the main two failure mechanisms are tensile cracking and compressive crushing of the concrete material. The evolution of the yield (or failure) surface is controlled by two hardening variables, ${\textstyle {\epsilon }_{t}^{pl}}$ and ${\textstyle {\epsilon }_{c}^{pl}}$, linked to failure mechanisms under tension and compression loading, respectively. In Lee and Fenves 1998 [22] model, for uniaxial compression and tension, the stress-strain relation under uniaxial loading in displayed in Figure 20.c, can be written as:

 ${\displaystyle {\sigma }_{c}\,=\left(1-{d}_{c}\right)\,{E}_{0\,}({\epsilon }_{c}-{\epsilon }_{c}^{pl})}$
(21)
 ${\displaystyle {\sigma }_{t}\,=\left(1-{d}_{t}\right)\,{E}_{0\,}({\epsilon }_{t}-{\epsilon }_{t}^{pl})}$
(22)

where subindexes c and t refer to compression and tension, respectively.

 style="text-align: center;" ${\displaystyle E=\left(1-d\right)\,{E}_{0}\,}$
(23)

In equation (23), E is the reduced tangent stiffness and d is a scalar degradation variable, which is a function of stress state and of compression and tension damage variables (dc and dt, respectively):

 ${\displaystyle 1-d=\left(1-{s}_{t}{d}_{c}\right)\,(1-{s}_{c}{d}_{t})\,}$
(24)

In equation (24), sc and st are dimensionless coefficients accounting for stress state and stiffness recovery effects, being given by

 ${\displaystyle {s}_{c}=\,1-{h}_{c}\,(1-{r}^{\ast }\left({\sigma }_{11}\right))\,}$
(25)
 ${\displaystyle {s}_{t}=\,1-{h}_{t}\,{r}^{\ast }\left({\sigma }_{11}\right)\,}$
(26)

In equations (25) and (26), ${\textstyle {\sigma }_{11}}$ is the first principal uniaxial stress (positive for tension), r* is a stress state parameter being ${\textstyle {r}^{\ast }\left({\sigma }_{11}\right)}$  = 1 for tension and ${\textstyle {r}^{\ast }\left({\sigma }_{11}\right)}$  = 0 for compression, and ${\textstyle {h}_{c}}$ and ${\textstyle {h}_{t}\,}$ are weighting factors ranging between 0 and 1. Factor ${\textstyle {h}_{c}}$ accounts for re-closing of cracks after tension-compression reversal; ${\textstyle {h}_{t}\,}$ represents recovery of crushed concrete after compression-tension reversal. If ℎc = 0.9 and ℎt = 0 for example; this means that 90% of the cracks close upon tension-compression reversal and the crushed concrete does not experience any recovery. Equations (25) and (26) show that ${\textstyle {s}_{c}}$ and ${\textstyle {s}_{t}}$ also range between 0 and 1.

For a better understanding of the effect of sc and st coefficients, Figure 24 displays plots of uniaxial stress-strain loading-unloading behavior. The initial elastic branch with slope E0 reaches the descending branch (Figure 27.b) at peak point 1, then cracking begins; later, unloading starts at point 2. At that point, there is no compression damage and dc = 0, r* = 1, and sc = 1; therefore, equation (24) shows that d = dt. Consequently, the linear unloading branch has slope (1 dt) E0. In the way to stress reversing point 3, cracks begin closing. After point 3, r* = 0, sc = 1 hc, st = 1, dc = 0, and equation (24) shows that d = (1 hc) dt. Therefore, the slope of the ongoing compression segment of the branch depends on parameter hc; three options are plotted in Figure 24: (i) hc = 0 (no crack is closed) with slope (1 dt) E0, (ii) hc = 0.5 (half of the cracks are closed) with slope (1 0.5 dt) E0, (iii) hc = 1 (all cracks are closed) with slope E0. Noticeably, in the third option (hc = 1), there is no compressive strength reduction. At point 4, an unloading branch arises; at that point, r* = 0, sc = 1 hc, st = 1, and equation (24) shows that 1 d = (1 dc) [1 (1 hc) dt] = 1 dc. Point 5 correspond again to stress reversal; after it, provided that ht = 0, the slope of the ongoing branch is equal to (1 dt) (1 dc) E0. Point 6 is the peak for the reduced tensile strength; after it, cracking reinitiates and a new descending branch is generated.

In this work and for cyclic loading simulation, ℎc = 0.9 is assumed in the main concrete bodies (90% of the stiffness in compression is recovered once cracks close), and lower values are taken in the zones where widely-opened cracks and bond-slip effects concentrate. ℎt = 0.0 is assumed for all the sections (referrer to subsection ‎6.1).

For multiaxial condition, the stress-strain relationship is given by:

 ${\displaystyle {\boldsymbol {\sigma }}=\left(1-d\right)\,{\boldsymbol {D}}_{0}^{el}\,:({\boldsymbol {\epsilon }}-}$${\displaystyle {\boldsymbol {\epsilon }}^{pl})}$
(27)

In equation (27), ${\textstyle {D}_{0}^{el}}$ is the elastic stiffness tensor, and ${\textstyle {\boldsymbol {\sigma }}}$ and ${\textstyle {\boldsymbol {\epsilon }}}$ are the stress and strain tensors, respectively. Scalar damage variable d keeps same meaning than for uniaxial condition, although replacing scalar factor ${\textstyle {r}^{\ast }}$ with a multiaxial one [22].

Regarding plasticity model, yield condition is based on the loading (yield) function F proposed in [30] with modifications suggested by [22] to account for different tension and compression strength evolution.

 ${\displaystyle F=\,{\frac {1}{1-\alpha }}\left(q-3\,\alpha \,p+\,\beta \,\left\langle {\sigma }_{max}\right\rangle -\right.}$${\displaystyle \left.\gamma \,\left\langle -{\sigma }_{max}\right\rangle \right)-\,{\overline {{\sigma }_{c}}}=}$${\displaystyle 0}$
(28)
 ${\displaystyle \alpha ={\frac {\left({f}_{b0}/{f}_{c0}\right)-1}{2\,\left({f}_{b0}/{f}_{c0}\right)-1}};\,\beta =}$${\displaystyle {\frac {\overline {{\sigma }_{c}}}{\overline {{\sigma }_{t}}}}\,\left(1-\alpha \right)-}$${\displaystyle \left(1+\alpha \right);\,\gamma ={\frac {3\left(1-{K}_{c}\right)}{2\,{K}_{c}-1}}}$
(29)

In equations (28) and (29), ${\textstyle \left\langle \cdot \right\rangle \,}$ is the Macaulay bracket, ${\textstyle p}$ is the hydrostatic pressure stress, ${\textstyle q}$ is the Von Mises-equivalent effective stress (effective stress accounts for stress divided by 1 d), and fb0 and fc0 are the biaxial and uniaxial compressive yield strengths, respectively; since fb0 ≥ fc0, α ranges between 0 (fb0 = fc0) and 0.5 (fb0 ${\textstyle \gg }$ fc0). ${\textstyle {\sigma }_{max}}$ is the maximum principal effective stress, and ${\textstyle {\overline {{\sigma }_{c}}}}$ and ${\textstyle {\overline {{\sigma }_{t}}}}$ are the effective compressive and tensile cohesion stress, respectively. ${\textstyle {\overline {{\sigma }_{c}}}}$ and ${\textstyle {\overline {{\sigma }_{t}}}}$ are defined as ${\textstyle {{\overline {{\sigma }_{c}}}=\sigma }_{c}\,/\left(1-\right.}$${\displaystyle \left.{d}_{c}\right)}$ and ${\textstyle {{\overline {{\sigma }_{t}}}=\sigma }_{t}\,/\left(1-\right.}$${\displaystyle \left.{d}_{t}\right)}$ . Kc is the ratio of second stress invariants on tensile and compressive meridians.

The plasticity model assumes non-associated potential plastic flow. The flow potential G is the Drucker-Prager hyperbolic function given by:

 ${\displaystyle G={\sqrt {{\left(\epsilon \,{\sigma }_{t0}\mathrm {tan} \,\psi \,\right)}^{2}+{q}^{2}}}\,-}$${\displaystyle \,p\mathrm {tan} \,\psi }$
(30)

In equation (30), σt0 is the uniaxial tensile stress at failure, ϵ is the eccentricity of plastic potential surface, and ѱ is the dilatancy angle measured in ${\textstyle p}$- ${\textstyle q}$ (deviatory) plan at high confining pressure.

As discussed previously, Kc is the ratio between the magnitudes of deviatoric stress in uniaxial tension and compression; Kc ranges between 0.5 (Rankine yield surface) and 1 (Von Mises). In this study, Kc is obtained from the Mohr-Coulomb yield surface function in cylindrical coordinates [57]:

 ${\displaystyle H\left(\rho ,\,x,\theta ,\phi ,c\right)=\,{\sqrt {2}}\,x\,sin\,\phi +{\sqrt {3}}\,\rho \,cos\,\theta -}$${\displaystyle \rho \mathrm {sin} \,\theta \,sin\phi -{\sqrt {6}}\,c\,cos\phi =0}$
(31)

In equation (31), ρ is the octahedral radius, ξ is the distance from the origin of stress space to the stress plan, θ is the Lode similarity angle, ϕ is the friction angle, and c is the cohesion. In equation (31), ${\textstyle r=}$${\displaystyle {\sqrt {2\,{J}_{2}}}}$, ξ = I1 / ${\textstyle {\sqrt {3}}}$, and ${\textstyle \mathrm {sin} \,q=}$${\displaystyle {\frac {3{\sqrt {3}}{J}_{3}}{2{J}_{2}^{3/2}}}}$; I1 is the first invariant of stress tensor and J2 and J3 are the second and third invariants of deviatoric stress tensor, respectively.

For ${\textstyle x}$ = 0 and θ = ± π / 6 (negative/positive for tension/compression meridian plans), the magnitudes of deviatory stress in uniaxial compression and tension at yield (ρc0 and ρt0) and Kc are

 ${\displaystyle \,{\rho }_{c0}=\,{\frac {2\,c\,{\sqrt {6}}\,cos\phi \,}{3-sin\phi \,}}}$ ${\displaystyle {\rho }_{t0}=\,{\frac {2\,c\,{\sqrt {6}}\,cos\phi }{3+sin\phi }}}$ ${\displaystyle {K}_{c}\,={\frac {{\rho }_{t0}}{{\rho }_{c0}}}={\frac {3-sin\phi }{3+sin\phi }}}$
(32)

By assuming that ϕ = 32° [57], last equation in (32) shows that Kc = 0.7. Figure 25 represents the yield surface in deviatory plan for several values of Kc ranging from 0.5 to 1. CM and TM account for compression/tension meridian plans, respectively.

 Figure 25. Yield surface in the deviatory plan for several values of Kc

Equations (30), (28) and (29) show that the concrete behavior depends on four constitutive parameters Kc, ѱ, fb0 / fc0, ϵ; it can be assumed that ѱ =13° [82]. Table 3 describes the values used in this study.

Table 3. Parameters of CPDM
 Kc Ѱ (º) fb0 / fc0 ϵ 0.7 13 1.16 0.1

Since this model uses the concept of effective stresses as it’s indicated in equations (27), (28) and (29), as well, the cyclic uniaxial and multiaxial behavior is dominated by the same variables. All these direct influences of damage variables make the approach to calculate them an essential issue. Defining these variables is a hard task, as no clear guide is available to calculate and implement such variables in the CPDM of Abaqus software. Next section presents a new approach for calculating and implementing these variables in CPDM, the proposed approach can be uses for the model already implemented in Abaqus, as well as any model requires the damage evolution to be defined.

## 3. Proposed methodology for calculating damage variables of CPDM

### 3.1 Introduction

As described before, damage evolution law plays a significant role in any damage model, particularly when this law is imposed. An important number of researchers have proposed different damage evolutions laws. Most of them are based on splitting damage into compressive and tensile parts and each one is determined separately by its evolution law; total damage is calculated with some combination rules e.g. [68], [83], [80], [84]. Few evolution laws are based on general formulae for calculating damage in compression and tension [85]. In most of the available damage evaluating formulae the parameters need to be calibrated experimentally. Following some approaches that can be used to calculate the damage evolution for CPDM:

• Mazars, Pijaudier-Cabot 1989 [34]. Defined a formula to calculate the total scalar damage as the sum of the uniaxial compressive and tensile damage evolution. The model is originally developed by Mazars 1984-1986 ( [33] [86]), and formulated in the next formula
 ${\textstyle d={\alpha }_{t}{d}_{t}+\,{\alpha }_{c}{d}_{c}}$
(33)

dt and dc are the damage evolution for uniaxial tension and compression respectively, αt and αc are weight factor and can be calculated according to [33]. Uniaxial damage evolution is presented in the following formulation, "i" takes the same meaning for "c" compression and "t" tension.

 ${\displaystyle {d}_{i}=\,1-\,{\frac {\left({1-A}_{i}\right){K}_{0}}{\overline {\epsilon }}}-{\frac {{A}_{i}}{exp\left[{B}_{i}\left({\overline {\epsilon }}-{K}_{0}\right)\right]\,}}\,}$
(34)

${\textstyle {\overline {\epsilon }}}$ is the equivalent strain, and the parameters Ac, At, Bc, Bt and K0 are derived from compression test on cylinders and bending test on beams.

• Birtel, Mark 2006 [83]. Tension and compression damage variables "damage evolution" are calculated by assuming constant ratios between compressive plastic and crushing strains, and between tensile plastic and cracking strains, respectively. The uniaxial damage evolution is presented in the next formulation, "i" takes the same meaning for "c" compression and "t" tension.
 ${\displaystyle {d}_{i}=\,1-\,{\frac {{\sigma }_{i}{E}_{0}^{-1}}{{\epsilon }_{i}^{pl}\left({\frac {1}{{b}_{i}}}-1\right)+{\sigma }_{i}{E}_{0}^{-1}}}\,}$
(35)

E0 is the original concrete stiffness modulus, εpl is the uniaxial plastic strain, bi is a parameter that needs to be calibrated from cyclic uniaxial tension and compression tests.

• Häussler-Combe, Hartig 2008 [85]. Defined a scalar damage evolution function with three material parameters. Damage function is given in the next formulation:
 ${\displaystyle D=\,1-{e}^{-{\left[{\frac {\left({k}_{d}-{e}_{d0}\right)}{{e}_{d}}}\right]}^{{g}_{d}}}\,\,}$
(36)

kd is the equivalent damage strain variable defined in [85], ed0, ed and gd are material parameters to be calibrated with experimental results.

• Yu et al. 2010 [87]. This work focuses on confined concrete. A linear expression of damage variable vs. stress is proposed for the compressive descending branch (compressive damage) as follow:
 ${\displaystyle d=\,1-\,{\frac {{\sigma }_{c}}{{f}_{c0}^{'}}}\,}$
(37)

Where ${\textstyle {f}_{c0}^{'}}$ is the maximum strength of concrete. For constant confinement pressure, the damage function becomes:

 ${\displaystyle d=\,1-\,{\frac {{\sigma }_{c}-\,{\frac {1+C+2A}{1-A}}{\sigma }_{1}}{{f'}_{cc}^{\ast }-{\frac {1+C+2A}{1-A}}{\sigma }_{1}}}}$
(38)

${\textstyle {f'}_{cc}^{\ast }}$ the maximum concrete strength under constant confining pressure, σ1 is the confining pressure, C and A are material parameters need to be calibrated from experiments.

• Zheng et al. 2012 [88]. Modified the function developed by Mazars, Pijaudier-Cabot 1989 [34], uniaxial damage evolution function is presented as follow:
 ${\displaystyle {d}_{i}=\,1-\,{\frac {{1-A}_{i}}{1+{\frac {{k}_{i}^{d}}{{\epsilon }_{i}^{f}}}}}-{\frac {{A}_{i}}{exp\left[{B}_{i}{\frac {{k}_{i}^{d}}{{\epsilon }_{i}^{f}}}\right]\,}}\,}$
(39)

${\textstyle {\epsilon }_{i}^{f}}$ is defined as the pick value of strain when damage occurs, ${\textstyle {k}_{c}^{d}}$ and ${\textstyle {k}_{t}^{d}}$ are calculated according to [Zheng et al. 2012], Ac, At, Bc, Bt shall be calibrated form experiments.

• López-Almansa et al. 2014 [84]. Damage variables are determined by an iterative empirical procedure aiming to fit experimental results. The damage evolution function for compression and tension is presented as follow:
 ${\displaystyle {d}_{i}=\,1-\,{e}^{\left(-{a}_{i}{\epsilon }_{i}^{pl}\right)}\,}$
(40)

Parameters ac, at needs to be calibrated from the uniaxial tension and compression stress-strain behavior.

After these investigations, this section proposes a new approach for obtaining damage variables. The need of proposing a new methodology arises from their advantages:

• Is based on the formulation by Lubliner and Lee/Fenves, which is the base for the ABAQUS plastic damage model for concrete. The proposed approach modifies this formulation and obtains closed-form expressions of the damage variables in terms of the corresponding strains; these expressions are derived after integration of concrete fracture and crushing energy.
• No calibration with experimental results is required.
• Implementation is particularly easy.
• Results are insensitive of mesh size, since a strategy aiming to avoid mesh-dependency is incorporated.

All these advantages make this methodology well suited for practical applications. Noticeably, the default ABAQUS model based on the Lubliner/Lee/Fenves approach requires that the uniaxial damage variables (i.e. the damage evolution) are provided by the user, and no guidelines are provided. Therefore, defining the damage evolution is a big concern for any user of the Concrete Damage Plastic Model. In the proposed model, constant mesh size is required in the elements corresponding to the same material. If elements with different size are employed for a single material, their parameters should be defined individually for each size following the proposed methodology.

This approach can consider any concrete constitutive law, either empirical (e.g. like formulations commonly recommended by design codes) or directly based on particular experiments; noticeably, any such constitutive law should account for mesh-sensitivity. In this paper, a particular algorithm that uses laws based on European recommendations is derived; in this case, the only input parameters are concrete compressive strength and mesh size. This algorithm is implemented in the software package ABAQUS [89]; could be also implemented in other computer codes that contain Plastic Damage Model and that require values of damage variables.

Next subsections describe the bases of the proposed approach as well as its implementation for CPDM. Mesh-insensitivity is validated in a simple tension example, accuracy and reliability are also verified by simulating a cyclic experiment on plain concrete specimens. Furthermore, three laboratory experiments are described with the proposed approach to investigate its capability to reproduce the behavior of reinforced concrete (RC) framed structures under monotonic and cyclic loading.

### 3.2 General description

The proposed approach for calculating the damage variables starts from definition of compressive and tensile variables as the portion of normalized energy dissipated by damage:

 ${\displaystyle {\,d}_{c}\,={\frac {1}{{\,g}_{c}}}\int _{0}^{{\epsilon }_{c}^{ch}}{\sigma }_{c}\,d{\epsilon }_{c}^{ch}}$ ${\displaystyle {\,d}_{t}\,={\frac {1}{{\,g}_{t}}}\int _{0}^{{\epsilon }_{t}^{ck}}{\sigma }_{t}\,d{\epsilon }_{t}^{ck}}$
(41)

In equation (41), ${\textstyle {\epsilon }_{c}^{ch}}$ and ${\textstyle {\epsilon }_{t}^{ck}}$ are the crushing and cracking strains respectively, see Figure 27. Normalization coefficients gc and gt represent the energies per unit volume dissipated by damage along the entire deterioration process:

 ${\displaystyle {\,g}_{c}\,=\int _{0}^{\infty }{\sigma }_{c}\,d{\epsilon }_{c}^{ch}}$ ${\displaystyle {\,g}_{t}\,=\int _{0}^{\infty }{\sigma }_{t}\,d{\epsilon }_{t}^{ck}}$
(42)

Equations (41) and (42) show that dc and dt range between 0 (no damage) and 1 (destruction). Figure 26 describes the meaning of gc and gt.

 a) Compression (gc) b) Tension (gt) Figure 26. Parts of energy dissipated by damage

Noticeably, the energies per unit area and per unit volume are related by ${\textstyle {\,g}_{c}\,=}$${\displaystyle {\,G}_{ch\,}/\,{l}_{eq}}$ and ${\textstyle {\,g}_{t}\,={\,G}_{F}\,/{\,l}_{eq}}$; ${\textstyle {\,G}_{ch}}$ and ${\textstyle {\,G}_{F}}$ are material parameters defined as crushing and fracture energies and ${\textstyle {\,l}_{eq}}$ is the characteristic length of the element (subsection ‎3.3).

Relation between compressive and tensile stress and, respectively, crushing and cracking strain, is established, according to [30], as:

 ${\displaystyle {\sigma }_{c}={f}_{c0}\left[\left(1+{a}_{c}\right)exp\left(-{b}_{c\,}{\epsilon }_{c}^{ch}\right)-\right.}$${\displaystyle \left.{a}_{c}exp\left(-{2b}_{c\,}{\epsilon }_{c}^{ch}\right)\right]}$
(43)
 ${\displaystyle {\sigma }_{t}={f}_{t0}\left[\left(1+{a}_{t}\right)exp\left(-{b}_{t\,}{\epsilon }_{t}^{ck}\right)-\right.}$${\displaystyle \left.{a}_{t}exp\left(-{2b}_{t\,}{\epsilon }_{t}^{ck}\right)\right]}$
(44)

In equations (43) and (44), fc0 and ft0 and are the compressive and tensile stresses that correspond to zero crushing ( ${\textstyle {\epsilon }_{c}^{ch}=}$${\displaystyle 0}$) and to onset of cracking ( ${\textstyle {\epsilon }_{t}^{ck}=}$${\displaystyle 0}$), respectively; see Figure 26. As well, ac, at, bc and bt are dimensionless coefficients to be determined. Replacing equations (43) and (44) in equation (42), provides the following relations among gc and gt and such coefficients:

 ${\displaystyle {g}_{c}={\frac {{f}_{c0}}{{b}_{c\,}}}\left(1+{\frac {{a}_{c\,}}{2}}\right)}$ ${\displaystyle {g}_{t}={\frac {{f}_{t0}}{{b}_{t\,}}}\left(1+{\frac {{a}_{t\,}}{2}}\right)}$
(45)

Coefficients bc and bt can be obtained by replacing ${\textstyle {g}_{c}\,=}$${\displaystyle {\,G}_{ch}\,/\,{l}_{eq}}$ and ${\textstyle {\,g}_{t}\,={\,G}_{F}\,/\,{l}_{eq}\,}$ in equation (45):

 ${\displaystyle {b}_{c}={\frac {{f}_{c0\,}{l}_{eq}}{{G}_{ch\,}}}\,\left(1+{\frac {{a}_{c\,}}{2}}\right)}$ ${\displaystyle {b}_{t}={\frac {{f}_{t0\,}{l}_{eq}}{{G}_{F\,}}}\,\left(1+{\frac {{a}_{t\,}}{2}}\right)}$
(46)

By substituting results (43) through (45) in equation (41), the proposed tensile and compressive damage functions are derived:

 ${\displaystyle {\,d}_{c}\,=1-{\frac {1}{2+{a}_{c}}}\left[2\,\left(1+{a}_{c}\right)\,exp\left(-\right.\right.}$${\displaystyle \left.\left.{b}_{c\,}{\epsilon }_{c}^{ch}\right)-{a}_{c}\,exp\left(-{2\,b}_{c\,}{\epsilon }_{c}^{ch}\right)\right]}$
(47)
 ${\displaystyle {\,d}_{t}\,=1-{\frac {1}{2+{a}_{t}}}\left[2\,\left(1+{a}_{t}\right)\,exp\left(-\right.\right.}$${\displaystyle \left.\left.{b}_{t\,}{\epsilon }_{t}^{ck}\right)-{a}_{t}\,exp\left(-{2\,b}_{t\,}{\epsilon }_{t}^{ck}\right)\right]}$
(48)

Therefore, provided that coefficients ac and at are not zero:

 ${\displaystyle exp\left(-{b}_{c\,}{\epsilon }_{c}^{ch}\right)={\frac {1}{{a}_{c}}}\left[1+{a}_{c}-\right.}$${\displaystyle \left.{\sqrt {1+{a}_{c}\,\left(2+{a}_{c}\right)\,{d}_{c}}}\right]}$
(49)
 ${\displaystyle exp\left(-{b}_{t\,}{\epsilon }_{t}^{ck}\right)={\frac {1}{{a}_{t}}}\left[1+{a}_{t}-\right.}$${\displaystyle \left.{\sqrt {1+{a}_{t}\,\left(2+{a}_{t}\right)\,{d}_{t}}}\right]}$
(50)

By zeroing derivatives of σc and σt (equations (43) and (44), respectively) with respect to, respectively, crushing and cracking strain, maximum values fcm and ftm (Figure 27) are obtained:

 ${\displaystyle {f}_{cm}={\frac {{f}_{c0}{\left(1+{a}_{c}\right)}^{2}}{4{\,a}_{c}}}}$ ${\displaystyle {f}_{tm}={\frac {{f}_{t0}{\left(1+{a}_{t}\right)}^{2}}{4{\,a}_{t}}}}$
(51)

Equations (51) provide:

 ${\displaystyle {a}_{c}=2\,\left({\frac {{f}_{cm}}{\,{f}_{c0}}}\right)-1+2\,{\sqrt {{\left({\frac {{f}_{cm}}{{\,f}_{c0}}}\right)}^{2}-\left({\frac {{f}_{cm}}{\,{f}_{c0}}}\right)}}}$
(52)
 ${\displaystyle {a}_{t}=2\,\left({\frac {{f}_{tm}}{{\,f}_{t0}}}\right)-1+2\,{\sqrt {{\left({\frac {{f}_{tm}}{\,{f}_{t0}}}\right)}^{2}-\left({\frac {{f}_{tm}}{\,{f}_{t0}}}\right)}}}$
(53)

These developments complete the description of the proposed methodology. Parameters ac, at, bc and bt can be determined from equations (52), (53) and (46) in terms of fcm, fc0, ftm, ft0, leq, Gch and GF. Then, the damage functions can be calculated according equations (47) and (48). An implementation of the proposed approach using the concrete behavior described in next subsection is presented in subsection ‎3.4.

### 3.3 Uniaxial concrete behavior

This subsection describes the concrete uniaxial stress-strain law that is selected for implementation. Figure 27.a and Figure 27.b explain the compressive and tensile models, respectively; noticeably, tensile model is smeared, i.e. strain includes both crack opening and actual tensile strain between cracks. Figure 27 displays the constitutive laws (thick solid lines) and the unloading / reloading branches (thin dashed lines). Figure 27 corresponds to the damage-plasticity behavior depicted in Figure 20.c. The ascending compressive segments in Figure 27.a follow the Model Code recommendations [90] and the descending segment is engendered as [28]. Tensile stress-strain relation consists of an initial linear segment and a nonlinear descending branch [91] [92], as shown in Figure 27.b. Both compressive and tensile descending branches are generated to ensure nearly mesh-independency; the regularization approach is based on selecting the softening branches of concrete constitutive laws depending on mesh size. In Figure 27, fcm and ftm represent compressive and tensile stress strength, respectively; corresponding strains are εcm and εtm, respectively. According to [90], it is assumed that εcm = 0.0022, fcm = fck + 8 (fck is the characteristic value of concrete compressive strength), and ${\textstyle \,{f}_{tm}\,=}$${\displaystyle 0.3016\,{f}_{ck}^{2/3}}$; these stresses and the deformation modulus are expressed in MPa. In Figure 27.a, ${\textstyle {\epsilon }_{c}^{ch}}$ and ${\textstyle {\epsilon }_{0c}^{el}}$ are the crushing and elastic undamaged components of strain; ${\textstyle {\epsilon }_{c}^{pl}}$ and ${\textstyle {\epsilon }_{c}^{el}}$ are the plastic and elastic damaged components. In Figure 27.b, ${\textstyle {\epsilon }_{c}^{ck}}$ and ${\textstyle {\epsilon }_{0t}^{el}}$ are the cracking and elastic undamaged strain components; ${\textstyle {\epsilon }_{t}^{pl}}$ and ${\textstyle {\epsilon }_{t}^{el}}$ are the plastic and elastic damaged components.

 a) Compression b) Tension Figure 27. Assumed uniaxial model of concrete behavior

First segment in Figure 27.a is linear, ${\textstyle {\,\sigma }_{c(1)}=}$${\displaystyle {E}_{0}\,{\epsilon }_{c}}$, reaching 0.4 fcm; second (ascending) segment (in between 0.4 fcm and fcm) is quadratic [90]:

 ${\displaystyle {\sigma }_{c(2)}=\,{\frac {{E}_{ci\,}\,{\frac {{\epsilon }_{c}}{{f}_{cm}}}-\,{\left({\frac {{\epsilon }_{c}}{{\epsilon }_{cm}}}\right)}^{2}}{1+\left({E}_{ci\,}{\frac {{\epsilon }_{cm}}{{f}_{cm}}}-2\right)\,{\frac {{\epsilon }_{c}}{{\epsilon }_{cm}}}}}{f}_{cm}}$
(54)

In equation (54), Eci is the modulus of deformation of concrete for zero stress, given by ${\textstyle {\,E}_{ci}=}$${\displaystyle 10000\,{{f}_{cm}}^{1/3}}$ and E0 = (0.8 + 0.2 fcm / 88) Eci (in MPa) [90]. In the initial linear branch, E0 is the secant modulus that corresponds to 0.4 fcm stress.

Third (descending) segment is given by:

 ${\displaystyle {\sigma }_{c(3)}=\,{\left({\frac {2+\,{\gamma }_{c}\,{f}_{cm}{\,\epsilon }_{cm}}{2\,{f}_{cm}}}-{\gamma }_{c}{\,\epsilon }_{c}+{\frac {{\,{\epsilon }_{c}^{2}\,\gamma }_{c\,}}{2\,{\epsilon }_{cm}}}\right)}^{-1}}$
(55)
 ${\displaystyle {\gamma }_{c}=\,{\frac {{\pi }^{2}{\,f}_{cm\,}{\epsilon }_{cm\,}\,}{2{\left[{\frac {{G}_{ch}}{{l}_{eq}}}-0.5\,{f}_{cm\,}\left({\epsilon }_{cm\,}\left(1-b\right)+b\,{\frac {{f}_{cm\,}}{{E}_{0\,}}}\right)\right]}^{2}}}}$ ${\displaystyle \qquad b={\frac {{\epsilon }_{c}^{pl}\,}{{\epsilon }_{c}^{ch}}}}$
(56)

In equations (55) and (56), Gch is the crushing energy per unit area [28], and leq is the characteristic length, which depends on the mesh size, the type of finite element and the crack direction [93] [28] [94]. Assuming idealized behavior of single band of cracks, the characteristic length can be determined after the mesh size; [95] relate the crack width with the square root of the finite element area for 2D elements. In this work, brick solid elements are utilized; the characteristic length is taken as average length of the brick element.

Based on experimental observations, b = 0.9 (equation (56)) can be initially assumed. After calculating the damage variables, the average value of b along the relevant strain range is obtained; iterative calculations are performed until reaching convergence. The final value of b affects the softening branch of the compressive stress-strain relation (equations (55) and (56)); therefore, the dissipated crushing energy will be changed. Anyway, this effect is not very intense.

Equation (55) shows that the descending branch approaches asymptotically zero; therefore, a fictitious maximum strain shall be selected for calculation purposes. The maximum strain value shall fulfill that the crushing energy in equation (60) is equal to the area under the corresponding compressive stress-strain law multiplied by the characteristic length.

Regarding tensile behavior, the ratio between tensile stress σt(w) (for crack width w) and maximum tensile strength ftm, is given by [96]:

 ${\displaystyle \,{\frac {{s}_{t}(w)}{{f}_{tm}}}=\,\left[1+\,{\left({c}_{1}{\frac {w}{{w}_{c}}}\right)}^{3}\right]{e}^{-{c}_{2}\,{\frac {w}{{w}_{c}}}}\,-}$${\displaystyle {\frac {w}{{w}_{c}}}\left(1+{c}_{1}^{3}\right)\,{e}^{-{c}_{2}}}$
(57)

In equation (57), c1= 3, c2= 6.93 [96], and wc is the critical crack opening. Equation (57) shows that σt(0) = ftm and σt(wc) = 0. Therefore, wc can be considered as the fracture crack opening. Equation (58) [96] relates wc with the tensile strength and fracture energy GF per unit area:

 ${\displaystyle {\,w}_{c}=5.14\,{G}_{F}\,/{f}_{tm}\,}$
(58)

According to [90], GF (N/mm) can be calculated as

 ${\displaystyle {G}_{F}=0.073\,{f}_{cm}^{0.18}}$
(59)

In equation (59), fcm is expressed in MPa. The ratio between crushing and fracture energies can be assumed proportional to square of the ratio between compressive and tensile strengths [94]:

 ${\displaystyle {G}_{ch}={\left({\frac {{f}_{cm}}{{f}_{tm}}}\right)}^{2}{G}_{F}}$
(60)

In this work, the actual crack spacing is not studied, it has been assumed that there is a single crack per element. This supposition is suitable for global-purpose simulation. After this assumption, in the descending segment of the tensile stress-strain curve (Figure 27.b), the strain can be obtained in terms of the crack opening from the following kinematic relation:

 ${\displaystyle {\,\epsilon }_{t}\,={\epsilon }_{tm}+w\,/\,{l}_{eq}}$
(61)

This subsection describes the concrete uniaxial concrete behavior; next subsection describes their implementation in the proposed approach.

### 3.4 Implementation

Following the formulation described in the previous subsection, coefficients ac, at, bc and bt (subsection ‎3.2) can be determined. Coefficient ac is obtained by replacing fc0 = 0.4 fcm (Figure 27.a) in equation (52). Figure 27.b shows that ftm = ft0, then equation (53) shows that at = 1. Coefficient bc can be determined from equation (46) by replacing fc0 = 0.4 fcm = 0.4 (fck + 8) MPa. Coefficient bt is determined by replacing at = 1 and ${\textstyle {f}_{t0}\,=}$${\displaystyle {f}_{tm}\,=0.3016\,{f}_{ck}^{2/3}}$ [90] in equation (46). Therefore:

 ${\displaystyle {a}_{c}=7.873}$ ${\displaystyle {a}_{t}=1}$ ${\displaystyle {b}_{c}={\frac {1.97({f}_{ck}\,+\,8)}{{G}_{ch}}}{l}_{eq}}$ ${\displaystyle {b}_{t}={\frac {0.453{f}_{ck}^{2/3}}{{\boldsymbol {\,}}{\,G}_{F}}}{l}_{eq}}$
(62)

In the last two expressions in equation (62), fck is expressed in MPa. After parameters ac, at, bc and bt, damage variables dc and dt are determined by equations (47) and (48) in terms of crushing and cracking strain, respectively.

A particular implementation of the proposed methodology by using the previously described concrete constitutive law is described next. All stress values are in MPa.

1. The input data are the concrete compressive strength fck, the parameters in Table 3, the mesh size ${\textstyle {l}_{eq}}$, and the ratio b (eqn. (56)). Initial assumption is b = 0.9

2. Calculate the compressive / tensile stress strength fcm = fck + 8 / ${\textstyle {f}_{tm}\,=}$${\displaystyle 0.3016\,{f}_{ck}^{2/3}}$

3. State the strain at compressive stress strength as εcm = 0.0022

4. Calculate the initial tangent modulus of deformation of concrete ${\textstyle {\,E}_{ci}=}$${\displaystyle 10000\,{{f}_{cm}}^{\frac {1}{3}}}$ and the undamaged modulus of deformation ${\textstyle {\,E}_{0}=}$${\displaystyle {\,E}_{ci}(0.8+0.2{\frac {{f}_{cm}}{88}})}$

5. Calculate the fracture / crushing energy (N/mm) ${\textstyle {G}_{F}=}$${\displaystyle 0.073\,{f}_{cm}^{0.18}}$ / ${\textstyle {G}_{ch}={\left({\frac {{f}_{cm}}{{f}_{tm}}}\right)}^{2}{G}_{F}}$

6. Calculate the critical crack opening ${\textstyle {\,w}_{c}=}$${\displaystyle 5.14\,{G}_{F}\,/\,{f}_{tm}}$

7. Build the first / second / third segments of the concrete uniaxial compressive law: ${\textstyle {\sigma }_{c(1)}=}$${\displaystyle {E}_{0}\,{\epsilon }_{c}}$ / eqn. (54) / eqn. (55). In eqn. (55), strain is bounded; the selected upper bound should fulfill the condition that the crushing energy Gch (equation (60)) is reached

8. Build the first / second segment of the concrete uniaxial tensile law: ${\textstyle {\sigma }_{t(1)}=}$${\displaystyle {E}_{0}\,{\epsilon }_{t}}$ / eqns. (57) and (61)

9. Calculate the damage parameters according equation (62): ${\textstyle {a}_{c}=}$${\displaystyle 7.873}$; ${\textstyle {a}_{t}=1}$; ${\textstyle {b}_{c}=}$${\displaystyle {\frac {1.97({f}_{ck}\,+\,8)}{{G}_{ch}}}{l}_{eq}}$; ${\textstyle {b}_{t}=}$${\displaystyle {\frac {0.453{f}_{ck}^{2/3}}{{\boldsymbol {\,}}{\,G}_{F}}}{l}_{eq}}$

10. Calculate the compressive / tensile damage variables (damage evolution) according to eqns. (47) / (48)

11. Calculate the compressive and tensile plastic strains as indicated in Figure 27: ${\textstyle {\epsilon }_{c}^{pl}\,=}$${\displaystyle {\epsilon }_{c}^{ch}-{\sigma }_{c}{d}_{c}/\left(1-{d}_{c}\right)\,{E}_{0\,}}$ , ${\textstyle {\epsilon }_{t}^{pl}\,=}$${\displaystyle {\epsilon }_{t}^{ck}-{\sigma }_{t}{d}_{t}/\left(1-{d}_{t}\right)\,{E}_{0\,}}$

12. Calculate the average value of ratio ${\textstyle b={\frac {{\epsilon }_{c}^{pl}\,}{{\epsilon }_{c}^{ch}}}}$ and compare with the assumption in step ‎1. Repeat until reaching convergence.

Once convergence is reached, the major output are curves of compressive / tensile stress and damage variables vs. crushing / cracking strain, respectively. These plots constitute the input of the implementation in any model describing the global structural behavior. This algorithm is suited for software package ABAQUS [89] [97] [98]. Noticeably, the proposed approach can be also fitted for other computer codes and for other concrete formulations, either empirical (code-type) or directly based on experiments. A flow chart of the methodology implementation is displayed in Figure 28.

 Figure 28. Flow chart of the proposed methodology implementation

### 3.5 Mesh insensitivity verification

When the material exhibits softening, finite element size influences significantly the entire model behavior due to localization since the dissipated energy decreases upon mesh refinement. This can be solved by the so-called mesh regularization techniques. One of the simplest remedy is the Crack Band Method; it uses energy-based scaling of the softening part of the stress-strain relation [99]. This technique can be considered as a simple and effective method for practical engineering analysis and has been implemented in many concrete plastic damage models, among others [21] [78]. More advanced techniques are the so-called Non-Local approaches; they are based on introducing non-locality in the constitutive model. This non-locality can be incorporated into an integral format [100] [100]; this approach has been implemented in different plastic damage models [101] [26]. Another strategy for incorporating non-locality is including higher-order deformation gradients in the model [102]. This has been implemented in several plastic damage models [103] [104].

As the model of [30] is the base of this work and this model uses the fracture energy-based regularization to describe the softening in tension and compression, the crack band technique is considered. Then both fracture and crushing energies are scaled in relation to the finite element size. The main assumption in this technique is that damage is localized in a single raw of elements; this is true for the physical failure mechanism of concrete in tension but is not exactly right for compression. Uniaxial compression tests show a great relation between boundary conditions, concrete strength, size, and failure mechanism [105] [106]. It varies between cone failure and vertical splitting modes. Cone failure occurs when friction at top and bottom sections restraints lateral expansion; this mode is characterized by inclined slip bands at each corner and crushing in the middle zone for very stocky specimens. Vertical splitting mode can happen for highly slender specimens when the effect of friction is neglectable. Intermediate scenarios can occur for different degrees of friction and slenderness; they are characterized mainly by inclined shear band zone. Regarding damage localization under compression, experimental observation by [105] [92] [107] showed strain localization; analytical investigation by [108] showed that the concept of fracture energy in tension holds true for compression. This gives the evidence to use the fracture energy-based regularization in describing the softening of concrete under compression [109].

In this subsection the capacity of proposed methodology to overcome mesh insensitivity is verified for a uniaxial tension problem. Given the above considerations, it is apparent that the uniaxial compression case cannot be verified in the same way. For practical engineering applications, behavior is mainly controlled by bending and there will be mixed stresses states with high stress gradients; therefore, the possible inaccuracies will have little effect on overall results.

This subsection describes the application of the proposed methodology to the problem of uniaxial tension of a 200 mm × 200 mm × 200 mm cube, as described in Figure 29. The cube is discretized with a uniform mesh of 3D 8-node hexahedron solid finite elements (C3D8R). Figure 29 displays three sketches of the cube discretized with a coarse mesh (leq = 200 mm, one element), a medium mesh (leq = 50 mm, 64 elements) and a fine mesh (leq = 25 mm, 512 elements), respectively. The objective of this analysis is to verify the allegedly low sensitivity to mesh size.

This problem is analyzed with the algorithm described in subsection ‎3.4. Table 4 displays the values of the parameters for each mesh size. Any of the values of parameter b in Table 4 is the result of an independent iterative process starting from 0.9. For consistency, the final value is independent on this initial assumption.

The algorithm is implemented in ABAQUS code [89]. The displacement is applied incrementally; for coarse / medium / fine meshes, the maximum number of iterations is 4, 6 and 5, respectively.

 (a) Coarse mesh (200 mm) (b) Medium mesh (50 mm) (c) Fine mesh (25 mm) Figure 29. Uniaxial tension example

Tabla 4. Parameters for the mesh-insensitivity verification example
 Mesh leq fck fcm ftm Gch GF b ac at bc bt (mm) (MPa) (MPa) (MPa) (N/mm) (N/mm) Coarse 200 25 33 2.58 22.43 0.137 0.6 7.873 1 581 5648 Medium 50 25 33 2.58 22.43 0.137 0.914 7.873 1 145.2 1412 Fine 25 25 33 2.58 22.43 0.137 0.967 7.873 1 72.6 706

Figure 30 presents the plots that constitute the major inputs for this implementation. Figure 30.a and Figure 30.b display plots of compressive/tensile stress vs. crushing/cracking strain, respectively; Figure 30.c and Figure 30.d display plots of compressive/tensile damage variable vs. crushing/cracking strain, respectively. Figure 30 shows that the inputs are strongly dependent on the mesh size; conversely, outputs are expected to be almost independent on it. Noticeably, compressive plots are displayed only for information.

 (a) Compressive stress vs. crushing strain (b) Tensile stress vs. cracking strain (c) Compressive damage variable vs. crushing strain (d) Tensile damage variable vs. cracking strain Figure 30. Results of the proposed algorithm for the uniaxial tension example using coarse, medium and fine meshes

Given that the problem under consideration is extremely simple, a closed-form solution can be provided. This solution is derived from the uniaxial tensile constitutive law given by equation (57) (Figure 27.b), with the critical crack opening wc obtained from equation (58) (wc = 0.273 mm).

After the implementation in ABAQUS of the plots displayed in Figure 30, Figure 32 presents the force- displacement plots that are obtained for each mesh size. For comparison purposes, Figure 32 contains also the aforementioned closed-form solution. Figure 32 shows a satisfactory agreement among the plots for the three mesh sizes; this observation corroborates the mesh-insensitiveness of the proposed methodology. As well, numerical results match analytical ones, thus confirming the accuracy. Noticeably, the best match is provided by the coarse mesh; this circumstance can be explained because a single element is the closest representation of the closed-form solution. On the other hand, the result of the fine mesh is better that the one of the medium mesh; certainly, further refinements would lead to higher accuracy.

 (a) Coarse mesh (200 mm) (b) Medium mesh (50 mm) (c) Fine mesh (25 mm) Figure 31. Uniaxial damage localization

 Figure 32. Force-deflection plots for the uniaxial tension example (coarse, medium and fine meshes)

The area under the plots in Figure 32 (between the peak and the displacement for zero stress in the closed-form solution) is equal to 5.5149, 5.6358 and 5.3965 Nm for coarse, medium and fine meshes, respectively. These quantities represent the fracture work; it can be converted into fracture energy per unit area (GF) dividing by the area perpendicular to the applied displacement following the regularization method [99]; results are: 0.1379, 0.1409 and 0.1349 N/mm for coarse, medium and fine meshes, respectively. Comparison with the values indicated in Table 4 shows a satisfactory agreement.

### 3.6 Experimental validation of damage evolution

Accuracy and reliability of the proposed methodology are verified by simulating one experiment with the particular algorithm described in subsection ‎3.4 and with the formulation to determine the damage variables described in subsection‎3.1. The experiment [110] consisted in imposing a cyclic force law to a number of plain concrete cylinders. Figure 33 displays the imposed force law (in terms of stress).

 Figure 33. Input force law for the plain concrete test

This experiment had been previously simulated [80] [111]. Analog to the mesh sensitivity verification example, the proposed approach has been implemented in ABAQUS FEM program.

The simulated test refers to a cube element with 125 mm side; it is discretized with a single element. The maximum number of iterations is 2, and loading increment is 4 × 10-4 mm. Table 5 displays the values of the parameters.

Tabla 5. Parameters for the simulation of the plain concrete test
 fck (MPa) fcm (MPa) ftm (MPa) Gch (N/mm) GF (N/mm) b ac at leq (mm) bc bt 18 26 2.07 20.7 0.1312 0.8 7.873 1 125 310.48 2960

Figure 34 displays major outputs of the simulations, compared with the corresponding experimental results. Figure 34.a exhibits stress-strain plots and Figure 34.b shows the evolution of damage variable (d) in terms of strain. In points “▲” in Figure 34.b, the abscissa corresponds to measured values of strain and the ordinate corresponds to damage determined applying equation (23) to the experimental plots in Figure 34.a, as average for each pair of unloading and reloading branches. In the numerical plots in Figure 34.b, the horizontal segments belong to constant damage along unloading and reloading branches.

 (a) Stress-strain (b) Damage-strain Figure 34. Comparison between experimental and simulated results for the plain concrete test

Figure 34.a shows a satisfactory agreement between numerical and experimental results, particularly in terms of envelope curve; noticeably, envelope results from both models are almost identical (because the same uniaxial constitutive laws have been considered). Figure 34.b shows an adequate agreement between the experimental damage (at the measured strains) and the numerical damage calculated with the proposed approach, given that points “▲” are highly close to the numerical curve.

## 4 Reinforcing bars modeling

### 4.1 Introduction

The behavior of reinforcing steel may control the response of reinforced concrete structural elements subjected to earthquake loading. Thus, it is necessary to use an appropriate model that predicts the fundamental characteristics of steel within a range of loading.

This section presents the analytical model used in this work to simulate the linear and nonlinear behavior of reinforcing bars. 2D truss finite element is used in this study to model the reinforcing bars, the model assigned to such element is used to perform the simulation of different reinforced concrete elements subjected to monotonic and cyclic loading.

A rate-independent elasto-plasticity model with isotropic hardening has been used for the simulation of monatomic loading, in the other hand a combined isotropic-kinematic hardening has been used for the simulation of cyclic loading. The available steel plastic models in ABAQUS [89] has been used to simulate the behavior of steel reinforcement. In these models the total stain is written in terms of elastic and plastic strain rates as:

 ${\displaystyle {\overset {\cdot }{\boldsymbol {\varepsilon }}}\,={\overset {\cdot }{{\boldsymbol {\varepsilon }}^{e}}}+}$${\displaystyle {\overset {\cdot }{{\boldsymbol {\varepsilon }}^{pl}}}}$
(63)

The elastic behavior follows Hook's law and is written as:

 ${\displaystyle {\boldsymbol {\quad \quad \,\,}}{\overset {\cdot }{\boldsymbol {\sigma }}}=\,{\boldsymbol {D}}_{0}^{el}\,:{\overset {\cdot }{{\boldsymbol {\varepsilon }}^{e}}}=}$${\displaystyle {\boldsymbol {D}}_{0}^{el}\,:\left({\overset {\cdot }{\boldsymbol {\varepsilon }}}-{\overset {\cdot }{{\boldsymbol {\varepsilon }}^{pl}}}\right)}$
(64)

${\textstyle {\boldsymbol {D}}_{0}^{el}}$ is the elastic stiffness tensor, ${\textstyle {\boldsymbol {\sigma }}}$ and ${\textstyle {\boldsymbol {\varepsilon }}}$ are the stress and strain tensors.

These models use the Von Mises yield condition with an associated plastic flow. Hence, the yield function and plastic potential are defined by the same function. The pressure-independent yield function is defined by the function

 ${\displaystyle F({\boldsymbol {\sigma }},{\boldsymbol {\alpha }})={\sqrt {{\frac {3}{2}}\left({\boldsymbol {\sigma }}^{'}-{\boldsymbol {\alpha }}^{'}\right):{\left({\boldsymbol {\sigma }}^{\boldsymbol {'}}-{\boldsymbol {\alpha }}^{'}\right)}}}\,-}$${\displaystyle \,{\sigma }_{y}}$
(65)

where σ' and α' are the deviatoric parts of the stress and backstress tensors σ and α respectively, σy is the uniaxial yield stress.

As it's mentioned above, an associated plastic flow is used. Therefore, as the material yields, the inelastic deformation rate is in the direction of the normal to the yield surface and there's no volumetric plastic strains. This assumption is generally acceptable for application of reinforced concrete cyclic simulation as long as microscopic details, such as localization of plastic flow occurring as a metal component ruptures due to cyclic fatigue loads, are not of interest.

The plastic flow is given by:

 ${\displaystyle {\boldsymbol {\quad \quad \,\,}}{\overset {\cdot }{{\boldsymbol {\varepsilon }}^{pl}}}=}$${\displaystyle \,{\overset {\cdot }{\overline {\varepsilon }}}^{pl}{\frac {\partial F}{\partial {\boldsymbol {\sigma }}}}}$
(66)

where ${\textstyle {\overset {\cdot }{\overline {\epsilon }}}^{pl}}$ is the equivalent plastic strain rate.

The nonlinear isotropic-kinematic hardening model is based on the work of Lemaitre and Chaboche (1990) [69]. In this model, the evolution law consists of two components; Nonlinear kinematic hardening describes the translation of the yield surface in stress space through the backstress tensor α, the isotropic hardening describes the change of the equivalent stress that defines the size of the yield surface σ0 as a function of plastic deformation. The nonlinear kinematic hardening component is defined to be two terms; Linear kinematic hardening and relaxation term which introduces the nonlinearity, nonlinear kinematic hardening law is defined in the next formulation:

 ${\displaystyle {\overset {\cdot }{\mathit {\boldsymbol {\alpha }}}}_{}={C}_{m\,}{\frac {1}{{\sigma }_{y}}}\left({\boldsymbol {\sigma }}-\right.}$${\displaystyle \left.{\boldsymbol {\alpha }}\right){\overset {\cdot }{\overline {\epsilon }}}^{pl}-}$${\displaystyle {\gamma }_{m\,}{\mathit {\boldsymbol {\alpha }}}_{k\,}{\overset {\cdot }{\overline {\varepsilon }}}^{pl}}$
(67)

where ${\textstyle {C}_{m\,}and\,{\gamma }_{m\,}}$are material parameters that can be calibrated form cycle test data, ${\textstyle {\overset {\cdot }{\overline {\epsilon }}}^{pl}}$ is the equivalent plastic strain rate. When ${\textstyle {\gamma }_{m\,}}$is zero, the model recovers to Linear kinematic hardening.

The isotropic hardening component of the model defines the evolution of the yield surface size, this evolution can be defined by a simple exponential law

 ${\displaystyle {\sigma }^{0}={\sigma }_{y}+{Q}_{\infty }\left(1-{e}^{-{b}_{m}{\overline {\varepsilon }}^{pl}}\right)}$
(68)

where ${\textstyle {Q}_{\infty }\,and\,{b}_{m}\,}$ are material parameters and can be calibrated from test data.

The required data for the kinematic hardening component have been obtained from the results of stabilized cycle simulation. In this work the FEM package seismostruct [112] has been used to obtain the cyclic behavior of reinforcing bars subjected to symmetric strain cycles. A stabilized cycle is obtained by cycling a specimen over fixed strain range Δε until the stress strain curve is no longer changes from one cycle to the next. Figure 35 shows the strain stress curve for a stabilized cycle according to [89], each pair of ( ${\textstyle {\sigma }_{i},{\epsilon }_{i}^{pl})}$ is specified after shifting the strain axes to ${\textstyle {\epsilon }_{p}^{0}}$.

 Figure 35. Stress-strain data for a stabilized cycle [89]

The shifted data are entered in tabular form in ABAQUS [89] as stress-strain, more detailed information about the calibration and the calculation of the back stresses can be found in [89].

Another possibility to simulate the nonlinear behavior of steel bars under cyclic loading is to use a linear-kinematic hardening model instead of combined hardening "isotropic and kinematic"; such model requires the hardening parameter to be defined. Hardening parameter is defined as the difference between the yield stress and higher stress value divided over the plastic strain "strain difference between the total strain correspond to the higher stress and the yield strain value".

## 5 Bond-slip effect modeling

### 5.1 Introduction

This section provides a general overview of the fundamental mechanisms governing the bond behavior between reinforcing bars and the surrounding concrete. A summary about the recent and past studies and models is also provided. Finally, a new implementation scheme is presented, this scheme is particular appropriate to implement the interface bond model developed by [20] in 3D FEM simulation.

The stress transfer capacity between concrete and reinforcing bars is generally referred as the bond of reinforcing bars. This phenomenon is an essential mechanism that engages both material generating the composite action of reinforced concrete structures. Bond of deformed bars (not plain) has been studied and characterized by large number of researchers, which has led to be well documented in design and modeling codes [113] and [114]. Transfer of forces from the reinforcement to the surrounding concrete for a deformed bar shown in Figure 36 occurs by :

• Chemical adhesion between the bar and the concrete;
• Frictional forces arising from the roughness of the interface, forces transverse to the bar surface, and relative slip between the bar and the surrounding concrete.
• Mechanical anchorage or bearing of the ribs against the concrete surface.
 Figure 36. Source of bond resistance in a deformed bar [114]

When deformed bar moves with respect to the surrounding concrete (slip), the chemical surface adhesion is lost while bearing and friction forces are mobilized; friction forces are mobilized due to bar roughness and bearing forces at the ribs caused by the wedding action against the concrete. The pressure exerted by ribs onto concrete creates micro crack referred as Goto cracks [115] starting at the tip of the ribs and propagating transversally away from the bar as shown in Figure 37. With increasing the demand and the slip, the wedging action of the ribs tends to introduce a radial expansion at the interface, which activates the passive confinement in the concrete. Radial expansion produces a hoop expansion in the concrete, which causes splitting cracks to develop at the surface in contact with the bar and propagate radially, as shown in Figure 37. This hoop expansion is restrained by the undamaged outer concrete ring aswell as the confining reinforcement if any. For low confinement conditions, splitting cracks propagate radially through the concrete cover and the bond fails abruptly. This type of failure is referred as splitting failure and can be represented by vertical drop of the bond stress at point A in Figure 38.a, the horizontal axes in Figure 38 the is referred to the relative displacement "slip".

 Figure 37. Goto cracks due to bond slip [114]

When sufficient confinement and cover are provided, the bond resistance will be increased due to introducing large radial stresses at the contact between the concrete and the bars. Further slip occurs after crushing the concrete front of the ribs, crushed particles of concrete contributes to increase the radial component of the bearing forces and lateral shear crack will be generated in as shown in Figure 38.b. At this point the slope of bond-slip decrease rapidly (approximately at point B in Figure 38.b). If the confinement of the concrete is not sufficient, splitting failure can happen, conversely if the concrete is well confined, higher bond strength can be achieved till reaching the maximum bond-strength at point C in Figure 38.c, at this state a part of or the total concrete keys will be crashed and sheared. Bond stress starts to decrease with increasing the slip and more concrete will be crushed, failure occurs due to loss of interlocking action caused by crushing and shearing the concrete keys between ribs. Eventually, the bar is pulled out from the concrete and only frictional resistance remains at point E in Figure 38.c, this type of failure is referred as pull-out failure. Figure 38 shows a representation of the damage evolution for pull-put failure type under monotonic loading.

 Figure 38. Monotonic bond-slip behavior for pull-out failure [116]

Bond-slip behavior and damage mechanism under cyclic loading are even more complex, the hysteretic behavior is well studied and theorized by [116] as shown in Figure 39.

 (a) (b) (c) Figure 39. Cyclic bond-slip behavior for pull-out failure [116]

Cyclic behavior for low demand reversal before developing shear crack is presented in Figure 39.a. After unloading at point A along path AF, the gap between the right side of the ribs and the adjacent concrete, after concrete crushing on the left side of the ribs, remains open with a width equal to the residual slip at point F. Only a small fraction of the slip is recovered by the elastic unloading of the concrete. When the slip is reversed along path GH, some frictional resistance is built up. At H, the ribs are in contact again with the concrete but a gap has opened on the left side of the ribs. Because of a resumed contact with the concrete, a sharp increase in stiffness occurs along path HI. With increasing load, the old inclined cracks close, allowing the transfer of compressive stresses across them with no noticeable reduction in stiffness. Inclined cracks perpendicular to the old ones appear as the stress increases in this direction. At point I, a gap equal to the distance between points F and I has opened. When reversing the slip, the path IKL is similar to AFH. However, the bond resistance starts to increase again at L, when the ribs start to press broken pieces of concrete against the previous bearing face. With further movement, the transverse cracks previously closed are opened and the cracks previously opened are closed. At M, the ribs and the concrete are in full contact and the monotonic loading curve is recovered.

If the slip reversal occurs after horizontal shear cracks have initiated, a different behavior is obtained as shown Figure 39.b. When loading in the opposite direction (along path HI), the ribs press against the concrete in between whose resistance has been lowered by the shear cracks. Therefore, the bond resistance is lowered compared to the monotonic curve. When reversing the slip again (along path IKLMN), the resistance is further lowered compared to that at point I because of the additional shearing damage in the concrete. When a large slip is imposed during the first cycle, almost all the concrete between the ribs can be sheared off and the behavior will be like the one shown in Figure 39.c. When moving the bar back (along path GH), the frictional resistance is higher than that for the previous cases, in which the slip in the first cycle is smaller, because the concrete surface along the shear crack is rougher. When reloading in the opposite direction, the peak resistance (point I) is lowered. When reversing the slip again, the frictional resistance is lowered because the surface has been smoothened (path KL).

From the simulation point of view, a number of numerical models have been proposed to simulate bond-slip behavior; they can be classified depending on their scale [117].

Rib scale. The detailed geometry of the interface, including bar ribs, is modelled, concrete and steel are discretized with continuum elements [118] [119]. These models can provide very high accuracy, although with a high computational cost.

Bar scale. Ribs are not directly included in the simulation; the contact surface is idealized as smooth. Concrete and steel can be discretized with different types of elements (e.g. solid for concrete and truss or beam for steel). The interaction forces (cohesion, friction and bearing) are represented by tangential (bond) and normal stress; particular attention is usually paid to the tangential component. Bond-slip is represented through relation between bond stress and bar slip implemented in zero-thickness interface element. [120] developed a four-node zero-thickness bond-slip element to be used for two-dimensional modelling; other studies using this model have been reported. Broadly speaking, bar scale models have good balance between computational cost and ability to reproduce accurately bond behavior.

Structural Element scale. The bond-slip effects are either incorporated in the formulation of the element or taken into consideration through zero-length springs. [121] introduced the beam-column model with bond-slip proposed by [122] into the force-based fiber-section element developed by [123]. In regards to zero-length approach; [124] proposed a law to relate bar stress and slip at end of anchorage in footing-column or beam-column connections; this law has been calibrated with experimental results and used as a constitutive relation for steel fibers in a zero-length fiber-section element to simulate end rotation of RC columns. These models have low computational cost and, therefore, are suitable for simulation of full-scale structures; conversely, they do not study explicitly the bond behavior, merely reflect the additional flexibility provided by bond-slip.

Rib scale models are too computationally expensive for the practical simulation of full-scale actual structures, and element scale models can be an attractive option for large-scale simulation but cannot capture strength reduction due to bond degradation (because of slip, concrete spalling or steel yielding) and cannot be implemented in continuum elements. Bar scale models are the most appropriate for implementation in continuum elements under seismic excitation, as they represent efficiently interaction between two materials and predict different deterioration modes. A number of models of this type have been proposed:

• [Gan 2000] [125]. A double-node zero-thickness contact element is used to connect steel bar elements with the surrounding concrete. The monotonic and cyclic behavior of this element is governed by a bond-slip law originally proposed by [116] and later modified by [125] according to experimental observations. This contact element was implemented in Finite Element Code TRIX99 and used to perform 2-D cyclic analysis for shear walls.
• [Ožbolt et al. 2002] [126]. These researchers proposed a monotonic bond-slip law based on experimental results [116] and [127]. This law was implemented in the finite element code MASA using a zero-width element. Concrete and steel are discretized with solid and truss elements, respectively. This model was been verified with pull-out tests and simply-supported RC beams.
• [Lowes et al. 2004] [120]. A four-node zero-thickness bond-slip element is used for two-dimensional monotonic and cyclic modelling. The model is characterized by a normalized bond stress-slip law and a relation between the maximum bond strength and stress-strain state of concrete and steel.
• [Rabczuk et al. 2005] [128]. A 2-D two-double-node interface element, bond behavior is described in terms of radial stress-strain relation. This element was implemented in FE package ABAQUS using user subroutine; it was used in the simulation of a prestressed RC beam subjected to static load.
• [Murcia-Delso, Shing 2015a] [129]. These researchers developed a general 3-D bond-slip model accounting for bond deterioration due to cyclic slip reversals, concrete splitting, and bar yielding in tension. This model is an extension, to represent different bond deteriorations, of the semi-empirical cyclic bond-slip law proposed by [20] and [130] which is similar to that proposed by [116]. This model was implemented in a four-node zero-thickness interface element with linear shape functions and two integration points.

### 5.2 Spring model

This subsection describes a simple and practical model that can be used only for monotonic loading simulation. The model can be implemented in any FEM package that provides nonlinear 3D spring element with ability to define the constitutive force displacement law.

The relation between concrete and steel bars is represented by four nonlinear spring that connect the nodes of truss elements representing reinforcement bars from one side, and the nodes of 3D solid elements representing concrete material form the other side. Only the longitudinal component of bond is assigned to the springs, while the normal and rotation stiffness are ignored. Figure 40.a.b show respectively a 3D and 2D representation of the model. In Figure 40, A and B are the start and end nodes of the truss element, "Spr" refers to spring. Under well confined conditions it can be assumed to consider average "local bond" versus "local slip" relation between reinforcement and concrete, this relation can be described by bond-slip law. [90] provides detailed monotonic bond-slip law for pull-out and splitting failure as a function of bond stress τ2 and the relative displacement S2 as shown in Figure 41. Once the bond stress-slip relation is defined, it can be converted to force-slip law and attached to each spring as its constitutive law.

 (a) (b) Figure 40. Implementation of Spring model for bond-slip [129]

In Figure 41; τf is the friction bond resistance, SR is the distance between the bars ribs.

Bond-slip law is converted to force displacement relation and assigned as a constitutive law for the spring elements.

 Figure 41. Monotonic bon-slip law [90]

For uniform mesh size discretization for truss and solid elements as shown in Figure 40, the force in one spring can be calculated according to Eq (69).

 ${\displaystyle {F}_{bs}=\,1.36\,{D}_{b}L{\tau }_{2}}$
(69)

Fbs is the equitant bond force in the spring, Db is the bar diameter, L is the mesh size.

What makes this model useful only for monotonic loading simulation is the hardening and plasticity evolution of spring element. As was described in the general explanation of bond-slip in subsection ‎5.1, the cyclic behavior of bond stress-slip is not possible to be described by standard isotropic hardening model This might require special treatment and developing of special constitutive law to cover the nonlinear behavior of this element.

### 5.3 Interface model

The 3-D interface model proposed by [129] is used in this research to represent the bond relation between the longitudinal reinforcement bars and the surrounding concrete. This model is assigned to a 2-D interface element connecting steel and concrete representing a fraction of the perimeter of the bar. Their interaction is represented in the longitudinal direction by the relation between an equivalent bond stress and the bar slip. The stress-slip constitutive law uses the phenomenological law proposed by [130] which is based on concepts originally developed by [116]. This law was extensively verified by pull-out tests on bars embedded in well confined concrete.

Figure 42 displays the three components of interaction stress (σ1, τ2, τ3) and relative displacement (S1, S2, S3) during sliding.

 Figure 42. Stress and relative displacement components at bar-concrete interface [129]

#### 5.3.1 Bond-stress vs slip law

The bond resistance τ2 in the model proposed by [129] is decomposed into bearing and friction components:

 ${\displaystyle {\tau }_{2}=\,{\rho }_{n}\left({\rho }_{b,s}\,\,{\rho }_{b,c\,}{\tau }_{b}+\right.}$${\displaystyle \left.{\rho }_{f,s}\,{\rho }_{f,c\,}{\tau }_{f}\right)}$
(70)

In equation (70), τb and τf are the full bearing and friction resistance of an elastic bar under monotonic pull-out action. Such resistances are multiplied by reductions factors. ρb,s and ρf,saccount for the influence of yielding of bar in tension,ρb,c and ρf,c represent the effect of slip history, and ρn takes into consideration splitting cracks and is a function of S1.

Figure 43 shows the considered monotonic and cyclic bond stress versus bar slip curves. These curves are defined piecewise using polynomial functions [129] in terms of three parameters: the peak bond strength τmax, the slip at the pick bond strength Speak and the clear spacing between the bar ribs SR.

 Figure 43. Bond stress-versus-slip law: (a) monotonic response; (b) cyclic response [129]

One of the main features of this model is the consideration of the reduction of bearing and friction resistance because of the diminution of the bar contact area due to yielding in tension [131] and [132]. Reduction factors that account for steel yielding in tension are expressed as follows:

 ${\displaystyle {\rho }_{b,s}({\epsilon }_{s})=\left\{{\begin{matrix}1\quad \quad \quad {\rm {for}}\quad {\epsilon }_{s}\leq {\epsilon }_{y}\\{\frac {{\epsilon }_{sh}-{\epsilon }_{s}}{{\epsilon }_{sh}-{\epsilon }_{y}}}\quad \quad \,\,{\rm {for}}\quad {\epsilon }_{y}<{\epsilon }_{s}<{\epsilon }_{sh}\,\,\\0\quad \quad \quad {\rm {for}}\quad {\epsilon }_{s}>{\epsilon }_{sh}\end{matrix}}\right\}}$
(71)
 ${\displaystyle {\rho }_{f,s}({\epsilon }_{s})=\left\{{\begin{matrix}1\quad \quad \quad \quad \,\,{\rm {for}}\quad {\epsilon }_{s}\leq {\epsilon }_{sh}\\{\frac {{\epsilon }_{u}-{\epsilon }_{s}}{{\epsilon }_{u}-{\epsilon }_{sh}}}\quad \quad {\rm {for}}\quad {\epsilon }_{s}>{\epsilon }_{sh}\,\,\end{matrix}}\right\}}$
(72)

In equations (71) and (72), εy, εsh andεu are yielding, hardening and ultimate strains, respectively.

Reduction factors due to the slip history are given by:

 ${\displaystyle {\rho }_{b,c}=1.2{\,e}^{-2.7{\left({\frac {{\overline {S}}_{max}}{{S}_{R}}}\right)}^{0.8}}\leq 1\,}$
(73)
 ${\displaystyle {\rho }_{f,c}=1-min\left({\frac {{S}_{max}^{+}+{S}_{max}^{-}}{{S}_{R}}},1\right)\left(1-\right.}$${\displaystyle \left.{e}^{-0.45{\left({\frac {{S}_{cum}}{{S}_{R}}}\right)}^{0.75}}\right)}$
(74)

In equations (73) and (74), ${\textstyle {S}_{max}^{+}}$ and ${\textstyle {S}_{max}^{-}}$ are the maximum absolute values of slip in positive and negative directions, respectively. ${\textstyle {S}_{cum}}$ is the cumulated slip after slip exceeds Speak for the first time. ${\textstyle {\overline {S}}_{max}}$ is equal to:

 ${\displaystyle {\overline {S}}_{max}=0.75\mathrm {max} \,\left({S}_{max}^{+},{S}_{max}^{-}\right)+}$${\displaystyle 0.25\,\left({S}_{max}^{+}+{S}_{max}^{-}\right)}$
(75)

In this work, the normal displacement S1 is assumed to be very small due to the well confined conditions in the region where it is expected that slip occurs. Therefore, the reduction factor that accounts for splitting cracks is equal to 1: ρn = 1.

The interface element has two additional components perpendicular to the bar longitudinal axis, one normal and one transverse tangential. The stress-displacement relations in the normal direction represent the splitting stresses introduced by the wedging action of the bar ribs. Assuming that the resultant bond force has a fixed angle of inclination of 60º as proposed in [129], the normal stress is proportional to the bond stress. For the transverse tangential direction, a penalty stiffness is introduced to restrain the rotation of the bar about its longitudinal axis.

### 5.4 Proposed Implementation Approach

#### 5.4.1 General Consideration

Next subsection presents a new modeling scheme to implement the bond-slip model described in the previous subsection "developed by [129] and [133]" in a FEM analysis.

#### 5.4.2 Modeling Scheme for Bond-Slip

The bond-slip model was originally implemented in an interface element in the FEM package ABAQUS using the user element subroutine UEL. This interface element can connect truss or beam elements representing reinforcing bars with solid elements representing surrounding concrete, as shown in Figure 44a. The modeling scheme proposed by [129] required the same discretization for steel, interface and concrete elements. This technique is efficient when longitudinal reinforcement consists of parallel bars, but can complicate the mesh when there are crossing bars with non-parallel directions. Figure 44b describes the basis of the new modeling scheme proposed herein.

 (a) Modeling scheme used in [129] (b) Proposed modeling scheme Figure 44. Modeling schemes with bond-slip interface model

Figure 44.b shows that the proposed technique is based on assigning, for each interface element, a pair of tie-constraint conditions; this being defined as two surfaces having the same degrees of freedom in all the contact points. The softest and stiffest surfaces are slave and master, respectively. The first tie-constraint (1) connects concrete nodes (master surface) with the interface nodes on the concrete side (slave surface). The second constrain (2) ties bar nodes (master surface) and the interface nodes on the bar side (slave surface). Hence, the relative displacement between the concrete and the bar is equal to the relative displacement at the interface element.

Figure 45 illustrates the modeling scheme for two perpendicular crossing bars. Figure 45.b and Figure 45.c show, respectively, that Bar 1 (CD) and Bar 2 (AB) are connected to four interface elements each.

 (a) Two bars crossing within a solid element (b) Four inteface elements attached to Bar 1 (CD) (c) Four inteface elements attached to Bar 2 (AB) Figure 45. Application of the proposed modeling scheme for two orthogonal longitudinal crossing bars

By using tie constraints, the proposed technique does not require that interface and solid elements have the same discretization. This strategy aims to overcome meshing problems and convergence difficulties. It is particularly suitable for regions with non-parallel longitudinal reinforcement bars, such as beam-column and slab-column joints, structural walls, foundations, and connections among these elements.

## 6. Finite element analysis verification examples

### 6.1 Integrating CPDM and bond-slip

In the next two subsections the proposed model is verified with large scale experiments, this involves a verification of 2D RC frames and bridge column subjected to monotonic and cyclic loading actions. CPDM described in subsection 2.4 and the proposed methodology presented in section 3 are used for the simulation of monotonic loading tests, while perfect bond conditions are assumed for such type of loading. Bond-slip model described in subsection 5.3 and its implementation scheme proposed in subsection 5.4 are implemented in the simulation of cyclic loading tests.

Bond-slip effects and damage evolution can be considered as essential issues in RC structures when they are subjected to cyclic loading. In this context, a few researchers have attempted also to integrate concrete damage-plastic behavior and bond-slip of reinforcement for practical 3-D continuum simulation of RC structures subjected to seismic loads:

• [Car et al. 2002] [60]. These researchers proposed a constitutive model based on assuming that stress transfer from concrete to steel depends on concrete damage and steel strength. Slip is represented by an irrecoverable inelastic steel strain.
• [Deaton 2013] [134]. Proposed a general model to study the 3-D cyclic behavior of non-seismic RC beam-column joints. Concrete behavior was represented by a damage model originally developed by [135] and implemented in the FE analysis program DIANA. Bond-slip was described with an interface model “line-to-solid” that follows a built-in stress-slip hysteretic law. Verifications with experimental results showed the inability of the bond model to represent the actual hysteretic behavior; concrete model showed good capacity to represent the capacity envelope but not the dissipated energy.
• [Ali et al. 2013] [136]. Studied the cyclic behavior of composed concrete–steel shear walls. The original CDPM that is implemented in ABAQUS was used to simulate concrete behavior; bond-slip is described with contact conditions. Comparison with experimental results, pointed out the lack of capturing the pinching effect (due to bond degradation and cracks opening-closing) and the strength degradation.
• [Murcia-Delso 2013] [20]. The 3-D bond-slip model was used together with the original CDPM that is implemented in ABAQUS. It is calibrated with pull-out tests and utilized to simulate 3-D cyclic behavior of RC columns supported on enlarged pile shaft foundations; satisfactory agreement was obtained. Discrete cracks were introduced at column base to simulate opening and closing of wide flexural cracks at this location and circumvent limitations of the CDPM in ABAQUS to handle the large stiffness degradation required to accurately represent cyclic response of cracks.

New smeared strategy to simulate the reclosing of cracks after tension-compression reversals is proposed in this work. As discussed in subsection 2.4 the reclosing of cracks after tension-compression reversals is governed by the parameter hc, representing the percentage of compression stiffness recovery in the reclosed cracks. Low values of ℎc reduce the concrete strength after reversing; subsequently, the relative displacement between concrete and steel will be higher, thus simulating bond degradation and poor quality concrete.

In actual structures, highly uneven behavior is expected, due to the irregular distribution of bond-slip effects and to the poorer performance in the casting joints. For both reasons, in any framed structure, significantly less reclosing of cracks is expected at the member ends and the nearby segments. To enhance the opening-closing simulation of widely opened cracks in the casting joint region, [20] introduced discrete cracks at the member ends. In the current study, an alternative solution is proposed by using smaller values of hc near the member ends. Figure 46 displays a frame modelled with different values of hc; in the darker shadowed zones, little or no stiffness is expected upon cracks reclosing.

 Figure 46. Zones of framed structures with different cracks opening-closing behavior

In practical applications, the darker shadowed areas in Figure 46 can correspond to the first and last rows of finite elements in the discretization of each column; the lighter shadowed segments can comprise the zones where bond-slip is modelled.

#### 6.2.1 One floor RC frame Pires 1990

This experiment [137] is a quasi-static test consisting of imposing a displacement law to a laboratory single-span, single-story 2-D RC frame. Figure 47.b describes the tested frame, and Figure 47.a and Figure 47.c display beam and columns sections, respectively. Figure 47.b shows that both columns were loaded with constant forces and that displacement was imposed to the top left joint. Noticeably, as in the first frame test, given the absence of distributed loads on beams, there was no initial cracking.

     (a) Beam section   (b) Tested frame (c) Column section Figure 47. [Pires 1990] frame experiment [138] [139]

Mechanical parameters of materials are based on nominal values. The characteristic value of the concrete compressive strength is 20 MPa (C20/25, [140], and the steel yield point is 400 MPa for the longitudinal reinforcement and 500 MPa for the stirrups [141]. As described in Figure 47.c, in the critical end segments (“confinement sections”), closer stirrup spacing was used; the lengths of these segments are 40 cm in beam and 30 cm in columns. This frame had been previously simulated by [138] [139] by using concentrated and distributed plasticity models.

This test is simulated implementing the methodology of calculating the damage variables (section ‎3) in ABAQUS code [89]. Time integration follows an implicit formulation, and the global algorithm is generated by imposing energy balance. In this work, analyses are conducted for large displacements, although not for large strains. The maximum number of iterations is 10, and loading increment ranges between 0.0001 mm and 0.01 mm. Table 6 displays the selected values of the parameters. Time integration follows an implicit formulation, and the global algorithm is generated by imposing energy balance. In this work, analyses are conducted for large displacements, although not for large strains.

Tabla 6. Parameters for monotonic simulation of [Pires 1990] frame experiment
 fck (MPa) fcm (MPa) ftm (MPa) Gch (N/mm) GF (N/mm) b ac at leq (mm) bc bt 20 28 2.222 21.12 0.133 0.9 7.873 1 25 65.48 626.67

Figure 48 displays the finite element mesh; the right part depicts steel discretization with 2-node truss elements (T3D2) and the left part describes concrete discretization with 3D 8-node hexahedron solid elements (C3D8R). Embedded element technique [142] has been used to connect reinforcing bars with the surrounding concrete assuming perfect bond conditions. Isotropic hardening has been assigned to steel material for both longitudinal and transverse bars.

 a) Initial 2D elastic domain b) Uniaxial cyclic behavior Figure 48. Finite element discretization of [Pires 1990] RC frame experiment

Figure 49 presents the plots that constitute the major inputs of the proposed methodology, Figure 49.a and Figure 49.b display plots of compressive/tensile stress vs. crushing/cracking strain, respectively; Figure 49.c and Figure 49.d display plots of compressive/tensile damage variable. These plots correspond to the only mesh size in the model as presented in Table 6.

For comparison results are plotted in Figure 50 together with the experimental performance. Plots from Figure 50 show that the proposed methodology captures the initial stiffness, the inception of overall yielding, maximum strength capacity, ductility phase, sequence of damage progression and the final state.

 (a) Compressive stress vs. crushing strain (b) Tensile stress vs. cracking strain (c) Compressive damage variable vs. crushing strain (d) Tensile damage variable vs. cracking strain Figure 49. Major inputs of the damage variables methodology for [Pires 1990] RC frame experiment

 Figure 50. Experimental and simulated capacity curves for [Pires 1990] RC frame

Figure 51 shows the final damage state at the base of the right column, describing the final state and the evolution from the undamaged state. Figure 51.a and Figure 51.c represent the distribution of the compressive and tensile damage variables, respectively. Figure 51.b and Figure 51.d display, for selected finite elements, plots of such variables vs. beam displacement, respectively. Figure 51.a and Figure 51.c show that there is hinging, since both damage variables attain values close to 1. Figure 51.b and Figure 51.d show that cracking occurs for smaller displacement (approximately 20 mm) than crushing (more than 100 mm).

 (a) Distribution of dc at the final state (b) Evolution of dc (c) Distribution of dt at the final state (d) Evolution of dt Figure 51. Damage at the right column base for [Pires 1990] frame experiment

Figure 51 shows that the observed phenomena are adequately reproduced by the proposed methodology, since the obtained damage distributions fit the expected results. Results from Figure 51.b.d show the rapid damage evolution in terms of top floor displacement when concrete is subjected to mainly tensile stresses in comparison to compressed side, this also fits with expected behavior and the inputs shown in Figure 49. First bars yielding under tension is detect at the top-right edge of the beam, subsequently both column at the bases, ending by yielding at top-right, bottom left edges of beam and at column bases. Compressed bars at the bases of the columns reached as well the yield stress. Figure 52.a shows the distribution of plastic stain in the longitudinal bars at the base of the right column. Figure 52.b and Figure 52.c show respectively the tension and compression stress versus the lateral displacement at the beam. The stresses values in Figure 52.b and Figure 52.c are taken as the average stress value in the highly stresses elements in the three bars shown in Figure 52.a at each side.

 (a) Plastic strain (b) Tensile stress vs. displacement (c) Compressive stress vs. displacement Figure 52. Plastic strain in reinforcing bars are the right base column of [Pires 1990] experiment.

Final damaged state of the frame is demonstrated in Figure 53, compression, tension and scalar damage are presented respectively in Figure 53.a.b.c. Final stress state (KN/m2) in reinforcing bars is demonstrated in Figure 53.d

 (a) (b) (c) (d) Figure 53. Final damaged state of [Pires 1990] RC frame experiment

#### 6.2.2 Two floors RC frame [Vecchio and Emara 1990]

This experiment [143] is a quasi-static test consisting of pushing monotonically until failure a laboratory, single-span, two-story, planar RC frame. Figure 54.b and Figure 54.c display front and side views of the tested frame, respectively. Figure 54.a and Figure 54.d exhibit cross sections of columns and beams, respectively. In Figure 54 dimensions are in mm. Figure 54.b shows that both columns were loaded with constant forces, and that pushing consisted in imposing a displacement law to the top left joint. Noticeably, since there were no distributed forces acting on the beams, there was no cracking prior to the lateral pushing. The tested frame was widely instrumented, thus providing extensive information of the damage progression.

 Figure 54. Frame experiment [Vecchio and Emara 1990]

Concrete mechanical parameters were determined from standard cylinder tests. The characteristic value of compressive strength is 30 MPa and the average value of secant deformation modulus is 23.67 GPa. Poisson ratio and shear deformation modulus were estimated as 0.2 and 9.86 GPa, respectively. The steel parameters were obtained from coupon tests; the stress-strain plots were approximately trilinear: an initial linear elastic branch, a horizontal yielding plateau and a plastic hardening branch. Table 7 displays the most relevant figures regarding reinforcement steel. In Table 7, bars No. 20 and 10 correspond to longitudinal and transverse reinforcement, respectively (Figure 54.a and Figure 54.d). Db is the bar diameter, fy / fu are the yield point / ultimate stress, Es is the steel modulus of elasticity, Esh is the slope of the hardening branch, εsh is the strain that corresponds to the onset of hardening, and εu is the ultimate strain.

Table 7. Reinforcement steel parameters for the [Vecchio and Emara 1990] frame experiment
 Bar No. Db (mm) fy (MPa) fu (MPa) Es (GPa) Esh (MPa) εsh εu 20 19.5 418 596 192.5 3100 0.0095 0.0669 10 11.3 454 640 200 3100 0.0095 0.0695

Figure 55 summarizes the results of the experiment. Figure 55.a displays plots of pushing force vs. top level displacement. The first observed damage was flexural cracking at the end sections of the first level beam, for force 52.5 kN and displacement 2 mm; this instant can be considered as the overall yield point, since, after that point, the capacity curve became nonlinear. Under force 145 kN and displacement 9.3 mm, flexural cracks were perceived at the columns bases, and shear cracks were simultaneously detected at the first story beam. When force reached 264 kN and displacement 26.4 mm, first steel yielding was perceived in the bottom tensioned longitudinal reinforcement of the left end section of the first story beam; for force 287 kN and displacement 31.6 mm, the top tensioned longitudinal reinforcement of the right end section of the first story beam also yielded. For force 323 kN and displacement 52.5 mm, the longitudinal reinforcements of the column bases yielded as well, and hinges at the ends of the first story beam failed; this failure involved yielding of longitudinal reinforcement, and crushing of compressed concrete. Then, for force 329 kN and displacement 74.7 mm, similar failure was apparent at the column bases. Almost simultaneously, same failure affected at ends of second story beam. Afterwards, lateral stiffness was almost non-existent; therefore, collapse mechanism consisted in formation of six hinges. The experiment was terminated, for pushing force 332 kN and lateral displacement 150 mm, due to stroke limitations of the actuator, see Figure 55.a. Figure 55.b and Figure 55.c display images of the damaged bottom section of the right column, and the top left connection, respectively.

 (a) Force-displacement plot (b) Right column base (c) Top left connection Figure 55. Results of [Vecchio and Emara 1990] frame experiment [143]

The main objective of this experiment was to investigate the influence of shear-related effects in the overall structural behavior; the results showed that approximately 20% of the nonlinear lateral displacement was due to shear effects. Noticeably, at failure, around 12% of the overturning moment was due to P-Δ effects. Supplementary information regarding this experiment is available in [144].

This test had been previously simulated by some of the authors of this work [145] [146] [147], and by other researchers [148]. [145], used a viscous damage model, which was implemented in a fiber model with Timoshenko frame elements. [146] [147], implemented a damage plasticity formulation in a finite element model with planar 2D elements. [148], employed commercial software packages and an ad-hoc fiber model developed at University of Toronto.

Analogously to the simulated experiment in the previous subsection ‎6.2.1, this test is simulated implementing the damage evolution algorithm in ABAQUS code [89]. The maximum number of iterations is 10, and loading increment ranges between 10-11 mm and 0.01 mm. Table 8 displays the selected values of parameters. Figure 56 displays the finite element mesh; the right side depicts steel discretization with 2-node truss elements (T3D2) and the left side describes concrete discretization with 3D 8-node hexahedron solid elements (C3D8R).Perfect bond is also assumed between concrete and reinforcing bars.

 Figure 56. Finite element discretization of [Vecchio and Emara 1990] RC frame experiment
Table 8. Parameters for simulation [Vecchio and Emara 1990] frame experiment
 fck (MPa) fcm (MPa) ftm (MPa) Gch (N/mm) GF (N/mm) b ac at leq (mm) bc bt 30 38 2.912 23.93 0.1405 0.75 7.873 1 50 156.83 1554.54

Figure 57.a displays plots of ratio ${\textstyle {\epsilon }_{c}^{pl}\,/\,{\epsilon }_{c}^{ch}}$ vs. ${\textstyle {\,\epsilon }_{c}^{ch}}$ (refer to section ‎3); Figure 57.b displays analogous plots for the tensile behavior. Figure 57 highlights that, in the proposed methodology, damage and plastic energy absorptions are not related, since ratios ${\textstyle {\epsilon }_{c}^{pl}\,/{\,\epsilon }_{c}^{ch}}$ and ${\textstyle {\epsilon }_{t}^{pl}\,/{\,\epsilon }_{t}^{ck}}$ do not approach zero when damage variables dc and dt are close to 1. Figure 57.a shows that the average value of ratio b (equation (56)) is approximately 0.75 (Table 8).

 (a) Compression (b) Tension Figure 57. Variation of ratio between plastic and crushing / cracking strains for [Vecchio and Emara 1990] experiment

Figure 58 displays experimental results (Figure 55.a) plotted together with numerical results obtained with the proposed methodology. Descriptions of observed damage states are also displayed.

 Figure 58. Experimental and simulated capacity curves for [Vecchio and Emara 1990] frame experiment

Plots from Figure 58 show the superior ability of the proposed methodology to reproduce the experimental results along the whole displacement range. It captures the initial stiffness, the onset of overall yielding, ductility phase and the sequence of damage progression and the final state. Noticeably, “overall yielding” does not refer to steel yielding but to inception of overall nonlinear behavior, due to cracking of tensioned concrete.

 (a) Distribution of dt at the first story right beam-column connection (b) Distribution of dt at the frame Figure 59. Tensile damage variable for force 52.5 kN and displacement 2 mm for [Vecchio and Emara 1990] frame experiment

To further highlight the capacity of the proposed methodology to capture damage progression, Figure 59 through Figure 61 display the damage predicted for some of the previously described stages. Figure 59 displays the distribution of the tensile damage variable for force 52.5 kN and top level displacement 2 mm; it corresponds to the first cracking at end sections of 1st story beam. Figure 59.a and Figure 59.b refer to the first story right beam-column connection and to the overall frame, respectively. Figure 59.a shows that cracking (indicated with lighter gray) actually occurred in the top part of beam, since the tensile damage variable reaches values close to one. Figure 59.b shows that the overall distribution of cracking fits the expected pattern according to structural analysis principles, with onset of cracking in the bottom part of the first story right beam-column connection. Figure 60 and Figure 61 refer to the final damaged state, distributions of the scalar damage variable (d) for the right column base and the top left beam-column connection are shown respectively in Figure 60.a.b.

 (a) Right column base (b) Top left connection Figure 60. Distribution of d at the final state of [Vecchio and Emara 1990] frame experiment

 (a) Distribution of dc at the final state (b) Evolution of dc (c) Distribution of dt at the final state (d) Evolution of dt Figure 61. Damage at the right column base [Vecchio and Emara 1990] frame experiment

Figure 61 refers to the right column base, describing both the final state and the evolution from the undamaged state. Figure 61.a and Figure 61.c represent the distribution of the compressive and tensile damage variables, respectively; Figure 61.b and Figure 61.d display, for selected finite elements, plots of such variables vs. top level displacement, respectively. In all the images, higher damage corresponds to lighter grey. Figure 60 and Figure 61 show that the observed phenomena are adequately reproduced by the proposed methodology, since the obtained damage distributions fit the expected results. Comparison between Figure 60.a and the observed damage in Figure 55.b shows a satisfactory fit, since cracking and crushing are detected by the obtained higher values of d. Comparison between Figure 60.b and the observed damage in Figure 55.c shows also a satisfactory match. Figure 61.a and Figure 61.c show that there is hinging, since both damage variables attain values close to 1; this circumstance is observed in Figure 55.b. Figure 61.b and Figure 61.d show that cracking occurs for smaller displacement (approximately 10 mm) than crushing (more than 150 mm).

Final damaged state of the frame is demonstrated in Figure 62, compression, tension and scalar damage are presented respectively in Figure 62.a.b.c. Final stress state in reinforcing (KN/m2) bars is demonstrated in Figure 62.d

 (a) (b) (c) (d) Figure 62. Final damaged state of [Vecchio and Emara 1990] frame experiment

#### 6.3.1 Bridge column [Tanaka 1990]

Several RC columns have been tested by [Tanaka 1990] [149] to study the effect of lateral confinement on their ductility. Main differences are the type of transverse reinforcement, its anchorage detailing, the axial force, and the mechanical and geometric parameters. Column unit 6 [149] has been analyzed with the proposed model. Figure 63 describes the tested column and the experimental mockup.

 (b) Vertical section (a) Loading arrangements (c) Horizontal section Figure 63. Experiment on [Tanaka 1990] bridge column [98]

Figure 63.b and Figure 63.c display the geometry of the column, dimensions are in mm. The column was loaded with constant force 0.1 f’c/Ag and a cyclic lateral displacement protocol was imposed by hydraulic jacks as shown in Figure 63.a. Concrete mechanical parameters were determined from standard cylinder tests. The characteristic value of compressive strength is 32 MPa and the average value of secant deformation modulus is 27.65 GPa; Poisson ratio and shear deformation modulus were estimated as 0.2 and 11.52 GPa, respectively. Steel parameters were obtained from coupon tests. Table 9 displays the most relevant information regarding reinforcement steel. In Table 9, bars HD20 and D12 correspond to longitudinal and transverse reinforcement, respectively. Db is bar diameter, fy / fu are yield point / ultimate stress, Es is modulus of elasticity, εsh is the strain corresponding to onset of hardening, and εu is ultimate strain.

Tabla 9. Reinforcement steel parameters for [Tanaka 1990] bridge column
 Bar Id. Db (mm) fy (MPa) fu (MPa) Es (GPa) εsh εu HD20 20 511 675 200 0.0165 0.14 D12 12 325 429 200 0.015 0.14

Figure 64 displays the imposed displacement law.

 Figure 64. Imposed displacement law for [Tanaka 1990] bridge column experiment

The results of the experiment showed stable hysteretic behavior, good energy dissipation and limited strength reduction up to the final stage. Substantial crushing of compressed cover concrete (spalling) was first observed near peak load during the first cycle. During the final stage, the visible damage was concrete crushing and slight buckling of longitudinal reinforcement bars.

Figure 65 depicts the element discretization. Figure 65.a describes the concrete discretization with 3-D 8-node hexahedron solid elements (C3-D8R). Figure 65.b and Figure 65.c represent the steel and interface elements; steel is discretized with 2-node truss elements (T3-D2). Bond-slip effect is considered for the column longitudinal bars at the base of the column and development region in the footing, as shown in Figure 65.b; full bond conditions are assumed for the other segments of longitudinal bars as well as for transverse reinforcement because minimal bar slip is expected for these bars.

 (a) Column discretization (b) Steel and interface elements (c) Plan view of steel and interface elements Figure 65. Finite element discretization of [Tanaka 1990] bridge column

Steel behavior is described with a classical plastic model; isotropic hardening is used for stirrups, and nonlinear isotropic/kinematic hardening for longitudinal bars. Concrete behavior is simulated by implementing the methodology described section ‎3 in ABAQUS code [89].

Figure 66.a displays plots of ratio ${\textstyle {\epsilon }_{c}^{pl}\,/\,{\epsilon }_{c}^{ch}}$ vs. ${\textstyle {\,\epsilon }_{c}^{ch}}$, where ${\textstyle {\epsilon }_{c}^{pl}}$ refers to plastic compressive strain. Figure 66.b displays analogous plots for the tensile behavior ( ${\textstyle {\epsilon }_{t}^{pl}\,/\,{\epsilon }_{t}^{ck}}$ vs. ${\textstyle {\,\epsilon }_{t}^{ck}}$).

 (a) Compression (b) Tension Figure 66. Variation of ratios plastic and crushing / cracking strains for [Tanaka 1990] bridge column experiment

Figure 66.a points out that the average of ${\textstyle {\epsilon }_{c}^{pl}\,/\,{\epsilon }_{c}^{ch}}$ is approximately equal to 0.88; this magnitude is termed b and is a required parameter in the implementation of the algorithm [150]. Table 10 displays the considered values of the parameters.

Tabla 10. Límites mínimos y máximos de riesgos de las funciones de trasformación utilizadas
 fck (MPa) fcm (MPa) ftm (MPa) Gch (N/mm) GF (N/mm) b ac at leq (mm) bc bt 24 32 2.51 22.14 0.136 0.889 7.873 1 50 142.65 1381.74

Figure 67 displays comparisons between experimental and numerical hysteresis loops calculated with the proposed model. In Figure 67.a, the numerical results are obtained by using hc = 0.9 in all the body, and generating a discrete crack at the footing-column interface; normal friction contact condition is assumed at that interface. In Figure 67.b, the approach to represent concrete discontinuity described in subsection ‎6.1 is applied by assigning ℎc = 0.01 to the first row of elements above the footing, ℎc = 0.1 to the elements where bond-slip is considered (Figure 65.b), and ℎc = 0.9 for the rest of the model (Figure 46). These values have been selected to provide the best fit with the experiment.

 (a) Experiment and proposed formulation with a discrete crack (b) Experiment and proposed formulation with reduced ℎc (subsection 6.1) Figure 67. Experimental and numerical force-displacement response of [Tanaka 1990] bridge column experiment

Figure 67 shows that the proposed model provides a satisfactory agreement with test results, capturing initial stiffness, strength and stiffness degradation, and pinching. Plots from Figure 67.a show that the model with discrete crack is able to reproduce the main aspects of cyclic behavior; the analysis ended prematurely due to the large nonlinearity and the important separation between both blocks. The model with reduced ℎc (Figure 67.b) captured suitably the dissipated energy from the experiment with difference 1.66 %; conversely, the model with discrete crack (Figure 67.a) overestimated that energy by 22.88 %.

Figure 68 displays comparisons between pairs of plots obtained by numerical simulation with the proposed model and with a similar model assuming perfect bond between the reinforcement and the surrounding concrete. In Figure 68.a a discrete crack is introduced in the proposed model (as in Figure 67.a), and in Figure 68.b reduced values of ℎc are considered (as in Figure 67.b).

 (a) Proposed formulation with a discrete crack vs. perfect bond (b) Proposed formulation with reduced ℎc (subsection 6.1) Figure 68. Numerical force-displacement response of [Tanaka 1990] bridge column experiment

Figure 68 highlights the relevance of bond-slip. Comparison between Figure 68.a and Figure 68.b shows that differences are more significant in the smeared model. The discrete crack model response in Figure 68.a exhibits some pinching due to the added flexibility at the column bottom; conversely, the smeared model response in Figure 68.b shows higher energy dissipation.

Figure 69 shows the experimental results plotted together with the ones obtained from the models assuming perfect bond, plots show that models with perfect bond reproduce satisfactorily the envelope behavior but lack to capture the unloading and reloading branches.

 (a) Experiment vs. Perfect bond proposed formulation with a discrete crack (b) Experiment vs. Perfect bond proposed formulation model with reduced ℎc (subsection 6.1) Figure 69. Numerical force-displacement response of [Tanaka 1990] bridge column

Experimental observations revealed buckling in the longitudinal reinforcement bars at column bottom [149]. Figure 70 presents numerical (Figure 70.a and Figure 70.b) and experimental (Figure 70.c) representations of this failure. Figure 70.a displays the buckling that is detected at the last cycle, Figure 70.b shows the final damage state and Figure 70.c depicts the observed damage. Comparison between numerical and experimental results highlights the accuracy of the proposed formulation.

 (a) Buckling in the longitudinal bars (b) Final damage state (c) Final state [149] Figure 70. Final state of the bridge column

The average stress-vertical displacement behavior of the bottom-left corner longitudinal bars presented in Figure 65.c is shown in Figure 71. The average stress-displacement along embedded part of the bar inside the footing is shown in Figure 71.a, while the average stress- displacement along the region where interface elements are assigning (Figure 65) is shown in Figure 71.b. Plots from Figure 71 shown the significant influnce of introducing the bond-slip elemnts in comparision with perfect bond model on the hesteritic behavior of stress-displacment in longitudinal bars. The higher negative displacement in the hysteretic behavior of the model with bond-element shown in Figure 65, can be understood as a result of penetration due to bond degradation.

 (a) Proposed formulation with reduced ℎc (subsection 6.1) vs. perfect bond along 55 cm (b) Proposed formulation with reduced ℎc (subsection 6.1) vs. perfect bond along 130 cm Figure 71. Numerical results of stress-slip relation in the reinforcing bars for [Tanaka 1990] bridge column experiment

The average cyclic damage evolution at the bottom of the column with depth 50 cm (See Figure 73) is presented in Figure 72. Plots in Figure 72 are obtained from the proposed model with reduced hc (subsection ‎6.1). Plots from Figure 72 show the evolution of the compressive, tensile and scalar damage along the entire range of loading. It can be seen that the scalar damage follows mainly the tensile damage at a certain value, subsequent, the value of the scalar damage exceeds the tensile one and rises up as long as the compressive damage is increasing. This results holds true for a low value of hc in the considered zone, as well low value of r* as it’s explained in subsection ‎2.4, equations (24) (25) (26).

 Figure 72. Cyclic damage evolution of [Tanaka 1990] bridge column experiment, proposed formulation with reduced ℎc (subsection 6.1)

Same plots from Figure 72 are plotted versus top lateral displacement in Figure 73.b.c.d for compressive, tensile and scalar damage respectively. Plots from Figure 73 confirm the expected results, in which the damage variables, as well the scalar damage increase along increasing the displacement, as well as, when at the last two cycles for constant value of displacement (see Figure 64)

 (a) Cinsidered zone for damage evolution (b) Compressive damage evolution (c) Tensile damage evolution (d) Scalar damage evolution Figure 73. Cyclic damage evolution vs. top displacement of [Tanaka 1990] bridge column experiment, proposed formulation with reduced ℎc (subsection 6.1)

#### 6.3.2 Building frame [Pires 1990]

The experiment presented subsection ‎6.2.1 and initially was simulated for monotonic loading, is further studied and simulated again under cycle loading. The RC frame is subjected to cyclic displacement law in Figure 74, this displacement is imposed at the beam side. Mechanical parameters of materials as well as the geometry of the frame are presented in subsection ‎6.2.1.

 Figure 74. Imposed cyclic displacement law for [Pires 1990] frame experiment

Analogously to bridge column experiment (subsection ‎6.2.1), this test is simulated by implementing the damage variables algorithm [129] in ABAQUS code and the interface bond-slip model described in ‎5.3, modeling schemes described in subsections ‎5.4 and ‎6.1 are also incorporated. The model has been refined to implement the bond-slip elements in the foundation. Table 6 displays the calculated values for the parameters of the damage variables methodology.

Table 11. Parameters for monotonic simulation of [Pires 1990] frame experiment
 fck (MPa) fcm (MPa) ftm (MPa) Gch (N/mm) GF (N/mm) b ac at leq (mm) bc bt 20 28 2.222 21.12 0.133 0.97 7.873 1 21.6 56.57 541.45 20 28 2.222 21.12 0.133 0.94 7.873 1 37 96.91 927.48

This experiment is simulated using bond-slip elements where bars are expected to slip. Figure 75 displays the frame discretization; Figure 75.a shows the overall discretization, Figure 75.b presents a general view of steel elements, Figure 76.a shows the discretization of interface element "bond-slip element" attached to longitudinal bars, Figure 76.b depicts a more detailed view of the right beam-column joint steel bars, and Figure 76.c shows the same of Figure 76.b with bond-slip element, Figure 76.d depicts a closer view of the right column base.

 (a) Frame discretization (b) Steel elements discretization Figure 75. Finite element discretization of [Pires 1990] frame experiment

As shown in Figure 76.a, bond-slip elements are assigned to the longitudinal bars along the column foundation, inside beam-column joints, at the column top and bottom ends, and at beam-ends. Noticeably, the proposed modeling scheme (subsection ‎5.4) allows that bar elements intersect in the same solid element without need to generate perpendicular holes nor to use double nodes. Regarding the approach described in subsection ‎6.1, similarly to Figure 46, ℎc = 0.01 for the bottom and top row of elements inside the columns; then, ℎc = 0.1 for the elements where bond-slip is considered (Figure 76), and ℎc = 0.9 for the rest of the model. The results from this approach are compared with those after introducing discrete cracks right above the footings.

 (a) Steel and interface elements (b) Steel element at right beam-column joint (c) Interface elements at the right beam-column joint (d) Interfce elements at the right column base Figure 76. Bond-slip elements discretization of [Pires 1990] frame experiment

Analogously to Figure 67, Figure 77 displays comparisons between experimental and numerical hysteresis loops obtained with the proposed model. Plots from Figure 77 provide similar conclusions than Figure 67.

 (a) Experiment and proposed formulation with a discrete crack (b) Experiment and proposed formulation with reduced ℎc (subsection 6.1) Figure 77. Experimental and numerical force-displacement response of [Pires 1990] frame experiment

Similarly to Figure 68, Figure 78 displays comparisons between numerical hysteresis loops obtained with the proposed model and assuming perfect bond. Figure 78 allows deriving parallel observations than Figure 68.

 (a) Proposed formulation with a discrete crack vs. perfect bond (b) Proposed formulation with reduced ℎc (subsection 6.1) vs. perfect bond Figure 78. Numerical force-displacement response of the frame experiment

 (a) Experiment vs. Perfect bond proposed formulation with a discrete crack (b) Experiment vs. Perfect bond proposed formulation model with reduced ℎc (subsection 6.1) Figure 79. Numerical force-displacement response of [Pires 1990] frame experiment

Figure 99 display similar comparison as shown in Figure 69, derived observation are also similar. Final state in terms of tensile, compressive, and total damage is represented in Figure 80, results from Figure 80 show that the damage is mainly localized at the top and bottom of the column, as well in the joints. This can be explained by the small sections of column in comparison with the beam section, also since the beam is not loaded by any vertical load before applying the lateral displacement conversely to the columns.

 (a) Compressive damage (b) Tensile damage (c) Scalar damage Figure 80. Final damage state, proposed formulation with reduced ℎc (subsection 6.1) of [Pires 1990] frame experiment

Similarly to Figure 72, the average cyclic damage evolution at the bottom of the column with depth 30 cm is presented in Figure 81. Plots in Figure 72 are obtained from the proposed model with reduced hc (subsection ‎6.1). Plots from Figure 81 show the evolution of the compressive, tensile and scalar damage along the entire range of loading. Similar results to Figure 72 can be derived from plots in Figure 81

 Figure 81. Cyclic damage evolution of [Pires 1990] frame experiment, proposed formulation with reduced ℎc (subsection 6.1), bottom of column

Same plots from Figure 81 are plotted versus top lateral displacement in Figure 82.b.c.d for compressive, tensile and scalar damage respectively. Plots from Figure 82confirm again the expected results, in which the damage increases along the last five cycles at a constant displacement value (see Figure 74).

 (a) Cinsidered zone for damage evolution (b) Compressive damage evolution (c) Tensile damage evolution (d) Scalar damage evolution Figure 82. Cyclic damage evolution vs. top displacement of [Pires 1990] frame experiment, proposed formulation with reduced ℎc (subsection 6.1), bottom of column

Similar plots to the ones in Figure 81 and Figure 82 are presented in the next figures for the average damage evolution at top part of the column, and at the beam column joint. Same results can be derived.

 Figure 83. Cyclic damage evolution of [Pires 1990] frame experiment, proposed formulation with reduced ℎc (subsection 6.1), top of column

 (a) Cinsidered zone for damage evolution (b) Compressive damage evolution (c) Tensile damage evolution (d) Scalar damage evolution Figure 84. Cyclic damage evolution vs. top displacement of [Pires 1990] frame experiment, proposed formulation with reduced ℎc (subsection 6.1), top of colum

 Figure 85. Cyclic damage evolution of [Pires 1990] frame experiment, proposed formulation with reduced ℎc (subsection 6.1), beam column joint

 (a) Cinsidered zone for damage evolution (b) Compressive damage evolution (c) Tensile damage evolution (d) Scalar damage evolution Figure 86. Cyclic damage evolution vs. top displacement of [Pires 1990] frame experiment, proposed formulation with reduced ℎc (subsection 6.1), beam column joint

## 7 Continuum mechanics-based approaches concluding remarks

A new methodology to calculate damage variables "damage evolution" for Concrete Plastic Damage Models has been proposed, the methodology is suitable to simulate the monotonic and cyclic behavior of reinforced concrete structures.

The proposed methodology is mesh-insensitive, is based on a sound continuum mechanics formulation, and does not require calibration with experimental results. Implementation is straightforward. A particular algorithm is presented to implement the methodology in the CPDM of FEM software ABAQUS. Mesh-insensitivity is validated in a simple tension example. Accuracy and reliability of the proposed methodology are verified by simulating a cyclic loading experiment on plain concrete specimens, the methodology has been verified as well by simulating large scale tests consist in pushing monotonically and cyclically until failure different RC structures. Following the main verification conclusions:

Tension example. This example confirms the mesh-insensitiveness of the proposed methodology. As well, numerical results match analytical ones, thus ratifying their accuracy.

Cyclic experiment on plain concrete. The simulation of the cyclic experiment on plain concrete corroborates the capacity of the proposed approach to reproduce damage evolution and stress-strain behavior, both in terms of envelope curve and stiffness degradation at each unloading / reloading cycle. Agreement between experimental and numerical damage at measured strain levels is satisfactory.

Validation with Experiments-Monotonic loading: Two RC frames have been simulated, results have proved the ability of the proposed methodology to reproduce accurately the experimental force-displacement results along the whole displacement range. Captures initial stiffness, onset of the overall yielding, sequence of the damage progression, collapse mechanism, and final state. All failure and degradation modes are adequately simulated: tensioned concrete cracking, compressed concrete crushing, and reinforcement steel yielding.

Aforementioned work has been extended to study the cyclic behavior of RC structures. In this context, an integrated continuum FEM model is proposed for monotonic and cyclic simulation of RC frames. The model combines the methodology proposed for calculating the damage variables of CPDM and the 3-D interface bond-slip model developed by [20]. A new technique to integrate the interface model in a continuum FEM model of regions with crossing bars is presented. The proposed model is used to simulate two experiments consisting in imposing cyclic displacement laws to an RC column and a frame. These experiments are also simulated with the same model although assuming perfect bond conditions.

Validation with Experiments-Cyclic loading: Obtained results of simulating the cyclic behavior of one RC bridge column and one RC building frame, show that the proposed model is able to predict the actual behavior of highly damaged RC elements and frames, capturing strength reduction, stiffness degradation and pinching. Noticeably, this model uses only a reduced number of parameters and does not require any calibration with experimental results. Comparison with the results obtained assuming perfect bond points out the high relevance of bond-slip in the hysteretic behavior and the energy dissipation capacity.

The proposed integrated model can constitute a practical tool for accurate simulation of highly damaged medium-size RC structures, both in research and in conventional analyses. It can be also used to calibrate more simplified models.

# PART II-STRUCTURAL COMPONENT-BASED APPROACHES

## 1 Introduction and motivation

The purpose of developing Part II is; First, investigate the capacity of the simplified numerical models and commonly used in earthquake engineering. Second, although the advanced model developed in Part I showed satisfactory performance, its high computational cost raises the need to develop a model that can capture the overall response with lower cost.

As discussed in the general introduction of this dissertation, numerical models can be divided into two main groups according to their formulation to describe the nonlinear behavior of structures. Part I dealt with complex models that are based on continuum mechanics-based approaches, this implied that structural components (e.g. columns, beams, walls) are discretized to many micro finite elements. Reinforced concrete elements are simulated by combining two materials, concrete simulated by solid elements while steel bars are described by truss bars, each is associated to a constitutive material model, relation between concrete and steel elements is described as well by a continuum model. Part II will deal with more simple models that are based on structural approaches. In other words, the models discussed in this part are derived from simplified theories of mechanics of material, structural components can be discretized with simple finite elements models (e.g. flexural bar element for column and beams). The composite material behavior is not described explicitly, rather the used element in these models aims to describe the response of the complete structural component. Different sources of degradation are also described, either by the element directly or by attaching specific element in the locations where the damage is expected to occur. These models are referred in this work as structural component-based models.

It can be distinguished between two type of structural component-based models according to the distribution of plasticity. These models are named as lumped and distributed models. Lumped plasticity models assume that the plastic deformation will be concentrated at one point along the length of the element, while distributed plasticity models provide more accuracy by capturing the plastic deformations along a specific plastic hinge length of at different locations along the element length and at different points at the cross section of the element. The lumped plasticity approach is advantageous due to its computational efficiency and for its analogy to simple elasto-plastic analysis procedures used in design practice. The distributed plasticity approach, however, is convenient for performance based design since it has the ability to capture local behaviors at intermediate element lengths and takes the spread of plasticity along an element length into account. Both approaches are discussed in the following sections, their capacity is verified by simulating the experiments described in the previous part of this dissertation.

After verifying the commonly used structural component-based models in earthquake engineering, a numerical study is conducted on the relation among the non-simulated deterioration modes of the elements of non-ductile RC frames and their final capacity.

The main contribution of this part is developing an advanced structural component-based model for simulating the nonlinear dynamic behavior of old reinforced structures, accounting for flexure, shear and axial deterioration modes. The capacity of the developed model is verified by simulating the nonlinear dynamic behavior of an existing non-ductile building and the prototype building.

## 2 Models for nonlinear analysis

### 2.1 Lumped plasticity models

The lumped plasticity approach utilizes the simplicity of the plastic hinge concept by separating a line element into elastic and inelastic components. A fixed point "zero length element" or predetermined length (typically the plastic hinge length), is determined to be the region in which all inelastic action is concentrated, while elastic properties are assigned to the remainder of the element. The inelastic behavior in the plastic hinge is determined from a sectional analysis of the critical section, which has detailed description of the geometry and material properties (uniaxial constitutive models) assigned to sections within this region. Uniaxial constitutive material models for both steel and concrete is defined and assigned to the corresponding components of a discretized section. This section is then assigned to the plastic hinge(s) of the element. Determining the section’s behavior usually requires a sectional analysis to be completed first, and then the moment-curvature relationship is assigned to the plastic hinge(s) to determine the element response under external loads. Other methods define a backbone curve (capacity curve), such methods are acceptable when the element behavior is previously well known and does not require the initial sectional analysis.

The first inelastic concentrated plasticity model was proposed by Clought et al in 1965. [151] shown in Figure 87.a. In this model, known as the two-component model, a bilinear elastic-strain hardening moment-curvature relationship is assumed along the length of the element. The beam model consists of two components acting in parallel: one which is linear elastic and one which is elastic-perfectly plastic with the plastic deformations concentrated in plastic hinges at the ends of the element. The elastic modulus of the linear component is equal to the strain hardening modulus p ⋅ EI of the moment-curvature relation, where EI is the pre-yield section stiffness. The elastic modulus of the elasto-plastic component is equal to q ⋅ EI where q = 1− p. One of the shortcomings of this model is the difficulty of accounting for the stiffness deterioration of RC elements during cyclic load reversals.

To overcome the problem of stiffness deterioration, Giberson proposed another model in 1967 [152] . This model is known as the one-component model shown in Figure 87.b . It consists of two nonlinear rotational springs which are attached at the ends of a perfectly elastic element representing the girder. All nonlinear deformations of the girder element are lumped in the two rotational springs. This is a simplification of experimental evidence which shows that inelastic deformations spread over a finite region at the ends of the girder. Giberson's model has the advantage that any kind of hysteretic law can be assigned to the nonlinear springs. This fact along with the simplicity of the model accounts for its wide use in analytical studies to date.

To describe the cyclic behavior of nonlinear spring, a hysteretic law is needed. The first model was proposed by [151] and a more refined model was Takeda et.al 1970 [153]. After this, Takeda's model has been widely used in the description of the hysteretic moment curvature or moment-rotation relation of RC members.

The concentrated plasticity concept was generalized to incorporate axial-moment interaction at the element ends, the first model with interaction between normal loads and bending moment was proposed by Lai et.al in 1985 [154]. Much later, these concepts were used to develop a simplified macro-element capable of accounting for inelastic lateral-torsional coupling in asymmetric buildings. This model was developed considering a perfect elasto-plastic interaction between the story-shears and torque [155]. More recently, other researchers have included effects such as hardening into bi-dimensional hinges by introducing a hardening rule to the yield surface [156].

 Two-component model [151] (a) One-component model [152] Figure 87. Nonlinear two and one component models for a nonlinear beam element

A complete model for the analysis of seismic response of RC structures was proposed by Banon et al.1981 [157]. The one-component model in its original form describes the nonlinear behavior of the girder. The hysteretic moment-rotation relation is based on a modified Takeda model. In order to reproduce the "pinching" effect due to shear and bond deterioration, a nonlinear rotational spring is inserted at each member end. The hysteretic model of the nonlinear springs is based on a bilinear skeleton curve with strength decay under large deformations and includes the effect of "pinching" during reloading.

The method of defining backbones curves is well approved as the most practical way in performance based design, such backbones or capacity curves can be represented by moment-rotation or force-deformation relation for ductile and brittle behavior respectively [158] [16].

Figure 88 shows a general representation of lumped plasticity model with assigned capacity curves for plastic hinges. Backbone curve can be obtained from experimental results (Figure 88.b) as the envelope of the hysteretic behavior, they can be classified according to the level of the ductility into three levels [16] as shown in Figure 88.c. For cyclic loading a hysteretic response must be assigned to the plastic hinge, the hysteretic response must reflect the main characteristics of the cyclic behavior (e.g. strength reduction, stiffness degradation, buckling in compressed reinforcing bars, cracks opening and closing and bond-slip degradation).

 (a)) Lumped plasticity model (b) Backbone representation of hysteretic behavior [16] (c) Backbone curve types [16] Figure 88. Nonlinear behavior of lumped plasticity model

As mentioned before, Takeda model has been widely used for reinforced concrete simulation, however, many calibrated models have been developed recently. The most recent model is Ibarra et al.2005 [2]that is based on modifying three of the basic hysteretic models used in seismic demand evaluation to include deterioration properties: bilinear, peak-oriented, and pinching. The modified models include most of the sources of deterioration: i.e. various modes of cyclic deterioration and softening of the post-yielding stiffness, and also account for a residual strength after deterioration. The models incorporate an energy-based deterioration parameter that controls four cyclic deterioration modes: basic strength, post-capping strength, unloading stiffness, and accelerated reloading stiffness deterioration. Figure 89 shows the main hysteretic features of Ibarra model

Some types of concentrated hinge models employ axial load-moment (P-M) yield surfaces. Whereas these models generally do a good job at tracking the initiation of yielding under axial load and bending, they may not capture accurately the post-yield and degrading response. On the other hand, some hinge elements with detailed moment-rotation hysteresis models may not capture P-M interaction, except to the extent that the moment-rotation response is defined based on average values of axial load and shear that are assumed to be present in the hinge.

 (a) Basic strength deterioration (b) Post-capping strength deterioration (c) Unloading stiffness deterioration (d) Accelerated reloading stiffness deterioration (e) Basic pinching model rules (f) Modification in pinching according to reloading Figure 89. Hysteretic bases of Ibarra et al. 2005 model [2]

### 2.2 Distributed plasticity models

From the physical behavior point of view, the damage in reinforced concrete structures is not lumped in a specific section, hence, spreading the plasticity along a specific length gives more accuracy to represent the nonlinear behavior of RC structure.

The first model which accounts for the spread of inelastic deformations into the member was introduced by Soleimani et al.1979 [159]. In this model a zone of inelastic deformations gradually spreads from the beam-column interface into the member as a function of loading history. The rest of the beam remains elastic. The fixed-end rotations at the beam-column interface are modeled through point hinges which are inserted at the ends of the member. These are related to the curvature at the corresponding end section through an "effective length" factor which remains constant during the entire response history. The model of Soleimani was further improved by Meyer et al. 1983 [160], by proposing different way to calculate the stiffness of the plastic zone during reloading and incorporating Takeda's law to describe the hysteretic behavior. This model was used to simulate the nonlinear behavior of beam and column with counting for the axial load. Further development of Meyer's model to incorporate the effect of shear and axial load has been done by Roufaiel and Meyer 1987 [161]. Another approaches were proposed by Filippou and Issa 1988 [162], in this model the element is divided into a finite number of short subelements. Each subelement describes a single effect, such as inelastic behavior due to bending, shear behavior at the interface or bond-slip behavior at the beam column joint.

Distributed plasticity models were stated at the beginning into displacement method (stiffness method). In this approach the unknowns are the displacement and cubic interpolation functions are used to approximate the deformation along the beam-column element. The first displacement-based model used the cubic Hermitian interpolation functions to describe the deformation of the element along its length. Relation between the element deformation and nodal displacement of the beam-column element represented in Figure 90 is writing as follow:

 ${\displaystyle \left\{\left.d(x)\right\}\right.=\left[N\right]\left\{\left.\delta \right\}\,\right.}$
(76)

Where ${\textstyle d(x)}$ is the transverse and longitudinal displacement along the longitudinal axes (x), ${\textstyle \left[N\right]}$ is the matrix of the interpolation functions and ${\textstyle \left\{\left.\delta \right\}\,\right.}$ is the vector of the nodal displacements presented in Figure 90. Equation (76) can be written as :

 ${\displaystyle \left\{\left.{\begin{matrix}u(x)\\v(x)\end{matrix}}\right\}\right.=\left[{\begin{matrix}{N}_{1}&0&0\\0&{N}_{3}&{N}_{4}\end{matrix}}\,\,{\begin{matrix}{N}_{2}&0&0\\0&{N}_{5}&{N}_{6}\end{matrix}}\right]\,\left\{\left.{\begin{matrix}{\begin{matrix}{u}_{i}\\{\nu }_{i}\\{\theta }_{i}\end{matrix}}\\{u}_{j}\\{\begin{matrix}{\nu }_{j}\\{\theta }_{j}\end{matrix}}\end{matrix}}\right\}\right.}$
(77)

where

 ${\displaystyle {\begin{array}{ll}\displaystyle {N}_{1}=1-{\frac {x}{L}}&\displaystyle {N}_{4}={\frac {{x}^{3}}{{L}^{2}}}-2{\frac {{x}^{2}}{L}}+x\\\displaystyle {N}_{2}={\frac {x}{L}}&\displaystyle {N}_{5}=-2{\frac {{x}^{3}}{{L}^{3}}}+3{\frac {{x}^{2}}{{L}^{2}}}\\\displaystyle {N}_{3}=2{\frac {{x}^{3}}{{L}^{3}}}-3{\frac {{x}^{2}}{{L}^{2}}}+1&\displaystyle {N}_{6}={\frac {{x}^{3}}{{L}^{2}}}-{\frac {{x}^{2}}{L}}\end{array}}}$
(78)

 Figure 90. Nodal displacement of beam-column element in 2D [163]

The generalized deformations of the problem are the axial strain ε(x) and the curvature about the z-axis φz (x). Under the assumptions that displacements are small and plane sections remain plane the section deformations "generalized deformations" a(x) are related to the nodal displacements by:

 ${\displaystyle \left\{\left.a(x)\right\}\right.=\left\{\left.{\begin{matrix}\epsilon (x)\\\varphi z(x)\end{matrix}}\right\}\right.=}$${\displaystyle \left\{\left.{\begin{matrix}{u}^{'}(x)\\{v}^{'}(x)\end{matrix}}\right\}\right.=}$${\displaystyle \left[B\right]\left\{\left.\delta \right\}\,\right.}$
(79)

Where ${\textstyle \left[B\right]=\,\left[dN\right]}$ , the derivative of the interpolation function matrix. Using the principal of virtual displacement of the principal of minimum potential energy, the element stiffness matrix can be writing as the integral of the section stiffness matrix ${\textstyle \left[k\right]}$ :

 ${\displaystyle \left[K\right]=\int _{0}^{L}{\left[B\right]}^{T}\left[k\right]\left[B\right]{d}_{x}}$
(80)

The section stiffness matrix ${\textstyle \left[k\right]}$ relates the section forces with the corresponding deformations:

 ${\displaystyle \left\{\left.D(x)\right\}\right.=\left[k\right]\left\{\left.a(x)\right\}\right.}$
(81)

Where ${\textstyle \left\{\left.D(x)\right\}\right.}$ is the vector of the section forces "the generalized stresses"; the axial force N(x) and the bending moment Mz(x) at section x.

The main shortcoming of this classical finite element displacement-based approach is their inability to describe the behavior of the member near its ultimate resistance and after the onset of strain softening, this is because of assumption of cubic interpolation functions, which result in a linear curvature distribution along the element. This assumption leads to satisfactory results under linear or nearly linear response. However, when the reinforced concrete member undergoes significant yielding at the ends, the curvature distribution becomes highly nonlinear in the inelastic region. This requires the use of a very fine discretization in the inelastic regions of stiffness based elements. This shortcoming has been demonstrated in the interesting study developed by Zeris and Mahin 1988 and 1991 [164], [165]; by studying the softening behavior of a cantilever beam it has been observed that assumed curvature distribution deviates significantly from the actual distribution during element softening, as the sharp jump in the curvature value near the fixed end. Results of this study are shown in Figure 91 as presented by [166]

 Figure 91. softening behavior for cantilever beam [166], (A) member and loading, (B) Moment distribution, (C) Curvature distribution, (D) Moment-Curvature relation.

It has been agreed that a more convenient format for distributed plasticity models is in terms of force method (flexibility method) where the unknowns are the forces, which implies assuming force interpolation functions. This approach uses the generalized deformation of the element instead of the displacement at the joints. The generalized deformation for the beam-column element represented in Figure 92 are: end rotations relative to the chord and the axial differential displacement, these deformations are demonstrated in

Figure 92.a. Accompaniment element forces "called as the generalized forces" are demonstrated in Figure92.b

 (a) Generalized deformations (b) Generalized forces Figure 92. Beam-column element in force-based method

Assuming that the bending moment distribution inside the element is linear and that the axial force distribution is constant, the section forces "generalized stresses" can be related to the element force in a similar concept of equation (76) as follow:

 ${\displaystyle \left\{\left.D(x)\right\}\right.=\left[b\right]\left\{\left.Q\right\}\,\right.}$
(82)

${\textstyle \left[b\right]}$ is the matrix containing the force interpolation functions, ${\textstyle \left\{\left.Q\right\}\,\right.}$ is the vector of the generalized force demonstrated in

Figure92.b. Equation (92) can be written as follow:

 ${\displaystyle \left\{\left.{\begin{matrix}N(x)\\{M}_{Z}(x)\end{matrix}}\right\}\right.=\left[{\begin{matrix}0&0&1\\{\frac {x}{L}}-1&{\frac {x}{L}}&0\end{matrix}}\right]\,\left\{\left.{\begin{matrix}{Q}_{1}\\{Q}_{2}\\{Q}_{3}\end{matrix}}\right\}\right.}$
(83)

The application of the virtual force principle yields the element flexibility

 ${\displaystyle \left[F\right]=\int _{0}^{L}{\left[b\right]}^{T}\left[f\right]\left[b\right]{d}_{x}}$
(84)

Where ${\textstyle \left[f\right]}$ is the section flexibility, equation (81) can be rewritten as:

 ${\displaystyle \left\{\left.a(x)\right\}\right.=\left[f\right]\left\{\left.D(x)\right\}\right.\,}$
(85)

The advantage of this formulation is satisfying the element equilibrium in a strict sense. Whatever the material nonlinearity take place at the section level and even as the element starts softening when deformed beyond its ultimate resistance, the assumed internal force distributions are exact [166].

In1984 the concept of fiber element was introduced by kaba and Mahin [167]. In this model the element is discretized into longitudinal sections called control sections and each section is divided into several fibers. The constitutive relation of the section is not explicitly specified, but is derived by integrating the response of the fibers, which follow the uniaxial stress-strain relation of the particular material as shown in Figure 93.

 Figure 93. Fiber element model [166]

This concept has been well accepted and approved as a good approach to simulate the nonlinear behavior of reinforced concrete beams and columns. Two main models have been developed according to the concept of fiber element. Displacement-Based and Force-Based models termed in this research as DB and FB respectively. Both approaches allow plastic hinges to form at any location and account for axial moment interaction by integrating the force-deformation response at sections along the element length. The number of sections and their location is determined by the numerical quadrature rule, such as those based on Gauss quadrature are used to integrate the element force-deformation relationship.

DB follow the standard finite element approach, in which the element displacement field is expressed as a function of the nodal displacement using displacement interpolation functions. The displacement field is approximate; thus several displacement-based elements are required along the length of a frame member to represent the deformations in a plastic hinge region [168] [169]. In contrast, FB approach interpolates the section forces in terms of the basic forces, satisfying equilibrium even in the range of nonlinear material response [123] [170]. The primary advantage of FB formulation over the DB is the ability to use one force-based element to simulate the material nonlinear response of a frame member, compared with several displacement-based elements, thereby keeping the number of degrees of freedom in the structural model to a minimum.

While distributed plasticity formulations model variations of the stress and strain through the section and along the member in more detail, important local behaviors, such as strength degradation due to local buckling of steel reinforcing bars, or the nonlinear interaction of flexural and shear, or the degradation due to bond-slip degradation are difficult to capture without sophisticated and numerically intensive models. On the other hand, concentrated hinge/spring models, may be better suited to capturing the overall response, through calibration using member test data on phenomenological moment-rotations and hysteresis curves as described in the previous subsection.

## 3 Structural component-based vs continuum mechanics-based approaches

### 3.1 Introduction

This section presents a comparison between a models derived from structural component-based approaches and the more complex models that are based on continuum mechanics formulation. Models involved in this section are the simple ones based on concentrated and distributed plasticity approaches described in the previous section, and the continuum-based model developed in Part I of this dissertation. The monotonic and cyclic experiments described in section ‎6 of Part I are simulated with both models, obtained results are compared with the ones from the experiments.

The one floor RC frame experiment [Pires 1990] described in subsection ‎6.2.1 is simulated with lumped and distributed plasticity models. The concentrated plasticity model is implemented in software package SAP-2000 v.16 [171] by using the displacement formulation [172], without accounting for shear-generated deformation. Second-order effects are taken into account. Zero-length flexural and axial-flexural hinges are assigned to beams and columns, respectively. As well, shear hinges are considered for both types of members. Hinges are located at faces of joints. In flexural and axial-flexural hinges, parameters of initial and further branches of the assigned backbones curves are obtained according [173] and [12], respectively. In shear hinges, parameters of initial branch are obtained according to [173], and after-yield behavior is assumed to be totally brittle. Neither shear-flexural interaction nor shear-axial-flexural interaction is considered. Figure 94 shows the capacity curves obtained by the aforementioned concentrated plasticity model with the experimental one. Results from Figure 94 shows that model is able to capture only the initial stiffness of the frame. The rest of the behavior is not captured by the model, degradation due to flexural cracking is not represented, nor the maximum capacity. As a first observation, the concentrated plasticity model with the parameters from [173] and [12] provides poor representation of the real behavior of this one floor RC frame.

 Figure 94. Capacity curves for [Pires 1990] RC frame. Experimental vs. concentrated plasticity model.

[Pires 1990] RC frame is also simulated by the distributed plasticity models. This involves the fiber element with both DB and FB approaches. The FEM package SeismoStruct V6.5 [112] is used to perform this analysis. In the DB models, each member is discretized with four 2-node finite elements. In the FB models, a single element with five integration points represents each member. For both FB and DB models, sections are discretized into 250 fibers. Second-order effects are accounted for. As in the concentrated plasticity model, the interaction between shear and flexure is not taken into consideration. Each fiber along the cross section is associated with a uniaxial stress-strain relation; concrete fibers are assigned to Mander et al. 1988 model [174] considering the confinement effect for the fibers in the inside the confinement zone, steel fibers are assigned to Menegotto and Pinto 1973 model [175] with isotropic hardening.

 Figure 95. Capacity curves for [Pires 1990] RC frame. Experimental vs. distributed plasticity models.

Obtained results by the DB and FB models are plotted with the experimental one in Figure 95. Plots from Figure 95 show a good agreement between the numerical models and experimental behavior, initial stiffness and subsequent degradation due to flexural cracking are well represented till reaching the maximum capacity. After reaching the maximum capacity, both DB and FB models started a rapid strength reduction, while the real behavior showed some ductility. Comparing between FD and DB response, the FB approach showed better agreement in representing the maximum strength.

For comparison with the results obtained by the continuum-base model in subsection ‎6.2.1, Figure 96 displayed the plots obtained in Figure 50, Figure 94 and Figure 95. Plots from Figure 96 show the superiority of the continuum-based model to describe the response of the RC frame, results of the structural component-based model show that the concentrated plasticity model describes satisfactorily the initial slope, but fails to predict the cracking and, therefore, the onset of the overall yielding; force and displacement ductility are underestimated because of the conservative assumptions in the predefined plastic hinges. Distributed plasticity models perform better, particularly the one based on FB formulations.

 Figure 96. Capacity curves for [Pires 1990] RC frame. Experimental vs. structural and advanced models.

The two floor RC frame experiment [Vecchio and Emara 1990] described in subsection ‎6.2.2 is also simulated similarly to the previous experiment. Concentrated and distributed plasticity models are used to describe the behavior of the two floor RC frame.

Obtained capacity curve by the contracted plasticity model is displayed with the experimental behavior in Figure 97. Plots from Figure 97 drive similar observation to Figure 94. The concentrated model describes satisfactorily the initial slope, but fails to predict the cracking and, therefore, the onset of overall yielding. The maximum capacity in terms of force is underestimated because of the conservative assumptions in the predefined plastic hinges. The final failure is earlier because the actual ductility of members is also underestimated by the assumed moment-rotation laws

 Figure 97. Capacity curves for [Vecchio and Emara 1990] RC frame. Experimental vs. concentrated plasticity model.

Similarly, to the previous experiments, the two floor RC frame is simulate as well using the fiber element with DB and FB approaches. Figure 98 shows the obtained capacity curves and experimental ones. In general, better agreement is obtained by the distributed plasticity models, particularly FB model captured the maximum strength and the ductility phase, however, earlier degradation is predicted by DB and FB models.

 Figure 98. Capacity curves for [Vecchio and Emara 1990] RC frame. Experimental vs. distributed plasticity models.

For comparison with continuum-based model, Figure 104 collects the plots of Figure 58, Figure 97 and Figure 98. Plots from Figure 99 show again the superior ability of the proposed continuum-based model to reproduce the experimental results along the whole displacement range. Regarding the lumped plasticity model, accuracy can be considered satisfactory, given the important simplifications involved in that model. Regarding the distributed plasticity models, Figure 99 shows that they perform better, particularly FB. DB model exhibits less accuracy, because the mesh is too coarse. However, these models cannot capture adequately the gradual progression of the global softening after the initial overall yielding, due to lack of consideration of the concrete tensile strength. Noticeably, the negative slope of the final branches is due to second-order effects.

 Figure 99. Experimental and simulated capacity curves for [Vecchio and Emara 1990] RC frame. structural and advanced models

An overall comparison between the experimental observation and the ones obtained by the structural and continuum-based models are demonstrated in Table 12.

Table 12. Experimental vs. numerical results [Vecchio and Emara 1990] frame. structural and advanced models
 SAP2000 SeismoStruct FB SeismoStruct DB ABAQUS Experiment Ultimate Force (KN) 302 322 361 333 332 First Plastic Hinge Developed at Beam 1st Story Beam 1st Story Beam 1st Story Beam 1st Story Beam 1st Story Failure Mechanism 6 Plastic Hinges 6 Plastic Hinges 6 Plastic Hinges 6 Plastic Hinges 6 Plastic Hinges Joints Damage N/A N/A N/A YES YES Force at First Crack in Beam First Story N/A N/A N/A 61 52.5 Force at First Crack in the column at the base N/A N/A N/A 135 145 Force at first yielding in the column's reinforcement N/A 295 327 305 323 N/A Not Available

This subsection presents similarly compression to the previous subsection. The cyclic loading experiments described in subsection ‎6.3 of Part I are simulated with structural component-based models. Results are compared with the real behavior and the one previously captured by the proposed continuum mechanics-based model.

As an observation from the previous subsection, the fiber element models showed higher capacity in representing the real behavior. Hence, the concentrate plasticity models are excluded from the cyclic loading simulation and just fiber models are implemented. As well, the advanced FEM code OpenSEES [176] is used to implement the fiber model. The open source OpenSEES provides wide range or materials and elements models, as well different algorithms are available to perform nonlinear static and dynamic analysis with powerful convergence controlling.

Since the cyclic behavior of RC structures involves new degradation modes, additional elements are incorporated in the simulation. Standard fiber models nor the concentrated plastic models are able to capture the bond-slip degradation and the results represented in the so called pinching effect. To overcome this shortcoming, the strain penetration model developed by Zhao and S. Sritharan 2007 [124] is assigned into a zero length element at the elements edges where bond-slip effects are expected to occur. The strain penetration model [124] is the most common strategy to capture bond-slip effect in structural component-based models; it is assumed that bars are well anchored and slip is entirely due to strain penetration in the anchorage zone. Then, rotations at member ends are generated. A constitutive hysteretic stress-slip law for steel bars based on pull-out tests is generated; it is used for zero-length elements situated at column bottom and is integrated into a global fiber model. Figure 100.a shows the hysteretic stress-slip law for steel bars, while Figure 100.b displays the implementation of the zero length in fiber model. This model is implemented in OpenSEES [176].

 (a) Influence of pinching parameter on bar stress vs slip (b) Integrating fiber element with strain penetration model Figure 100. Strain penetration model for bond-slip effect Zhao and S. Sritharan 2007 [124]

The RC column experiment Tanaka 1990 [98] described in subsection ‎6.3.1 is discretized with one 2-node finite elements using the force-based formulation, and one zero-length element to represent bond-slip. Sections are discretized in fibres using Concrete01 material for both confined and unconfined concrete. while Steel02 material is used for reinforcement. Figure. 101 shows the uniaxial cyclic behavior for both concrete and steel materials model OpenSEES code.

 (a) Concrete01 model (b) Steel02 model Figure 101. Uniaxial cyclic model for concrete and steel OpenSEES

Concrete confinement is taken into consideration using the modified kent and park 1971 [177] model, this model has been developed by Scott et al. 1982 [178] by incorporating the confinement effects on the maximum strength and strain. Stress-strain response for confined and unconfined concrete according to Scott model is shown in Figure 104.a.

 (c) Scott et al. 1982 concrete model [178] (d) Kent and Park 1971 [177] Figure 102. Confinement effect on Uniaxial concrete stress-strain relation

The following formulations are used to calculate Concrete01 model parameters (MPa)

 ${\displaystyle {\rm {For\,}}{\epsilon }_{c}\,\leq \,0.002K\qquad f_{c}=Kf_{c}^{'}\left[{\frac {2{\epsilon }_{c}}{0.002K}}-\left({\frac {{\epsilon }_{c}}{0.002K}}\right)^{2}\right]}$ ${\displaystyle {\rm {For\,}}{\epsilon }_{c}>\,0.002K\qquad f_{c}=Kf_{c}^{'}\left[1-{z}_{m}\left({\epsilon }_{c}-0.002K\right)\right]\,\geq 0.{2Kf}_{c}^{'}}$ ${\displaystyle {z}_{m}={\frac {0.625}{{\epsilon }_{50h}+{\epsilon }_{50u\,}-0.002K}}\qquad \qquad K=1.25\left[1+{\frac {{\rho }^{''}{f}_{yh}}{{f}_{c}^{'}}}\right]}$ ${\displaystyle {\epsilon }_{50h}={\frac {3}{4}}{\rho }^{''}{\sqrt {\frac {{b}^{''}}{S}}}\qquad {\epsilon }_{50u}={\frac {3+0.29{f}_{c}^{'}}{145{f}_{c}^{'}-1000}}\qquad {\rho }^{''}={\frac {2\left({b}^{''}+{d}^{''}\right){{A}_{s}}^{''}}{{b}^{''}{d}^{''}S}}}$
(86)

In equations (86), f'c and fyh are concrete compressive strength and steel yield strength respectively of transverse bars, ρ´´ is volumetric ratio of confining hoops to volume of concrete core measured to the outside of the perimeter hoops, where b´´ and d´´ are the width and depth of the confined core respectively, As´´ is the cross-sectional area of the hoop bar and S is the center to center spacing of the hoops. ε50u and ε50c are the strain corresponds to 50% of concrete maximum compressive strength for unconfined and confined concrete respectively, while ε50h is the differencebetween the previous two values as shown in Figure 104.b. Implementing concrete model Concrete01 requires the strain at the maximum strength to be calculated ε0 as 2f'c/Ec where Ec is concrete modulus of elasticity, as well the ultimate compressive strain shall to be defined, in Figure 104.b this corresponds to ε520c.In this work εcu is used to express the ultimate compressive strain, it’s calculated according to the formulation proposed by [178] as :

 ${\displaystyle {\epsilon }_{cu}=0.004+0.9\,\left({\frac {{f}_{yh}{\rho }_{sh}}{300}}\right)}$\qquad ${\displaystyle {\rho }_{sh}={\frac {{A}_{sh}}{Sd}}}$
(87)

ρsh and fyh are the transverse reinforcement ratio and yield stress respectively Ash is the total area of transverse bars acting in the considered direction of the section, d is the effective depth of the section.

Steel02 material model is based on the isotropic hardening model for steel developed by Menegotto and Pinto 1973 model [175]. Parameters of this model are obtained from OpenSEES software.

The zero-length element is discretized as any section along the length of the element, Steel02 material is replaced with the strain penetration material. This model follows the phenomenological hysteretic law proposed displayed in Figure 100.a; one of the parameters that governs reloading branch is the pinching factor Rc shown in Figure 100.a. This parameter modifies completely the response of the zero-length element; consequently, this will alter significantly column behavior. Rc = 1 represents perfect bond situation and, thus, no pinching effect; lower values of Rc correspond to significant pinching. For the simulation of Tanaka column, Rc = 0.6 has provided the best results.

Figure 103 displays comparisons between experimental and numerical hysteresis loops obtained with the abovementioned structural component-based model.

 Figure 103. Hysteretic force displacement response of [Tanaka 1990] column. Experimental vs. structural component-based model

Plots from Figure 103 show that the structural component-based model is able to capture acceptably the envelope of column response, however, the maximum strength is not well captured nor the degradation upon unloading and reversing the loading direction as well the pinching effect. For comparison with results obtained previously in the Part I by the continuum-based model, Figure 104 collects the plots from Figure 67.b and Figure 103

 Figure 104. Hysteretic force displacement response of [Tanaka 1990] column. Experimental vs. structural and advanced model

Results from Figure 104 show that the simplified model reveals less accurate in comparison with the continuum-based model.

Similarly to the simulation of Tanaka column, the one floor RC frame described in subsection ‎6.3.2 of Part I is analyzed with the aforementioned structural component-based model. Herein the zero-length element is assigned in three places where bond-slip is expected to be significant: columns base, columns top, and beam ends.

As discussed previously, this model is particularly sensitive to pinching factor Rc; As a preliminary attempt, the calibrated value from the simulation of the bridge column Rc = 0.6 is used in the three locations. Analogously to Figure 103 comparison between experimental and numerical hysteresis loops obtained by the structural component-based is demonstrated in Figure 105.

 Figure 105. Hysteretic force displacement response of [Pires 1990] frame. Experimental vs. structural component-based model (constant Rc)

Plots from Figure 105 show that the response obtained by the structural component-based model is shifted from the experimental response, this can be explained by the high flexibility generated by the low values of Rc in the three locations.

To stabilize the numerical response, the value of Rc must be calibrated in each location individually. After several attempts, better results are obtained by assigning the values 0.8, 0.9 and 0.4 to the columns base, columns top, and beam ends respectively. Analogously to Figure 105, Figure 106 shows a comparison between the experimental response and the one obtained with the aforementioned Rc values.

 Figure 106. Hysteretic force displacement response of [Pires 1990] frame. Experimental vs. structural component-based model (varying Rc)

Plots from Figure 106 show that the structural component-based model shows very good agreement in terms of capturing the unloading stiffness and pinching effect, however the maximum strength is not well captured, as well as the strength reduction.

Analogously to Figure 104, Figure 107 displays comparisons between experimental and numerical hysteresis loops obtained with the proposed continuum-based and the simplified one with the abovementioned values of Rc. Figure 107 collects the plots from Figure 77.b and Figure 106. Results from Figure 107show that the structural component-based model has better capacity in representing the pinching effect, maximum capacity and strength reduction are better represented by the continuum-based model.

 Figure 107. Hysteretic force displacement response of [Pires 1990] frame. Experimental vs. structural component-based (varying Rc) and advanced model