An inelastic beam element with hysteretic damping

In this work, an inelastic beam macro-element that incorporates hysteretic damping is presented. Based on classical theory of plasticity, a Bouc–Wen type model is utilized that simulates the hysteretic behavior of an inelastic spring element. Using this model, an inelastic nonlinear beam element is formulated based on the appropriate combination of two coupled nonlinear spring elements. The equations of motion are determined and are cast in a state-space form for the vector of the end displacements, velocities and hysteretic forces. The system is solved by employing a Runge-Kutta type of algorithm. The proposed inelastic beam model is then employed to simulate the experimental dynamic behavior of steel beams. The model parameters are estimated with the aid of a nonlinear system identification algorithm using existing experimental data. The proposed element approximates the inelastic behavior of steel beams adequately within plastic regions that do not undergo substantial stiffness degradation, or other relevant phenomena. Finally, the hysteretic damping features of the model are demonstrated.


Introduction
Advances in the nonlinear analysis of structures have revealed the particular features of structural behavior under dynamic loading, offering the ability to incorporate further physical phenomena and propose more appropriate means for their modeling.
Structures under dynamic loading usually exhibit inelastic deformations.In the case of cyclic loading the nonlinear restoring force of structures, as a function of displacement, generates hysteretic loops, which are memory-dependent, not only on the instantaneous deformation, but also on the history of deformation.In 1927, Kimball and Lovell [9] demonstrated that many engineering materials exhibit a type of internal or hysteretic damping in which the energy loss per cycle is proportional to the square of the strain amplitude and independent of the frequency of the applied sinusoidal strain.This result was confirmed further by experimental results and strengthened the notion of hysteretic damping.Due to its mathematical simplicity though, the rate dependent viscous damping is the one that dominates the formulation and solution of the general structural dynamic problem.This is usually selected to represent all forms of damping by establishing an equivalent viscous damping coefficient equating the dissipated energies.One way to incorporate hysteretic damping is to use the complex stiffness matrix formulation [6], which from an engineering and mathematical point of view leads to complications.Another widely used model for describing hysteretic behavior is that of Bouc-Wen [17].The hysteretic force is defined as a function of the hysteretic variable, the evolution of which is determined by a nonlinear differential equation.In this work, a generic beam macro element that takes into account the rate independency of the hysteretic damping forces, suitable for the inelastic dynamic analysis of structures, is presented.The proposed model is capable of simulating the dynamic behavior of inelastic systems with practically constant stiffness in loading and unloading, either with hardening or softening behavior.Important contribution in this field has been made, among others, by Cherng and Wen [1], Iwan [7,8] and Song et al. [16].The influence of stiffness degradation, and/or other deteriorating phenomena, as well as of geometric non-linearities are not taken into account, but can be incorporated easily according to Sivaselvan, and Reinhorn [15].Moreover, Gaul and Bohlen [5] have suggested similar type of elements for structural connections.The main features of the proposed element are related to the internal force formulation and the incorporation of rate independent hysteretic damping, during oscillations.
The element is formed by coupling two inelastic springs, the behavior of which is controlled by their end rotations and transverse displacements as nodal degrees of freedom.The behavior of each spring is determined by a Bouc-Wen type model, which is based on the classical theory of plasticity [12].For the elastic case, the beam element reduces to the classical Euler-Bernoulli beam.In the present formulation only nodal loads are considered, while the assemblage of different elements that form a structure follows the steps of the direct stiffness method to derive the mass and stiffness matrices.In addition, two evolution first order nonlinear differential equations are deduced for every element.The entire set of governing equations is cast into a state-space form, which is solved numerically.The validity of the proposed element is demonstrated with a cantilever beam excited at its free edge [13].
The response of the model is matched with experimental data and the parameters of hysteretic model are identified.From a mathematical standpoint, identification of the parameters of the model constitutes of a non-linear multivariate optimization problem, which is solved using Levenberg-Marquardt [3] gradient optimization algorithm.Finally, the hysteretic damping properties of the model are investigated and the rate independency is revealed, as well as the amplitude dependency.

Hysteretic model
The inelastic response of a single degree-of-freedom (SDOF) elastoplastic oscillator with kinematic hardening can be modeled using the system shown in Fig. 1.The restoring force can be expressed as the sum of the forces acting on the springs.Assuming that u is the displacement of the oscillator and using the theory of classical plasticity [12,15], it is proven that the equation of motion of this system can be described by the following system of ordinary differential equations: where dot denotes differentiation with respect to time, k represents the linear stiffness coefficient, a is the post-to pre-yielding stiffness ratio, m denotes the mass attached to the oscillator, F ext the external force and non-dimensional parameters β, γ, n control the size and shape of the hysteretic loop.The restoring force resulting from the hysteretic spring is denoted as F p , whereas its ultimate value, after which the friction slider starts to move, is denoted as F Y .The resulting hysteretic model is very similar to the Bouc-Wen model [17].In addition, the following relations, between the parameter values, must hold, so that the model is consistent with classical plasticity and also thermodynamically admissible [4]: The proposed model has been extensively used for modeling a wide range of inelastic behavior under dynamic excitation by selecting different parameter values.Recently it was shown that it also incorporates rate independent damping in its autonomous linear response [2].

Stiffness matrix formulation
Based on the dynamic properties of the spring element presented above, as well as previous considerations [7,8,16], one can assemble an inelastic beam element that also incorporates hysteretic damping.The proposed beam element is depicted in Fig. 2. The formulation follows an internal force approach, in which the end forces and displacements are related to the internal forces and displacements.
Consider the internal forces in the vertical and horizontal springs denoted as F 1 and F 2 .The end forces of the element i.e. the external shearing forces and bending moments are related to the internal ones via equilibrium relations as follows: On the other hand, the extensional deformations on the two springs can be expressed as: which constitutes a congruent relation, as the transformations are defined by transpose matrices.
The internal forces on the springs can also be expressed as functions of their deformations using the hysteretic model formulation: where the stiffness coefficients of the two springs can be selected to follow classical beam theory: The stiffness matrix for the proposed beam element results by substituting Eqs (4), ( 5) and ( 6) into relation (3) as follows: Notice that, in the linearly elastic case where a 1 = a 2 = 1, the first term of Eq. ( 7) reduces to the stiffness matrix of the standard Euler-Bernoulli beam element, while the second term vanishes as implied by relation (1).

Equation of motion
Assembling a number of elements to form a frame structure and neglecting axial deformations, the equations of motion in global coordinates are of the form: where F external denotes the external loads applied to the system.The global mass and stiffness matrices are established using the following relations: where N is the number of elements, B is a matrix incorporating the boundary conditions, containing unities for the restrained degrees of freedom and zeros for the free, A is a connectivity matrix that collects the contribution of each member stiffness coefficient and T is the transformation matrix from local to global coordinates [14].The consistent mass matrix is used for the element formulation, whereas the local stiffness matrix is computed with the aid of relation (7).Finally, the global matrix H of the hysteretic forces is established using the following relation: A T H HL HL = T T H L (10) where the local matrix H L of the hysteretic forces is computed using Eq. ( 7): In the Appendix A the M, K and H matrices that appear in the first part of the equations of motion (8), are determined for a cantilever beam consisting of two elements.
The second part of the equations of motion consists of 2N nonlinear differential equations that describe the rate of the hysteretic forces based on the Bouc-Wen model presented above.In the case of a single element, the rate of the two hysteretic forces is expressed in the following form by substituting relation (4) into Eq.(1): where the quantities F p 1Y , F p 2Y denote the ultimate values of the hysteretic forces.The system of differential Eq. ( 8) that constitutes the linear equations of motion and the nonlinear evolution of the hysteretic forces is transformed into two groups of state-space equations.The solution is obtained by solving successively at every time step for the rate of the displacement vector ẇ1 , θ1 , ẇ2 , θ2 T and then substituting these values to the rhs of the nonlinear evolution equations to determine the rate of the hysteretic forces . A Runge Kutta type of integration algorithm is implemented that uses a variable time step to achieve the required accuracy.

Description of experimental data
The validity of the proposed beam model is tested against experimental data regarding the inelastic behavior of steel beams.Two experiments related to the elastoplastic behavior of cantilever steel beams were used [10,11].These experiments were conducted at the University of Berkeley and the results are summarized in the UCB/EERC-70/03 Report [13].In Figs 14-16 some details of the test specimens and the test arrangement are shown.
In the first experiment, the inelastic behavior of a cantilever steel beam with a fully welded fixed end was studied.In the second experiment, a steel beam of same dimensions, but with a bolted fixed end, was tested.The cross section of both specimens was a typical Wide Flange US section 24WF76 with Area A = 0.0452 m 2 and area moment of inertia I = 8.7409 × 10 −4 m 4 .Their length was 2.1336 m and their nominal yield stress 248 Mpa.The tests were conducted using pseudo-static displacement-control loading.All measurements and actions were at the tip of the cantilever.

Parameter identification method
From Eqs (7) and (12) it can be observed that, except the cross section and material properties that are readily available, an additional vector of model parameters need to be determined.For the single element representation of the cantilever structure, the partial fixity at the end is incorporated in the stiffness of the element in both experiments (E = 120000 Mpa for the first experiment and 110000Mpa for the second experiment).As far as the post-to preyielding stiffness ratios (a 1 and a 2 ) are concerned, these can be obtained graphically from the experimental forcedisplacement diagrams.The ultimate values of the hysteretic forces (F p 1Y and F p 2Y ) were calculated assuming the formation of plastic hinges at the fixed ends of the steel beams according to Eqs ( 13) and ( 14), where M P represents the plastic moment of the cross-section.
Finally, the hysteretic loop parameter values ( β i , γ i , n i , i = 1, 2 ) were obtained using a nonlinear system identification algorithm [2].The unknown loop parameter values were identified, using the Levenberg -Marquardt algorithm, by matching the numerical results with the experimental data.
The identified parameter values for both specimens are summarized in Table 1.Using these parameter values and the proposed model, the experimental behavior of both specimens was simulated.The comparison between the experimental and simulated data is shown in Figs 3 and 4 from which it becomes evident that the results are in good agreement.

Features of the hysteretic damping
For a forced excitation of the model, the area enclosed by the hysteresis loop is equal to the energy dissipated over a cycle [6].On the other hand, in the case of viscous damping, the energy dissipated over a cycle E D is equal to: where c is the viscous damping coefficient, Ω the excitation frequency and X 0 denotes the steady-state amplitude of the response.
In the particular case of resonance, the energy dissipation can be expressed as a function of the damping ratio ζ and the stiffness k as: By equating the energy dissipated by a viscously damped system to that of a hysteretic system at resonance, one can define the equivalent viscous damping ratio as: where E Dhyst denotes the energy dissipated over a cycle of a system with hysteretic damping at resonance and is equal to the area enclosed by the hysteresis loop.The notion of equivalent damping ratio was applied to the beam element models of the two specimens.Numerical, displacement -controlled experiments were conducted by solving the equations of motion for each specimen using the identified model parameters, as shown in Eqs ( 18) and ( 19).
Initially, each model is excited at various frequencies at constant amplitude.By calculating numerically the area of the resulting steady-state hysteretic loops, an estimation of the dissipated energy per cycle of oscillation is obtained.The evolution of the dissipated energy as a function of the excitation frequency for the two models is presented in Fig. 5.By examining this figure, one can observe that the dissipated energy is independent of the excitation frequency (any fluctuations shown are due to the numerical error in estimating the loop area); therefore, the rate independency of the proposed model is verified numerically.Moreover, the hysteretic behavior of the model is examined with respect to the amplitude dependency of the dissipated energy.This is accomplished by another series of numerical experiments conducted at fixed frequencies and over a range of excitation amplitudes.Each model is excited at its first natural frequency.The amplitude range is chosen in such a way that none of the two models enters into the plastic regime.The energy dissipation for each case is calculated numerically and its variation with respect to the amplitude is presented in Fig. 6.
In addition, if an equivalent viscous damping formulation is required then using Eq. ( 15), based on the fact that the systems were excited at resonance (17), an estimation of the equivalent damping ratio can be deduced.

Discretization with more macro-elements
The accuracy of the predictions of different models with one, four and eight macro elements is presented for the same experiment on the cantilever with a sinusoidal tip force of 300 kN at a frequency of 90 Hz.The identified parameters of the welded specimen are used for all the different models.In Fig. 7 the tip displacement with respect of time is presented from which a slight variation in the amplitude is observed.This variation is plotted for a short period of time in detail.In Fig. 8, the force-displacement diagrams of the horizontal and vertical hysteretic springs are plotted for the discretization with eight macro-elements.It is observed that the inelastic behavior is attributed mainly to the horizontal spring at the clamped end, whereas the three others behave almost elastically, exhibiting small hysteresis, as demonstrated in the detailed plot in Fig. 9.
The first six natural frequencies and mode shapes of the cantilever are plotted in Fig. 10.They are in good agreement with the periodogram of Fig. 11, which is derived from a numerical experiment conducted for the autonomous case and for the discretization with eight elements.
From the above comparisons it becomes evident that the proposed macro-element is capable of capturing the overall inelastic behavior of the structure.The use of more elements though, is important in the sense that more elements can reveal unknown plastic regions depending on the loading, especially in indetermined structures.
The experimental validation of the model, previously presented, was based on the experimental measurements at the tip of the cantilever.As demonstrated clearly by the experiments, the different types of connection i.e. welded or bolted, absorb different amounts of energy.Therefore, a realistic modeling of the structure requires at least one connection element and a number of elastoplastic beam elements.If measurements were available at sections just after the connection and other sections along the beam, a proper identification of the parameters of both the connection and the beam elements would have been possible.This more realistic modeling would have offered representative stiffness properties for any type of loading, whereas with only the tip measurements, average properties of the connection and the beam are implemented that correspond to these particular experiments.
In addition, a different modeling with one connection inelastic element and seven elastic beam elements with nominal properties, three inelastic elements and four elastic ones and eight inelastic elements is considered.All models are excited with a sinusoidal load at a frequency of 90 Hz, close to the first natural frequency of the elastic beam.In Fig. 12 the tip displacement of the three models is plotted for a small time window.From this it is observed that the stiffer model, with the one inelastic element, deviates from the response of the other two models which are in good agreement between them.In Fig. 13 the hysteretic loops of the horizontal spring at the clamped end of the three models is presented.It is observed that most of the energy is dissipated at the clamped end.Moreover, the one plastic element model dissipated all the energy at the clamped end, whereas the other two models dissipated most of the energy at the clamped end and the rest in the remaining plastic elements.It is observed that less energy is dissipated at the clamped end in the models with increased number of plastic elements.
However, as it is shown in this numerical example, even with one set of parameters for the whole beam, the model can predict where plasticity occurs and simulate it adequately.

5.
In this work, a generic inelastic beam element with rate independent damping proposed.For the inelastic behavior of steel beams, aid of nonlinear identification methods, it was shown that proposed model is capable to simulate experimental data adequately.In addition, numerical results revealed that the energy dissipation is rate independent and amplitude dependent.the proposed model simulates the beam dynamic behavior both in linear and nonlinear case, while the stiffness of the element reduces to that of Euler-Bernoulli beam for the elastic case.
For E = 1.965e +11 Pa (modulus of elasticity), I = 8.7409e −4 (area moment of inertia), a 1 = a 2 = 0.1 and L = 1.0668 m (length of each element) the local stiffness matrix is computed as: The stifness matrix KL is equal K L since the current example: The stifness matrix K G is computed with the aid of the connectivity matrix A. procedure is depicted in relation (23): The stifness matrix K is computed with the aid of the boundary conditions matrix B: 0 0 0 0 0 0 0 0 1 0 0 0 0 1 0 0 0 0 1 0 0 0 0 1 In the diagonal of the B matrix, zeros are used for the fixed and ones for the free degrees of freedom of the structure.The final stifness matrix K of the supported structure is: The same procedure is followed for the mass matrix M of the model.The local mass matrix is given as: where ρ is the material density and A the cross section area.The local matrix of hysteretic forces is described by Eq. (11).For L = 1.0668 m (length of each element) and h = 0.60706 m (height of each element) the local matrix of hysteretic forces is computed: The matrix of hysteretic forces HL is equal to H L in the current example.The matrix of hysteretic forces H G is computed with the aid of the connectivity matrix A H .This procedure is determined in relation ( 28):

Fig. 5 .
Fig. 5. Dissipated Energy per cycle of oscillation with respect to excitation frequency.

Fig. 10 .Fig. 13 .
Fig. 10.Natural frequencies and mode shapes of the specimen of the first experiment

Table 1
Identified model parameter values for the two experiments