“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
El comportamiento estructural de edificios y otras construcciones bajo severas excitaciones sísmicas es muy complejo e implica temas como, la interacción sueloestructura, grandes esfuerzos y desplazamientos, daños, plasticidad y el comportamiento de la estructura cerca del colapso. Por otra parte, en estructuras de hormigón armado, existen varios modos de fallo y de degradación: agrietamiento, aplastamiento y desprendimiento del hormigón, plastificación y extracción de las armaduras traccionadas y plastificación y pandeo de las armaduras comprimidas. Además, otras circunstancias hacen que la situación sea más alarmante: dada la creciente conciencia y preocupación por el enorme riesgo sísmico mundial, la ingeniería sísmica ha experimentado en los últimos años avances sustanciales, para lo cual se han propuesto nuevas estrategias de análisis y diseño, lo que conduce a desarrollos relevantes. Estos desarrollos se basan en pruebas y simulaciones numéricas basadas principalmente en modelos simplificados referidos en este trabajo como modelo basados en la estructura, resultando un costo computacional moderado. Por lo tanto, existen una gran necesidad de verificar la fiabilidad de los nuevos desarrollos en comparación con los análisis realizados utilizando herramientas de simulación más avanzadas y con ensayos. Esta trabajo se organiza en dos partes; en la primera se describe un modelo preciso basado en la mecánica del medio continuo y en la segunda se presenta otro modelo más simplificado basado en los componentes de la estructura.
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 3D 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.
Segunda parte. El modelo desarrollado en el Primer Parte ha sido comparado con modelos simplificados basados en los componentes estructurales de uso común en la ingeniería sísmica, los resultados mostraron la superioridad del modelo propuesto para predecir el comportamiento real de los elementos y pórticos RC altamente dañados, capturando la reducción de la resistencia, la degradación de la rigidez y el efecto pinzamiento (“pinching”). Sin embargo, algunos de los modelos basados en componentes estructurales han mostrado un desempeño aceptable teniendo en cuenta el costo computacional de la ley en comparación con el modelo avanzado basado en la mecánica del medio continuo. Con de esta conclusión, este parte de este trabajo presenta un estudio numérico sobre la relación entre los modos de deterioro nosimulados de pórticos de hormigón sin ductilidad y su capacidad última. Se ha desarrollado un modelo avanzado basado en los componentes de la estructura para simular el comportamiento dinámico no lineal de las estructuras sin ductilidad, teniendo en cuenta los modos de deterioro de flexión, corte y axial. El modelo desarrollado es numéricamente eficiente, siendo pues adecuado para el uso profesional en ingeniería sísmica. La capacidad del modelo desarrollado se verifica mediante la simulación del comportamiento dinámico no lineal de un edificio no dúctil existente y del edificio prototipo. Los resultados obtenidos muestran que el desarrollado, a pesar de su coste computacional moderado, detecta y reproduce con precisión el comportamiento dinámico no lineal de estructuras RC no dúctiles, capturando también los modos de deterioro que no contemplan los modelos más simplificados. La comparación con los resultados de modelos más simplificados resalta la importancia de los modos de falla no considerados en el comportamiento de cada elemento y en los mecanismos generales de colapso. Se aborda la relación entre los modos de fallo no simulados y el fenómeno de la “Resurrección Estructural”.
Under severe seismic excitation, structural behavior of buildings and other constructions is highly complex. It involves, among other issues, soilstructure interaction, large strains and displacements, damage, plasticity, and nearcollapse behavior. Moreover, in reinforced concrete structures, there are several coupled degradation and failure modes: cracking, crushing and spalling of concrete, yielding and pullout 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 componentbased 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 mechanicsbased 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 3D interface bondslip 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. Meshinsensitivity, 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 componentbased 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 componentbased models have shown an acceptable performance considering the law computationally cost in comparison with the advanced continuum mechanicsbased model. After this conclusion, this part presents a numerical study on the relation among the nonsimulated deterioration modes of the elements in nonductile RC frames and their final capacity. A structural componentbased model has been developed for simulating the nonlinear dynamic behavior of nonductile 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 nonductile 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 nonductile 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 nonsimulated failure modes and the socalled "Structural Resurrection" is addressed.
a_{c} / a_{t} / b_{c} / b_{t}: dimensionless coefficients in equations
b: ratio
c / c_{1} / c_{2}: cohesion / coefficients in the uniaxial tensile behavior of concrete
d / d_{c} / d_{t}: damage variable / compression damage variable / tension damage variable
b´´/d´´: are the width and depth of the confined core respectively
f / f_{cm} / f_{tm} / f_{c0} / f_{t0} / f_{ck}: 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 / f_{b0} / f_{c0 }/ f_{cm} / f_{tm} / f_{c0} / f_{t0} / f_{ck}/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.
f_{y} / f_{u}: longitudinal reinforcing bars stresses yield / ultimate
f_{yh} / f_{u}: transverse reinforcing bars yield stress
g_{c} / g_{t}: compressive / tensile energies per unit volume dissipated by damage along entire deterioration process
/ : weighting factors accounting for stiffness recovery
l_{eq}: mesh size (finite element characteristic length)
r*: stress state; for uniaxial stress r*(σ_{11}) = 1 for tension and r*(σ_{11}) = 0 for compression
s_{c} / s_{t}: coefficients accounting for stress state and stiffness recovery effects
p: hydrostatic pressure stress
q: Von Misesequivalent effective stress
r*: stress state; for uniaxial stress r*(σ_{11}) = 1 for tension and r*(σ_{11}) = 0 for compression
s_{c} / s_{t}: coefficients accounting for stress state and stiffness recovery effects
w / w_{c}: crack opening / crack opening at fracture
D_{b}: reinforcement bar diameter
E / E_{0} / E_{ci}: modulus of deformation / undamaged modulus of deformation / tangent modulus of deformation of concrete for zero stress
E_{s} / E_{sh}: steel modulus of elasticity / slope of hardening branch
EI_{flex} : column effective flexural rigidity
F: loading (yield) function
G / G_{ch} / G_{F}: flow potential / crushing energy per unit area / fracture energy per unit area
H: MohrCoulomb yield surface function
I_{1}: first invariant of stress tensor
J_{2} / J_{3}: second / third invariants of deviatoric stress tensor
K_{c}: ratio of second stress invariants on tensile and compressive meridians
K_{deg}: shear degradation stiffness (Elwood 2005 mode [1])
C_{slip}: Stiffness of linear slip spring
L: column height
L: mesh size for the spring bondslip model.
M_{SP}: moment at spalling of concrete
F_{bs} : is the equitant bond force
S_{1 }/ S_{2 }/ S_{3}: relative normal displacement / longitudinalslip / transverse slip
S_{peak }/ : peak slip / maximum positive slip / maximum negative slip / cumulated slip
S_{R}:clear spacing between bar ribs
S: is the center to center spacing of the hoops
A_{s}´´: is the crosssectional area of the hoop bar
A_{g}/A_{eff}: are the gross section area/ the effective shear section area.
R_{c} : is the pinching parameter of [124] model.
M: fourthorder tensor for anisotropic damage model
I identity tensor
ε / ε_{c} / ε_{t} / ε^{el} / ε^{pl} / ε_{cm} / ε_{tm}: strain / compression strain / tensile strain / elastic strain / plastic strain / strain at compressive strength / strain at tensile strength
/ / / / / / / : 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 }/ ε_{y} /ε_{sh }/ε_{u}: axial steel strain / yielding steel strain / hardening steel strain / ultimate steel strain
/ : crushing compressive strain / cracking tensile strain
/ : 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 postcapping stiffness (Ibarra et al model [2])
λ: is the energy dissipation parameter Ibarra et al model [2])
ρ: octahedral radius
σ / σ_{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)
/ : 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
: Steel material parameters
Φ_{y}: Yield curvature
Δ_{y}/Δ_{flex}/Δ_{slip}/Δ_{shear}: 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 flexureshear and flexureshearaxial deterioration (Elwood 2005 mode [1])
This section presents a brief introduction and a concise historical review of earthquakeresistant 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 socalled 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 earthquakeresistant 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 earthquakeresistant 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 overresistance 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 timehistory responses; then, the maximum values will be selected, they would represent the design demands. This formulation is often referred to as earthquakeresistant 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.
Earthquakeresistant 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 nonstructural 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 singledegreeoffreedom system and of the input acceleration. Figure 2 shows an elastic model of a singledegreeof freedom system undergoing a horizontal ground motion z_{g}.
Figure 2. Elastic singledegreeoffreedom 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 (degreeoffreedom) and z_{g} 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

(1) 
By dividing both sides by m, (1) becomes

(2) 
In this relationship, is the undamped natural frequency of the system and ζ is the critical damping factor. These coefficients are given by

(3) 
The damped natural frequency is related to and to ζ by

(4) 
It is remarkable that, unless the damping ζ takes extremely high values, and are nearly coincident.
The acceleration, velocity and displacement spectra are obtained, for each input z_{g}(t), as the maximum values of the absolute acceleration (where ), relative velocity 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]:

(5)  

(6)  

(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 (EW 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 ICA2EW 
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 nonzero 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:

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 [NCSE02 2002].
Figure 4. Design acceleration spectrum [NCSE02 2002]. 
The spectrum shown in Figure 4 consists of three branches: a linearly increasing one, a constant one and a hyperbolically decreasing one. Periods T_{A} and T_{B} 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) shortperiod 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 T_{A} and T_{B}: for stiff soil the range of building periods whose motion is highly amplified (in between T_{A} and T_{B}) 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 S_{a} are dimensionless).
Performancebased earthquakeresistant 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 nonstructural) 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 nonstructural damage, the following four levels of performance (“Performance States”) are defined according to SEAOC 1995 [9]:
More recently, another similar classification is considered according to ATC40 [10] and FEMA 356 [12];
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.
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 socalled MCE (“Maximum Considered Earthquake”) corresponds to a return period of about 2475 years. The relationship between the return period T and the probability p_{n} of being exceeded n years is given by the expression ; it is often used to indicate the severity of an earthquake by the probability p_{50} to be exceeded in 50 years, for example, in the case of MCE
is and in the case of an earthquake “Rare” is .
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.
Nonlinear static analysis (Pushover)
The method of earthquakeresistant 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 ATC40 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 pushover. Figure 6 shows an example of a capacity curve obtained from a pushover analysis.
Figure 6. Capacity curve obtained from pushover analyses ATC40 1996 [10] 
In the pushover 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 pushover 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 PushOver 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 S_{a} (vertical axis) vs. the relative displacement spectrum S_{d} (horizontal axis). This type of representation is commonly known as “AccelerationDisplacement Response Spectra” (ADRS).
The methods mostly used to obtain the target displacements are: Capacity Spectrum Method ATC40 1996 [10], Displacement Coefficient Method FEMA356 2000 [12], Equivalent linearization method FEMA440 2005 [16], Modified displacement coefficient method FEMA440 2005 [16], Modified Capacity Spectrum ATC40 1996 [10].
The analysis pushover 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 pushover 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 (interstory 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 nonstructural 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 (pushover) 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, secondorder 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 earthquakeresistant 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 pushover analysis cannot take into account the accumulated plastic strain, the socalled 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 secondorder 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 socalled IDA curves. These representations consist of capacity curves similar to the result of the pushover 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 S_{a}(T_{1}, 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 (interstory 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 performancebased 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 nondamaged 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 performancebased 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 bondslip 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 beamcolumn elements, same concepts can be applied for other element (e.g. braces and shear walls).
Figure 10. BeamColumn nonlinear models [19] 
Concentrated Plasticity Models are the simplest formulation, representation of such models are shown in Figure 10a,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 zerolength or spring element at location of the inelastic deformation, the behavior of this nonlinear element is governed by the socalled “backbone” or capacity curve that can be represented by momentrotation or forcedeformation 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 bondslip 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 10c 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 momentcurvature relationships or explicit fibersection integrations that enforce the assumption that plane sections remain plane. The inelastic hinge length may be fixed or variable, as determined from the momentcurvature 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 10d. 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 stressstrain characteristics in the cross sections, each fiber is assigned to a material (e.g. steel, confined and unconfined concrete). The planesectionsremainplane assumption is enforced, where uniaxial material “fibers” are numerically integrated over the cross section to obtain stress resultants (axial force and moments) and incremental momentcurvature and axial forcestrain 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 10e. 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 “bondslip” 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 ComponentBased Formulation represented by the models in Figure 10a,b,c,d, and the most complex and general formulation is the ContinuumBased Formulation represented by Figure 10e. 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.
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 mechanicsbased formulations, and the simplified structural componentbased models. This research will focus on RC frames behavior.
This global objective will be pursued through the following particular objectives:
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 mechanicsbased 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 bondslip model developed by [20] and implemented in ABAQUS. Combining the methodology for calculating damage evolution in CPDM, as well the modeling scheme of the bondslip model, results in an integrated continuumbased 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 plasticdamage 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 bondslip 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 componentbased 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 componentbased model for nonlinear dynamic simulation of nonductile RC framed structures. A numerical study is conducted using the developed model on the relation among the nonsimulated deterioration modes of the elements of nonductile 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 componentbased 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 mechanicsbased model developed in Part I. Chapter 4 presents the developed model to simulate the nonlinear dynamic behavior of nonductile reinforced concrete structures. Chapter 5 verifies the capacity of the presented model in chapter 4 by performing numerical simulation of a two nonductile 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.
Under severe seismic excitation, structural behavior of buildings and other constructions is highly complex. It involves, among other issues, soilstructure interaction, large strains and displacements, damage, plasticity, and nearcollapse behavior. Moreover, in reinforced concrete structures, there are several coupled degradation and failure modes: cracking, crushing and spalling of concrete, yielding and pullout 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 componentbased 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 mechanicsbased 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 damageplastic model. Regarding steel bars, the basic plasticity model is described with its implementation. The complex relation between concrete and steel bars "bondslip" 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 bondslip model in 3D FEM simulation described in subsection 5.4. As a final contribution, a new approach to integrate CPDM and bondslip model in simulating the cyclic and monotonic behavior of RC structures is presented in subsection 6.1
Quasibrittle materials, as concrete, exhibit nonlinear stressstrain response mainly because of microcracking. Cracks are oriented as the stress field and generate the failure modes. In tension, failure is localized in a narrow band; stressstrain 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; stressstrain 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]. Damagebased 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.
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.
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 microcracks 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.
Uniaxial compression under cyclic loading shows important characteristics as shown in Figure 11.a; the developed stress follows the path of monotonic stress response, initial stiffness in each unloading and reloading stiffness faces a progressive deterioration with increasing the demand. The phenomena of stiffness degradation upon increasing the demand is referred as damage. Figure 11.b; shows a plot a data from [41] of normalized unloading stiffness as a function of normalized plastic compressive deformation. Plastic deformation is defined as the deformation that is not recovered upon unloading to zero compressive stress and this deformation is normalized with respect to the deformation at approximately zero compressive strength.
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 aggregatemortar 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.
Figure 12. Uniaxial concrete response under cyclic tension loading [42] 
Since plain concrete in RC elements is usually subjected to multidimensional 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.
Tow envelopes must be considered to describe the behavior of concrete; elastic limit surface defining the elastic region and the failure surface characterizing the maximumstrength 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; I_{1} the first invariant of the stress tensor, J_{2} and J_{3 }the second and third invariants of deviatoric stress tensor respectively. HaighWestergaard 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.

(9) 

(10) 

(11) 

(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.
In principle, it is desired that the abovementioned 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.
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 plasticitybased 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

(13) 
Elastic constitutive relationship follows Hook's law

(14) 
is the elastic stiffness tensor, and 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 nonassociate. 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 yieldfailure 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 MohrColumb and DrukerPrager [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 MohrCoulomb criterion is defined as follows:

(15)  

(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 plane as shown in Figure 15.
Figure 15. MohrCoulomb and DruckerPrager Failure Criteria [42] 
The DruckerPrager criterion represents moderately well the response of plain concrete
subjected to multiaxial compression and provides a smooth yield surface Figure 15. This criterion is defined as follows:

(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. DruckerPrager 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 DruckerPrager criterion with experimental data shows that while
the criterion may be used to represent the response of concrete subjected to multiaxial
compression, the model overestimates the capacity of concrete subjected to compressiontension or tensiontension type loading. Variation in concrete response under various load regimes has been addressed by a number of researcher through the use of multisurface 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 threedimensions by Chen and Chen 1975 [31] and Lubliner et al. 1989 [30].
The main shortcoming of MohrColumb and DrukerPrager surfaces is that they assume a linear relationship between J_{2} and I_{1 }(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 J_{2} and I_{1} 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 nonassociated 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.
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 and its abscissa x along the direction .
Figure 18. Definition of damage variable [61] 
The damage value D(M, , x) at point M in the direction and at abscissa x is defined as:

(18) 
δS is the area of the intersection considered plan and the RVE, and 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 . If microcracks and microcavities are uniformly distributed in the RVE, it's adequate to assume isotropic damage and damage variable D(M, , 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 S – SD. 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 / (S − SD) = σ (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:

(19) 
Mis a fourthorder 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 = (1d) I, where I is an identity tensor. The equation (19) can be written for scalardamage model:

(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 nondamaged 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 nondamage 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 damagebased model is controlled by damage variable. The evolution law of damage, which plays a very important role in any damagebased 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], ).
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. [2123,2529]; these models have shown good performance in capturing concrete behavior in tests on fullscale 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,7180].
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 pressuredependent 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 multiaxial 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 damageplasticity model by Yazdani and Schreyer 1990 [81] 
Compared to other coupled damageplasticity 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 yieldfailure surface is defined and evolves with damage variables, for which the evolution laws are postulated. The yield surface in twodimensional 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 DruckerPrager 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 damageplasticity 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.
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 nonassociated multihardening 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.
Going back to Figure 20, it displays uniaxial stressstrain plots typical of plasticity, damage and damage plasticity models; loading branches are represented with solid thick lines and unloading / reloading branches are plotted with dashed thin lines. E_{0} is the initial (undamaged) elastic stiffness (deformation modulus), and ε^{el} and ε^{pl} are the elastic (recoverable) and plastic (irrecoverable) strain, respectively. Figure 20 shows that damage generates stiffness degradation [57] since the slope of unloading / reloading branch is (1 d) E_{0} where d is a damage variable ranging between 0 (no damage) and 1 (destruction).
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, and , linked to failure mechanisms under tension and compression loading, respectively. In Lee and Fenves 1998 [22] model, for uniaxial compression and tension, the stressstrain relation under uniaxial loading in displayed in Figure 20.c, can be written as:

(21) 

(22) 
where subindexes c and t refer to compression and tension, respectively.
For uniaxial cyclic loadingunloading conditions, Lee and Fenves 1998 [22] model assumes that the degradation in the elastic stiffness is given by

(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 (d_{c} and d_{t}, respectively):

(24) 
In equation (24), s_{c} and s_{t} are dimensionless coefficients accounting for stress state and stiffness recovery effects, being given by

(25)  

(26) 
In equations (25) and (26), is the first principal uniaxial stress (positive for tension), r* is a stress state parameter being = 1 for tension and = 0 for compression, and and are weighting factors ranging between 0 and 1. Factor accounts for reclosing of cracks after tensioncompression reversal; represents recovery of crushed concrete after compressiontension reversal. If ℎ_{c} = 0.9 and ℎ_{t} = 0 for example; this means that 90% of the cracks close upon tensioncompression reversal and the crushed concrete does not experience any recovery. Equations (25) and (26) show that and also range between 0 and 1.
For a better understanding of the effect of s_{c} and s_{t} coefficients, Figure 24 displays plots of uniaxial stressstrain loadingunloading behavior. The initial elastic branch with slope E_{0} 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 d_{c} = 0, r* = 1, and s_{c} = 1; therefore, equation (24) shows that d = d_{t}. Consequently, the linear unloading branch has slope (1 d_{t}) E_{0}. In the way to stress reversing point 3, cracks begin closing. After point 3, r* = 0, s_{c} = 1 h_{c}, s_{t} = 1, d_{c} = 0, and equation (24) shows that d = (1 h_{c}) d_{t}. Therefore, the slope of the ongoing compression segment of the branch depends on parameter h_{c}; three options are plotted in Figure 24: (i) h_{c} = 0 (no crack is closed) with slope (1 d_{t}) E_{0}, (ii) h_{c} = 0.5 (half of the cracks are closed) with slope (1 0.5 d_{t}) E_{0}, (iii) h_{c} = 1 (all cracks are closed) with slope E_{0}. Noticeably, in the third option (h_{c} = 1), there is no compressive strength reduction. At point 4, an unloading branch arises; at that point, r* = 0, s_{c} = 1 h_{c}, s_{t} = 1, and equation (24) shows that 1 d = (1 d_{c}) [1 (1 h_{c}) d_{t}] = 1 d_{c}. Point 5 correspond again to stress reversal; after it, provided that h_{t} = 0, the slope of the ongoing branch is equal to (1 d_{t}) (1 d_{c}) E_{0}. Point 6 is the peak for the reduced tensile strength; after it, cracking reinitiates and a new descending branch is generated.
Figure 24. Uniaxial loadingunloading law 
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 widelyopened cracks and bondslip effects concentrate. ℎ_{t} = 0.0 is assumed for all the sections (referrer to subsection 6.1).
For multiaxial condition, the stressstrain relationship is given by:

(27) 
In equation (27), is the elastic stiffness tensor, and and are the stress and strain tensors, respectively. Scalar damage variable d keeps same meaning than for uniaxial condition, although replacing scalar factor 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.

(28)  

(29) 
In equations (28) and (29), is the Macaulay bracket, is the hydrostatic pressure stress, is the Von Misesequivalent effective stress (effective stress accounts for stress divided by 1 d), and f_{b0} and f_{c0} are the biaxial and uniaxial compressive yield strengths, respectively; since f_{b0} ≥ f_{c0}, α ranges between 0 (f_{b0} = fc_{0}) and 0.5 (f_{b0} fc_{0}). is the maximum principal effective stress, and and are the effective compressive and tensile cohesion stress, respectively. and are defined as and . K_{c} is the ratio of second stress invariants on tensile and compressive meridians.
The plasticity model assumes nonassociated potential plastic flow. The flow potential G is the DruckerPrager hyperbolic function given by:

(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  (deviatory) plan at high confining pressure.
As discussed previously, K_{c} is the ratio between the magnitudes of deviatoric stress in uniaxial tension and compression; K_{c} ranges between 0.5 (Rankine yield surface) and 1 (Von Mises). In this study, K_{c} is obtained from the MohrCoulomb yield surface function in cylindrical coordinates [57]:

(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), , ξ = I_{1} / , and ; I_{1} is the first invariant of stress tensor and J_{2} and J_{3} are the second and third invariants of deviatoric stress tensor, respectively.
For = 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 K_{c} are

(32) 
By assuming that ϕ = 32° [57], last equation in (32) shows that K_{c} = 0.7. Figure 25 represents the yield surface in deviatory plan for several values of K_{c} 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 K_{c} 
Equations (30), (28) and (29) show that the concrete behavior depends on four constitutive parameters K_{c}, ѱ, f_{b0} / f_{c0}, ϵ; it can be assumed that ѱ =13° [82]. Table 3 describes the values used in this study.
K_{c}  Ѱ (º)  f_{b0} / f_{c0}  ϵ 
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.
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:

(33) 
d_{t} and d_{c} 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.

(34) 
is the equivalent strain, and the parameters A_{c, }A_{t,} B_{c,} B_{t }and K_{0} are derived from compression test on cylinders and bending test on beams.

(35) 
E_{0} is the original concrete stiffness modulus, ε^{pl} is the uniaxial plastic strain, b_{i} is a parameter that needs to be calibrated from cyclic uniaxial tension and compression tests.

(36) 
k_{d} is the equivalent damage strain variable defined in [85], e_{d0}, e_{d} and g_{d} are material parameters to be calibrated with experimental results.

(37) 
Where is the maximum strength of concrete. For constant confinement pressure, the damage function becomes:

(38) 
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.

(39) 
is defined as the pick value of strain when damage occurs, and are calculated according to [Zheng et al. 2012], A_{c, }A_{t,} B_{c,} B_{t} shall be calibrated form experiments.

(40) 
Parameters a_{c}, a_{t} needs to be calibrated from the uniaxial tension and compression stressstrain 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:
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 meshsensitivity. 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. Meshinsensitivity 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.
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:

(41) 
In equation (41), and are the crushing and cracking strains respectively, see Figure 27. Normalization coefficients g_{c} and g_{t} represent the energies per unit volume dissipated by damage along the entire deterioration process:

(42) 
Equations (41) and (42) show that d_{c} and d_{t} range between 0 (no damage) and 1 (destruction). Figure 26 describes the meaning of g_{c} and g_{t}.
a) Compression (g_{c})  b) Tension (g_{t}) 
Figure 26. Parts of energy dissipated by damage 
Noticeably, the energies per unit area and per unit volume are related by and ; and are material parameters defined as crushing and fracture energies and 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:

(43)  

(44) 
In equations (43) and (44), f_{c0} and f_{t0} and are the compressive and tensile stresses that correspond to zero crushing ( ) and to onset of cracking ( ), respectively; see Figure 26. As well, a_{c}, a_{t}, b_{c} and b_{t} are dimensionless coefficients to be determined. Replacing equations (43) and (44) in equation (42), provides the following relations among g_{c} and g_{t} and such coefficients:

(45) 
Coefficients b_{c} and b_{t} can be obtained by replacing and in equation (45):

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

(47)  

(48) 
Therefore, provided that coefficients a_{c} and a_{t} are not zero:

(49)  

(50) 
By zeroing derivatives of σ_{c} and σ_{t} (equations (43) and (44), respectively) with respect to, respectively, crushing and cracking strain, maximum values f_{cm} and f_{tm} (Figure 27) are obtained:

(51) 
Equations (51) provide:

(52)  

(53) 
These developments complete the description of the proposed methodology. Parameters a_{c}, a_{t}, b_{c} and b_{t} can be determined from equations (52), (53) and (46) in terms of f_{cm}, f_{c0}, f_{tm}, f_{t0}, l_{eq}, G_{ch} and G_{F}. 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.
This subsection describes the concrete uniaxial stressstrain 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 damageplasticity 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 stressstrain 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 meshindependency; the regularization approach is based on selecting the softening branches of concrete constitutive laws depending on mesh size. In Figure 27, f_{cm} and f_{tm} represent compressive and tensile stress strength, respectively; corresponding strains are ε_{cm} and ε_{tm}, respectively. According to [90], it is assumed that ε_{cm} = 0.0022, f_{cm} = f_{ck} + 8 (f_{ck} is the characteristic value of concrete compressive strength), and ; these stresses and the deformation modulus are expressed in MPa. In Figure 27.a, and are the crushing and elastic undamaged components of strain; and are the plastic and elastic damaged components. In Figure 27.b, and are the cracking and elastic undamaged strain components; and 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, , reaching 0.4 f_{cm}; second (ascending) segment (in between 0.4 f_{cm} and f_{cm}) is quadratic [90]:

(54) 
In equation (54), E_{ci} is the modulus of deformation of concrete for zero stress, given by and E_{0} = (0.8 + 0.2 f_{cm} / 88) E_{ci} (in MPa) [90]. In the initial linear branch, E_{0} is the secant modulus that corresponds to 0.4 f_{cm} stress.
Third (descending) segment is given by:

(55)  

(56) 
In equations (55) and (56), G_{ch} is the crushing energy per unit area [28], and l_{eq} 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 stressstrain 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 stressstrain law multiplied by the characteristic length.
Regarding tensile behavior, the ratio between tensile stress σ_{t}(w) (for crack width w) and maximum tensile strength f_{tm}, is given by [96]:

(57) 
In equation (57), c_{1}= 3, c_{2}= 6.93 [96], and w_{c} is the critical crack opening. Equation (57) shows that σ_{t}(0) = f_{tm} and σ_{t}(w_{c}) = 0. Therefore, w_{c} can be considered as the fracture crack opening. Equation (58) [96] relates w_{c} with the tensile strength and fracture energy G_{F} per unit area:

(58) 
According to [90], G_{F} (N/mm) can be calculated as

(59) 
In equation (59), f_{cm} 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]:

(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 globalpurpose simulation. After this assumption, in the descending segment of the tensile stressstrain curve (Figure 27.b), the strain can be obtained in terms of the crack opening from the following kinematic relation:

(61) 
This subsection describes the concrete uniaxial concrete behavior; next subsection describes their implementation in the proposed approach.
Following the formulation described in the previous subsection, coefficients a_{c}, a_{t}, b_{c} and b_{t} (subsection 3.2) can be determined. Coefficient a_{c} is obtained by replacing f_{c0} = 0.4 f_{cm} (Figure 27.a) in equation (52). Figure 27.b shows that f_{tm} = f_{t0}, then equation (53) shows that a_{t} = 1. Coefficient b_{c} can be determined from equation (46) by replacing f_{c0} = 0.4 f_{cm} = 0.4 (f_{ck} + 8) MPa. Coefficient b_{t} is determined by replacing a_{t} = 1 and [90] in equation (46). Therefore:

(62) 
In the last two expressions in equation (62), f_{ck} is expressed in MPa. After parameters a_{c}, a_{t}, b_{c} and b_{t}, damage variables d_{c} and d_{t} 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 f_{ck}, the parameters in Table 3, the mesh size , and the ratio b (eqn. (56)). Initial assumption is b = 0.9
2. Calculate the compressive / tensile stress strength f_{cm} = f_{ck} + 8 /
3. State the strain at compressive stress strength as ε_{cm} = 0.0022
4. Calculate the initial tangent modulus of deformation of concrete and the undamaged modulus of deformation
5. Calculate the fracture / crushing energy (N/mm) /
6. Calculate the critical crack opening
7. Build the first / second / third segments of the concrete uniaxial compressive law: / eqn. (54) / eqn. (55). In eqn. (55), strain is bounded; the selected upper bound should fulfill the condition that the crushing energy G_{ch} (equation (60)) is reached
8. Build the first / second segment of the concrete uniaxial tensile law: / eqns. (57) and (61)
9. Calculate the damage parameters according equation (62): ; ; ;
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: ,
12. Calculate the average value of ratio 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 (codetype) 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 
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 socalled mesh regularization techniques. One of the simplest remedy is the Crack Band Method; it uses energybased scaling of the softening part of the stressstrain 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 socalled NonLocal approaches; they are based on introducing nonlocality in the constitutive model. This nonlocality 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 nonlocality is including higherorder 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 energybased 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 energybased 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 8node hexahedron solid finite elements (C3D8R). Figure 29 displays three sketches of the cube discretized with a coarse mesh (l_{eq} = 200 mm, one element), a medium mesh (l_{eq} = 50 mm, 64 elements) and a fine mesh (l_{eq} = 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 
Mesh  l_{eq}  f_{ck}  f_{cm}  f_{tm}  G_{ch}  G_{F}  b  a_{c}  a_{t}  b_{c}  b_{t} 
(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.
Given that the problem under consideration is extremely simple, a closedform 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 w_{c} obtained from equation (58) (w_{c} = 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 closedform solution. Figure 32 shows a satisfactory agreement among the plots for the three mesh sizes; this observation corroborates the meshinsensitiveness 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 closedform 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. Forcedeflection 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 closedform 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 (G_{F}) 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.
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 subsection3.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.
f_{ck}
(MPa) 
f_{cm}
(MPa) 
f_{tm}
(MPa) 
G_{ch}
(N/mm) 
G_{F}
(N/mm) 
b  a_{c}  a_{t}  l_{eq}
(mm) 
b_{c}  b_{t} 
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 stressstrain 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) Stressstrain  
(b) Damagestrain  
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.
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 rateindependent elastoplasticity model with isotropic hardening has been used for the simulation of monatomic loading, in the other hand a combined isotropickinematic 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:

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

(64) 
is the elastic stiffness tensor, and 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 pressureindependent yield function is defined by the function

(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:

(66) 
where is the equivalent plastic strain rate.
The nonlinear isotropickinematic 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:

(67) 
where are material parameters that can be calibrated form cycle test data, is the equivalent plastic strain rate. When 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

(68) 
where 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 ( is specified after shifting the strain axes to .
Figure 35. Stressstrain data for a stabilized cycle [89] 
The shifted data are entered in tabular form in ABAQUS [89] as stressstrain, 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 linearkinematic 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".
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 :
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 bondslip 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 bondstrength 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 pullout failure. Figure 38 shows a representation of the damage evolution for pullput failure type under monotonic loading.
Figure 38. Monotonic bondslip behavior for pullout failure [116] 
Bondslip 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 bondslip behavior for pullout 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 bondslip 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. Bondslip is represented through relation between bond stress and bar slip implemented in zerothickness interface element. [120] developed a fournode zerothickness bondslip element to be used for twodimensional 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 bondslip effects are either incorporated in the formulation of the element or taken into consideration through zerolength springs. [121] introduced the beamcolumn model with bondslip proposed by [122] into the forcebased fibersection element developed by [123]. In regards to zerolength approach; [124] proposed a law to relate bar stress and slip at end of anchorage in footingcolumn or beamcolumn connections; this law has been calibrated with experimental results and used as a constitutive relation for steel fibers in a zerolength fibersection element to simulate end rotation of RC columns. These models have low computational cost and, therefore, are suitable for simulation of fullscale structures; conversely, they do not study explicitly the bond behavior, merely reflect the additional flexibility provided by bondslip.
Rib scale models are too computationally expensive for the practical simulation of fullscale actual structures, and element scale models can be an attractive option for largescale 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:
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 bondslip law. [90] provides detailed monotonic bondslip law for pullout and splitting failure as a function of bond stress τ_{2} and the relative displacement S_{2} as shown in Figure 41. Once the bond stressslip relation is defined, it can be converted to forceslip law and attached to each spring as its constitutive law.
(a)  (b) 
Figure 40. Implementation of Spring model for bondslip [129] 
In Figure 41; τ_{f } is the friction bond resistance, S_{R }is the distance between the bars ribs.
Bondslip law is converted to force displacement relation and assigned as a constitutive law for the spring elements.
Figure 41. Monotonic bonslip 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).

(69) 
F_{bs} is the equitant bond force in the spring, D_{b} 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 bondslip in subsection 5.1, the cyclic behavior of bond stressslip 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.
The 3D 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 2D 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 stressslip constitutive law uses the phenomenological law proposed by [130] which is based on concepts originally developed by [116]. This law was extensively verified by pullout tests on bars embedded in well confined concrete.
Figure 42 displays the three components of interaction stress (σ_{1}, τ_{2}, τ_{3}) and relative displacement (S_{1}, S_{2}, S_{3}) during sliding.
Figure 42. Stress and relative displacement components at barconcrete interface [129] 
The bond resistance τ_{2} in the model proposed by [129] is decomposed into bearing and friction components:

(70) 
In equation (70), τ_{b} and τ_{f} are the full bearing and friction resistance of an elastic bar under monotonic pullout action. Such resistances are multiplied by reductions factors. ρ_{b,s} and ρ_{f,s}account 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 S_{1}.
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 S_{peak }and the clear spacing between the bar ribs S_{R}.
Figure 43. Bond stressversusslip 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:

(71)  

(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:

(73)  

(74) 
In equations (73) and (74), and are the maximum absolute values of slip in positive and negative directions, respectively. is the cumulated slip after slip exceeds S_{peak} for the first time. is equal to:

(75) 
In this work, the normal displacement S_{1} 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 stressdisplacement 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.
Next subsection presents a new modeling scheme to implement the bondslip model described in the previous subsection "developed by [129] and [133]" in a FEM analysis.
The bondslip 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 nonparallel 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 bondslip interface model 
Figure 44.b shows that the proposed technique is based on assigning, for each interface element, a pair of tieconstraint 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 tieconstraint (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.
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 nonparallel longitudinal reinforcement bars, such as beamcolumn and slabcolumn joints, structural walls, foundations, and connections among these elements.
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. Bondslip model described in subsection 5.3 and its implementation scheme proposed in subsection 5.4 are implemented in the simulation of cyclic loading tests.
Bondslip 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 damageplastic behavior and bondslip of reinforcement for practical 3D continuum simulation of RC structures subjected to seismic loads:
New smeared strategy to simulate the reclosing of cracks after tensioncompression reversals is proposed in this work. As discussed in subsection 2.4 the reclosing of cracks after tensioncompression reversals is governed by the parameter h_{c}, 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 bondslip 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 openingclosing 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 h_{c} near the member ends. Figure 46 displays a frame modelled with different values of h_{c}; in the darker shadowed zones, little or no stiffness is expected upon cracks reclosing.
Figure 46. Zones of framed structures with different cracks openingclosing 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 bondslip is modelled.
This experiment [137] is a quasistatic test consisting of imposing a displacement law to a laboratory singlespan, singlestory 2D 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.
f_{ck}
(MPa) 
f_{cm}
(MPa) 
f_{tm}
(MPa) 
G_{ch}
(N/mm) 
G_{F}
(N/mm) 
b  a_{c}  a_{t}  l_{eq}
(mm) 
b_{c}  b_{t} 
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 2node truss elements (T3D2) and the left part describes concrete discretization with 3D 8node 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.
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).
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 topright edge of the beam, subsequently both column at the bases, ending by yielding at topright, 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/m^{2}) in reinforcing bars is demonstrated in Figure 53.d
(a)  (b) 
(c)  (d) 
Figure 53. Final damaged state of [Pires 1990] RC frame experiment 
This experiment [143] is a quasistatic test consisting of pushing monotonically until failure a laboratory, singlespan, twostory, 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 stressstrain 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). D_{b} is the bar diameter, f_{y} / f_{u} are the yield point / ultimate stress, E_{s} is the steel modulus of elasticity, E_{sh} is the slope of the hardening branch, ε_{sh} is the strain that corresponds to the onset of hardening, and ε_{u} is the ultimate strain.
Bar No.  D_{b}
(mm) 
f_{y}
(MPa) 
f_{u}
(MPa) 
E_{s}
(GPa) 
E_{sh}
(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 nonexistent; 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) Forcedisplacement 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 shearrelated 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 adhoc 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 2node truss elements (T3D2) and the left side describes concrete discretization with 3D 8node 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 
f_{ck}
(MPa) 
f_{cm}
(MPa) 
f_{tm}
(MPa) 
G_{ch}
(N/mm) 
G_{F}
(N/mm) 
b  a_{c}  a_{t}  l_{eq}
(mm) 
b_{c}  b_{t} 
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 vs. (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 and do not approach zero when damage variables d_{c} and d_{t} 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.
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 1^{st} story beam. Figure 59.a and Figure 59.b refer to the first story right beamcolumn 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 beamcolumn 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 beamcolumn 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 
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/m^{2}) bars is demonstrated in Figure 62.d
(a)  (b) 
(c)  (d) 
Figure 62. Final damaged state of [Vecchio and Emara 1990] frame experiment 
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}/A_{g} 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. D_{b} is bar diameter, f_{y} / f_{u} are yield point / ultimate stress, E_{s} is modulus of elasticity, ε_{sh} is the strain corresponding to onset of hardening, and ε_{u} is ultimate strain.
Bar Id.  D_{b} (mm)  f_{y} (MPa)  f_{u} (MPa)  E_{s} (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 3D 8node hexahedron solid elements (C3D8R). Figure 65.b and Figure 65.c represent the steel and interface elements; steel is discretized with 2node truss elements (T3D2). Bondslip 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 vs. , where refers to plastic compressive strain. Figure 66.b displays analogous plots for the tensile behavior ( vs. ).
(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 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.
f_{ck}
(MPa) 
f_{cm}
(MPa) 
f_{tm}
(MPa) 
G_{ch}
(N/mm) 
G_{F}
(N/mm) 
b  a_{c}  a_{t}  l_{eq}
(mm) 
b_{c}  b_{t} 
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 h_{c} = 0.9 in all the body, and generating a discrete crack at the footingcolumn 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 bondslip 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 forcedisplacement 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).
Figure 68 highlights the relevance of bondslip. 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.
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 stressvertical displacement behavior of the bottomleft corner longitudinal bars presented in Figure 65.c is shown in Figure 71. The average stressdisplacement 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 bondslip elemnts in comparision with perfect bond model on the hesteritic behavior of stressdisplacment in longitudinal bars. The higher negative displacement in the hysteretic behavior of the model with bondelement shown in Figure 65, can be understood as a result of penetration due to bond degradation.
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 h_{c} (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 h_{c} 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)
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 bondslip 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 bondslip elements in the foundation. Table 6 displays the calculated values for the parameters of the damage variables methodology.
f_{ck}
(MPa) 
f_{cm}
(MPa) 
f_{tm}
(MPa) 
G_{ch}
(N/mm) 
G_{F}
(N/mm) 
b  a_{c}  a_{t}  l_{eq}
(mm) 
b_{c}  b_{t} 
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 bondslip 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 "bondslip element" attached to longitudinal bars, Figure 76.b depicts a more detailed view of the right beamcolumn joint steel bars, and Figure 76.c shows the same of Figure 76.b with bondslip 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, bondslip elements are assigned to the longitudinal bars along the column foundation, inside beamcolumn joints, at the column top and bottom ends, and at beamends. 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 bondslip 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.
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.
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.
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 h_{c} (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).
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 
Figure 85. Cyclic damage evolution of [Pires 1990] frame experiment, proposed formulation with reduced ℎ_{c} (subsection 6.1), beam column joint 
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 meshinsensitive, 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. Meshinsensitivity 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 meshinsensitiveness 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 stressstrain 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 ExperimentsMonotonic loading: Two RC frames have been simulated, results have proved the ability of the proposed methodology to reproduce accurately the experimental forcedisplacement 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 3D interface bondslip 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 ExperimentsCyclic 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 bondslip in the hysteretic behavior and the energy dissipation capacity.
The proposed integrated model can constitute a practical tool for accurate simulation of highly damaged mediumsize RC structures, both in research and in conventional analyses. It can be also used to calibrate more simplified models.
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 mechanicsbased 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 componentbased models.
It can be distinguished between two type of structural componentbased 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 elastoplastic 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 componentbased models in earthquake engineering, a numerical study is conducted on the relation among the nonsimulated deterioration modes of the elements of nonductile RC frames and their final capacity.
The main contribution of this part is developing an advanced structural componentbased 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 nonductile building and the prototype building.
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 momentcurvature 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 twocomponent model, a bilinear elasticstrain hardening momentcurvature 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 elasticperfectly 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 momentcurvature relation, where EI is the preyield section stiffness. The elastic modulus of the elastoplastic 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 onecomponent 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 momentrotation relation of RC members.
The concentrated plasticity concept was generalized to incorporate axialmoment 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 macroelement capable of accounting for inelastic lateraltorsional coupling in asymmetric buildings. This model was developed considering a perfect elastoplastic interaction between the storyshears and torque [155]. More recently, other researchers have included effects such as hardening into bidimensional hinges by introducing a hardening rule to the yield surface [156].
Twocomponent model [151] 
(a) Onecomponent 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 onecomponent model in its original form describes the nonlinear behavior of the girder. The hysteretic momentrotation 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 momentrotation or forcedeformation 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 bondslip 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, peakoriented, and pinching. The modified models include most of the sources of deterioration: i.e. various modes of cyclic deterioration and softening of the postyielding stiffness, and also account for a residual strength after deterioration. The models incorporate an energybased deterioration parameter that controls four cyclic deterioration modes: basic strength, postcapping 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 loadmoment (PM) 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 postyield and degrading response. On the other hand, some hinge elements with detailed momentrotation hysteresis models may not capture PM interaction, except to the extent that the momentrotation response is defined based on average values of axial load and shear that are assumed to be present in the hinge.
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 beamcolumn interface into the member as a function of loading history. The rest of the beam remains elastic. The fixedend rotations at the beamcolumn 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 bondslip 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 beamcolumn element. The first displacementbased 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 beamcolumn element represented in Figure 90 is writing as follow:

(76) 
Where is the transverse and longitudinal displacement along the longitudinal axes (x), is the matrix of the interpolation functions and is the vector of the nodal displacements presented in Figure 90. Equation (76) can be written as :

(77) 
where

(78) 
Figure 90. Nodal displacement of beamcolumn element in 2D [163] 
The generalized deformations of the problem are the axial strain ε(x) and the curvature about the zaxis φ_{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:

(79) 
Where , 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 :

(80) 
The section stiffness matrix relates the section forces with the corresponding deformations:

(81) 
Where is the vector of the section forces "the generalized stresses"; the axial force N(x) and the bending moment M_{z}(x) at section x.
The main shortcoming of this classical finite element displacementbased 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) MomentCurvature 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 beamcolumn 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. Beamcolumn element in forcebased 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:

(82) 
is the matrix containing the force interpolation functions, is the vector of the generalized force demonstrated in
Figure92.b. Equation (92) can be written as follow:

(83) 
The application of the virtual force principle yields the element flexibility

(84) 
Where is the section flexibility, equation (81) can be rewritten as:

(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 stressstrain 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. DisplacementBased and ForceBased 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 forcedeformation 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 forcedeformation 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 displacementbased 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 forcebased element to simulate the material nonlinear response of a frame member, compared with several displacementbased 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 bondslip 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 momentrotations and hysteresis curves as described in the previous subsection.
This section presents a comparison between a models derived from structural componentbased 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 continuumbased 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 SAP2000 v.16 [171] by using the displacement formulation [172], without accounting for sheargenerated deformation. Secondorder effects are taken into account. Zerolength flexural and axialflexural 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 axialflexural 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 afteryield behavior is assumed to be totally brittle. Neither shearflexural interaction nor shearaxialflexural 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 2node 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. Secondorder 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 stressstrain 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 continuumbase 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 continuumbased model to describe the response of the RC frame, results of the structural componentbased 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 momentrotation 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 continuumbased 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 continuumbased 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 secondorder 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 continuumbased models are demonstrated in Table 12.
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 componentbased models. Results are compared with the real behavior and the one previously captured by the proposed continuum mechanicsbased 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 bondslip 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 bondslip effects are expected to occur. The strain penetration model [124] is the most common strategy to capture bondslip effect in structural componentbased 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 stressslip law for steel bars based on pullout tests is generated; it is used for zerolength elements situated at column bottom and is integrated into a global fiber model. Figure 100.a shows the hysteretic stressslip 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].
The RC column experiment Tanaka 1990 [98] described in subsection 6.3.1 is discretized with one 2node finite elements using the forcebased formulation, and one zerolength element to represent bondslip. 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. Stressstrain 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 stressstrain relation 
The following formulations are used to calculate Concrete01 model parameters (MPa)

(86) 
In equations (86), f^{'}_{c} and f_{yh} 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, A_{s}´´ is the crosssectional 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}/E_{c} where E_{c} 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 :

(87) 
ρ_{sh} and f_{yh} are the transverse reinforcement ratio and yield stress respectively A_{sh} 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 zerolength 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 R_{c } shown in Figure 100.a. This parameter modifies completely the response of the zerolength element; consequently, this will alter significantly column behavior. R_{c} = 1 represents perfect bond situation and, thus, no pinching effect; lower values of R_{c} correspond to significant pinching. For the simulation of Tanaka column, R_{c} = 0.6 has provided the best results.
Figure 103 displays comparisons between experimental and numerical hysteresis loops obtained with the abovementioned structural componentbased model.
Figure 103. Hysteretic force displacement response of [Tanaka 1990] column. Experimental vs. structural componentbased model 
Plots from Figure 103 show that the structural componentbased 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 continuumbased 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 continuumbased 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 componentbased model. Herein the zerolength element is assigned in three places where bondslip is expected to be significant: columns base, columns top, and beam ends.
As discussed previously, this model is particularly sensitive to pinching factor R_{c}; As a preliminary attempt, the calibrated value from the simulation of the bridge column R_{c} = 0.6 is used in the three locations. Analogously to Figure 103 comparison between experimental and numerical hysteresis loops obtained by the structural componentbased is demonstrated in Figure 105.
Figure 105. Hysteretic force displacement response of [Pires 1990] frame. Experimental vs. structural componentbased model (constant R_{c}) 
Plots from Figure 105 show that the response obtained by the structural componentbased model is shifted from the experimental response, this can be explained by the high flexibility generated by the low values of R_{c} in the three locations.
To stabilize the numerical response, the value of R_{c} 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 R_{c }values.
Figure 106. Hysteretic force displacement response of [Pires 1990] frame. Experimental vs. structural componentbased model (varying R_{c}) 
Plots from Figure 106 show that the structural componentbased 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 continuumbased and the simplified one with the abovementioned values of R_{c}. Figure 107 collects the plots from Figure 77.b and Figure 106. Results from Figure 107show that the structural componentbased model has better capacity in representing the pinching effect, maximum capacity and strength reduction are better represented by the continuumbased model.
Figure 107. Hysteretic force displacement response of [Pires 1990] frame. Experimental vs. structural componentbased (varying Rc) and advanced model 