A degrading Bouc–Wen model for the hysteresis of reinforced concrete structural elements

Abstract This paper presents a smooth hysteresis model for reinforced concrete (RC) structural elements based on the differential equation of the Bouc–Wen model. Stiffness degradation and strength degradation are defined by introducing a damage index that includes both dissipated energy and maximum displacement. The pinching effect acts directly on the stiffness of the system and is controlled by an activation energy. The degrading functions are connected to the actual processes with which the damage occurs, thereby giving each parameter a physical meaning. The simple formulation of the model allows a straightforward identification of the best-fitting parameters and an easy interpretation of the results. Applications to the cyclic behaviour of RC structural elements demonstrate that the model is well capable of describing complex hysteretic behaviours involving degradation and pinching effects.


Introduction
The hysteresis phenomenon is observed in various mechanical systems. For instance, engineering structures often exhibit hysteretic behaviour under severe cyclic loads caused by earthquakes (Roufaiel & Meyer, 1987). Repeated cyclic deformations frequently cause a deterioration of the characteristics of reinforced concrete (RC) mechanical systems. In fact, building structures may experience opening and closing of cracks, postyielding and buckling of metallic elements (Kashani, Salami, Goda, & Alexander, 2018), strength and stiffness deterioration, and other local inelastic behaviours (Sivaselvan & Reinhorn, 2000). All these phenomena contribute to the structural damage. An accurate prediction of the behaviour resulting from a cyclic loading must take into account the damage of the structure. Therefore, constitutive models capable of representing non-linear deterioration are required.
The evaluation of the seismic damage of structures needs, firstly, a definition of a suitable measure for the damage itself. Most of the damage indicators available in literature involve two quantities: the energy dissipated during the inelastic cycles of hysteresis and the displacement ductility, defined as the ratio between maximum displacement and yielding displacement. For instance, Park and Ang (1985), Kunnath, Reinhorn, and Park (1990), Banon and Veneziano (1982) and Singhal and Kiremidjian (1996) considered both ductility and dissipated energy to assess structural damage phenomena. Greco and Marano (2015) and Marano and Greco (2006) employed a damage index involving displacement ductility and hysteretic energy to perform a stochastic analysis of the cumulative damage in structural systems subjected to ground motions. Elenas (2000) and Bassam, Iranmanesh, and Ansari (2011) performed, respectively, an evaluation of the damage of RC bridges and RC frame structures based on ductility and energy dissipated by the system. Consequently, mathematical models able to represent hysteretic cycles with degradation should involve both displacement ductility and dissipated energy.
The BWBN model is an extension of the Bouc-Wen model, developed through the incorporation of functions that describe degradation and pinching. However, as pointed out by Marano et al. (2017), the formulation of the BWBN model is complicated and not all the parameters involved have a direct physical meaning. Furthermore, degradation and pinching are described only on the basis of the dissipated energy.
In this work, a smooth hysteresis model for degrading RC structural elements is presented. The model is based on the differential equation of the Bouc-Wen model, to which appropriate degrading functions are added in order to depict the cyclic damage. Stiffness degradation and strength degradation are defined through a damage index that involves both maximum displacement and dissipated energy. Moreover, the pinching effect is depicted along the same lines of the simple formulation proposed by Pelliciari et al. (2018), incorporating modifications that correct its physical deficiencies.
The mathematical description of the degradation of the system is based on the physical processes with which the structural damage takes place. Thereby, each parameter of the proposed model is physically based and related to a certain damage phenomenon, allowing an easy control and interpretation of the results. The formulation of the present model is straightforward and oriented towards practical engineering applications. Indeed, the purpose of this study is to provide a simple and reliable model able to predict the hysteresis of RC structural elements with pinching and degradation, including at once a physically based description of the damage.
The paper is organised as follows. In Section 2, an outline of the Bouc-Wen model is provided. The damage index is then introduced and the functions that account for the mechanical degradation of the system are described, including the pinching effect. Section 3 provides an overview of the parameter optimisation strategy, which is based on a genetic algorithm (GA). Consequently, in Section 4, the effectiveness of the proposed model is verified by its application to the cyclic behaviour of a circular RC bridge pier and two rectangular RC columns. Such structural systems exhibit complex cyclic behaviours, involving significant damage. Finally, conclusions are drawn in Section 5.

Hysteresis model for degrading systems
In this section, the Bouc-Wen model of hysteresis is briefly introduced. At the same time, each improvement that allows to account for degrading phenomena is discussed.

Preliminaries on the Bouc-Wen model
The equation of motion for a single-degree-of-freedom system is: where, x is the relative displacement of the mass of the system m with respect to the ground, c is the linear viscous damping coefficient, F s ½xðtÞ, zðtÞ, t is the non-damping restoring force, z(t) is the hysteresis displacement and F(t) is the external excitation. The overdots indicate the derivative with respect to the time, thus _ x and € x represent velocity and acceleration, respectively.
The Bouc-Wen model gives the following expression for the restoring force: where, k is the elastic stiffness of the system and a is the ratio between the final tangent stiffness k f and the elastic stiffness (0 a 1). Equation (2) is composed of two contributions: the linear elastic component akxðtÞ and the hysteresis component ð1ÀaÞkzðtÞ, which depends on the past history of stresses and strains. The hysteresis displacement z(t) is given by the differential equation: with the initial condition zð0Þ ¼ 0: The parameters b, c and n control the shape of the hysteresis cycles, while A determines the tangent stiffness. As is well-known in literature (Charalampakis & Koumousis, 2008b;Ma, Zhang, Bockstedte, Foliente, & Paevere, 2004), the aforementioned parameters are redundant. Thus, A is commonly set to unity in order to eliminate this redundancy. However, since in this work the parameter A will be used in the mathematical developments for the description of the degradation, this position will be performed later (see Section 2.3). Note also that Ismail, Ikhouane, and Rodellar (2009) gave the following conditions for the Bouc-Wen model parameters in order to satisfy the thermodynamic admissibility, the accordance with Drucker and Il'iushin stability postulates and the uniqueness of the solution: The hysteresis energy dissipated by the system is defined as the area under the hysteresis restoring force F h ½zðtÞ, t ¼ ð1ÀaÞkzðtÞ along the total displacement x(t). The hysteresis energy can be normalised with respect to the mass, as follows (Foliente, 1995;Ortiz, Alvarez, & Bedoya-Ru ıZ, 2013): Some authors found out that for many structural systems the two parameters b and c are related to each other. For instance, Sues, Mau, and Wen (1988) suggested that b ¼ c for steel structures and b ¼ À3 c for RC structures, whereas Sengupta and Li (2013) applied a Bouc-Wen model to RC beam-column joints, obtaining b ¼ À5 c: Ye and Wang (2007) estimated the Bouc-Wen model parameters considering both the cases of softening and hardening behaviour. The results showed that for a system with softening b ffi 2 c, while for a system with hardening b ffi À c: In view of the above, it appears reasonable to establish a linear proportionality between the parameters b and c c ¼ g 0 b: Given this position, the conditions (4) reduce to: b > 0, À1 < g 0 1, n ! 1: Note that, with the position (6), the constraints on the variability of the parameters are now decoupled, as expressed by Equation (7). This simplifies the mathematical formulation of the model and represents a great advantage during the identification of the optimal parameters. Substituting Equation (6) into Equation (3), the differential equation for the hysteretic displacement becomes:

Damage index
As already pointed out, the most popular damage indicators involve both hysteretic dissipated energy and maximum displacement (Park & Ang, 1985). Hence, in this work, the authors propose a model that connects the degradation phenomena to the following dimensionless damage index: where, x u is the ultimate displacement under a monotonic loading, f y ¼ F y =m is the yielding force normalised with respect to the mass and x max ðtÞ is the maximum displacement of the system until the time t. Equation (9) brings two new parameters into the hysteresis model: q and q x . Their purpose is to weigh the influence of the degradation due to the energy dissipation and the maximum displacement, respectively. Note that the magnitude of both ðtÞ and x max ðtÞ is not known a priori. On the contrary, it depends on the excitation to which the mechanical system is subjected. Therefore, it is not possible to normalise Equation (9) in a way that q and q x vary in the range between 0 and 1.

Stiffness and strength degradation
For systems with softening behaviour, the following derivatives can be computed (Baber & Wen, 1980): where, the ultimate value of the hysteretic displacement, z u , is expressed by: The ultimate hysteretic restoring force component can be thus defined as (Cunha, 1994): The stiffness degradation can be introduced into the model by expressing k i as a function that degrades with the damage index d i ðtÞ: Nevertheless, the elastic stiffness k should not be subjected to any degradation. Thus, looking at Equation (10), one can conclude that the degradation of the stiffness must be entirely included into the parameter A, by substituting this parameter with a function of d i ðtÞ: In order to decouple stiffness and strength degradation, the ultimate hysteretic restoring force component must remain constant while the degradation of the function Aðd i Þ occurs. Therefore, the parameter b 0 in Equation (13) is replaced by a function b k ðd i Þ, in such a way that the following relation holds: It is often assumed that the damage of a structural system evolves with the dissipation of energy following a negative exponential law (see, among the others, Ba zant, Pan, & Pijaudier-Cabot, 1987;Lanzoni & Tarantino, 2014Miehe, 1995;Mazars, 1986;Tarantino, 2014). Thus, the following form for the function Aðd i Þ is adopted: where, d k is a parameter that controls the amount of stiffness degradation. Obviously, for d k ¼ 0, the model does not exhibit stiffness degradation. Note that the redundancy of the Bouc-Wen model described in Section 2.1 is now removed by turning the parameter A into the function Aðd i Þ, where the only unknown parameter is d k . At this point, given (15), Equation (14 This procedure allows incorporating the whole description of the stiffness degradation into the function b k ðd i Þ: Figure 1 shows the effect of variations of d k on the hysteresis cycles. The higher the value of d k , the more the degradation of the stiffness is pronounced. It is worth noting that the maximum value of the force acting on the system is not influenced by the stiffness degradation. In the following, the strength degradation is described separately from the stiffness degradation and the two formulations will be combined at the end of the present section. Similarly to what has been done above, the degradation of the strength is introduced in a way that it is decoupled from the degradation of the stiffness. This is done by defining the ultimate hysteretic restoring force F u h as a degrading function of d i ðtÞ: Hence, the initial stiffness k i must remain constant. This means that, looking at Equation (10), the parameter A must be constant. Therefore, the degradation of the strength is introduced by replacing b 0 with b f ðd i Þ into Equation (13), with b f ðd i Þ defined by the following exponential law: where, d f is a parameter that controls the amount of strength degradation. Note that for d f ¼ 0, there is no strength degradation. Substituting Equation (17) into Equation (13), the following relation is obtained: which states that the degradation of the strength follows a negative exponential law with the damage index d i ðtÞ: Note that when strength and stiffness degradation are combined, the parameter A is replaced by the function Aðd i Þ expressed by Equation (15), where the redundancy of the Bouc-Wen model has already been removed. In addition, it should be pointed out that the description of the strength degradation is entirely included in the function b f ðd i Þ, in analogy with b k ðd i Þ for the stiffness degradation. Figure 2 shows the effect of variations of d f on the hysteresis cycles. Increasing the value of d f , the degradation of the strength increases. It is noted that the initial slope of the cycles, which represents the initial stiffness k i , remains unchanged.
Stiffness and strength degradation have been described above by establishing proper degrading functions of the damage index. However, as already pointed out by Pelliciari et al. (2018), the rate of degradation of the stiffness of a system might not always maintain the same proportionality to the energy (in this case to the damage index). For instance, the degradation rate could be stronger at the beginning of the inelastic phenomena. Then, when the system is severely damaged, its proportionality to the damage index could become less. Thus, the increasing of stiffness degradation rate is here controlled introducing a negative exponential function p k ðd i Þ into Equations (15) and (16): where, w is a parameter that controls the rising of stiffness degradation. Specifically, for w ¼ 0 the stiffness degradation rate follows the increasing of d i ðtÞ through the exponential law expressed by Equation (15). Instead, as w increases, the stiffness degradation rate decreases its proportionality to d i ðtÞ as the event proceeds. The limit situation is w ! 1, for which Aðd i Þ ! 1: Therefore, in this last case, no stiffness degradation occurs. In order to clarify the role of w, its effect on the hysteresis cycles is displayed in Figure 3. At this point, the effects of both strength and stiffness degradation are combined by defining a single function that includes both degradation phenomena: and the differential equation for the degrading hysteresis model assumes the form: where, the function p k ðd i Þ is expressed by Equation (21). For the sake of clarity, the dependence on the time t has been omitted.  Equation (23) governs the hysteresis model including the degradation of stiffness (controlled by the parameters d k and w) and the degradation of strength (controlled by the parameter d f ). In the following, the pinching effect will be introduced.

Pinching effect
A phenomenological and straightforward way of simulating the pinching phenomenon is to consider it as a degradation process that acts on the stiffness of the system. More precisely, it can be represented as a reduction of the stiffness that occurs in a central range of the force-displacement curve, followed by an increase in the outside regions. It goes without saying that this is not a strict definition, being the nature of pinching accidental and extremely complex. Nevertheless, this definition was employed by Marano et al. (2017) and Pelliciari et al. (2018), obtaining promising results from the application on RC bridge piers.
In the aforementioned works, the authors described this effect through a Gaussian function. However, in their model, this effect starts at the beginning of the event. Hence, the pinching affects every hysteresis cycle, including the very first ones. This is obviously not consistent with the reality, or at least it does not happen for every kind of structural elements and for every displacement history. Therefore, in this work, the same analytical description of the pinching will be adopted, but with the introduction of a new parameter that acts as an activation energy for the pinching effect.
The pinching function is defined as a normalised Gaussian function with mean equal to zero: The higher the value of the exponent u, the steeper the stiffness change, while the standard deviation r controls the width of the Gaussian bell (see Pelliciari et al., 2018, for a more detailed description). The function g p ðxÞ is then introduced directly into the elastic stiffness k, that turns into a degrading stiffness through the following relation: where, q p 2 ½0, 1 and f ðÞ is defined as: The parameter q p defines the amount of reduction with respect to the elastic stiffness (q p ¼ 0 means no pinching), while the function f ðÞ has been introduced in order to make the pinching effect start at a certain activation energy, defined by the parameter p . The exponent in Equation (26) has been set to 8 in order to provide a fast activation of the pinching effect as soon as p is reached. This high exponent does not produce singularities or numerical problems. Note that, regarding the simulation of the pinching phenomenon, the introduction of f ðÞ represents the novelty of this work with respect to what was proposed by Pelliciari et al. (2018). Figure 4(a) presents the effect of p on the forcedisplacement response of the system. As p increases, the activation of the pinching is delayed. The corresponding shape of the function f ðÞ is shown in Figure 4(b). It is clearly visible that when the system reaches the activation energy p the value of f ðÞ goes rapidly from 0 to 1, which causes the starting of the pinching effect. It should be noted that the case p ¼ 0 cannot happen, because it would produce a singularity. However, since a real system never exhibits a pinching phenomenon that starts at the very beginning of the event, this case is not of interest.
The restoring force is then computed by substituting the elastic stiffness k with the degrading stiffness k p ðÞ, expressed by Equation (25), into Equation (2): F s ½xðtÞ, zðtÞ, t ¼ ak p ðÞxðtÞ þ ð1ÀaÞk p ðÞzðtÞ: In conclusion, the parameters of the model are the following: a, k, b 0 , g 0 , n, q x , q , d k , d f , w, r, u, q p , p : Table 1 summarises the role that each parameter has on the definition of the hysteresis cycles. Although the model is composed of a significant number of parameters, each one has a specific physical meaning. Moreover, the pinching effect acts directly on the stiffness of the system through the parameters r, u, q p and p , which can be easily identified and controlled. The relatively simple formulation of the degrading phenomena and the pinching effect brings benefit to the interpretation of the results and the use of the proposed model.

Parameters optimisation algorithm
The values of the model parameters must be identified on the basis of experimental data. The aim is to find the optimal set of parameters, which provides the simulated response of the system that best fits the experimental one.
The set of parameters is collected into a vector h, named parameter vector. The optimal parameter vector is denoted as h Ã and it is the one that minimises an objective function (OF), which is a measure of discrepancy between numerical and experimental response. For the model proposed in this work, the parameters variation is constrained by defining the lower bound h l and the upper bound h u : Since the model must provide the best simulation of the cyclic behaviour of the system during the whole event, the OF must be an integral measure over the whole set of experimental data. Therefore, the OF is defined as: where, x i and x f are the initial and final displacement records, F e ðxÞ is the force derived from the experimental data and F s ðh, xÞ is the one simulated by the model. Hence, as expressed by Equation (28), the objective function is here defined as a normalised integral of the difference between experimental force and simulated force.
It is clear that the problem of minimisation of the objective function described above is not simple. This is mainly due to the complexity of the differential equation that governs the hysteresis model (Equation (23)). Therefore, for such complex optimisation problems, an evolutionary approach is usually adopted (Greco & Vanzi, 2019;Quaranta, Marano, Greco, & Monti, 2014).
Many identification of the parameters of structural models via a genetic algorithm (GA) can be found in literature. Marano, Quaranta, and Monti (2011) developed a modified real-coded genetic algorithm to identify the parameters of structural systems subjected to dynamic loads. Sengupta and Li (2013) used a GA to calibrate the parameters of a modified BWBN model applied to RC beam-column joints. Charalampakis and Koumousis (2008a) developed a hybrid evolutionary algorithm based on a GA, which was employed for the prediction of the behaviour of a steel cantilever beam. Sireteanu, Giuclea, and Mitu (2010) used a GA to identify the Bouc-Wen model parameters for elastomeric isolators and buckling restrained dissipative braces, while Ha, Kung, Fung, and Hsien (2006) performed the identification of a Bouc-Wen model via a real-coded genetic algorithm (RGA), for the case of piezoelectric actuators.
In view of this, the genetic algorithms represent suitable tools for complex identification problems. Thus, in this work, the optimisation of the objective function is performed via a genetic algorithm.

Application to reinforced concrete structural elements
In the following, the model is applied to a RC bridge pier subjected to a cyclic test carried out in the lab of the Fuzhou University. First, the experimental test is presented. Then, the parameters of the model are calibrated using the experimental data. Furthermore, comparisons with the results obtained with a model deprived first of the contribution of the maximum displacement (q x ¼ 0) and second of the dissipated energy (q ¼ 0) are performed. This allows to assess the benefits of the introduction of both the energy and the maximum displacement into the damage index.
Two further applications are then reported in order to ensure the reliability of the proposed model.

Experimental test
The identification of the model parameters is performed on the basis of the experimental data of a RC bridge pier. The experiment was carried out in the lab of the Fuzhou University (Xue et al., 2018). The pier considered is named R16-1B and it comes out from a repairing process of the original pier P16-1B, which was previously damaged by a cyclic test. The pier P16-1B had a circular cross section and it was made of Chinese concrete C30, steel HRB335E for the longitudinal reinforcement and steel R235 for the transversal reinforcement (D62-2004, 2004). Its height and diameter were 1.17 m and 0.42 m, respectively. The repairing consisted in the restoration of the damaged portion of the longitudinal steel reinforcement by means of turned rebar parts and substitution of the stirrups and the damaged concrete parts with HPFRC (High Performance Fiber Reinforced Concrete). A detailed description of the repairing process can be found in the following studies (Albanesi, Lavorato, Nuti, & Santini, 2009;Lavorato, Nuti, Santini, Briseghella, & Xue, 2015;Savino, Lanzoni, Tarantino, & Viviani, 2018).
The cyclic test, which was carried out first on the pier P16-1B and consequently on the pier R16-1B, consisted in the imposition of a displacement history at the top of the pier. As Figure 5 shows, the displacement history is  composed of a first part and a second part. The first part corresponds to the response of the bridge column to the Tolmezzo accelerogram (Xue et al., 2018), referred to an earthquake occurred in Italy in 1976. The second part represents the response of the bridge column to the Tolmezzo accelerogram scaled to double. The Tolmezzo accelerogram has been selected because it is representative of a strong earthquake for the bridge analysed. The second part (scaled to double) has been added in order to further investigate the hysteresis behaviour of the pier when subject to severe damage. The test set-up is composed of a vertical load system and an horizontal actuator (Figure 6(a,b)). The vertical load system applies a constant load of 266 kN on the top of the pier while the horizontal actuator impresses the displacement history. The constant vertical load represents the bridge deck load, properly scaled. As a result, the force-displacement curve is obtained. More details on the mechanical characteristics, the repairing process and the test set-up are reported in Xue et al. (2018). Figure 7 presents a schematic representation of the cross section of the pier R16-1B.
Similar experimental loading systems have been employed by several authors in order to obtain the behaviour of RC columns under cyclic lateral loading and a constant axial load, which simulates the dead load deriving from the superstructure (see, among the others, Colomb, Tobbi, Ferrier, & Hamelin, 2008;He, Sneed, & Belarbi, 2013;Sun, Wang, Du, & Si, 2011).

Calibration of the model parameters
The experimental force-displacement data of the pier R16-1B are used to perform the identification of the optimal parameters of the model, which are gathered in the following parameter vector: As pointed out by Pelliciari et al. (2018), a suitable value of u for the description of the behaviour of RC piers is u ¼ 4. Thus, its value is fixed and the identification involves only the remaining 13 parameters. The values of ultimate displacement and horizontal normalised yielding force for the pier R16-1B are, respectively, x u ¼ 60 mm and f y ¼ 0:63 kN/kg, considering the presence of a constant vertical load of 266 kN applied on the top of the pier. The lower and upper bounds of the parameter vector are reported in Table  2 and the number of generations performed by the genetic algorithm is fixed to 100.
The identification process is performed in the following three cases: 1. the damage index includes both maximum displacement x max ðtÞ and dissipated energy ðtÞ, as expressed by Equation (9); 2. the damage index does not include the maximum displacement (q x ¼ 0), therefore the degradation phenomena are only function of the dissipated energy ðtÞ; 3. the damage index does not include the dissipated energy (q ¼ 0), namely the degradation phenomena are only function of the maximum displacement x max ðtÞ: Table 3 lists the optimal parameters and the final OF value for the three cases. As expected, the introduction of both maximum displacement and dissipated energy in the damage index (case 1) provides the best final simulated response, reaching an OF around 9:5 %: Instead, the simulations in cases 2 and 3 provide an OF that is approximately 5 % and 1 % higher, respectively. Figure 8(a-c) present the comparison between experimental data (F exp ) and simulated result (F s ) in terms of force-displacement curves in case 1, 2 and 3, respectively. The simulation in case 1 (Figure 8(a)) shows a good fitting of the experimental hysteresis curve during the whole test. Instead, in case 2 (Figure 8(b)), it can be observed that the simulated behaviour is not accurate for the last hysteresis cycles, i.e. when the damage is severe. This is confirmed by the plot of experimental energy exp and simulated energy s (Figure 9(b)), which shows that after a certain point of the test the simulated dissipation of energy is significantly lower than the experimental one. It follows that the description of the damage as a function of the sole hysteretic energy is not satisfactory for the complex behaviour analysed here, while case 1 provides a more sound and reliable result (Figure 9(a)).
Turning to the force-displacement results in case 3 (Figures 8(c) and 9(c)), it is noted that the simulation fits quite well the experimental behaviour. Although, the first hysteresis cycles (i.e. when the degradation is still not severe) are not adequately described. The surprisingly low value of the OF (see Table 3), almost as low as the simulation in case 1, is explained by the fact that the main discrepancy of the simulation is concentrated in the first part of the test, during which the restoring force of the system assumes relatively low values. Therefore, this discrepancy has a small influence on the magnitude of the OF, that is an integral measure over the whole event (Equation (28)). Indeed, the OF alone would denote that the simulation in case 3 is nearly as accurate as that in case 1. However, it follows from the analysis of the simulated response that its accuracy is significantly lower than in case 1.  Figure 10 shows the hysteretic response of the system divided into two parts. The dividing point is the instant of time when the displacement history switches from the first part to the second part ( Figure 5). As already discussed, the simulation in case 1 provides a well-fitting of the experimental behaviour throughout the event (Figure 10(a,b)). Instead, the simulated behaviour in case 2 is accurate in the first part, while there is a significant discrepancy in the second part (Figure 10(c,d)). Last, the simulation in case 3 is capable of describing the second part, whereas the first part is affected by a considerable error (Figure 10(e,f)).
In order to provide a comparison of the local accuracy of each simulation, the following relative error in terms of energy is computed: e ðh, tÞ ¼ j exp ðtÞ À s ðh, tÞj exp ðtÞ : The errors made in case 1 (e , 1 ), case 2 (e , 2 ) and case 3 (e , 3 ) are displayed in Figure 11. The plot confirms all the above observations. In fact, the error in case 3 is higher during the first part of the test, while the error in case 2 is higher during the second part of the test. The simulation in case 1 provides the lowest local error during most of the test. Therefore, it is concluded that the incorporation of both maximum displacement and energy in the damage index brings benefits to the capability of the model to simulate the hysteresis behaviour of the system. Finally, Figure 12 shows the value of the objective function with the progress of the generations in cases 1, 2 and 3. For each case, the simulation has already reached a good result after 20 generations. This indicates that the genetic algorithm      does not struggle too much to find the optimal solution, thanks to the relatively simple formulation of the model.

Further applications of the proposed model
In order to guarantee the reliability of the hysteresis model presented in this paper, two applications on rectangular RC columns are now presented. The experimental data were obtained from the PEER database for RC elements (Berry, Parrish, & Eberhard, 2004, https://nisee.berkeley.edu/spd).
Since the purpose of the following applications is to further investigate the reliability of the proposed model, the two cases with q x ¼ 0 and q ¼ 0 are not considered. The previous application already demonstrated that the introduction of a damage index containing both q x and q allows reaching a higher accuracy. The first application regards the column ORC1, whose cyclic response was experimentally investigated by Aboutaha and Machado (1999). The column ORC1 was composed of a concrete with strength 83 MPa, reinforced with eight 25 mm and grade 414 MPa longitudinal reinforcing bars. The size of the cross section was 305 Â 508 mm, while the height of the column was 1.829 m. The transversal reinforcement was formed by 4-legged shear stirrups with a diameter of 10 mm and steel grade 414 MPa, spaced at 75 mm in the potential plastic hinge region and 150 mm outside that region. The ultimate displacement and horizontal normalised yielding force for the specimen ORC1 are x u ¼ 109:74 mm and f y ¼ 0:29 kN/kg. The second application concerns the column C5-00, which was tested by Matamoros (2000). The column C5-00 was composed of a concrete with strength 37.9 MPa. The rectangular cross section had dimensions of 203 Â 203 mm and it was longitudinally reinforced with four bars with diameter 15.9 mm and yielding strength 572.3 MPa. The height of the column was 0.61 m and the transversal reinforcement was performed by introducing 2-legged shear stirrups with a diameter of 9.5 mm and yielding stress 513.7 MPa, at a spacing of 76.2 mm. The column C5-00 has an ultimate displacement x u ¼ 50:8 mm and an horizontal normalised yielding force f y ¼ 0:84 kN/kg.
Both the columns ORC1 and C5-00 were tested under the imposition of a cyclic horizontal displacement at the top of the system, while no axial load was applied. The section at the bottom of the columns was fixed, thus the test configuration was a cantilever. Figure 13(a,b) shows a schematic representation of the cross sections of the columns ORC1 and C5-00, respectively. The experimental forcedisplacement hysteresis curves of the two columns are used to perform the identification of the optimal parameters of the model for both of them. Once again, the parameter u is fixed to the value 4 (see Section 4.2) and the number of generations performed by the genetic algorithm is set to 100. The lower and upper bounds of the parameters are the same as the ones listed in Table 2, with the exceptions of the ranges of k, r and p . These three parameters are not dimensionless, and therefore their ranges must be adjusted basing on the particular system being examined. For this reason, both for the column ORC1 and C5-00, their ranges have been modified before running the optimisation process. Table 4 lists the optimal parameters and the final OF value for the columns ORC1 and C5-00. Figure 14(a,b) shows the experimental data and the simulated force- Figure 12. OF values along with the generations during the identification process for the pier R16-1B in case 1, case 2 and case 3.  . Relative error e r in terms of dissipated energy committed in case 1 (e , 1 ), case 2 (e , 2 ) and case 3 (e , 3 ). displacement curves for the columns ORC1 and C5-00, respectively. It can be observed that, for both systems, the hysteresis behaviour simulated by the model is well fitting the experimental one. In fact, the error in terms of OF is around 10 % in both cases.
It is worth noting that the pronounced strength degradation and pinching effect of the column C5-00 with respect to the column ORC1 is directly reflected on the values of the parameters d f and q p . In fact, as reported in Table 4, the optimal values of d f and q p for the column C5-00 are sensibly higher than those for the column ORC1.

Conclusions
This paper presents a smooth hysteresis model for RC structural elements with damage and pinching effects. Its formulation is based on the differential equation of the Bouc-Wen model, which has been extensively used for the description of the hysteresis of various mechanical systems.
The deterioration of the mechanical characteristics is introduced through a damage index that is defined as a linear combination of dissipated energy and maximum displacement. Stiffness and strength degradation are mathematically described as decoupled phenomena using negative exponential laws. Furthermore, the pinching effect is introduced in a similar way as what was done by Marano et al. (2017) and Pelliciari et al. (2018). However, the model proposed in the abovementioned works involves a pinching phenomenon that starts at the beginning of the event, affecting also the very first hysteresis cycles. In the present paper, this physical inconsistency is resolved by including an activation energy for the pinching effect.
The model is applied to the hysteresis loops of a RC bridge pier tested in the lab of the Fuzhou University. The following three cases are analysed: damage index involving both dissipated energy and maximum displacement, damage index involving only dissipated energy and damage index involving only maximum displacement. The results show that the description of the damage as a function of both dissipated energy and maximum displacement produces a simulation that well fits the experimental hysteresis curve during the whole test. Instead, modelling the damage as a function of the sole dissipated energy leads to an inadequate simulation of the behaviour of the pier when the damage is severe. On the contrary, a description of the degradation considering only the maximum displacement is not accurate during the first part of the test, where the damage is still slight.
Two further applications on rectangular RC columns are also presented, with the aim of ensuring the reliability of the proposed model. The results demonstrate that the model is capable of reproducing the complex hysteresis behaviour of the structural elements analysed, that involves stiffness degradation, strength degradation and pinching effect.
The simple formulation of the model allows to perform the identification of the best-fitting parameters with a small computational burden. Moreover, each parameter has a direct physical meaning, which helps in the interpretation of the results.
The hysteresis model presented in this work is oriented towards practical engineering applications, with a view to providing an instrument that can be easily managed and applied to simulate the hysteresis of RC structural elements. Although the applications presented provided  accurate results, the general applicability of the model has to be further investigated. In order to verify this, future developments will involve simulations on other types of RC elements, such as beam-column joints, frames, walls, etc.