$$\frac{\partial \Pi(u)}{\partial u\_i} = \int\_{\Gamma} \psi\_j(\mathbf{x}) \phi\_i(\mathbf{x}) p\_j \, \mathrm{d}x - \int\_{\Gamma\_N} \psi\_j(\mathbf{x}) \phi\_i(\mathbf{x}) p\_{Nj} \, \mathrm{d}x. \tag{46}$$

The integration in (46) does not need to be extended all over the boundary Ŵ since ψj and φi are local vector shape-functions whose support is given only by the boundary elements associated to the DOFs j and i, respectively. Taking also into account that, for the minimization problem under study, the relevant DOF are associated to ŴC and their support is a part of ŴC, supp φi ⊂ ŴC, Equation (46) can be further simplified as

$$\frac{\partial \Pi(u)}{\partial u\_i} = \int\_{\Gamma} \psi\_i(\mathbf{x}) \phi\_i(\mathbf{x}) p\_i \, \mathrm{d}x. \tag{47}$$

This expression gives the relevant numerically computed gradient of the potential energy in the present problem in terms of the numerically computed boundary tractions on ŴC.

### 9. A CRITICAL COMPARISON OF ALGORITHMS M1 AND M2

In **Tables 1**, **2** of section 7, two variants of the conjugate gradient method considered, M1 and M2, respectively, were described by pseudocodes, where M2 is the original version introduced by Dostál (2009). The main difference between these two variants, is that in M1 we do not explicitly calculate the matrix **A**, since instead of using the vector-matrix-vector product 1 2 **x** T**Ax** and the matrix-vector product **Ax**, we use the quadratic form A (**x**) and the derivative A ′ (**x**). In reference to M2, we initially need to compute the matrix **A**, as it is described in section 5 (actually it is the Step 0 in **Table 2**), which might be time consuming, but then the rest of the procedure, that is actually the conjugate gradient loop, showed to be fast enough. On the contrary, for M1 the initialization (Step 0) is direct and fast, however each time we need to compute the value of the quadratic form A (**x**) and its derivative A ′ (**x**), that is in each k th step of the conjugate gradient loop (Step 1 in **Table 1**), a BEM solution is required. Nevertheless, the influence BEM matrices need to be computed only once, for all the iterations and time steps. Additionally, having decomposed the matrix of the linear algebraic system of equations produced by BEM (e.g., by LU- or QR-decomposition), each time that the solution of this system is required, only a back-substitution takes place.

Details of the algorithms for constrained minimization problems can be found in Dostál (2009). In this work, for both variants of the conjugate algorithms, we employed Polyak's algorithm for bound constrained minimization (Algorithm 5.2 in Dostál, 2009), for which the number of iterations is bounded by N˜ = n 2 2 n , while for the standard conjugate gradient algorithm by N = n, with n being the number DOFs, i.e., the order of the matrix **A**. However, according to Dostál (2009), this bound for Polyak's algorithm may be considered too pessimistic.

It was found, in this work, that in M2 the usage of the symmetric matrix described in section 5, leads to exactly the same results for the energy calculation as in M1, i.e., A (**x**) = 1 2 **x** T**Ax**. However, the accuracy of the results for displacements and tractions was worse in M2 using the symmetric matrix than in M1, especially for points near a change of boundary condition type (that is from ŴC to ŴN). Previously, similar observations were made in BEM by de Paula and Telles (1989) and Tullberg and Bolteus (1982), when the symmetric stiffness matrices was chosen instead of the non-symmetric one. We conclude that the source of this lose of accuracy in our context is because the numerically computed derivatives (gradient) of the quadratic form used in M1 and M2 may be quite different, i.e., A ′ (**x**) 6= **Ax**, when the symmetric matrix **A** is used in M2 and if in M1, instead of computing explicitly the gradient of the quadratic form, we use its expression in terms of the boundary tractions according to section 8. In order to overcome this difficulty in M2 application we proceed to construct a consistent non-symmetric matrix **A**˜ for its usage in M2.

In Antes and Panagiotopoulos (1992), it is claimed that the stiffness matrix constructed by BEM should be symmetric due to Betti's theorem, and not because of the numerical method used to compute the unit displacement response. This might be equivalent with our consideration here, that the reciprocity statements in the context of a numerical implementation of BEM are only approximately valid. In order to construct the nonsymmetric matrix **A**˜ , we solve by BEM a problem, for each ith DOF of displacements on ŴC, such that each time this is the only non zero component and equal to the unit value. The ith column of the **A**˜ matrix consists of the derivatives ∂A ∂uj , with uj being a DOF of displacements on ŴC.

It can be shown, and also it was numerically verified, that the matrix **A** of section 5, is the symmetric part of matrix **A**˜ constructed here, i.e., **A** = **A**˜ T+**A**˜ 2 .

# 10. PROPERTIES OF THE SYMMETRIC AND NON-SYMMETRIC BEM MATRICES: A AND A˜

As mentioned in **Table 2**, the matrix **A** is assumed to be symmetric and positive definite. However, it is well known that it is possible to apply the conjugate gradient method of Hestenes and Stiefel (1952) also to non-symmetric systems after a minor modification of the method (e.g., solving **A**˜ T**A**˜ ξ = **A**˜ T**b** instead of **A**˜ ξ = **b**). The well established and reported fact of a convergence property reduction in the conjugate gradient method for nonsymmetric, Saad and Schultz (1985), and general inconsistent systems, Axelsson (1980), did not have any significant influence on the numerical examples studied here. The linear system matrix computed by the collocation BEM, using the procedure described in section 5, may not be positive definite even if it is symmetric. Actually, this is the case when the Dirichlet and potential contact boundary parts, ŴD and ŴC, are empty, whereas a non-empty Dirichlet part guarantees also the positive definiteness.

In order to numerically explore properties of the eigensystem of the BEM system matrix, we consider a circular elastic body of radius r = 5.0, whose Young's modulus is E = 1000.0 and Poisson's ratio ν = 0. Assuming the entire boundary of the body as a potential contact zone, i.e., ŴC = Ŵ, we construct the matrix **A** as described in section 5, and compute its eigenvalues and eigenvectors. It is observed that three small, possibly negative, eigenvalues exist which approximate the theoretically zero eigenvalues associated to the three rigid body motions of the circular domain. Two negative eigenvalues have almost equal value of −6.116E-6, and correspond to translations, whereas the third one is several orders of magnitude smaller, taking the value of 7.395E-9 and corresponds to a rotation. In **Figure 6** deformed shapes given by some eigenvectors of the matrix **A** are plotted, for the above defined circular body. Note that for n nodes we have 2n kinematic DOFs and eigenvectors.

In **Figure 7** (left) it is shown how the theoretically zero eigenvalues associated to rigid body translations are approximated by the two equal and small (possible negative) eigenvalues of the symmetric and non-symmetric BEM matrices **A** and **A**˜ , respectively, for a refinement of the boundary element mesh. The absolute values of the numerically computed eigenvalues are plotted in logarithmic scale. These eigenvalues are negative for the symmetric matrix **A**. Recall that the number of nodes is one half of the number of DOFs. We may observe that these two small eigenvalues are numerically zero regardless discretization in the case of the non-symmetric matrix **A**˜ . Thus, the symmetric matrix **A** describes the rigid body motions only approximately (an approximation which is improved with a mesh refinement), while in the case of the non-symmetric matrix **A**˜ these rigid body motions are well represented. Finally, in **Figure 7** (right) the largest and smallest positive eigenvalues for a mesh refinement are plotted, for both matrices **A** and **A**˜ , a very good agreement being observed.

In order to numerically study the influence of the application of **A** instead of **A**˜ in BEM calculations, we present the following trivial example whose solution is a rigid body translation. A quadrangular domain is considered with zero tractions on the right and left sides, thus they belong to ŴN. On the upper side, which belongs to ŴD, positive normal displacements (in the outward normal direction) are prescribed together with zero horizontal displacements. The bottom side is in frictionless contact with an obstacle, thus it belongs to ŴC. This problem is quite trivial, since the solution is a rigid body translation in the vertical direction with zero tractions along the whole boundary of the domain. Recall, that in the present computational implementation, both displacements and tractions along the bottom side are considered as unknowns. **Figure 8** (left) shows that numerically zero horizontal displacements are obtained, and similarly **Figure 8** (right) shows that the vertical displacements accurately represent the rigid body translation by applying the method M1, and also by applying the method M2 where however the non-symmetric matrix **A**˜ is used. On the other hand, when the symmetric matrix **A** is used in the method M2, both results for horizontal and vertical displacements are not so accurate although they improve with the mesh refinement. The main difference from the exact solution is observed near the corners where ŴC intersects with ŴN.

# 11. NUMERICAL EXAMPLES

A numerical implementation of EC-BEM for both cases of contact was accomplished using an in-house BEM software (Panagiotopoulos, 2017) developed, mainly by the first author, in Java programming language. Linear continuous elements with two nodes are used for the analysis of several examples assuming plane strain conditions. The results obtained by this code are compared with the results obtained by 2D BEM codes developed previously by Blázquez et al. (1998a,b, 2006), and Graciani et al. (2005) in which standard techniques of BEM for contact problems essentially based on the so-called displacement and load scaling techniques (recently also referred to as Sequential Linear Analysis, SLA) (París and Blázquez, 1994) are deployed, as well as by the commercial FEM code ANSYS.

The energetic implementation of the first case of contact described, the Signorini contact, is tested in section 11.1 by solving two problems: one of them with receding contact, in section 11.1.1, and another with advancing contact in section 11.1.2. The case of adhesive unilateral contact is tested in section 11.2 with a similar problem to that used as example of receding contact to test the numerical solution of the Signorini contact problem.

# 11.1. Numerical Examples to Test the Signorini Contact Implementation

#### 11.1.1. A Case of Receding Signorini Contact

The first problem studied is shown in **Figure 9**. The geometry is composed by two rectangular solids 1 and 2 in frictionless

FIGURE 7 | Left: Absolute values of the two equal and small eigenvalues (associated to rigid body translations) of the BEM matrices A and A˜ for a mesh refinement. Right: The largest and smallest positive eigenvalues of A and A˜ for a mesh refinement.

contact. The shorter rectangle 1, located below, has dimensions (L7 + L8 + L9) × h1, see **Table 3** for lenght values. Vertical and horizontal displacements are restrained, respectively, along a horizontal segment of length L1 at the lower edge and the lower-left corner of the rectangle. A larger rectangle 2 with dimensions P6 i=1 (Li) × h2 is located at the top. Its vertical and horizontal displacements are restrained, respectively, along an horizontal segment at the lower edge and the upper-right corner. In addition, a uniform vertical displacement uD is imposed along a segment with dimension L3 at the upper edge of this rectangle. A segment of length L8 along the common edge of both solids is defined as potential contact zone, while L7 and L9 are assumed to be traction free zones, since separation of the solids is expected, thus, no interpenetration of solids happens along L7 and L9. The material of both solids is assumed to be linear elastic and isotropic with properties of aluminum shown in **Table 3**.

In the case of the BEM models presented below, the mesh used is set by the number of elements ni with i = 1, 2, ..., 9

corresponding to the segment of length Li and nh1 and nh2 corresponding to the segments with lengths h1 and h2 respectively. Element length is constant along each segment. In the case of the FEM model, solids are meshed in such a way that the nodes along their edges coincide with those in BEM models. Values for the different meshes used here are presented in **Table 4**.

TABLE 4 | Values for the variables ni .

Along the potential contact zone defined between both solids, two distinct types of contact are considered: (a) unilateral frictionless Signorini contact analyzed in this section, as well as (b) adhesive unilateral frictionless contact shown in section 11.2.

Firstly, computation is carried out taking the model A with the mesh 1 (see **Tables 3**, **4** respectively). The boundary of the domain 1 is discretized by 340 elements, and that of 2 by 580 elements. The EC-BEM computation works with 160 interface elements along the segment denoted as L8, requiring the optimization of a quadratic problem of 322 DOFs. Note that n8 in **Table 4** is the same for both solids, so the mesh along the potential contact zone is conforming.

The model is solved by applying the above algorithm. **Figure 10** shows the deformed shape, with the scale factor of 500, predicted by the EC-BEM. Note that the potential contact zone is divided into three regions: two segments at both extremes of the initial contact zone where both solids are separated and an active contact zone in the central segment.

The numerical results along the potential contact zone are of most interest. **Figure 11** shows the normal tractions and relative normal displacements along this boundary part. The central region has compressive tractions and null opening, as expected. On the contrary, normal tractions are almost zero along the regions with a positive relative normal displacement. It is interesting to recall that, in any iteration, during the minimization in the EC-BEM procedure the initial potential contact zone is assumed as a Dirichlet boundary part. According to this, the simultaneity between compressive tractions and null opening or null tractions and positive opening is not explicitly imposed by the algorithm used, as in most contact algorithms. However, the results obtained show that these conditions are well predicted by the energy optimization procedure.

Results from the EC-BEM algorithm described here are compared in **Figure 12** with the results computed for the same problem with a classic point method algorithm implemented in BEM and the code ANSYS using an Augmented Lagrange algorithm. Note that the agreement between the predictions is accurate in spite of the strong qualitative difference between the different algorithms. The most significant difference is the presence of fictitious oscillations predicted by the EC-BEM at the extremes of contact zone, which do not appear in other numerical results. Also in **Figure 12**, a second solution obtained by the EC-BEM algorithm and a finer mesh (mesh 2 in **Table 4**), where additionally a reduced potential contact zone (model B) was assumed, is shown. As can be seen there, the oscillations are reduced, together with the maximum values of the fictitious tensional tractions.

Previous results were obtained using conforming meshes along the potential contact zone. Results obtained with nonconforming meshes are presented in the following. However, strategies presented here for non-conforming meshes are very simple, more advanced formulations may be combined with the present framework for contact problems using energetic principles.

In order to evaluate the influence of a non-conforming mesh on the results obtained by the present algorithm, a computation with the model B and the mesh 3 (see **Tables 3**, **4** respectively) is carried out. Mesh 3 as well as mesh 4, which have a reduced potential contact zone given by L8, are defined in view of results of the previous analysis. The number of elements along the potential contact zone of 2 (top) is twice that of 1 (bottom).

As can be seen in **Figure 13**-top, the normal tractions along the contact boundary of 2, which is the one with finer mesh, present some oscillation around the smooth solution of tractions of 1. This behavior does not depend on the choice of main and secondary domain, a distinction which actually is just auxiliary. Note that, this oscillatory behavior reduces with mesh refinement. In order to observe the reduction of these oscillation, a finer mesh, at the contact boundaries, is employed (mesh 4 in **Table 4**), maintaining the ratio of number of elements at both sides of the contact zone. **Figure 13**-bottom shows the comparison of tractions along the boundary of 2 for the coarser and finer mesh. Note that these oscillations are slightly reduced.

#### 11.1.2. A Case of Advancing Signorini Contact

The indentation of a half cylinder against an elastic foundation is considered. This is an advancing contact problem as the length of

FIGURE 12 | Comparison of tractions computed along the contact zone by the algorithm described here (EC-BEM) and the point method with BEM (PM-BEM) and the FEM code ANSYS using an Augmented Lagrange algorithm. The EC-BEM solution corresponds to model A with Mesh 1, while a finer one also included here, for model B with Mesh 2 (see Tables 3, 4, respectively).

the contact zone depends on the value of the prescribed loading. The problem geometry, shown in **Figure 14**, is composed by a rectangle with dimensions L1 × h1 and a half circle with radius

h2. The potential contact zone is defined by the length L2 and the angle φ for the rectangle and the circle respectively. Vertical uD and null horizontal displacements are imposed along the whole straight edge of the half circle.

Linear elastic and isotropic behavior is assumed for both solids, values of Young's moduli E and Poisson's ratios ν being given in **Table 5**, where also the geometrical dimensions can be found. The number of elements along each boundary part are detailed in **Table 6**. Frictionless Signorini contact

TABLE 5 | Geometry, elastic properties and imposed displacements.

TABLE 6 | Values for the variables ni defining the different meshes.

is considered on the potential contact zone. Due to the fact that this is an advancing contact problem, we are interested in tracking the evolution of the load, which is a non-linear function of the prescribed displacements. In the numerical solution the vertical displacement is introduced by increasing 10 times this displacement in a monotonic way leading to the value uD shown in **Table 5** at the step 10.

The problem is solved using the EC-BEM code described above with the mesh 1 (see **Table 6**), where n1, n2, nh1 and nφ correspond to the number of elements on the boundary part defined by L1, L2, h1 and φ respectively. n90◦−φ is the number of

elements along the arc of the half cylinder which is not a potential contact zone and nL1−L2 is the number of elements at the two parts of the upper edge of the rectangle which are not potential contact zones.

The deformed shape predicted by this code is represented in **Figure 15**. Consequently, as expected, the vertical displacement imposed on the half circle results in its indentation onto the rectangle. A part of the potential contact zone is in active contact.

In a similar way to the example described in section 11.1.1 and for comparison purposes, the problem is solved by the BEM code implementing a classic algorithm of contact by Blázquez et al. (1998a, 2006). The results of the deformed shape around the contact zone for both meshes in **Table 6** are represented in **Figure 16**. Results obtained by the codes using different meshes (see **Table 6**) and distinct lengths of the potential contact zone (see **Table 5**). The aim of this figure is to compare both algorithms, and, for that reason, results obtained with the same mesh are grouped. Macroscopically both solutions for EC-BEM and PM-BEM show a very good agreement. However, it is interesting to see in the detailed plot of solutions for mesh 1, plotted in **Figure 17**, that results for EC-BEM, even for this coarse mesh, show a very good agreement with the Reference solution.

A comparison of tractions along the potential contact zone predicted by EC-BEM, PM-BEM and the reference solution is represented in **Figure 18**. In this figure, unlike those described previously, results are grouped by the method, in order to evaluate the influence of the mesh on the results for each method. Notice that both codes present more accurate results for finer meshes, thus confirming the expected solution convergence with mesh refinement.

The resultant force along the upper part ŴD of the half cylinder can be estimated from the results of tractions predicted by both methods. In particular, we focus on the vertical component Fy of this force because it is the only one which produces work since horizontal displacements are zero by the

boundary conditions. This value is calculated by the integration of stresses over ŴD,

$$F\_{\mathcal{V}} = \int\_{\Gamma\_{\mathcal{D}}} \sigma\_{\mathcal{V}} \mathrm{d}x \tag{48}$$

where σyy represent the normal tractions at the nodes. The integral is calculated by a standard quadrature.

The values of the resultant force Fy, calculated for the 10 steps by both codes, are represented in **Figure 19** as a function of the vertical displacement imposed. Both codes predict that stiffness is increasing with the vertical displacement, as typical for the Hertz contact problems. The reason is that the contact zone between both solids is enlarging, with a consequence that the global system becomes stiffer.

FIGURE 17 | Detailed comparison of the deformed boundaries around the contact zone predicted by EC-BEM and PM-BEM for coarse meshes (Meshes 1 and 2 of Table 6) together with the Reference solution by PM-BEM (using Mesh 3 of Table 6).

### 11.2. Numerical Examples to Test the Adhesive Contact Implementation

This section aims to compare the results obtained by the energetic approach and commercial FEM code (ANSYS, 2010) in the case of the elastic contact described in section 2. To this effect, the problem of receding contact studied above is solved again, this time with the contact conditions of the adhesive elastic contact.

The problem represented in **Figure 9** is studied with the properties of model A, see **Table 3**, and the mesh 1, see **Table 4**. The stiffness parameters, assumed for this example, are κn=150 GPa/m and κt=0.

Results from the EC-BEM computations for displacements and tractions along the boundary between both solids are presented in **Figure 20**. As expected, zones with tensions correspond to positive relative displacements (external zones for example) and zones with compressions correspond to vanishing relative displacements (zones adjacent to the previous ones). In this case, a strongly complex zone is predicted at the center of the boundary, which will be studied below.

For comparison purposes, the problem is solved using the FEM code (ANSYS, 2010) setting the contact properties in the following way:

Results for tractions in the solids obtained by ANSYS and the EC-BEM code are compared in **Figure 21**. In the upper figure we can observe that the agreement along the four zones

described above is accurate. It is interesting to focus on the zone in the middle of the boundary where tractions looks like vanishing at the scale of the upper figure. Lower figure shows the complexity of this zone, which is composed by two lateral regions where both solids are separated and a central region where both solids are under contact (compression). The agreement between both results is reasonable in spite of the very small scale compared to the tractions obtained on the whole boundary.

A subsequent analysis is carried out by multiplying the displacement uD by 100 without modifying any other parameters in both codes. In this case, the results from the energetic approach presented here remain qualitatively similar, whereas the results from ANSYS present some problems to match the complex central zone. To the knowledge of authors, this could be solved by modifying some of the values of the keyoptions which configure the contact algorithm of ANSYS. This shows the robustness of the energetic approach compared to other contact algorithms, where convergence depends on the value of certain non-physical parameters of the numerical model.

#### 12. CONCLUDING REMARKS

A novel approach for contact problems of elastic bodies was presented. This method is based on the energy principles expressed by boundary integrals. The solution of contact problems is obtained by the minimization of the total potential energy. Appropriate algorithms to solve the quadratic problems were presented. The cases of conforming as well as non-conforming meshes for contact zone were addressed.

Numerical examples, showing the suitability of the present framework to solve typical contact problems were examined,

and also successfully compared to other well established BEM methodologies as well as a commercial FEM code.

The proposed framework can be easily incorporated to existing BEM or FEM codes. Possible extension to the case of elastodynamics, where instead of the principle of virtual

#### REFERENCES

work we may use the principle of virtual power in order to compute the total energy (kinetic plus elastic), using only boundary values, of tractions and velocities, is possible. For this it might be convenient to utilize a formulation such as that introduced in Panagiotopoulos and Manolis (2010, 2011).

Finally, the present framework is found to be very useful in the case where also dissipative mechanisms exist on the boundaries or common interfaces of elastic bodies. In such cases we may apply energetic approaches to solve the corresponding non-linear problems as was shown in previous publications by some of the present authors (e.g., Panagiotopoulos et al., 2013).

### AUTHOR CONTRIBUTIONS

All authors listed have made a substantial, direct and intellectual contribution to the work, and approved it for publication.

#### FUNDING

The authors acknowledge the support by the Junta de Andalucía and European Social Fund (Proyecto de Excelencia TEP-4051), the Spanish Ministry of Science and Innovation (Project MAT2009-14022), the Spanish Ministry of Economy and Competitiveness and European Regional Development Fund (Projects MAT2012-37387 and MAT2015-71036-P), and Spanish Ministry of Education (FPU grant AP2009-3968). Furthermore, CP acknowledges the financial support of the Stavros Niarchos Foundation within the framework of the project ARCHERS (Advancing Young Researchers' Human Capital in Cutting Edge Technologies in the Preservation of Cultural Heritage and the Tackling of Societal Challenges).

#### ACKNOWLEDGMENTS

The authors would like to express their gratitude to Prof. Tomáš Roubícek of the Charles University of Prague, as well as to Dr. ˇ Roman Vodicka from the Technical University of Košice, for ˇ fruitful discussions.

**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2018 Panagiotopoulos, Mantiˇc, García and Graciani. This is an openaccess article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Hybrid Micro-Modeling Approach for the Analysis of the Cyclic Behavior of RC Frames

#### Gianni Blasi 1,2, Flavia De Luca2 \* and Maria Antonietta Aiello1

1 Department of Innovation for Engineering, University of Salento, Lecce, Italy, 2 Department of Civil Engineering, University of Bristol, Bristol, United Kingdom

The present study is aimed at developing a hybrid approach to consider the effect of concrete cracking on the hysteretic response of RC frames. The mechanical behavior of the concrete is defined according to the smeared cracking approach, while discrete cracking surfaces are included in the geometrical model. The interface behavior of the discrete cracking surfaces is defined by the combination of contact and cohesive elements. The proposed approach is adopted in ABAQUS to simulate an experimental test on a double cantilever column for the calibration of the numerical model. Therefore, a test conducted on a RC portal and modeled numerically for the first time is simulated. Numerical and experimental results are compared in terms of hysteretic forcedisplacement behavior and cumulative dissipated energy, in order to assess the reliability of the proposed model in simulating the energy dissipation capacity of RC members subjected to lateral cyclic loading. The hybrid modeling approach proposed allows an accurate description of the stress distribution and a fairly satisfactory matching of the hysteretic behavior with a reasonable compromise in terms of computational effort.

Keywords: reinforced concrete frames, cyclic behavior, hysteretic response, micro-modeling, abaqus, hybrid micro-model

#### INTRODUCTION

Finite element (FE) analysis is widely adopted in earthquake engineering, particularly when dealing with systems with a high number of degrees of freedom (DOF), representing the structural configuration of buildings. On the other hand the simulation of the seismic response of RC structures through finite element analysis can be complex and challenging. This is mainly caused by the difference between the mechanical response of concrete and steel rebars and by their consequent complex interaction. In numerical studies, either micro- or macro-modeling approaches are usually adopted depending on the aim of the study, which can require the accurate simulation of several phenomena influencing the response of the frames. Macro-models generally adopt simplified formulation which implicitly consider the influence of the micro-mechanical interaction phenomena on the response of the elements/materials. This approach leads to a significant simplification of the numerical problem. In micro-modeling approaches, the micro-mechanical interaction phenomena are explicitly modeled numerically (e.g., Karavelic et al., 2017; Sinaie ´ et al., 2018) and generally high computational efforts are required. Therefore, when the target of the analysis is global behavior of RC structures, generally macro-models are preferred (e.g., Vamvatsikos and Fragiadakis, 2010; De Risi et al., 2017). This is mainly due to the substantial reduction of the computational effort required with respect to micro- and meso-models. In the

#### Edited by:

Vagelis Plevris, OsloMet–Oslo Metropolitan University, Norway

#### Reviewed by:

Michalis Fragiadakis, National Technical University of Athens, Greece Francesco Tornabene, Università degli Studi di Bologna, Italy Zeljana Nikolic, University of Split, Croatia

> \*Correspondence: Flavia De Luca<EMAIL_ADDRESS>
#### Specialty section:

This article was submitted to Computational Methods in Structural Engineering, a section of the journal Frontiers in Built Environment

> Received: 08 August 2018 Accepted: 22 November 2018 Published: 10 December 2018

#### Citation:

Blasi G, De Luca F and Aiello MA (2018) Hybrid Micro-Modeling Approach for the Analysis of the Cyclic Behavior of RC Frames. Front. Built Environ. 4:75. doi: 10.3389/fbuil.2018.00075 last decades, large numbers of macro-models have been developed, able to accurately reproduce the response of RC members. On the other hand, for the accurate reproduction of stress/strain distribution in the elements subjected to cyclic loads, micro- and meso-models still represent the best option. The latter approach represents the compromise between microand macro-modeling, featuring a simplification of the numerical problem aimed at accurately evaluate only the phenomena which are object of the analysis. For this reason, the difference between macro- and meso-models can't be unambiguously defined. In macro-models, RC frames are usually modeled through beam elements, while both micro- and meso-models generally feature solid or shell elements to represent the concrete region, whereas the steel reinforcement is simulated with embedded beam/truss elements (e.g., Viswanathan et al., 2014; Redmond et al., 2018).

Mechanical models for concrete and steel, available in most of the FE platforms, allow reproducing the stress-strain response under monotonic loading and the degradation of the strength and stiffness due to cyclic loading. The mechanical response of the steel is generally reproduced adopting a symmetric bi-linear behavior, characterized by kinematic or isotropic hardening laws; even if in some cases the compressive response is modified to take into account the buckling phenomenon on rebars.

Although the simulation of the concrete response is more challenging; as well-known, the compressive behavior is characterized by a post peak softening slope, influenced by different parameters, such as the confinement (Mander et al., 1988), while the tensile response is generally composed of an elastic branch up to the tensile strength of the concrete, followed by concrete cracking.

The post-cracking response of concrete can be simulated adopting the smeared cracking approach (Lotfi and Shing, 1991; Lee and Fenves, 1998), where cracks are not explicitly defined in the geometry of the model but they are implicitly considered in the mechanical response of the material.

The smeared cracking approach has been widely adopted in numerical studies to reproduce the monotonic response of RC frames (Sinaei et al., 2012; Sümer and Akta¸s, 2015; Ors et al., 2016); although this mechanical model is not always able to accurately reproduce the response of concrete under cyclic loading.

The discrete cracking approach is an alternative method to reproduce the crack opening in the concrete (Dolbow et al., 2001; Giner et al., 2009), it generally allows high simplifications of the mechanical model adopted for the concrete. Discrete cracks develop in the mesh once the tensile strength or the cracking energy is attained and they can be explicitly identified in the post processing.

The two methods discussed allow the accurate reproduction of the effects of crack development in concrete members, but, in many cases, they are time consuming and lead to convergence issues.

Furthermore, when dealing with the cyclic behavior of the concrete members, the effects of the cyclic opening and closure of the cracks should be considered in the numerical simulation; in the smeared cracking approaches, this leads to a significant increase of the complexity of the mechanical model (Koutromanos and Shing, 2012; Redmond et al., 2018), while the eXtended Finite Element Method (XFEM) requires the definition of a contact law along the cracking surfaces, which is not predefined in the geometrical model (Yu et al., 2016).

The present study is aimed at proposing a hybrid approach using ABAQUS (Dassault Systemes, 2016) which combines the discrete cracking approach and the smeared cracking approach. The cracking surfaces are pre-defined in the numerical model by splitting the geometry in different "macro blocks"; along the cracking surfaces, cohesive and contact laws are defined, in order to simulate the crack opening and closure, respectively.

This approach makes possible to simplify the definition of the mechanical response of the concrete, since the cyclic opening/closure is not implicitly taken into account in the tensile behavior definition. Furthermore, since the cracking paths are predefined in the geometry, the contact and cohesive laws are only defined for specific surfaces in the geometry, leading to a significant reduction of the computational efforts.

The proposed approach is aimed at reproducing the hysteretic response of RC members subjected to cyclic lateral loading. In order to assess the reliability of the proposed model, a test on a double cantilever column was simulated (Ohue et al., 1985) by proposing different modeling approaches. The monotonic and hysteretic response of the numerical simulations were compared to the experimental results with a particular focus on the influence of the adopted model on the accuracy of the reproduction of the shape of the hysteresis loops. Then, a recent experimental test on a RC portal conducted by Verderame et al. (2016) is modeled numerically for the first time in order to assess preliminarily the reliability of the proposed approach. Hysteretic force-displacement response is evaluated and the numerical and experimental results are compared, in order to assess the accuracy of the simulation of the cyclic loops.

# MICRO-MODELING OF RC FRAMES—ISSUES IN THE SIMULATION OF THE CYCLIC BEHAVIOR

Most of the mechanical models available on FE platforms to simulate the behavior of the concrete take into account the significant asymmetry of the response in compressive and tensile direction. The compressive stress-strain relationship of the concrete is represented through a non-linear formulation (**Figure 1A**) whose parameters depend on the reinforcement details, as the post-peak softening response, being highly influenced by the transverse reinforcement (Mander et al., 1988). On the other hand, the tensile behavior definition depends on the approach adopted to simulate the concrete cracking.

The response of RC members under cyclic lateral loading is indeed highly influenced by the cracking of concrete, which modifies the energy dissipation capacity. The pinching effect characterizing the flexural response of frames is caused by the cyclic closure/opening of cracks, leading to a significant reduction of the stiffness in the load-inversion phase, due to the closure of the cracks. Depending on the adopted approach, the cyclic loops characterizing the lateral response of a RC frame

can be reproduced in different ways, explained in the following section.

#### The Smeared Cracking Approach

In the smeared cracking approach, the development of cracks in the concrete is considered by defining a post-peak tensile response of the material. Thus, rather than tracking individual "macro" cracks in the geometrical model, the presence of cracks is considered implicitly by modifying the stress and the material stiffness associated with the integration point. After the attainment of the tensile strength of the concrete (fct), a post-peak softening branch (tension stiffening) simulates the progressive degradation due to the crack openings. Tension stiffening is reported in **Figure 1B** and can follow either linear or exponential laws.

In this approach, the simulation of the crack closure, requires an increase of the complexity of the material model (Koutromanos and Shing, 2012) and, as said before, highly influences the dissipation capacity of the frame members.

As reported in **Figure 1B**, once unloading takes place, after the crack has opened, the unloading branch stiffness is equal to the elastic stiffness, reduced by a factor δt that takes into account the degradation due to damage. As soon as the elastic displacement is restored, the crack closure stage begins and a significant reduction of the stiffness up to the complete closure of the crack is observed. At this point, the elastic stiffness is restored and compression takes place; the variation of the stiffness during the load inversion phase (pinching) is more marked in case of relevant crack patterns in RC members.

The analytical formulation to simulate the monotonic response of concrete usually requires the definition of two yielding surfaces; the first one is referred to the end of the elastic branch in compression, at which the linear stress-strain law turns into an exponential one. The second yielding surface is referred to the end of the elastic slope in tension, where tension stiffening takes place.

The model proposed by Oliveira and Lourenço (2004) was adopted in finite-element analysis of masonry structures to simulate the cyclic behavior of interface elements; although the same formulation can be used for concrete. The cyclic response is defined by the adoption of two auxiliary yielding surfaces moving within the monotonic yielding surface, representing the unloading/reloading phase.

The monotonic surfaces are active when either only monotonic loading has occurred or the stress point reaches the monotonic yielding surfaces during an unloading process. The proposed model was validated by comparisons with experimental results of uniaxial tests on concrete and was assessed to be highly reliable to simulate all the phenomena characterizing the cyclic response of asymmetric materials as concrete.

In some cases, mechanical models available in finiteelement software packages are unable to simulate the response under cyclic loads, particularly in case of large tensile strain/displacement (**Figure 2A**); compression stiffness recovery in the reloading phase can occur earlier than expected, depending on the criterion adopted in the model to discern the transition from tension. As observed by Nikaido et al. (2015), if the compression stiffness recovery is detected at zero stress point, the pinching effect simulation can be challenging, since the crack closure cannot be accurately reproduced in the numerical model (as shown in **Figure 2B**).

In other cases, the auxiliary yielding surfaces are not even defined and the stiffness degradation, representing the cracking closure, does not take place.

If the mechanical model is unable to produce a satisfactory simulation of the pinching effect, a significant mismatch is obtained in the numerical simulation of the flexural response of RC members in terms of energy dissipation, since cracking of the concrete leads to a noticeable decrease of the energy dissipated.

# The Discrete Cracking Approach

In the discrete cracking approach, the cracking process is ruled by cracks developed in the geometrical model; the most used modeling approaches are the extended finite element method and the adoption of cohesive elements to express the post-cracking response of the concrete.

The XFEM was developed by Moës and Belytschko (2002); it allows to significantly simplify the definition of the tensile response in the mechanical model of the material. The cracking path is generated by adding discontinuous basis functions to standard basis functions for nodes where cracking energy/displacement is computed, in order account for the displacements due to crack opening. The main advantage in the XFEM is that no update of the mesh is required to track the crack path, thus, the mesh-dependency of the results is significantly reduced with respect to a smeared cracking approach.

Since the XFEM allows to model the discontinuity in a displacement field along the crack path, the simulation of the crack growth requires an accurate definition of the fracture modes; generally, the cracking path is computed adopting the mixed mode in order to consider not only the response under uniaxial traction (Ballatore et al., 1990). Dolbow et al. (2001) also incorporated friction laws in order to simulate more accurately the crack growth under compression, which is an important feature when dealing with cyclic response of the frame members.

Since the crack path is not pre-defined, the accurate simulation of the cyclic response of reinforced concrete members requires the inclusion of contact laws in order to consider the effect of the crack closure. This assumption generally leads to high computational effort, and in many cases can cause convergence issues, particularly when dealing with RC members with embedded rebars, for which multiple crack paths are expected to develop under cyclic lateral loading.

Discrete cracking approaches adopting cohesive elements rather than XFEM have been also widely adopted in the literature (Bocca et al., 1991). The so-called discrete interelement cracks require remeshing in the model after the generation of the cracks, since the path is defined between the elements. This method is similar to the XFEM since the crack path is not predefined and can lead to high computational efforts as well as convergence issues. Notwithstanding the discussed issues, it was used specifically to simulate the lateral response of reinforced concrete members as well as crack growth in ductile materials (Tvergaard and Hutchinson, 1996; Gullerud et al., 2000).

The so-called intraelement approach features the incorporation of embedded discontinuities in the geometry (Belytschko et al., 1988), to simulate the post cracking response of the material. This approach allows the pre-definition of the crack paths, identified by a preliminary analysis of the stress distribution along the member, in order to define the regions where cracks are most likely to develop. Despite this approach reduces the complexity of the model, in some cases the preliminary definition of the crack paths may be difficult, as in the case of RC members subjected to cyclic loading.

In the study proposed by Yu et al. (2016), the simulation of the pinching effect characterizing the flexural cyclic response of a RC column, was conducted providing a novel modeling approach, developed through combination of the XFEM, to detect the crack paths in concrete, traction-separation cohesive laws to simulate the post-cracking response, and contact elements to take into account the clack closure process. A satisfactory simulation of the experimental response of RC frame members was obtained and the comparison between experimental and numerical hysteretic loops demonstrated the high reliability of the model in simulating the pinching effect.

An alternative discrete cracking approach for the simulation of the lateral response of RC members was proposed by Stavridis and Shing (2010). Triangular elements were connected with double-node interface elements, simulating the post cracking behavior of the concrete. The model was validated by comparing numerical results with experimental tests on infilled frames, showing a good match with the experimental backbone curve (Stavridis and Shing, 2010; Redmond et al., 2018).

A similar approach was recently developed from Zivaljic et al., 2013, combining triangular plane-stress elements with discrete cracking interfaces to reproduce concrete cracking and embedded one-dimensional elements simulating steel rebars. The bar slip was also considered by defining a pre- and post-cracking mechanical model for steel. The proposed model was adopted for the simulation of several reinforced concrete elements and was assessed to be very reliable in reproducing the crack paths in concrete (Zivaljic et al., 2013; Živaljic and Nikoli ´ c, 2014; Nikoli ´ c´ et al., 2017).

Also in this cases, a satisfactory simulation of the pinching required a complex mechanical model for interfaces elements, increasing the computational efforts.

# THE SIMULATION OF THE HYSTERETIC RESPONSE IN THE PROPOSED APPROACH

As discussed in section Micro-Modeling of RC Frames—Issues in the Simulation of the Cyclic Behavior, the main issues in the simulation of the cyclic response of RC members under lateral loading are related to the pinching effect; when the post cracking behavior is defined in the mechanical model of the material, the process of crack closure requires a significant increase of the complexity of the formulation, while in discrete cracking approaches, the cracking closure must be simulated by the inclusion of contact laws, which can lead to convergence issues. XFEM method was assessed to be reliable when simulating the lateral response of RC members under monotonic lateral loads even if the post-cracking behavior of the concrete is not defined in the material model (Zi and Belytschko, 2003). On the other hand, the adoption of this method requires an oversimplification of the response of the material in compression, which in many cases was assumed to be linear elastic (Dolbow et al., 2001; Liao and Huang, 2015). When analyzing the non-linear response of frame members experiencing displacements far beyond the yielding point, the simplified assumption on the compressional behavior of the material becomes questionable and the post-peak actual response of the element can be mismatched.

#### Description of the Numerical Model

In the numerical model developed using ABAQUS, an equivalent bi-dimensional system is proposed to simulate the cyclic lateral behavior of a RC portal. 4-node bilinear shell elements with reduced integration (CPE4R) are adopted for concrete, while both longitudinal and transverse rebars are modeled using linear truss elements (T2D2), embedded in the concrete geometry (**Figure 3A**).

The mechanical behavior of the concrete is defined adopting a Concrete smeared cracking model available in the software. As evidenced by Nikaido et al. (2015), the latter is unable to accurately simulate the crack closure process, since earlier compressive stiffness restore occurs in the unloading phase when the tensile deformation is relevant. Thus, a combination of the smeared and discrete cracking approach is adopted to simulate the pinching effect.

A bilinear kinematic hardening relationship was adopted for the definition of the mechanical behavior of the steel rebars, assuming a symmetric response in tension and compression.

In order to accurately define the compressive response of the concrete, the well-consolidated and widely adopted formulation proposed by Mander et al. (1988) was considered; which takes into account the transversal reinforcement ratio to define the ductility and strength capacity of the material. The tensile response was assumed to be elastic-perfectly-plastic with no tension stiffening, since the crack opening was modeled following a discrete approach by the adoption of cohesive laws. The model was implemented by direct input of the stress-strain backbone curve in the mechanical model available in the software (Concrete damaged plasticity); the parameters required for the definition of the plastic flow and the yield function are assumed considering values suggested from the literature (Lubliner et al., 1989; Lee and Fenves, 1998).

#### The Interface Interaction

According to different experimental and numerical studies (Sneed et al., 2016; Duan et al., 2017), flexural cracks are most likely to develop along the discontinuity surfaces, represented by the transversal rebars running through the concrete. Thus, in the proposed model, discrete crack paths are preliminarily defined in the geometry.

The path initiation position corresponds to the stirrups in the frame (**Figure 3B**), while the inclination is defined according to the variable strut inclination truss model (EN 1992-1-1, 2004; Fardis, 2009), assuming an angle of inclination of the compression stress field equal to 21.8◦ . Moreover, in order reduce the computational efforts, the paths are only included in the zones with the highest value of flexural moment. The geometry of concrete members is divided in different "macro-blocks," whose interaction is defined through cohesive and contact relationships.

The interface zone between the macro blocks is divided in three parts, as shown in **Figure 3B**; the external parts represent the crack paths, while the central part is a rough connection assumed as the center of the rotation of each block. The assumption of a central rough connection means that cracks in the members can only evolve up to the external boundaries of the central zone; as consequence, the latter must be adequately restricted, in order to assume that the collapse of the member is attained before cracks develop up to the central zone.

Depending on the direction of the load, the two external zones are characterized by compression and tension stress regime, respectively (**Figure 3B**); thus, either cohesive or contact laws are activated along the surfaces, on the basis of the stress regime.

The contact law is defined in order to simulate the complete closure of the crack, when the elastic compression stiffness is restored. The augmented Lagrange method (Powell, 1969) is adopted to compute the contact stress, while the contact stiffness is automatically defined depending on the stiffness of the elements experiencing the contact. This approach also allows to take into account the compressive stiffness reduction due to cyclic

loading, by including compressive degradation in the mechanical model of the concrete.

The stiffness factor for the contact surfaces is assumed to be unitary, both for external and central surfaces.

In the initial stage, no gap exists between the blocks, and all the contact surfaces are in a closed status; since for the external zones separation of the surfaces is expected, the contact detention is ruled by the definition of a minimum gap between the nodes, this assumption assures contact elements being active every time that surfaces cyclically experience contact.

Once traction takes places along the surface, the separation between the surfaces follows a stress-displacement law, obtained as a function of the post-cracking stress-strain relationship of concrete, reported in **Figure 2A**.

The first elastic branch is defined considering the stiffness of the shell elements, thus, it takes into account the stiffness degradation by including a degrading tensile law in the mechanical model of the concrete.

The stress limit at the end of the elastic branch is equal to the tensile strength of the concrete. The post-peak behavior is ruled by an exponential law, according to analytical models expressing the tensile response of concrete available in the literature (Pramono and Willam, 1989; Mehta and Monteiro, 2006). In this study, an exponent equal to three is assumed to define the softening branch after the attainment of the tensile stress. Once the displacement between the surfaces reaches the residual limit, the cohesive law is inactive, meaning that no stress-displacement relationship exists between the surfaces. At this stage, once the load is reversed, the closure of the crack is represented by a displacement with zero-stress variation along the elements, making possible to simulate the pinching effect in the global response of the frame. As said before, once the crack is completely closed, the elastic compressive stiffness is partially restored, according to the compressive degradation law defined in the mechanical model of the material.

By this procedure, the pinching effect is "lumped" in the cracking interfaces; beside this, the effect of the cracking of the concrete on the flexural response in terms of strength is still considered in the whole model, by the definition of a tensile strength limit for the concrete, as well as tension stiffening.

# NUMERICAL SIMULATION OF QUASI-STATIC CYCLIC TEST ON RC PORTAL

The numerical model developed was firstly validated considering the experimental test conducted on a double-cantilever column by Ohue et al. (1985). The column tested (named non-ductile column or NDC), reported in **Figure 4** featured a cross section equal to 200 × 200 mm and a height equal to 800 mm. The longitudinal reinforcement was composed of 2 + 2 rebars with 16 mm diameter, while transverse reinforcement was realized adopting 5.5 mm stirrups with spacing equal to 50 mm.

Axial load was applied on the top of the column to simulate the presence of the gravity loads. The cyclic lateral load was applied in displacement-control and it was composed of 9 load steps, each

of them characterized by two fully reversed cycles at the same displacement amplitude.

A mixed shear-flexural failure mode of the specimen was observed in the test, featuring major shear cracks at a displacement amplitude equal to 17.5 mm (drift equal to 2.2%).

In order to assess the reliability of the proposed model in predicting the lateral behavior of RC portals, a campaign available in literature and never modeled numerically, was considered to carry out a numerical simulation. The reinforced concrete frame (named ductile frame or DF) tested by Verderame et al. (2016) was designed according to provisions of the Italian seismic code NTC08 (2008), in order to simulate the behavior of seismic designed frames.

In DF (**Figure 5**), the longitudinal reinforcement of the beam is symmetric and made of three top and three bottom 10 mm bars. 12 mm bars were used for the column reinforcement, in order to follow the capacity design approach. Transverse reinforcement of the beam and the columns was composed of 6 mm stirrups, whose spacing was equal to 100 mm in all the three elements and 50 mm in the "critical zones," i.e., the end regions of columns and beams and in the mid span of the beam.

The experimental test was carried out by the application of cyclic lateral load at the mid-span of the beam; the loading protocol was composed of fully reversed displacement-controlled load cycles, whose displacement amplitude was increased at each load step up to the failure of the specimen. For each load step, three fully reversed cycles were executed.

The axial load was applied on top of the column during the quasi-static cyclic tests, to simulate gravity loads due to the presence of higher floors in a five-story building.

In the test, the first damage was characterized by flexural cracking at the beam ends. By increasing the displacement amplitude, flexural cracks developed in the column ends, too. The peak lateral strength was attained after the widening of the flexural cracks in beam ends and in the bottom of the columns. As expected for seismic designed frames, the collapse mechanism was caused by major cracking in beams, according to a strong column/weak beam mechanism behavior. The failure of the

specimen was attained at a drift 1f equal to 5.3% (see Verderame et al., 2016 for further details).

# Numerical Simulation of the Tests in ABAQUS

In the numerical simulation in ABAQUS, different configurations of the finite element model are analyzed, in order to evaluate the influence of the proposed modeling approach on the numerical results. A Finite-Element "standard" (FES) model is firstly developed, with the aim of assessing the accuracy of the material model available in ABAQUS in simulating the pinching effect. The only difference between FES model (reported in **Figure 6A**) and the Finite-Element model with macro Blocks (FEB) model proposed in section Description of the Numerical Model is the absence of the macro blocks in the geometry, meaning that a fully smeared cracking approach is adopted, without the definition of contact and cohesive laws to simulate the cyclic opening and closure of cracks, leading to a high simplification of the numerical problem.

In FES model, the size of shell elements is assumed equal to 50 mm in the whole geometry, while, in FEB model, it is gradually reduced from 50 to 10 mm in proximity to the interfaces between the macro blocks. Referring to the truss elements adopted to simulate longitudinal rebars, a length equal to the stirrups spacing is assumed.

Furthermore, referring to FEB model, two different configurations were developed (**Figures 6B,C**) for the test considered; in order to analyse the influence of the simulation of bond between concrete and steel rebars on the accuracy of the simulation of the hysteretic loops characterizing the lateral response of the portal.

In the first configuration (FEB1), the end zones of beams and columns are divided in two macro blocks, meaning that two discrete cracking surfaces are included. Bar slip is not considered and the truss elements are embedded in the concrete region. In the second configuration (FEB2), which features the same number of blocks adopted for FEB1, bar slip is modeled through linear elastic connectors between trusses representing longitudinal rebars and shell elements. The mechanical behavior of the connectors is defined according to CEB-FIP (2010).

It is worth mentioning that the same approach was adopted for each of the experimental tests simulated, even if **Figure 6** is referred to the model adopted for the portal tested by Verderame et al. (2016). The double-cantilever model for the simulation of the test by Ohue et al. (1985) featured an encastre restraint for both the bottom and the top nodes of the column.

A dynamic analysis adopting implicit method was carried out; monotonic and cyclic displacements were applied on top of the frame and the column as boundary condition to reproduce the loading protocol applied in the experimental tests, while the base shear at the fixed supports was computed during post-processing in order to plot monotonic and cyclic loaddisplacement relationships.

# Comparison Between Numerical and Experimental Results

The numerical simulation of the test on the NDC confirms the low reliability of the mechanical model of the concrete available in the software in simulating the cyclic response of reinforced concrete columns. Firstly a monotonic push-over analysis was carried out, in order to obtain the monotonic response of each model developed (**Figure 7**). No significant variation in observed in terms of yielding strength and pre-yielding stiffness for FES and FEB1, confirming the equivalence of the two models in the simulation of monotonic behavior.

On the other hand, FEB2 model featured a lower stiffness due to the bar slip. The variation of the yielding strength with respect to the experimental test is equal to +1.5, +.2%, and +1.1% for FES, FEB1, and FEB2, respectively. The elastic stiffness increased by 13% for both FES and FEB1 while a decrease of 36% was observed for FEB2. The FES and FEB1 featured a hardening response leading to an increase of the ultimate strength with respect of the experimental results. This feature suggests that the strength degradation obtained in the experimental test could be

an effect of the cyclic loading, therefore it couldn't be reproduced through a push-over analysis. On the other hand, the strength degradation obtained in FEB2 simulation, was caused by the bar slip, leading to a softening behavior and a higher accuracy of the simulation in terms of strength.

The comparison between the hysteretic response of FES and FEB1 and the experimental test is reported in **Figure 8**. The introduction of the macro-block in the model leads to a significant reduction of the energy dissipation and the elastic stiffness of the backbone curve with respect to the FES model. Furthermore, the pinching effect observed in the experimental response is only simulated in FEB1, confirming the limit of the original model for the concrete in accurately simulate the cracking closure process. Lastly, the isotropic hardening of the reinforcement steel observed in FES simulation led to a significant increase of the maximum strength with respect to the experimental test as well as to the corresponding maximum strength obtained from the push-over analysis. Referring to FES, the maximum strength obtained increased by 14% with respect to the experimental result (+7% with respect to the monotonic response), Also in FEB1 model the isotropic hardening affected the global response in the first cycles, as evidenced by the different shape of the hysteretic loops comparing numerical and experimental response for low displacement amplitude. The increase in maximum strength observed for FEB1 was 11% comparing to experimental test (+3% with respect to the monotonic response). Despite being affected by the isotropic hardening model of the reinforcement steel, the cyclic response of FEB1 featured pinching effect and a satisfactory simulation of the experimental response was obtained comparing to experimental curve.

The introduction of the bond between steel rebars and concrete led to a significant increase of the accuracy of the experimental response simulation (**Figure 9**). The pinching effect simulation caused by the bar slip significantly affected the shape of hysteretic loops and led to a reliable simulation of the total energy dissipated (−19% with respect to the experimental test).

A good correspondence in terms of maximum strength was also observed, with a difference of 4.1% between numerical and experimental results. On the other hand, the inclusion of the bond in the numerical model affected the elastic stiffness, which was 45% lower in the numerical simulation.

The total energy dissipated, ED is reported in **Table 1**. In the FES a significantly higher value is obtained with respect to the ED of the experimental test (+252%), due to the absence of pinching and the isotropic hardening effect obtained in the model. For FEB1 and FEB2, a higher accuracy in the simulation of the crack opening/closure due to the presence of the macroblocks led to a satisfactory matching of the experimental curve and, consequently, the ED. Furthermore, the bar slip simulation characterizing the FEB2 approach reduced the influence of the cyclic response of the reinforcement steel on the global behavior of the column.

The numerical simulation of the test on the DF specimen shows a significant variation of the results depending on the

TABLE 1 | Total energy dissipated for experimental test and numerical models referring to NDC.

modeling approach adopted. The hysteretic force-displacement (F-d) behavior obtained considering the FES model (**Figure 10A**) confirms the assessments made in previous sections. No pinching effect is observed in the cyclic response of the frame, leading to a significant increase of the energy dissipated in the numerical model with respect to the experimental test. Furthermore, isotropic hardening effect is observed in hysteretic loops obtained from the numerical simulation, caused by the influence of the cyclic behavior of steel rebars on the lateral response of the frame.

In FEB1 (**Figure 10B**), the bar slip is not simulated since embedded constraints were adopted between the nodes of shell and truss elements. For this reason, an increase of the maximum lateral stiffness is obtained in the numerical model compared to the experimental test. The absence of bar slip increases the influence of the rebar deformation on the strain developed in the shell elements, and the crack opening is not adequately simulated. On the other hand, higher quality of the matching of the experimental test can be observed in terms of strength and stiffness, even if the pinching effect is not simulated.

Comparing FES results to experimental ones, an increase of 199 and 24% is obtained in terms of elastic stiffness, Ek , and yielding strength, Fy, respectively. Considering FEB1, the Fy is 11% higher than the value obtained from the test, while Ek increases by 111%.

For FEB2 model (**Figure 11**), a better matching of the experimental results is obtained both in terms of strength and stiffness. The simulation of the bond leads to a significant reduction of the elastic stiffness with respect to FEB1 (−64%), while bar slip was found to have a sensible influence on the pinching effect.

Comparing experimental to numerical results, an increase of 10% is obtained for Fy considering FEB2, while the elastic stiffness decreases by 23%. It is worth saying that FEB2 simulation is the only one which features a decrease of the elastic stiffness with respect to the experimental test.

Lastly, the energy dissipation capacity is also computed for the analyzed frame. The cumulative dissipated energy (CDE), evaluated as the sum of the areas of the hysteretic loops in the F-d response, is reported in **Figure 12**. Considering FES and FEB1, a significant increase of energy dissipated is observed even for low displacement amplitude, while in case of FEB2, a very satisfactory match is obtained in the early stage.

This assessment is confirmed by the comparison of the hysteretic loops of **Figures 10**, **11**.

The pinching effect is more noticeable in the first loops in case of FEB2, while the more the displacement amplitude increases, the more the difference between the experimental and the numerical CDE increases.

This effect could be caused by the increase of the strain in shell elements adjacent to the cracking surface for high displacement amplitudes.

For high displacements, the interaction between the stress components in shell elements adjacent to the crack leads to a significant increase of the strain, causing a distortion of the geometry in the crack closure phase. For this reason, the cyclic

loops obtained for high displacement amplitude are significantly influenced by the hysteretic behavior of the steel, which is not characterized by pinching.

The total energy dissipated, ED, is reported in **Table 2**. As expected, ED in the FES and FEB1 models is significantly higher with respect to the ED of the experimental test, due

TABLE 2 | Total energy dissipated for experimental test and numerical models referring to DF.

to the absence of pinching and because of the higher stiffness obtained in these simulations. For FEB2, a better simulation is obtained, confirming the need of accounting for the bar slip in the simulation of the hysteretic response when dealing with reinforced concrete members.

The efficacy of the proposed model was also assessed by comparing CPU time of the numerical analysis (**Figure 13**). Since the time is significantly influenced by the processor clock rate, the comparison is provided in terms normalized CPU time (CPUn) referring to the FES model, which is the simplest configuration developed.

As expected, the absence of the discrete cracks in the FES model led to a significant reduction of the CPU time with respect of both FEB1 and FEB2. The presence of the macro blocks and, consequently, the inclusion of contact and cohesive laws increased by 307% the solver time, while by including the bond between concrete and steel rebars the increase was 405%.

Even if more complex and more computationally demanding, the FEB2 model is still an appealing and simply-to-model alternative for FE analyses with reduced effort with respect to more complex approaches (e.g., Stavridis and Shing, 2010; Yu et al., 2016).

# CONCLUSIONS

Different numerical studies focussed on the seismic response of reinforced concrete members pointed out the issues when attempting to accurately reproduce the effects of the cyclic opening and closure of the cracks in concrete when adopting meso- and micro-modeling approaches. In different finite element platforms, mechanical models simulating the concrete tensile response as well as extended finite element methods available allow to accurately reproduce the crack opening process but lead to a significant increase of the computational efforts as well as convergence issues when simulating the effects of the crack opening on the cyclic response of the frames. Furthermore, in some cases the results obtained adopting a smeared cracking approaches show a significant difference compared to the actual behavior observed in experimental tests.

The hybrid modeling approach proposed in the present study was assessed to significantly improve the accuracy of the simulation of the hysteretic loops characterizing the lateral behavior of reinforced concrete members. The smeared cracking approach adopted to define the mechanical response of the concrete was combined with a discrete cracking approach through the introduction of pre-defined cracking path in the regions subjected to higher values of tensile stress.

Contact and cohesive interaction laws were adopted to define the mechanical response along the cracking paths in order to simulate the cyclic crack opening and closure. Two different configurations of the proposed model were developed, in order to evaluate the influence of the simulation of the bond on the accuracy of the matching of the experimental hysteresis loops. The simulation of the test on the double cantilever column evidenced the major limits of the existing model of the

#### REFERENCES

Dassault Systemes (2016). Abaqus Theory Manual. Vélizy-Villacoublay.

concrete in simulating the pinching effect. A better result was obtained adopting the proposed modeling approaches, although the simulation of the bond significantly affected the elastic stiffness. The comparison with a portal frame experimental test available in the literature also showed that the adoption of the smeared cracking approach results in a numerical model unable to accurately reproduce the hysteretic response, while the introduction of discrete cracking paths leads to a better fit of the experimental results in terms of strength and stiffness.

Also in this case, again, the simulation of the bond was found to be fundamental to better reproduce the pinching effect and to obtain a fair matching in terms of energy dissipation capacity.

Despite the significant improvement in the simulation, results of the most sophisticated model with bond showed that for high drift amplitudes the stress developed in shell elements along the cracking surface led to a distortion of the geometry of the concrete region. Consequently, the hysteretic loops corresponding to higher imposed displacement are mostly influenced by the behavior of the steel and pinching becomes less pronounced.

The results obtained are useful to lay the basis for the development of hybrid models to be adopted for the numerical simulation of the cyclic behavior of reinforced concrete members having fairly sustainable computational efforts. Further studies will be conducted in order to increase the accuracy of the simulation even for higher displacement amplitudes.

# AUTHOR CONTRIBUTIONS

GB, FD, and MA conceived of the presented idea. All of the authors contributed to the development of the numerical model and carried out a literature review. All authors discussed the results and contributed to the final manuscript.

# FUNDING

This work was supported by the Royal Society International Exchange Scheme. Award Number: IE161016. Project: Cumulative damage of masonry infills caused by anthropogenic earthquakes.

Methods Appl. Mech. Eng. 190, 6825–6846. doi: 10.1016/S0045-7825(01) 00260-2

and numerical implementation. Comput. Methods Appl. Mech. Eng. doi: 10.1016/j.cma.2017.09.013. [Epub ahead of print].

**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2018 Blasi, De Luca and Aiello. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# A Comparison of Three Different Methods for the Identification of Hysterically Degrading Structures Using BWBN Model

#### Ying Zhao1 , Mohammad Noori 1,2 \*, Wael A. Altabey 1,3,4 and Taher Awad4

1 International Institute for Urban Systems Engineering, Southeast University, Nanjing, China, 2 Department of Mechanical Engineering, California Polytechnic State University, San Luis Obispo, CA, United States, 3 Nanjing Zhixing Information Technology Company, Nanjing, China, 4 Department of Mechanical Engineering, Faculty of Engineering, Alexandria University, Alexandria, Egypt

#### Edited by:

Vagelis Plevris, OsloMet – Oslo Metropolitan University, Norway

#### Reviewed by:

Aristotelis E. Charalampakis, National Technical University of Athens, Greece Eleni N. Chatzi, ETH Zürich, Switzerland

> \*Correspondence: Mohammad Noori<EMAIL_ADDRESS>
#### Specialty section:

This article was submitted to Computational Methods in Structural Engineering, a section of the journal Frontiers in Built Environment

> Received: 16 October 2018 Accepted: 11 December 2018 Published: 07 January 2019

#### Citation:

Zhao Y, Noori M, Altabey WA and Awad T (2019) A Comparison of Three Different Methods for the Identification of Hysterically Degrading Structures Using BWBN Model. Front. Built Environ. 4:80. doi: 10.3389/fbuil.2018.00080 Structural control and health monitoring scheme play key roles not only in enhancing the safety and reliability of infrastructure systems when they are subjected to natural disasters, such as earthquakes, high winds, and sea waves, but it also optimally minimize the life cycle cost and maximize the whole performance through the full life cycle design. In this scheme, system identification is regarded as a major technique to identify system states and related parameter variables, thus preventing degradation of structural or mechanical systems when unexpected disturbances occur. In this paper, three different strategies are proposed to identify general hysteretic behavior of a typical shear structure subjected to external excitations. Different case studies are presented to analyze the dynamic responses of a time varying shear structural system with the early version of Bouc-Wen-Baber-Noori (BWBN) hysteresis model. By incorporating a "Gray Box" strategy utilizing an Intelligent Parameter Varying (IPV) and Artificial Neural Network (ANN) approach, a Genetic algorithm (GA), and a Transitional Markov Chain Monte Carlo (TMCMC) based Bayesian Updating framework system identification schemes are developed to identify the hysteretic behavior of the structural system. Hysteresis characteristics, computational accuracy, and algorithm efficiency are further discussed by evaluating the system identification results. Results show that IPV performs superior computational efficiency and system identification accuracy over GA and TMCMC approaches.

Keywords: hysteretic behavior, BWBN model, Intelligent Parameter Varying (IPV), Genetic algorithm (GA), Transitional Markov Chain Monte Carlo simulation (TMCMC), Bayesian updating

# INTRODUCTION

In recent years, an increasing attention is witnessed to face the challenging issues of safety, serviceability, reliability, risk and life-cycle management, and performance improvement of structures and infrastructure due to changing and more frequently occurring natural and man-made hazards, infrastructure crisis, and sustainability issues. These disturbances are dealt with innovative technologies to enhance structural functionality and safety in various stages of research and development (Spencer, 2003; Altabey, 2017a). Several types of structures that employ control strategies for different application scopes can be found in Constantinou et al. (1998), Soong and Spencer (2002), and Altabey (2014, 2017b,c); Altabey (2018). Proper modeling of inherent non-linearity in vast majority of structural systems plays an important role in understanding structural response under hazardous loadings. System identification is an important approach in control strategy regarded as the interface between the mathematical world of control theory and the real world of application and model abstractions (Zadeh, 1956; Ljung, 2010; Altabey, 2016, 2017d,e; Altabey and Noori, 2017a, 2018; Zhao et al., 2018), and it handles a wide range of system dynamics problem without the prior knowledge of actual system physics. The schematic diagram of system identification process is depicted in **Figure 1**.

Hysteresis can be described as the hereditary and memory nature of a non-linear or inelastic system behavior where the restoring force is dependent on both instantaneous as well as past history of deformations. In general, under cyclic loading, mechanical and structural systems are capable of dissipating considerable energy and they exhibit appreciable hysteretic behavior with hysteresis loops. Each loop enclosing the area in the restoring force vs. displacement curve depicts the energy dissipated over a complete cycle resulting from internal friction within the structural system.

Various empirical hysteresis models have been proposed in the past few decades. A class of smoothly varying hysteresis models used in engineering fields are Bouc-Wen class of hysteresis models. Bouc suggested a smooth and versatile hysteresis model for non-linear systems and hysteretic systems (Bouc, 1967; Wen, 1975, 1976, 1980, 1986, 1989; Park et al., 1986; Wen and Yeh, 1989; Ikhouane and Rodellar, 2007; Ikhouane et al., 2007; Ikhouane and Gomis-Bellmunt, 2008). Baber and Wen extended the Bouc model to take the degradation in strength or stiffness of structural systems into account (Baber and Wen, 1980). Baber-Noori, and later Noori, further extended the capabilities of Bouc-Wen model by including pinching behavior and studied the response of these systems under random excitation (Noori, 1984; Baber and Noori, 1986). Baber-Noori and subsequently Noori-Baber's work on integrating the pinching phenomenon in hysteretic behavior and extending Bouc-Wen-Baber (BWB) model was the first work in developing a smooth hysteresis model capable of taking into account strength and stiffness degradation as well as shear pinching phenomenon (Baber and Noori, 1985). BWBN was incorporated in structural design software, OpenSees developed at the University of California Berkeley (Hossain, 1995). A toolbox for computing the parameters of BWBN hysteresis model using multi-objective optimization evolutionary algorithms was also developed by SourceForge, an Open Source community (Bouc Wen Baber Noori Model of Hysteresis, Source Forge). Foliente showed Bouc-Wen-Baber-Noori (BWBN) model could produce previously observed inelastic behavior of wood joints and structural systems using BWBN smooth hysteresis model (Foliente, 1995; Zhao et al., 2017a,b; Noori et al., 2018). Deb et al. developed a toolbox that identifies structural parameters of Bouc–Wen–Baber–Noori hysteresis model through a noval multi-objective optimization evolutionary algorithms (MOBEAs) (Deb et al., 2002; Deb, 2013). Ortiz et al. analyzed and identified BWBN model via a multi-objective optimization algorithm (Ortiz et al., 2013). Peng et al. utilized BWBN model for identifying the parameters of a magneto-rheological damper and depicts its force-lag phenomenon (Peng et al., 2014). Muller et al. investigated the application of BWBN in their work and conducted performancebased seismic design through a Search-Based Cost Optimization (Muller et al., 2012). Chan et al. made a prediction of the hysteretic behavior of passive control systems by applying BWBN in a nonlinear-autoregressive-exogenous model (Chan et al., 2015).

Traditional artificial neural networks technique shows its superiority in the identification, monitoring, and control of complicated and non-linear dynamic systems (Narendra and Parthasarathy, 1990; Masri et al., 1992; Lu and Basar, 1998; Abouelwafa et al., 2014; Altabey and Noori, 2017b). However, a priori knowledge of the characteristics of restoring force is necessary and important for traditional parametric system identification approaches, while the non-parametric methods do not need information beforehand, lacking direct association between system dynamics and system model. In order to overcome the limitations of conventional parametric and nonparametric approaches, an noval Intelligent Parameter Varying (IPV) method was proposed, which makes full use of the embedded radial basis function networks to make an estimation of the hysteretic and inelastic characteristics of restoring forces constitutively for a multi degree of freedom system. A scaled three story base excited structure was designed to experimentally verify the non-linearity and its associated hysteresis of a structure using a displacement controlled shaking table (Saadat et al., 2003, 2004a,b, 2007). Further, a data-driven identification strategy for non-linear and hysteretic behaviors of steel wire strands was compared and verified using polynomial basis functions and neural networks. The results showed that neural networks were found more promising for the prediction of slightly pinched, hardening hysteresis, strongly pinched, hardening hysteresis, and classical quasi-linear softening hysteresis (Brewick et al., 2016). Genetic algorithms have been used for system identification of non-linear and hysteretic systems. The application of Real Coded Genetic Algorithms (RCGA) was demonstrated and applied to fit curves of synthetic and experimentally obtained Bouc-Wen hysteresis loops for a sandwich composite material (Hornig and Flowers, 2005). Different real coded genetic algorithms and their related criteria for efficiently identifying non-linear systems are regards as non-classical and optimized identification techniques (Monti et al., 2009). A Bayesian probabilistic framework was proposed to detect damage of continuous monitored structures by incorporating load-dependent Ritz vectors as an alternative to modal vectors (Sohn, 1998). A large body of work was conducted to track, estimate and identify structural parameters, system status and hysteretic and degrading behavior of structures using Kalman filters, extended Kalman filters and unscented Kalman filters (Jeen-Shang and Yigong, 1994; Yang et al., 2006; Wu and Smyth, 2007, 2008; Chatzi and Smyth, 2009; Chatzi et al., 2010; Lei and Jiang, 2011; Mu et al., 2013; Kontoroupi and Smyth, 2017; Erazo and Nagarajaiah, 2018). Traditional Markov Chain Monte Carlo approach in conjunction with Bayesian updating method were applied for structural response predictions and performance reliability evaluation (Yuen and Katafygiotis, 2001; Zhang and Cho, 2001; Beck and Au, 2002). Later, a transitional Markov Chain Monte Carlo (TMCMC) approach was developed by designing optimized sampling strategy from a series of intermediate probability density functions (PDFs) that converge to the target PDF, thus avoiding sampling difficulties. The TMCMC theory and algorithm were verified and demonstrated through the performance of the developed sampling approach, different PDFs as well as higher dimensional problems (Ching and Chen, 2007; Muto, 2007; Muto and Beck, 2008; Worden and Hensman, 2012; Zheng and Yu, 2013; Behmanesh and Moaveni, 2014; Green, 2015; Green et al., 2015; Ortiz et al., 2015).

It is a major barrier to successfully design hysteretic structures against degradation under severe cyclic loading. Most structural systems degrade with significant hysteresis, for example wood structures, dams, highways, reinforced concrete towers, steel bridges are critical and key elements of our built environment. In spite of their obvious importance, and their huge rehabilitation and replacement costs, design, construction, and analysis of the majority of these structures requires overly simplistic or in some cases flawed assumptions regarding hysteretic evolution. Development of a practical structural degrading identification approaches is much deserving. A comparative study of online and offline identification strategies for UAVs were discussed, and it is found that online approach is more adaptive to changes but with lower prediction accuracy (Puttige and Anavatti, 2007). Therefore, offline learning is employed in this paper. Based on what was discussed in the introduction above, the main contributions of the research are to present a three story hysteretically degrading shear structure by incorporating BWBN slip lock hysteresis to represent the hysteretic restoring forces in this system. This BWBN model will be capable of producing all significant and prominent features of structural strength and stiffness degradations as well as slip lock behavior, and conduct a comparative study using three system identification approaches including an Intelligent Parameter Varying Artificial Neural Network developed in an earlier research work by a group that involved one of the authors (a "gray box" model that considers linear as well as non-linear parts of the dynamic system), genetic algorithm optimization method, and a novel TMCMC statistical approach.

The comparative study of the aforementioned approaches for system identification and their application in a structural system using BWBN MODEL is an original work. To the best of the authors' knowledge such comparative study has not been reported in the literature.

# STRUCTURAL SYSTEM MODELING AND SYSTEM IDENTIFICATION

# BWBN Hysteresis Model

The model employed herein is an earlier version of BWBN hysteresis degradation model, which incorporates the previous smooth system degrading element by Bouc as modified by Baber and Wen in series with a slip-lock element (a nonlinear hardening spring) developed by Baber and Noori. Under cyclic excitation, degradation manifests itself in the evolution of progressively varying hysteresis loops. A non-linear system governed by Equation (1) is given with the incorporation of BWBN model.

$$m\ddot{\mathbf{x}} + c\dot{\mathbf{x}} + R = F\left(t\right) \tag{1}$$

$$R = \alpha kx + (1 - \alpha) \, kz \tag{2}$$

$$\dot{z} = \frac{A\dot{x}\_1 - \nu \left[\beta \left|\dot{x}\_1\right| \left|z\right|^{n-1} z + \nu \dot{x}\_1 \left|z\right|^n\right]}{n} \tag{3}$$

$$
\dot{x}\_2 = \sqrt{\frac{2}{\pi}} \frac{s}{\sigma} \exp\left[-\frac{z^2}{2\sigma^2}\right] \dot{z} \tag{4}
$$

$$
\dot{\kappa} = (1 - \alpha) \, k \dot{\mathbf{x}} z \tag{5}
$$

$$
\delta s = \delta\_s s \tag{6}
$$

$$\mathbf{x} = \mathbf{x}\_1 + \mathbf{x}\_2 \tag{7}$$

where, parameters m, c, k are, respectively, the mass, damping, and stiffness coefficients, and parameters x¨, x˙, x are quantities that describe the system acceleration, velocity and displacement, and R is the restoring force and F (t) is the ambient excitation. Parameter α is the weighting value denoting the ratio of postelastic to initial stiffness. Parameters A, β, and γ are basic hysteresis shape control parameters. Parameter z is the hysteretic displacement, and n is the degree of the sharpness of yield. Strength and stiffness degradation coefficients are, respectively denoted by ν(δν ) and η(δη). Parameters x1 and x2 are Bouc-Wen hysteretic system displacement and the additional displacement that considers slip-lock behavior. Parameter ε is the measure of the combined effect of duration and severity of the energy dissipated through hysteresis, σ is a measure of the sharpness of the peak of the hysteresis, and δs measures the slip magnitude. All 10 parameters are essential to produce the common features of hysteretic behavior. It would be very helpful if a small number of unspecified parameters for system identification can be reduced in that large numbers of parameters increase the uncertainty of convergence for updating parameters in search space. It was proved that the redundancy of specific hysteresis parameters can be eliminated through mathematical transformations in the parameter space devised to freeze them without affecting the system response (Ma et al., 2004; Charalampakis and Koumousis, 2008a,b; Charalampakis and Dimou, 2010).

#### Structural System Modeling

Restoring force curves of reinforced or steel structures show complex hysteresis characteristics, revealing material nonlinearity, crack opening and closing, bond and slip between steel bars and concrete and low cycle fatigue that result in structural strength and stiffness degradation. The hysteresis model used herein considers the specific case appropriate for both strength and stiffness degradations, and slip-lock behavior of restoring hysteretic forces of the structure. All the restoring forces of the structure are assumed to follow the BWBN hysteresis model. The structural system considered in this study is a shear type model subjected to ambient sinusoidal wave and ElCentro seismic excitation as input signals x¨g . The ground excitation motion x¨g makes an integral transformation to be incorporated into the structural equations. The structural system mainly includes three lumped mass coupled subsystems, as shown in **Figure 2**. These masses are lumped at floor (floor mi), levels and these floors are assumed and constrained to only move laterally. The restoring force Ri , in conjunction with the stiffness between the adjacent floors, are represented by dampers and springs, with the corresponding coefficients ci and ki , respectively. The time varying system yields energy dissipation due to the hysteretic behavior of the inner structure.

For this three story shear structure, the equations of motion are represented by Equations (8–10):

$$m\_3\ddot{\mathbf{x}}\_3 + c\_3\left(\dot{\mathbf{x}}\_3 - \dot{\mathbf{x}}\_2\right) + R\_3 = \mathbf{0} \tag{8}$$

$$m\_2\ddot{\mathbf{x}}\_2 + c\_2\left(\dot{\mathbf{x}}\_2 - \dot{\mathbf{x}}\_1\right) + R\_2 - c\_3\left(\dot{\mathbf{x}}\_3 - \dot{\mathbf{x}}\_2\right) - R\_3 = 0\tag{9}$$

$$m\_1 \ddot{\mathbf{x}}\_1 + c\_1 \left(\dot{\mathbf{x}}\_1 - \dot{\mathbf{x}}\_\emptyset\right) + R\_1 - c\_2 \left(\dot{\mathbf{x}}\_2 - \dot{\mathbf{x}}\_1\right) - R\_2 = 0\tag{10}$$

#### System Identification Theory

A typical system identification framework mainly has two components including system itself and system identification model. By defining an equivalence criterion, the parameters of the system identification model are updated via a comparison with the original system, until the system identification model will be eventually equivalent to the original system. System identification approaches used in this paper include intelligent parameter varying based approach, GA based approach and transitional markov chain monte carlo (TMCMC) based approach. IPV approach employs ANNs and establishes a "Gray Box" system, where the original system parameters are replaced by the parameters of ANNs, and the system structure is replaced by ANNs. For GA and TMCMC methods, they are used to optimize and identify the parameters of system model, which is regarded as an approximate model of the original real system. The

theoretical background for the three approaches are presented in the following subsections.

#### Intelligent Parameter Varying Based System Identification

Artificial neural network (ANN) is a non-linear and adaptive information processing system composed of large numbers of neuron units. A radial basis function (RBF) neural network is employed to build the intelligent parameter varying model. In ANN architecture, the Euclidean distance between the clustering center and input vector is calculated and the result is activated to pass through the output layer. The activation function is Gaussian output function and is formulated as:

$$g\_{j} = \exp\left(-\frac{\sum\_{i=1}^{n} \left(\chi\_{i} - \varepsilon\_{j}\right)^{2}}{\sigma\_{j}^{2}}\right) \tag{11}$$

where gj is the output of the j th unit in the hidden layer, xi is the input data fed to the network, cj is the center of the j th unit in the input space, and σj is the width of the j th function. j = 1, 2, · · · , m. Parameter m is the number of the centers of neurons, and n represents the dimension of the input space. Linearly weighted summation of hidden layer node outputs produces the output nodes. Therefore, the output of the network is calculated by:

$$\mathcal{Y} = \sum\_{j=1}^{n} \omega\_{j} \mathbf{g}\_{j} \tag{12}$$

where wj is the weight of the j th node.

The objective is to find a series of weights that minimize the square of error between the actual and desired network outputs, i.e.,:

$$E\left(k\right) = \frac{1}{2} \sum\_{j=1}^{N} \left(d\left(k\right) - \wp\left(k\right)\right)^2\tag{13}$$

where d(k) is the desired output and y(k) the actual output of RBF and j = 1, 2, · · · , N, and N is the number of data sets. There are commonly two approaches for utilizing ANNs, i.e., supervised and unsupervised learning. In supervised approach the input and the referred output are usually known. The network then processes the input values and makes a comparison between its resulting output and the desired output. The network system makes the errors propagate back through different layers, causing the system to adjust tuned weights which control the network. This process repeats over and over until all the weights are continually tweaked in an appropriate way. During the training process of a network the connection weights are continually refined to a specific generalization level and a good network performance level. In unsupervised training, the system itself must then decide what features and how many features need to be extracted to group the input data, which is often referred to as self-organization or adaption. Herein, we use the non-supervised learning algorithm to acquire the centers and variance of radial basis function, and meanwhile the least mean squared error is acquired by using supervised learning algorithm. The response of the hysteretic system is used as the desired signal, and the error between the desired signal and the simulated signal is back propagated to modify the weights and threshold values of neural network model.

Parametric system identification approaches have been widely used but in most published literature in this area a priori knowledge of the characteristics depicting the behavior of restoring force is required. Non-parametric approaches generally do not need information beforehand but they typically lack direct associations between system dynamics and associated model. When ANNs are implemented using the "Black Box" approach, little of the system information might be obtained from the traditional techniques due to the fact that the "Black Box" only considers system input and output. Intelligent Parameter Varying (IPV) method preserves the benefits of both traditional parametric and non-parametric approaches, and utilizes the embedded radial basis function as the activation of neurons to estimate the constitutive characteristics of inelastic and hysteretic restoring forces for a multi degree of freedom structural system.

IPV technique, i.e., a gray box approach (**Figure 3**) that incorporates the advantages of both "White Box" and "Black Box" approaches, was developed in such a way that the model structure can be determined using the first principle (Equations 18–20), while non-linear and adaptive learning capabilities of ANNs can be used to identify the non-linear, time varying system's dynamics (Equations 24–26) that would be difficult to model and identify using the traditional "White Box" and "Black Box" (Saadat et al., 2003, 2004a,b, 2007).

A non-linear system with full state measurement represented by the Linear Parameter Varying (LPV) model structure is given by:

$$\dot{\boldsymbol{x}} = f\_1 \left( \mathbf{x}, \boldsymbol{u} \right) \cdot \boldsymbol{x} + f\_2 \left( \mathbf{x}, \boldsymbol{u} \right) \cdot \boldsymbol{u} \tag{14}$$

$$y = x \tag{15}$$

The IPV approach introduced herein would preserve the model structure inherent in Equation (14) without requiring a priori representations of non-linearities f1 (x, u) and f2 (x, u). Instead, these terms would be represented by separate artificial neural networks g1 (x, u,w1) and g2 (x, u,w2) as depicted in Equations (16, 17):

$$\dot{\mathbf{x}} = \operatorname{g}\_1(\mathbf{x}, \boldsymbol{\mu}, \boldsymbol{\nu}\_1) \cdot \boldsymbol{\varkappa} + \operatorname{g}\_2(\mathbf{x}, \boldsymbol{\mu}, \boldsymbol{\nu}\_2) \cdot \boldsymbol{\mu} \tag{16}$$

$$\mathbf{y} = \mathbf{x} \tag{17}$$

By modeling the non-linearities f1 (x, u) and f2 (x, u) via separate artificial neural networks g1 (x, u,w1) and g2 (x, u,w2), the model structure (Equations 16, 17) is preserved. Therefore, the relationship between the model structure and artificial neural network parameters is preserved. The structural model is preserved by incorporating ANN, preserving a portion of information of the structural model. The IPV approach preserves the direct association between the construction of ANN and the system dynamics, used for structural health monitoring for system identification.

Based on first principle, system dynamics can be transformed to the following form:

$$-\left.c\_3\left(\dot{\mathbf{x}}\_3 - \dot{\mathbf{x}}\_2\right) - R\_3 - m\_3\ddot{\mathbf{x}}\_{\mathcal{G}} = m\_3\ddot{\mu}\_3\tag{18}$$

$$-c\_2\left(\dot{\mathbf{x}}\_2 - \dot{\mathbf{x}}\_1\right) - R\_2 + c\_3\left(\dot{\mathbf{x}}\_3 - \dot{\mathbf{x}}\_2\right) + R\_3 - m\_2\ddot{\mathbf{x}}\_\% = m\_2\ddot{u}\_2 \tag{19}$$

$$-\left(\dot{\varkappa}\_1 \left(\dot{\varkappa}\_1 - \dot{\varkappa}\_\emptyset\right) - R\_1 + \varepsilon\_2 \left(\dot{\varkappa}\_2 - \dot{\varkappa}\_1\right) + R\_2 - m\_1 \ddot{\varkappa}\_\emptyset = m\_1 \ddot{u}\_1 \tag{20}$$

where u1, u2, u3, respectively represent the relative displacements of each floor, i.e., u1 = x1 − xg , u2 = x2 − xg , u3 = x3 − xg .

$$-R\_3 - m\_3 \ddot{\mathbf{x}}\_{\mathcal{G}} = m\_3 \ddot{\boldsymbol{u}}\_3 \tag{21}$$

$$-R\_2 + R\_3 - m\_2\ddot{\mathbf{x}}\_{\mathcal{g}} = m\_2\ddot{\mathbf{u}}\_2\tag{22}$$

$$-R\_1 + R\_2 - m\_1 \ddot{\mathbf{x}}\_{\mathcal{K}} = m\_1 \ddot{\boldsymbol{\mu}}\_1 \tag{23}$$

The stiffness and damping terms are lumped into restoring forces, since in the hysteresis models the restoring force R associates with the lateral relative displacement xr and the restoring displacement z, where z is expressed by the function of lateral relative velocity x˙r .

The modeling for the restoring forces using radial basis function based neural network is as follows:

$$
\hat{R}\_3 = \mathfrak{g}\_3 \left( \ddot{\mu}\_3, \ddot{\mu}\_{\mathfrak{g}} \right) \tag{24}
$$

$$
\hat{R}\_2 = g\_2 \left( \ddot{u}\_2, \ddot{x}\_\emptyset, \hat{R}\_3 \right) \tag{25}
$$

$$
\hat{R}\_1 = \mathcal{g}\_1 \left( \ddot{u}\_1, \ddot{x}\_\mathcal{g}, \hat{R}\_2 \right) \tag{26}
$$

where, Rˆ 1, Rˆ 2, Rˆ 3 represent, respectively, the identified restoring forces through the training of ANN as shown in **Figure 4**.

#### Genetic Algorithm Based System Identification

Genetic algorithm (GA) is an approach that searches the global optimal solution by simulating a natural evolution process. The objective function can be formulated as the normalized mean square error (MSE) of the predicted time history y˜ t|p as compared to the reference time history y (t). Herein, the acceleration response signal is employed as the time history series. The purpose of the following optimization approach is to minimize the difference (or the error) between predicted time history and the referred time history. The hysteretic structural system objective function is introduced below:

$$OF\left(p\right) = \frac{\sum\_{i=1}^{n} \left(\wp\left(t\_i\right) - \tilde{\wp}\left(t\_i|p\right)\right)^2}{N\sigma\_\mathcal{\mathcal{V}}^2} \tag{27}$$

where p is a parameter vector, σ 2 y the variance of the reference history, and N the number of points used. Sum of three acceleration response signal differences of the hysteretic structural system is used as the objective function. The optimization problem can be stated as the minimization of the objective function OF p when the parameter vector has the following side constraints as:

$$
\varkappa\_{LB} \le p \le \varkappa\_{UB} \tag{28}
$$

where xLB and xUB are vectors defining the lower and the upper values of the model parameters, respectively. The basic strategy for the parameter identification using GA is shown in **Figure 5**. GA operates starting from a population of the potential solutions to a representive problem, and one population is composed of numbers of individuals coded by genes. Each individual is chromosomes with the characteristics of entity. GA initializes on a population of individuals (coded candidate solutions to the problem) that are manipulated by some operators such as selection, crossover, and mutation. In short, the selection process drives the search direction toward the region of best individuals, and the cross operator combines individuals to generate offsprings.

If it is indispensable to make selection and crossover operators converge toward the optimum, mutation alters one or more gene values (individuals) in a chromosome from its initial state, thus, maintaining genetic diversity from one generation of a population to the next. In this way, a complete exploration of global search space is forced by algorithm within the search space. Each individual in the population is represented as chromosome, indicating the collection of parameters are supposed to identify. GA adopts an elitist strategy, which consists of the preservation of the most fit individuals obtained in the current generation. Population representation and initialization generate population and individuals, and the initialized value will be assigned to the parameter space for solving the hysteretic system model. The fitness function is established by using the prediction error between the simulated signal response and predicted signal response. After a series of successive mathematical operators for optimization, the next generational loop begins.

#### Transitional Markov Chain Monte Carlo Based System Identification

Markov Chain Monte Carlo (MCMC) is an analytic approach which replaces numerical integration through summation over numbers of samples generated from iteration. A Markov chain is a stochastic process where one state is transformed to another state after a sufficiently long sequence of transition procedure. The next state is conditionally based on the last state. A key property of Markov chain is that the starting state has no influence on the state of the chain via a series of sequential transitions. The chain reaches its steady state at a specific point

where it reflects sampling distribution from stationary status. The principle of Monte Carlo simulation is applied for the integration to approximate the expected complex distribution status of numbers of samples. By increasing the number of samples, the approximation accuracy can be measured and achieve a desired value, which mainly depends on the independence of the samples.

Markov chain involves a stochastic sequential process where a series of states can be sampled from some stationary distributions while Monte Carlo sampling can make an estimation of various characteristics of a specific distribution. The goal of MCMC is to design a Markov chain to meet the target distribution of the chain which is what we are interested in sampling from.

Transitional Markov Chain Monte Carlo Theory (TMCMC) is introduced to avoid the difficulty of sampling from complicated target probability distributions (e.g., multimodal PDFs, PDFs with flat manifold, and very peaked PDFs) but sampling from a series of intermediate PDFs that converge to the target PDF and are easier to sample.

Bayesian Inference describes a process of solving posterior density functions given the likelihood and prior probability. The target probabilistic model can be depicted by M, D is the data acquired from the system, and the uncertain parameters of the model are described as θ. Sampling from the posterior PDF of θ conditioned on D is the aim of the Baysian model updating, which is given as:

$$\begin{split} \, \_f f \left( \theta | M, D \right) &= \frac{\, \_f \left( D | M, \theta \right) \cdot f \left( \theta | M \right)}{\int \left( D | M \right)} \\ &= \frac{\, \_f f \left( D | M, \theta \right) \cdot f \left( \theta | M \right)}{\int \, \_f f \left( D | M, \theta \right) \cdot f \left( \theta | M \right) \cdot d \theta} \end{split} \tag{29}$$

where f (θ|M) is the prior PDF of θ, f (D|M, θ) is the likelihood of D given θ, and f (D|M) is the evidence of the model M.

Bayesian model updating generally employs simulation based methods in that it is effective to obtain samples from f (θ|M, θ), which can be estimated at a specific quantity of interest E g|M, D based on the Law of Large Number.

$$E\left(\mathbf{g}|M,D\right) \approx \frac{1}{N} \sum\_{k=1}^{N} \mathbf{g}\left(\theta\_k\right) \tag{30}$$

where θk : k = 1, 2, · · · , N represents a set of N samples from f (θ|M, D). Consider the equation as follow:

$$\int \left(\theta \vert M, D\right) \propto f\left(\theta \vert M\right) \cdot f\left(D\vert M, \theta\right) \tag{31}$$

It is usually difficult to sample from f (θ|M, D) using Importance Sampling (IS) and Metropolis–Hastings (MH) in that it is not so easy to understand the geometry of the likelihood f (D|M, θ). To converge to the target PDF f (θ|M, D) from the prior PDF f (θ|M), a series of intermediate PDFs are constructed as the following:

$$f\_j(\theta) \propto f\left(\theta|M\right) \cdot f\left(D|M,\ \theta\right)^{\mathfrak{P}j}, \ j = 0, 1, \cdots, m$$

$$0 = p\_0 < p\_1 < \cdots < p\_m = 1 \tag{32}$$

$$\text{Note that}\\
f\_0\left(\theta\right) = f\left(\theta|M\right),\\
f\_m\left(\theta\right) = f\left(\theta|M,D\right).$$

Frontiers in Built Environment | www.frontiersin.org

TABLE 1 |

Parameter assignment

 for different cases.

In GA approach, the population size is assigned 100, the elite count is 10, and crossover fraction is 0.6, and the generation is 30.

In TMCMC approach, the sampling points are assigned 5000, and the scaling factor is assigned 0.2.

where the index j denotes the stage number. Although the geometry changing from f (θ|M) to f (θ|M, D) is large, the status of two changing adjacent intermediate PDFs is small. It is efficient to sample from fj+1 (θ) according to the previous sample from f j (θ) through this small change.

The function fj (θ) is used to extract samples and make an estimation of the PDF itself as a kernel density function (KDF), a combination of weighted Gaussian functions centered at the samples. The kernel density function can be regarded as the proposal PDF of the MH method to sample from f j+1 (θ). This will subsequently and ultimately result in f (θ|M, D) samples. This approach is called adaptive Metropolis–Hastings (AMH) algorithm. Given that the proposal PDF (KDF) function is fixed, rendering the MH method is as similar as IS, not efficient in high dimension situation.

It is a totally different strategy for TMCMC to acquire fj+1 (θ) samples based on fj (θ) samples, KDF method is replaced by a resampling algorithm. It covers a battery of resampling stages, with each stage completing the following, given Nj samples from fj (θ), depicted by θj, k : k = 1, · · · , Nj , acquire samples from f j+1 (θ), depicted by θj+1, k : k = 1, · · · , Nj+1 . It can be calculated by the following in a more easier way, with the samples θj, k : k = 1, · · · , Nj from fj (θ). The "plausibility weights (w θj, k )" of these samples regarding fj+1 (θ) can be computed by:

$$\mathbb{W}\left(\theta\_{\hat{\mathbb{M}},k}\right) = \frac{f\left(\theta\_{\hat{\mathbb{M}},k}|M\right)f\left(D|M,\ \theta\_{\hat{\mathbb{M}},k}\right)^{p\_{\hat{\mathbb{M}}}}}{f\left(\theta\_{\hat{\mathbb{M}},k}|M\right)f\left(D|M,\ \theta\_{\hat{\mathbb{M}},k}\right)^{p\_{\hat{\mathbb{M}}}}} = f\left(D|M,\ \theta\_{\hat{\mathbb{M}},k}\right)^{p\_{\hat{\mathbb{M}}+1-p\_{\hat{\mathbb{M}}}}},$$

$$k = 1,\cdots,\ \ N\_{\hat{\mathbb{M}}}\tag{33}$$

Based on the normalized weights, the uncertain parameters can be resampled, i.e. let: θj+1, k = θj, l and w.p. w θj, l PNj l=1 w θj, l k = 1, · · · , Nj+1

where "with probability" is represented by w.p. and pacifier index is denoted as l. It is shown that if Nj and Nj+1 achieve a relatively large quantity, θj+1, k : k = 1, · · · , Nj+1 will be distributed as f j+1 (θ). Moreover, w θj, k is expected as the following value:

$$\begin{aligned} E\left[\boldsymbol{\w}\left(\theta\_{\boldsymbol{\flat},k}\right)\right] &= \int \boldsymbol{\nu}\left(\theta\right) \cdot \boldsymbol{f}\_{\boldsymbol{\flat}}\left(\theta\right) \cdot d\theta \\ &= \int \boldsymbol{f}\left(\boldsymbol{D}|\boldsymbol{M},\theta\right)^{p\_{j+1}-p\_{j}} \cdot \boldsymbol{f}\_{\boldsymbol{\flat}}\left(\theta\right) \cdot d\theta \end{aligned}$$

$$\begin{split} &= \int f\left(D|M,\,\theta\right)^{p\_{j+1}-p\_j} \cdot \frac{f\left(\theta|M\right)f\left(D|M,\,\theta\right)^{p\_j}}{f\left(\theta|M\right)f\left(D|M,\,\theta\right)^{p\_j}d\theta} \cdot d\theta \\ &= \frac{\int f\left(\theta|M\right)f\left(D|M,\,\theta\right)^{p\_{j+1}}d\theta}{\int f\left(\theta|M\right)f\left(D|M,\,\theta\right)^{p\_j}d\theta} \end{split} \tag{34}$$

Therefore, PNj k=1 w θj, k /Nj is the automatically unbiased estimation made for R f(θ|M)f(D|M, θ) pj+1 dθ R f(θ|M)f(D|M, θ) pj dθ . According to the results above, the following method is used to sample from f (θ|M, D) and make an estimation of f (D|M).

More precisely, with probability w θj, k / PNj l=1 w θj, l , by using a covariance matrix equal to the scaled version of the estimated covariance matrix of f j+1 (θ), a Markov chain sample

Frontiers in Built Environment | www.frontiersin.org

in the k th chain can be generated from a Gaussian proposal PDF centered at the current sample of the k th chain:

$$\sum\_{j} = \beta^{2} \sum\_{k=1}^{N\_{j}} \le \left\{ \theta\_{j,k} \right\} \left\{ \theta\_{j,k} - \left[ \frac{\sum\_{j=1}^{N\_{j}} \le \left(\theta\_{j,l}\right) \theta\_{j,l}}{\sum\_{j=1}^{N\_{j}} \le \left(\theta\_{j,l}\right)} \right] \right\}$$

$$\times \left\{ \theta\_{j,k} - \left[ \frac{\sum\_{j=1}^{N\_{j}} \le \left(\theta\_{j,l}\right) \theta\_{j,l}}{\sum\_{j=1}^{N\_{j}} \le \left(\theta\_{j,l}\right)} \right] \right\}^{T} \quad \text{(35)}$$

where β is the prescribed scaling factor, P j = product of β 2 and the estimated covariance of f j+1 (θ). The rejection rate is chosen as the value β, and MCMC may probably achieve a larger value accordingly. The value 0.2 is found to be more reasonable for the scaling parameter β.

It is essential to choose pj : j = 1, · · · , m − 1 . The larger value of p is desirable to make the transition between the intermediate adjacent PDFs smoother. The number of intermediate stages achieves a huge value if the increase of p-values has slow change rates.

The degree of uniformity of the plausibility weights w θj, k : k = 1, · · · , Nj can appropriately indicate how close fj+1 (θ) approaches fj (θ), so pj+1 should be chosen so that the coefficient of variation (COV) of the plausibility weights can be equivalent to a prescribed threshold. The Bayesian inference framework for system identification is established (**Figure 6**) for structural model updating, which is regarded as a determinant reason for choosing the most suitable model parameters related to the hysteretic behavior of the structures by minimizing the difference between the predicted structural response and the simulated structural response. For the hysteretic structural model, the uncertain model parameters are selected as the ones that need to be updated through Bayesian inference (Equation 31) by drawing samples of parameters from the posterior PDF of parameters.

The TMCMC based Bayesian Updating algorithm is coded and implemented to establish the computational environment via exchange of data between the model and the algorithm. The evolution of parameter updating process is that, the samples from the prior PDFs are approximately uniformly distributed in the model parameter space at the first stage p0 = 0 . Through applying Bayesian inference with TMCMC probabilistic simulations, the samples eventually populate well in the high probability region of the posterior PDFs close to the true model parameters at the last stage pm = 1 .

TMCMC approach is employed to make an identification of the parameters of Bouc-Wen class models, henceforth represented by the vector θ ≤ 2 ⊆ R d . The appropriate choice of θ reflects the corresponding non-linear and hysteretic behavior of structure. By applying Bayesian model updating, the major advantage is that the result gives a probability distribution expressing the likelihood probability distribution of different parameters rather than yielding a single value for θ. It is clear that the evolution of the model parameter variation represents the Bayesian inference process.

To make a restriction of the parameter space 2, two side constraint vectors θmin and θmax are defined such that:

$$
\theta\_{\min} \left( i \right) \le \theta \left( i \right) \le \theta\_{\max} \left( i \right) \qquad \qquad 1 \le i \le d \tag{36}
$$

The vector θ has specific constraints such that the generated initial samples can determine the feasible values that the parameters can take. By defining prior PDF of θ a uniform distribution between the likelihood PDF and side constraints are regarded as the prediction error assumed as Gaussian function distributed with unknown variance and zero mean. The prediction error is described as the error between the predicted system response and the simulated system response given by:

$$f\left(D|\theta\right) = \prod\_{i=1}^{l} \frac{1}{\sigma\_{acc}\sqrt{2\pi}} \exp\left[-\frac{1}{2\sigma\_{acc}^{2}} \left(\frac{\mathbf{x}\left(t\_{i}\right) - \hat{\mathbf{x}}\left(t\_{i}|\theta\right)}{S\_{acc}\left(t\_{i}\right)}\right)^{2}\right] \text{(37)}$$

where σ 2 acc is the variance of the prediction errors and Sacc is the weighting function used to normalize the acceleration response of the hysteretic system. To achieve computational convenience, the log-likelihood function lnf (D|θ) is employed in the actual implementation of the TMCMC algorithm as:

$$\ln f\left(D|\theta\right) = -\frac{1}{2} N\_t \ln\left(2\pi\right) - N\_t \ln\left(\sigma\_{\text{acc}}\right)$$

$$-\frac{1}{2\sigma\_{\text{acc}}^2} \left(\frac{\varkappa\left(t\_i\right) - \hat{\varkappa}\left(t\_i|\theta\right)}{S\_{\text{acc}}\left(t\_i\right)}\right)^2\tag{38}$$

The j th stage of parameter evolution process by correspondingly choosing the values of pj can be shown as the contours of PDF f j (θ).

# NUMERICAL ANALYSIS OF SYSTEM IDENTIFICATION

#### Parameter Settings

**Table 1** lists the parameter assignments of the structural system associated with BWBN hysteresis model for different cases. The mass coefficient of each floor is 1 kg, stiffness coefficient 10 N/m (sin) and 20 N/m (ElCentro), damping coefficient 0.05 N/(m/s) for each floor, respectively.

For hysteresis model, A = 1, α = 0.01, β = 2, γ = 2, n = 1 (IPV), n = 1.5 (GA and TMCMC), δν = 0, δν1 = 1.5 (sin) and 3 (ElCentro), δη = 0, δη1 = 1.5 (sin) and 3 (ElCentro), σ = 0.01, δs = 0.05. Damage occurrences are assumed at 30 s, 60 s, 100 s (sin) and 10 s, 20 s, 40 s (ElCentro) when both the strength degradation factor δυ and stiffness degradation factor δη change to δυ1and δη1, respectively. For GA and TMCMC approaches, parameter initial values are assigned between the corresponding lower and upper bounds before the optimization processes. Parameters are updated between the two bounds of parametric searching space, and eventurally achieve the exact values.

#### Numerical Analysis Results

**Figures 7**, **8** show the restoring force identification results when the structure is subjected to sinusoidal signal and

The running time of system identification programs is computed using Samsung computer [Win 7 Ultimate, Processor: Intel (R), Xeon (R), CPU E3-1231 v3 @3.4G Hz, RAM: 8G, System type: 64 operating system].

ElCentro signal, respectively. Black curve represents the restoring force of the original structural system while the red curve represents the identified restoring force using IPV, GA, and TMCMC approaches. For the case of sinusoidal excitation, the identification results show that the degradation phenomenon of restoring forces of 2nd−3rd floor and 1st−2nd floor are more evident compared with ground-1st floor by using all of the three approaches. For the case of ElCentro excitation, the identification results show that the degradation phenomenon of restoring forces existing in all adjacent floors by using IPV and GA methods is more evident compared with that using TMCMC method. The hysteresis curves rotate accordingly when damage occurs. More results in details are discussed in the next section.

### DISCUSSIONS

# Choice of Signal Excitation Type and Objective Function

Both sinusoidal excitation signal and ElCentro excitation signal are used as input signals to the system identification of hyestetic system. The frequency of sin signal is single and history curve is smooth and periodic variation. The ElCentro wave is stochastic, non-periodic, and non-steady state, and this represents the other type of excitation which is totally counter to harmonic wave such as sin signal. These two types of signals are employed to the response analysis and verification of the generalization capability of three different algorithms for identifying restoring forces of hysteretic structure. In addition, the establishment of objective functions employed in the GA and TMCMC are formulated through the combination of desired and simulated acceleration signal differences of all three degrees of freedom.

#### Parameter Choice for System Identification

**Table 2** shows the parameter identification error using GA and TMCMC for the case of noise free, SNR = 30 and SNR = 10 (SNR = signal noise ratio). These error represents the relative error between the "true value" and "identified value." All errors are no more than 10%. Generally for the cases with higher SNR they perform better parameter identification results since the disturbance of noise to signal is low. From **Figures 7**, **8**, and **Table 2**, for the noise free case, it is shown that during the degradation process (damage evolution), strength degradation factor δυ changes from 0 to 1.4695 (sin,GA) /3.0130 (ElCentro,GA) /1.5913 (sin,TMCMC) /2.9890 (ElCentro,TMCMC), and stiffness degradation factor δη changes from 0 to 1.5370 (sin,GA)/2.9142 (ElCentro,GA)/1.5400 (sin,TMCMC)/2.7357 (ElCentro,TMCMC).

The hysteresis loop rotates clockwise with a certain degree, indicating large non-linearity and considerable degradation. However, from the identification results, it is also found that in the sin excitation case, the hysteresis loop exhibits better slip-lock phenomenon, and can absorb more energy than the case of ElCentro excitation. Noise corrupted cases have similar parameter identification results.

**Table 3** shows the analysis of correlation coefficients for different cases. For the cases with small SNR, the identification effectiveness is not relatively good compared with noise free and

TABLE

4

Degradation

 parameter assignment

 for different cases.

higher SNR (R-t is the restoring force-time relationship, and R-x is the restoring force-relative displacement relationship).

In this paper, IPV is a radius basis neural network which is a completely data-driven approach, while GA and TMCMC are not fully data-driven in that they assign the identification model within specific initial parameter intervals and then update and search to optimize structural model.

IPV method performs better identification results due to its adaptive learning capability and anti noise property in the noise environment. From computational time, it can be concluded that the IPV approach is a more efficient identification method compared with GA and TMCMC methods. Concretely, for the case of noise free, the computational time for GA and TMCMC are 9.59 and 760.51% longer than IPV (sin), and 8.77 and 694.40% longer than IPV (ElCentro). Similar results are also shown in noise corrupted cases. This demonstrates IPV method has higher computational efficiency than GA and TMCMC approaches.

Strength and stiffness degradation parameters are both very important in the hysteresis model, which determine the hysteretic behavior of structural systems. **Table 4** lists five groups of parameter assignments regarding the variance of strength and stiffness degradation coefficients. The objective of setting these cases is studying the influence of change of degradation parameters on system identification accuracy.

To simulate damage occurrence, for IPV approach, strength degradation factor δυ and stiffness degradation factor δη change from 1.5 to 1.8, 2.1, 2.4, and 2.7, respectively (sin), while change from 3 to 3.3, 3.6, 3.9, and 4.2, respectively (ElCentro). For GA (sin) method, for the first three groups, the true value of strength degradation factor δυ and stiffness degradation factor δη keep fixed value 1.5 but the lower and upper bound change from 0.1–3 to 0.5–2.5 and 1–2, respectively. For the last two groups, the true value of strength degradation factor δυ and stiffness degradation factor δη are 1.8 and 2.1, respectively but the lower and upper bound correspondingly change from 0.1 to 3.6 and 0.1 to 4.2, respectively.

For GA (ElCentro) method, for the first three groups, the true value of strength degradation factor δυ and stiffness degradation factor δη keep fixed value 3 but the lower and upper bound