Parametric study on a Bouc-Wen model with degradation features for the study of cyclic behavior of a reinforced concrete frame

1.   Introduction

Nonlinearities of hysteresis are present in several phenomena of science and technology, such as physics (e.g. plasticity, friction, ferromagnetism, ferroelectricity, superconductivity, adsorption, and desorption, etc.), chemistry, biology, mechanics, even economics and experimental psychology, but particularly in the present scope, in civil engineering ^[1,2].

According to Visintin ^[2], hysteresis is regarded as a rate-independent memory effect. Rate independence requires that the couple (input and output) is invariant to any increasing time homeomorphism, meaning that at any instant of time, the output only depends on the range of the restriction of the input and on the order in which values have been attained, being independent of the derivatives of the input, which may even fail to exist. Hysteresis is then a phenomenological concept that is based on experimental data relating to the aforementioned couple.

In reality, the response of structural systems depends always on the rate of the applied load, and when this needs to be considered, viscosity will be introduced. Although generally in seismic response of structures viscosity and hysteresis contribute to the overall energy dissipated by the structure, if the load is applied at a very low rate, the viscosity phenomena can be neglected ^[3].

Hysteresis models represent constitutive relationships of a structural system or element. These represent a relationship between a couple of two physical quantities, that in the present scope can be resistance-deformation or stress-strain, shear-distortion, force-displacement, and moment-curvature. More realistic hysteresis models are essential to accurately characterize an inelastic response analysis of a structural system subjected to dynamic loads, such as seismic excitations.

Several hysteresis models were created and studied over the years. These models were developed based on how a certain structural system composed of a certain material behaves under lateral cyclic loading. The characterization of the mechanical non-linear behavior of structural components, e.g., structural elements and connections, and the assessment of the seismic response of structural systems (moment frames, braced frames, and shear and masonry walls) constitute the main purpose in the development of these hysteresis models ^[4].

The non-linear or inelastic behavior that a structural system exhibits when subjected to a dynamic excitation, can be simulated using a special kind of hysteresis model, an empirical hysteresis model. This kind of model, as the name suggests, is based on a phenomenological approach. The hysteresis characteristics observed in experimental tests of a certain structural system are idealized in a hysteresis model capable of expressing the resistance-deformation relations under any loading history with load reversals ^[5].

The early models in this context were mainly developed based on experimental observations of reinforced concrete (RC) and steel structural members' behavior. The later and recent hysteresis models are focused on more versatile models capable of simulating different structural systems and different materials under diverse loading conditions. Recent works have presented model formulations capable of accurately simulating complex hysteretic phenomena in rate-independent mechanical systems and materials ^[6,7], and the complex hysteretic response of elastomeric isolation bearings for seismic isolation ^[8]. Other studies have presented hysteretic models capable of reproducing damage of structural systems and components presenting stiffness and strength degradation ^[9,10].

Depending on the loading conditions of a structural member, its behavior can be dominated by flexure, shear, or axial deformation. The interaction between these behaviors can be complex and must be taken into account to obtain more reliable simulations.

These models can be generally divided into two types, polygonal hysteresis models (PHMs) and smooth hysteresis models (SHMs) ^[11].

The PHMs are based on piecewise linear behavior and driven by real behavioral stages of an element or structure, such as initial or elastic, cracking, yielding, stiffness and strength degradation stages, and crack and gap closures ^[11]. Many PHMs were developed over the years mostly to simulate the flexure-dominance behavior of structural members (e.g. elastic-perfectly-plastic bilinear model, strength-hardening bilinear model, Clough's model ^[12], degrading bilinear model ^[13], modified Clough's model ^[14], trilinear degrading model ^[15], Takeda's model ^[16,17]. Additionally, to account for shear interaction other models have been developed simulating the pinching effect (e.g. slip-type Takeda's model ^[18], Kabeyasawa et al. model ^[19], Three-parameter model ^[20], Costa and Costa model ^[21], Rodrigues et al. ^[22].

Smooth hysteresis models present continuous stiffness changes due to yielding but with sharp changes due to unloading and deteriorating behavior ^[11]. Different SHMs were developed over the years, although the most popular ones are the Bouc-Wen models, which are a set of models including the original and its further modifications. With proper mathematical manipulations, this constitutive model can represent many different empirical behaviors in a phenomenological way, i.e., softening/hardening and smoothly varying hysteresis curves, e.g., stiffness and strength degradation and the pinching effect.

Bouc-Wen models' can be considered as "semi-physical" models since they come from a combination of some physical understanding of the hysteresis system along with some kind of black-box modeling rather than a detailed analysis of the physical behavior of the systems with hysteresis ^[23]. This model essentially consists of a first-order non-linear differential equation relating an input displacement and an output restoring force in a hysteresis form. With the proper choice/tuning of parameters, it is possible to replicate the response of the model regarding real hysteresis behavior.

One example of the set of Bouc-Wen models is the Sivaselvan and Reinhorn smooth hysteresis model ^[11,24]. This model was developed based on previous models, viz., Bouc ^[25], Wen ^[26,27], Baber and Noori ^[28], Foliente ^[29], and Reinhorn et al. ^[30]. It is a versatile model that carries a physical interpretation that is appropriate to model the non-linear behavior of structural systems and structural components.

In this paper a parametric variation of a smooth hysteresis model based on the works of Sivaselvan and Reinhorn's ^[11,24] is performed, to evaluate the role of each parameter for further calibration, attempting to emulate the global cyclic behavior of an experimental RC frame structure using a macro model developed in MATLAB ^[31].

2.   Parametric study

2.1.   Characteristics of the experimental frame structure

To study the cyclic response of RC frame structures with or without infill masonry walls an experimental campaign was performed at Laboratório Nacional de Engenharia Civil (LNEC) ^[32].

The RC bare frame structure used in this experimental study is represented in Figure 1a. A constant vertical load of 100 kN is applied at each column, and the structure is submitted to an increasing cyclical load/displacement pattern, as represented in Figure 1b, at the beam's center level.

The materials used in the experimental frame structure were concrete of the class C20/25, longitudinal steel reinforcement of S400, and transverse reinforcement (stirrups) of S500.

At the critical regions, i.e., location of the expected plastic hinges, near the extremities of the structural members (Figure 1c), adequate concrete confinement is provided by tightening the spacing between the stirrups as detailed in Figure 1d.

The experimental model was cyclically loaded in a quasi-statical way by displacement control, according to the displacement law in Figure 1b. Hence, the equation of motion will only comprise the restoring force, which will be further divided into the elastic and inelastic components. Having the displacement that the structure is subjected to and the structure's stiffness, the restoring force can be determined.

The following parametric variation will make use of the structural properties and loading conditions of the described experimental frame structure and testing, presenting all or only specific hysteresis cycles according to the current parameter under study.

2.2.   Brief description of the smooth hysteresis model considered

Sivaselvan and Reinhorn developed a polygonal and a smooth hysteresis model ^[11,24] with degradation features. As mentioned, this modification to the original Bouc-Wen model offers a more physical understanding than preceding modifications since the parameters are selected in a way to have physical meaning. A brief description of the model used in this study is outlined below based on ^[11].

The model considered herein comprises three different springs (Figure 2): a post-yielding spring (Spring 1); a hysteresis spring (Spring 2), and a slip-lock spring (Spring 3).

MATLAB was used to develop an algorithm to apply the aforementioned model. The differential equation that computes the hysteretic force is given in rate-dependent form as follows ^[11],

in which F* and F*_Y are the hysteretic force and respective yielding force, x is the total displacement, k₀ is the initial elastic stiffness, k_H is the total non-linear stiffness, k_hyst is the hysteretic stiffness, a is the post-yielding to initial stiffness ratio, N is the parameter controlling the smoothness transition from the elastic to the post-yielding range, and η₁ and η₂ are parameters controlling the shape of the unloading path. R_K is related to stiffness degradation, having a positive decreasing value, ^[11]

where α controls the degree of stiffness degradation. In addition, stiffness degradation of the elastic spring, i.e., the post-yielding stiffness k_post-yield = ak₀, can also be considered, as used by other authors ^[3,34], by the following modification in the elastic spring

The slip lock stiffness, k_slip-lock, translated by Spring 3 is thus implemented in the model in series with the hysteresis spring to simulate pinching behavior. The combined stiffness can thus be obtained by the following expression ^[11]

where s is the slip length, R_s controls the slip length, σ controls the fraction of the yield force beyond which slip will not occur, and λ defines the part of yield force about which slip will occur.

The deterioration of the strength capacity can be simulated by using a rule based on continuous energy degradation and backbone degradation due to the exceeded maximum deformation. This rule can be given by the following expression ^[11]

where F_Y^+/– and F_Y0^+/– are the yield force and the initial yield force, respectively, x_max^+/– is the maximum displacement, x_ult^+/– = x_Y × µ_ult^+/– is the ultimate displacement (product between the yield displacement and the ultimate ductility factor), and β₁ and β₂ are parameters based on ductility and energy demands, respectively. The hysteresis energy can be given in incremental form as follows ^[11]:

and the ultimate hysteretic energy under monotonic loading until the ultimate deformation without degradation can be computed by the following expression

Asymmetric yielding can also be considered by using the following expression ^[11]:

The experimental model was loaded in a quasi-static cyclic way. Hence, Eq 1 is rewritten for quasi-statically loaded systems, eliminating dt and noting that sgn($\dot x$) = sgn(dx) ^[30], considering now the motion independent of time,

and the non-linear equation is solved using a numerical solver such as functions, fzero, fsolve, or vpasolve ^[31].

2.3.   Model parameters variation

The following parametric variation investigates the meaning and the influence of each parameter in the above-presented model. The first graph of Figures 3–13 corresponds to the 18^th cycle for different parameter values, and the second plot is the complete hysteretic response for specific parameters' values.

2.3.1.   Parameter N

Figure 3 presents the variation of parameter N, controlling the smoothness transition from pre- to the post-yielding range. As can be observed the model becomes approximately bilinear from N = 15, showing very little variation between this value and the last one considered (N = 250).

2.3.2.   Parameters η₁ and η₂

These two parameters can be reduced to one, η, if one considers the following operation: η₁ + η₂ = 1, leading to ^[11]:

This is the result of prior studies ^[35,36], that suggest parameter A (in the original model ^[26]) is equal to one (Eq 1 and 10) due to redundancies in the model and for mathematical consistency, restoring the physical significance of the initial stiffness; and the previous operation (η₁ + η₂ = 1) is considered for compatibility with plasticity, i.e., returns the physical meaning of the yield force. It should be referred that these values correspond, respectively, to β and γ in the original model ^[26].

Variation of parameter η is presented in Figure 4 as combinations of η₁ and η₂. It can be verified that this parameter or parameters control the shape of the unloading path, i.e., the variation of the unloading stiffness, as well as the size of the hysteresis loop, in terms of area and shape. Parameter η will take only positive real values, which corresponds to positive real values for η₁ and positive or negative real values for η₂. When η < 0.5 the discharge path is linear. For η < 0.5 the unloading path is non-linear, compressing the loop area, and thus reducing the energy dissipation. For η > 0.5 the unloading path is also non-linear, although the loop area increases by slightly enlarging the cycle.

Other parametric studies ^[36,37] revealed that if parameters η₁ and η₂ are allowed to vary independently in the context of the original Bouc-Wen model, softening and hardening hysteresis behavior can be obtained for larger and smaller values of η₁, respectively.

2.3.3.   Parameters a and η_d

Figures 5 and 6 show, respectively, the variation of the values of parameters a and η_d, while fixing the remaining parameters.

As mentioned, parameter a accounts for the value of stiffness after yielding, i.e., in the plastic range. This parameter may be seen as a percentage of the initial elastic stiffness.

Recent studies ^[38,39] proposed the use of negative values for this parameter, resulting in a strain-softening behavior after yielding. Hence, this variation also accounts for negative values of a.

The effect of a negative value of parameter a leads to a clockwise rotation of the hysteresis loops, which allied to the softening behavior results in an increase of strength capacity at the end of the reloading phase, since this model is a Masing type hysteresis model.

The degradation of the post-yielding stiffness can also be implemented in the model ^[3,34] using the degradation parameter η_d, whose effect presented in Figure 6 reveals higher degradation for larger values of the maximum displacement and parameter η_d.

2.3.4.   Parameter α

Hysteresis stiffness degradation is regulated by parameter α, which takes positive real values. This parameter is implemented in the pivot rule ^[20], leading to a positive decreasing function, R_k (Eq 2), with the unity as its initial and maximum value.

From Figure 7 it is evident that for larger values of α the less stiffness degradation occurs. For values of α between 50 and 250, stiffness degradation is almost inexistent. Despite the influence in the stiffness of the unloading path, an obvious correlation with the hysteresis energy dissipation is seen by observing the reduced area of the hysteresis loops.

Severe stiffness degradation can be verified for values below α = 5. Values of α below 0.5 may not represent results with physical significance.

2.3.5.   Parameters β₁ and β₂

Figures 8 and 9 show the variation of the parameters controlling the strength deterioration in terms of ductility and energy demands, respectively, β₁ and β₂.

Parameter β₁ depends on the maximum displacement achieved as can be verified by the second term of Eq 5. Hence, whenever a maximum displacement is attained, degradation of the strength capacity occurs by reducing the backbone curve's ordinate. Figure 8 shows severe strength degradation for values between 0.4 and 0.6, and minor strength deterioration for values below 0.1.

Different from β₁, parameter β₂ exhibits continuous resistance degradation based on the incremental energy dissipation. This leads to greater degrees of deterioration for smaller values of β₂ compared with β₁. As can be observed severe resistance degradation can be verified for values of β₂ above 0.15. Values above 0.3 led to numerical instabilities for this case under study.

2.3.6.   Parameters R_s, σ, and λ

In this subsection variation of the parameters' values related to the pinching effect is carried out, viz., parameters R_s, σ and λ.

The slip length s (Figure 2) is regulated by parameter R_s that possess positive real values. The higher the value of R_s the more pronounced the pinching effect is as can be observed by Figure 10. The severe pinching effect can be verified for values bigger than R_s = 0.4, and minor pinching effect for very small values of this parameter.

The effect of parameter σ is clear by the observation of Figure 11, focusing the pinching effect in a smaller region near the center for smaller values of σ. Slip or the pinching effect will not occur beyond the fraction defined by σ. A value of σ near the unity leads to an almost uniform distribution of the pinching effect in the loading/reloading and unloading branches. Numerical instabilities occur for values of σ very close to zero.

Parameter λ defines the location around which slip will occur. It can take positive or negative real values, and has the positive or negative unity as its maximum or minimum value, respectively. A value different from zero introduces an asymmetric pinching effect in the numerical model. A null value for this parameter means that the structural element or the structure is symmetric. Figure 12 shows the effect of this parameter, which becomes more pronounced for values close to positive or negative one.

2.3.7.   Asymmetric yielding

Asymmetric yielding can be considered by the use of Eq 8. Figure 13 shows this feature by varying the yield force in the negative direction to have lower and higher yielding levels than in the positive direction. This feature is important since structural systems are not perfect or ideal, and yielding will possibly have different values in opposite directions.

3.   Calibration of model parameters

The cyclic response of the RC experimental bare frame is presented in Figure 14a. It can be verified a smooth evolution with a subsequent maximum before concrete cracking at the top and bottom ends of the columns. Then, a gradual stiffness decrease occurs without collapse, although with substantial damage, strength deterioration, and inelastic hinge spread in the columns.

Based on the previous parametric study that investigated the influence and effect of each parameter in the hysteresis response, an optimization procedure is carried out to search for the parameters' values that provide the best fit to the experimental data.

The optimal parameters' values defining the smooth hysteresis model that provide the best fitting to the experimental data are obtained by following the procedure illustrated in Figure 15 through a flowchart. The procedure starts with a manual calibration to obtain a starting point or initial guess of the hysteretic parameters for the implementation of function lsqcurvefit in MATLAB ^[31]. This function aims to minimize the sum of squares error between the numerical and experimental results of the restoring force. The optimized results with this function are then used as an initial guess for the use of function, fminsearchbnd ^[40], based on the fminsearch function algorithm ^[31] that follows the Nelder-Mead simplex algorithm, although with bound constraints applied to the parameters. This function aims to minimize the root mean square error (RMSE) of the same results mentioned, with emphasis given to the points belonging to the cyclic backbone curve through the use of weighting factors. This error was further normalized (NRMSE) with the range of the experimental restoring force and converted to a percentage. The stopping criteria considered for a converged solution were an NRMSE less than 5% and a tolerance less than 1 × 10⁻⁴ for the difference between two subsequent iterations on the calculation of the NRMSE.

The initial stiffness and the yielding force in both directions were also considered in the optimization procedure, verifying that the optimum values are relatively close to the actual values. Figure 14b presents the outcome of the optimization with the following optimized parameters: N = 2.07, a = −0.061, η_d = 0, η = 0.81, α = 1.33, β₁ = 0.03, β₂ = 0.005, R_s = 0.17, σ = 0.40, λ = 0.004, k₀ = 3.567 × 10⁶ N/m, F_Y0⁺ = 4.8203 × 10⁴ N, F_Y0^- = 4.9907 × 10⁴ N, µ_ult⁺ = 7.85, µ_ult^- = 14.85. The model can reasonably estimate the cyclic backbone curve as well as the global hysteretic cycles. The model presents some setbacks in the estimation of the experimental unloading branches and transitions from elastic to plastic range. In addition, further improvements may be applied in the model for the present purpose, viz., in the strength degradation parameters, whose influence and numerical values seem to not correspond to the level of deterioration in the experimental model, taking very small values. However, the small values of these parameters might be explained by the high degree of stiffness degradation and pinching effect.

The functions used for optimization gave satisfactory results, though a large number of variables were used, exceeding their limit of good performance. As a future development, an optimization procedure should be created to verify if a better solution exists to fit the experimental data.

4.   Conclusions

The parametric study performed herein highlighted the importance of each parameter in the hysteresis model investigated. The knowledge of the effect of each parameter in the model provides the necessary insight into the possibilities of the model to estimate real behavior, in the present scope, the non-linear cyclic behavior of RC frame structures.

Based on the parametric study carried out, the calibration of the model's parameters was performed using an optimization procedure to replicate the cyclic behavior of an experimental RC bare frame structure. Satisfactory results were obtained, though a large number of variables were used in the process of optimization, which may compromise the good performance of the implemented functions. Hence, a different optimization approach should be considered to support the current and future solutions.

Further improvements may be carried out for the present purpose, namely, to the strength deterioration formulation to better predict the experimental response, since the values obtained for the corresponding parameters are very small, which may not represent the actual behavior of deterioration. Nevertheless, the level of stiffness degradation and the pinching effect may also explain such low values of the strength deterioration parameters.

Future investigations may comprise validation of the studied hysteresis model with more experimental data considering different hysteretic configurations.

Acknowledgments

This paper is within the scope of the first author's Ph.D. degree in progress, financially supported by the Portuguese Foundation for Science and Technology (FCT) through the PhD grant reference SFRH/BD/139570/2018 under the programme POCH (N2020-P2020) and subsidized by the European Social Fund (FSE) and national funds from MCTES. This work was financially supported by: Base Funding-UIDB/04708/2020 of the CONSTRUCT-Instituto de I & D em Estruturas e Construções-funded by national funds through the FCT/MCTES (PIDDAC).

Conflict of interest

There is no conflict of interests with the authors and the publication of this paper.