# Nonlinear Resonance Analysis of Slender Portal Frames under Base Excitation

Authors: Luis Fernando Paullo Muñoz; Paulo B. Gonçalves; Ricardo A. M. Silveira; Andréa Silva
Journal: Shock and Vibration
 (2017)

---

## Abstract

The dynamic nonlinear response and stability of slender structures in the main resonance regions are a topic of importance in structural analysis. In complex problems, the determination of the response in the frequency domain indirectly obtained through analyses in time domain can lead to huge computational effort in large systems. In nonlinear cases, the response in the frequency domain becomes even more cumbersome because of the possibility of multiple solutions for certain forcing frequencies. Those solutions can be stable and unstable, in particular saddle-node bifurcation at the turning points along the resonance curves. In this work, an incremental technique for direct calculation of the nonlinear response in frequency domain of plane frames subjected to base excitation is proposed. The transformation of equations of motion to the frequency domain is made through the harmonic balance method in conjunction with the Galerkin method. The resulting system of nonlinear equations in terms of the modal amplitudes and forcing frequency is solved by the Newton-Raphson method together with an arc-length procedure to obtain the nonlinear resonance curves. Suitable examples are presented, and the influence of the frame geometric parameters and base motion on the nonlinear resonance curves is investigated.


---

## Body

## 1. Introduction

Lightweight framed structures are extensively used in engineering applications such as buildings, industrial constructions, off-shore platforms, and aerospace structures. In civil engineering, framed metal structures such as gabled frames are efficient structural forms to withstand various loads such as live, snow, wind, earthquake, and crane loads. Linear models are usually used for the design of these structures. However, as its slenderness increases, the degree of geometric nonlinearity increases as well as the occurrence of instability phenomena. Under dynamic loads, slender frames may be subjected to unwanted (or unsafe) large vibration amplitudes and accelerations and dynamic jumps in the main resonance regions due to hardening or softening frequency-amplitude relation and modal couplings. Therefore, an accurate nonlinear dynamic analysis is required. Recent earthquakes have increased the interest in the development of reliable analytical methods of assessing the dynamic response of existing structures and, if necessary, developing effective retrofit strategies. So, understanding the nonlinear dynamics of framed structures under base excitation constitutes an essential step for their safe design.The seminal works by Koiter [1] and Roorda [2] have shown the complex postbucking paths exhibited by slender frames; subsequently, several authors have studied the nonlinear behavior of reticulated structures liable to buckling under prescribed loading conditions. A compilation of these first results can be found in Bazant and Cedolin [3] who summarize many topics about the stability of slender frames. Design prescriptions are found, for example, in Galambos [4]. More recently, Silvestre and Camotim [5] presented a semianalytical method to analyze the geometrically nonlinear behavior of plane frames, including the influence of eventual buckling mode interaction phenomena. Later, they presented the results concerning the elastic in-plane stability and second-order behavior of unbraced single-bay pitched-roof steel frames and proposed a methodology to design this type of commonly used structure. In particular, they showed that, due to the rafter slope, the geometrically nonlinear behaviors of orthogonal beam-and-column and pitched-roof frames are qualitatively different [6]. The dynamics of such frames will be investigated in this work. More recently, Galvão et al. [7] assessed the effect of geometric imperfections and boundary conditions on the postbuckling response and imperfection-sensitivity of L-frames.The dynamic response of slender frames in the main resonance regions is an important topic in the analysis of these structures under time varying loads. Several nonlinear phenomena have been studied in the last decades. Simitses [8] studied the effect of static preloading and suddenly applied loads on the nonlinear vibrations and instabilities of slender structures. Chan and Ho [9] conducted the nonlinear vibration analysis of steel frames with flexible connections. Mazzilli and Brasil [10] presented an analytical study of the nonlinear vibrations in a three-time redundant portal frame, considering the effect of axial forces upon the natural frequencies. The axial forces play an important role in tuning the sway mode and the first symmetrical mode into a 1 : 2 internal resonance. Harmonic support excitations resonant with the first symmetrical mode are then introduced and the amplitudes of nonlinear steady states are computed based upon a multiple scales solution. Later, Soares and Mazzilli [11] described the implementation of a computer program to calculate nonlinear normal modes of structural systems and generate individual modes of planar framed structures exhibiting geometrically nonlinear behavior. The procedure follows the invariant manifold approach, adapted to handle equations of motion of systems discretized by finite element techniques. Chan and Chui [12] examined, based on their previous works, the nonlinear static and cyclic behavior of steel frames with semirigid connections. McEwan et al. [13] proposed a method for modeling the large deflection beam response involving multiple vibration modes, while Ribeiro [14] analyzed the geometrically nonlinear vibrations of beams and plane frames by the hierarchical finite element method and examined the suitability of the proposed formulation for time domain nonlinear analyses. Da Silva et al. [15] studied the nonlinear dynamics of a low-rise portal frame using the ANSYS finite element software. The results show the influence of semirigid joints and geometrical nonlinearity on the steel frames dynamics. More recently, Galvão et al. [16] investigated the effect of semirigid connections on the nonlinear vibrations of slender frames and obtained the nonlinear relation between the natural frequencies and static preloading. Su and Cesnik [17] introduced a strain-based geometrically nonlinear beam formulation for structural and aeroelastic modeling and analysis of slender wings of very flexible aircraft. Solutions of different beam configurations under static loads and forced dynamic excitations are compared against ones from other geometrically nonlinear beam formulations. Gonçalves et al. [18] published an experimental analysis of the nonlinear vibrations of a slender metal column under self-weight. Masoodi and Moghaddam [19] studied the nonlinear dynamics and natural frequencies of gabled frames having flexible restraints and connections. To control unwanted nonlinear phenomena and vibrations of frames, several control techniques are proposed in literature. For example, Palacios Felix et al. [20] examined the nonlinear control method based on the saturation phenomenon. The interaction of a nonideal source with a portal frame leads to the occurrence of interesting nonlinear phenomena during the passage through resonance when the nonideal excitation frequency is near the portal frame natural frequency and when the nonideal system frequency is approximately twice the controller frequency (two-to-one internal resonance). They also suggest the use of the portal frame for energy harvesting.The determination of a structure’s response in the frequency domain, characterized by the resonance curves and bifurcation diagrams, is, in some cases, indirectly obtained through a series of analyses in time domain, varying step by step a selected control parameter (see Parker and Chua [21]). In this brute force approach, for each control parameter value, the equations of motion must be integrated long enough to reach the steady-state regime. This leads to a huge computational effort when analyzing structures with a large number of degrees of freedom [16]. In nonlinear cases, the response in the frequency domain becomes even more cumbersome because of the possibility of multiple solutions for certain forcing frequencies. Those solutions can be stable and unstable, and bifurcations can appear, in particular saddle-node bifurcation at the turning points along the resonance curves. To predict the nonlinear response in frequency domain, some methods are used, such as the harmonic balance method (HBM) and other integral transformations. Based on the HBM, Nakhla and Valch [22] proposed a method for the solution of nonlinear periodic networks, avoiding the time domain solution of the dynamic equations. Lau et al. [23] formulated an incremental variational process based on the HBM to study the nonlinear vibrations of elastic systems. Further, Cheung and Chen [24] used the formulation of the incremental harmonic balance method (IHB) to solve a system of differential equations with cubic nonlinearity, which governs a wide range of engineering problems such as large-amplitude vibration of beams or plates. An incremental arc-length method combined with a cubic extrapolation technique is adopted to trace the response curves. Chen et al. [25] generalized this technique to use the incremental HBM in a finite element context, being able to use the HBM in the study of systems with several degrees of freedom. Due to its simple formulation and relatively easy implementation, the HBM is one of the most popular methods to study the nonlinear vibrations of structures [26]. The algebraic system of nonlinear equations, which results from the HBM, requires a numerical procedure to obtain the solutions. Frequently, Newton’s method is adopted to solve the nonlinear system. However, in many cases, the resonance curves exhibit limit points and, in such cases, the use of continuation methods is necessary. For example, Ribeiro and Petyt [27] studied the free and steady-state forced vibrations of thin plates. Symbolic computation is employed in the derivation of the model and the equations of motion are solved by the Newton and continuation methods. Von Groll and Ewins [28] described a numerical algorithm based on the harmonic balance method to study rotor/stator interaction dynamics under periodic excitation. The algorithm also calculates turning points and follows solution branches via arc-length continuation. More recently, Ferreira and Serpa [29] described the application of the arc-length method to solve a system of nonlinear equations obtaining as a result the nonlinear frequency response. In a recent work, Londoño et al. [30] used continuation methods in frequency domain to obtain the backbone curves of several nonlinear systems, while Renson et al. [31] presented numerical methodologies for the computation of nonlinear normal modes in mechanical system.Stoykov and Ribeiro [32, 33] investigated the geometrically nonlinear free vibrations of beams using a p-version finite element method. The variation of the amplitude of vibration with the frequency of vibration is determined and presented in the form of backbone curves. Coupling between modes is investigated, internal resonances are found, and the ensuing multimodal oscillations are described. Formica et al. [34] present a computational framework to perform parameter continuation of periodic solutions of nonlinear, distributed-parameter systems represented by partial differential equations with time-dependent coefficients and excitations. As a case study, they consider a nonlinear beam subject to a harmonic excitation.An important design concern is the dynamics of slender frames subjected to seismic loads. In such cases, the response of the structure in frequency domain is important, since the vulnerability of structure during an earthquake is related to the relation between the natural vibration frequencies of the structure and the frequency content of the seismic load, as studied in Paullo Muñoz et al. [35] who assess the frequency domain response of slender frames under seismic excitation. In seismic events, the ground acceleration magnitude is a crucial parameter, since the effect of the excitation on the structure is directly related to acceleration magnitude. Most recorded events have a magnitude lower than 1g. However, peak acceleration greater than 1g has been registered. For example, an acceleration magnitude of 1.7g was registered for the Los Angeles earthquake, and a peak acceleration of 2.99g was registered for the 2011 Tohoku earthquake [36]. Douglas [37] presented a review of equations for the estimation of peak ground acceleration using strong-motion records, while Katsanos et al. [38] reviewed alternative selection procedures for incorporating strong ground motion records within the framework of seismic design of structures. There are many other research areas where the influence of base acceleration on structural components is important, including subsurface explosions, rotating machinery, and launching vehicles.In this work, an incremental technique for the direct calculation of the nonlinear dynamic response in frequency domain of nonlinear plane frames discretized by the finite element method and subjected to a base excitation is proposed. The transformation of discretized equations of motion, in the finite element context, to the frequency domain is accomplished here through the classical harmonic balance method (HBM). For the nonlinear analysis, a particular adaptation of the HBM-Galerkin methodology presented by Cheung and Chen [24] and generalized for the use in FEM context by Chen et al. [25] is proposed here. The resulting system of nonlinear equations in terms of the modal amplitudes and forcing frequency is solved by the Newton-Raphson method together with an arc-length procedure to obtain the nonlinear resonance curves. Examples of commonly used frame geometries are presented and the influence of vertical and horizontal base motions on the frame nonlinear response as a function of the frame geometric parameters is analyzed.


## 2. Formulation

### 2.1. Equations of Motion

For a structural system under the action of harmonic forcing with a prescribed forcing frequencyΩ, the nonlinear equations of motion can be written as(1)Mu¨t+Cu˙t+Fiut=Fcos⁡Ωtf,where M, C, and Fi are the mass matrix, damping matrix, and the vector of nonlinear elastic forces, respectively, with Fi being dependent on the nodal displacement vector u(t) and on the geometric nonlinear formulation; f is the vector of external applied loads at the nodal points, which gives the amount of forces at each DOF, and F and Ω are, respectively, the magnitude and frequency of the harmonic excitation.


### 2.2. Harmonic Balance Method (HBM) for Linear Analysis

The HBM is one of the most popular methods for the nonlinear analysis of dynamical systems. In this method, the displacement field is approximated by a finite Fourier series of the form(2)ut=∑n=0NAncos⁡nΩt+Bnsin⁡nΩt,where An and Bn are the modal amplitude vectors corresponding to the nth harmonic and N is the number of terms considered in the approximation. For the linear damped case, it is only necessary to retain the second term in (2). Thus, the nodal displacement vector takes the form(3)ut=A1cos⁡Ωt+B1sin⁡Ωt.Introducing (3) into (1) results in(4)ΩCA1-Ω2MB1cos⁡Ωt-Ω2MA1+ΩCB1sin⁡Ωt=Ffcos⁡Ωt.TakingFf=F0 and Fi=Kut, where K is the linear stiffness matrix, and applying the HBM to equation (4), the following system is obtained:(5)C-ΩM+K-Ω2M+K-ΩCA1B1=F00.Equation (5) can be expressed in a more compact form as(6)K-ΩD--F-=0,where(7)K-Ω=ΩC-Ω2M+K-Ω2M+K-ΩC,D-=A1B1,F-=F00.As a result of the transformation of the system of equations to the frequency domain and the use of continuation techniques together with an arc-length constraint, the forcing frequencyΩ and the modal amplitude vector D- are the unknowns in (7) [29, 39, 40]. Even in the linear case, the transformation of (1) to the frequency domain leads to a nonlinear system of algebraic equations, with the use of nonlinear techniques such as the Newton-Raphson method being necessary.


### 2.3. HBM-Galerkin Methodology for Nonlinear Analysis

For the nonlinear analysis, considering that the steady-state response is periodic, it is convenient to use a periodic nondimensional time variable defined as(8)τ=Ωt.Substituting (8) into (1), the equation of motion can be rewritten as (9)Ω2Mu¨τ+ΩCu˙τ+Fiuτ=Fcos⁡τfand each component of the displacement vector can be approximated by(10)uiτ=∑n=0Nai,ncos⁡τ+bi,nsin⁡τ.For a damped system with quadratic nonlinearity, at least the first two terms of the Fourier series must be considered [39]:(11)uiτ=ai,0+bi,1cos⁡τ+ci,1sin⁡τ.For a damped system with only cubic nonlinearity, it is necessary to consider at least the following terms [39]:(12)uiτ=ai,1cos⁡τ+bi,1sin⁡τ.Plane frames exhibit quadratic and cubic nonlinearities.Define(13)ui=Cτd-iwith(14)Cτ=1,sinτ,cos⁡τ,d-i=ai,0,bi,1,ci,1T.Consider the following relation:(15)uτ=SτD-,where(16)Sτ=Cτ000⋯00Cτ00⋯000Cτ0⋯0000Cτ⋯⋮⋮⋮⋮⋮⋱0000⋯0CτN×3N,D-=d-1,d-2,…,d-N1×3NT.N is total number of degrees of freedom. Introducing (15) into (9), the following matrix equation is obtained:(17)Ω2MS¨τD-+ΩCS˙τD-+FiSτD-=F0cos⁡τ.Considering that the solution is periodic, multiplying both sides of (17) by S(τ)T and integrating the resulting expression over a period, one obtains [40](18)∫02πSτTΩ2MS¨τD-+ΩCS˙τD-+FiSτD-dτ=∫02πSτTF0cos⁡τdτ.In a compact form, (18) can be expressed as(19)Ω2M-D-+ΩC-D-+Fi-D-=F-,where(20)M-=∫02πSτTMS¨τ·dτ,C-=∫02πSτTCS˙τ·dτF-=∫02πSτTF0cos⁡τdτ,Fi-D-=∫02πSτTFiSτ·D-dτ.Equation (19) is a nonlinear system of algebraic equations, in which the unknowns are the forcing frequency Ω and the modal amplitude vector D-. The internal force vector in frequency domain Fi-D- is obtained from the transformation of Fiu(τ). In this context, an explicit expression for Fiu(τ) is necessary. The next step is the solution of the nonlinear system of algebraic equations.


### 2.4. Solution of the Nonlinear System of Equations

The transformation of the equations of motion from time to frequency domain results in a nonlinear system of algebraic equations both for the linear and nonlinear cases, as shown by (6) and (19), with the unknown variables being the forcing frequency and the modal amplitudes. In this section, the technique to solve the nonlinear system of equations must consider the possibility of frequency and amplitude limit points. In this context, (6) can be defined in a general form as(21)RΩ,D-=0.Through an incremental analysis, the dynamic equilibrium in theith step is given by(22)RΩi,D-i=0.The unknown variables in theith step are obtained by the incremental analysis as(23)Ωi=Ωi-1+ΔΩi,D-i=D-i-1+ΔD-i.Introducing (22) into (21) results in(24)RΩi-1+ΔΩi,D-i-1+ΔD-i=0.The frequency and amplitude incrementsΔΩi and ΔA-i can be calculated iteratively as(25)ΔΩik=ΔΩik-1+δΩik,ΔD-ik=ΔD-ik-1+δD-ik,where δΩik and δD-ik are the correctors which can be obtained through the first variation of (24) as(26)δRk=∂Rk-1∂D-iδD-ik+∂Rk-1∂ΩiδΩikor, in a more compact form, as(27)Km-k-1δD-ik+δΩikf-k-1=δRk.To solve (27), an additional constraint equation is necessary. In this work, a spherical arc-length constraint is used [41, 42]:(28)ΔD-iTΔD-i+ΔΩik2f-k-1Tf-k-1-Δli2=0.From (24), (25), (27), and (28), it is possible to obtain the iterative frequency corrector as(29)δΩik=±Δli2vTv+f-k-1Tf-k-1fork=1-rkTvkvkTvkfork>1,where(30)vk=Km-k-1-1f-k-1,rk=Km-k-1-1Rk-1.The signal of the first frequency corrector, also called frequency predictor, can be obtained using the positive work criterion as(31)sign⁡δΩi1=sign⁡vkT·f-k-1.For linear systems,Km- and f- can be defined by the following expressions:(32)Km-=ΩC-Ω2M+K-Ω2M+K-ΩC,f-=C-2ωM-2ωM-CD-.Analogously, for the nonlinear case,Km- and f- can be obtained from the following relations:(33)Km-=Ω2M-+ΩC-+∂Fi-D-∂D-,f-=2ΩM-+C-·D-.


### 2.5. Geometric Nonlinear Consideration

An explicit expression forFiu(τ) is necessary. Considering an incremental analysis, the nodal displacement vector can be expressed as(34)uti=uti-1+Δut.Analogously, the internal forces vector can be calculated as(35)Fiti=Fiti-1+ΔFit.F
i can be derived once the kinematic relations and constitutive law are defined. These relations depend on the adopted nonlinear strain measure formulation [7, 16]. In this paper, only the geometric nonlinearity is considered. In this context, the internal force increment vector can be approximated as(36)ΔFint=KL+KNLΔu,where KL is the linear elastic stiffness matrix and KNL is the nonlinear stiffness matrix. Considering a linearized formulation, the nonlinear stiffness matrix for an Euler-Bernoulli beam-column plane element in a local system of reference can be defined as [7, 16](37)KNL=PLM1+M2L20-PL-M1+M2L20M1+M2L26P5LP10-M1+M2L2-6P5LP100P10PL150-P10-PL30-PL-M1+M2L20PLM1+M2L20-M1+M2L2-6P5L-P10M1+M2L26P5L-P100P10-PL300-P10PL15,where P is the average internal axial force and M1 and M2 are the internal bending moments in the initial and the end nodes, respectively. For a prismatic bean-column element, the internal forces can be expressed as a function of the axial and transversal displacements, respectively, ux and vx, as(38)P=EA∂ux∂x+1L∫0L∂vx∂x2·dxM=EI∂2vx∂x2,where E is Young’s modulus, A is the cross-sectional area, and I is the moment of inertia. The displacement fields can be expressed by interpolation of nodal displacement components as(39)ux=∑i=12Hixui,vx=∑i=36Hixui,where Hi(x) is the shape function used for the finite element discretization. Analogous to (15), the increment of the nodal displacement vector can by calculated as(40)Δut=SτΔD-.Solving (15), (35), (36), and (40), the nodal internal force vector Fi at the ith increment can be obtained as(41)Fii=KLSτD-+Fii-1+KNLSτΔD-.Finally, the nodal internal force vector in frequency domain defined in (20) can be calculated as(42)Fi-iD-=K-LD-i+K-NLΔD-+Fi-i-1D-,where(43)K-L=∫02πSτTKLSτdτK-NL=∫02πSτTKNLSτdτFi-i-1D-=∫02πSτTFNLi-1dτ.In linear case,Fi-D-=K-LD-, and the resultant nonlinear system of equations defined by (19) is the same as the system obtained through the classical HBM (see (5)).


### 2.6. Equivalent Rotation Matrix

In order to assemble the global system of equations, the matricesM, C, KL, and KNL must be defined in global coordinates. The matrices M, C, and KL are the same for any reference system, since they are constant. On the other hand, KNL, defined by (37), is formulated in a local reference frame. Thus, its rotation to the global reference system is necessary. It is defined in global reference frame by(44)K-NL=∫02πSτTTTKNLTSτdτ,where T is a rotation matrix. The rotation matrix for a plane beam-column element can be found in [15]. In terms of the computational implementation, it is not convenient to rotate the stiffness matrix KNL before frequency domain transformation, since the elements of KNL have large expression which depend on nodal displacement, even in the local reference frame. So, it is convenient to rotate KNL after the transformation of the equations to the frequency domain, using an equivalent rotation matrix. In this case, it takes the form (45)K-NL=T′T∫02πSτTKNLSτdτT′,where T′ is the equivalent rotation matrix, which satisfies the following condition:(46)T′TSτT=SτTTT.Taking into account the block-diagonal characteristic of Sτ and since T is an orthogonal matrix, it is possible to find an explicit expression for T′ as a function of T through the following expression:(47)T′=t1,1′…t1,NGL′⋮⋱⋮tNGL,1′…tNGL,NGL′, where NGL is the number of degrees of freedom and ti,j′ is a submatrix given by(48)ti,j′=Ti,jI3,where I3 is an identity matrix of order three.


## 2.1. Equations of Motion

For a structural system under the action of harmonic forcing with a prescribed forcing frequencyΩ, the nonlinear equations of motion can be written as(1)Mu¨t+Cu˙t+Fiut=Fcos⁡Ωtf,where M, C, and Fi are the mass matrix, damping matrix, and the vector of nonlinear elastic forces, respectively, with Fi being dependent on the nodal displacement vector u(t) and on the geometric nonlinear formulation; f is the vector of external applied loads at the nodal points, which gives the amount of forces at each DOF, and F and Ω are, respectively, the magnitude and frequency of the harmonic excitation.


## 2.2. Harmonic Balance Method (HBM) for Linear Analysis

The HBM is one of the most popular methods for the nonlinear analysis of dynamical systems. In this method, the displacement field is approximated by a finite Fourier series of the form(2)ut=∑n=0NAncos⁡nΩt+Bnsin⁡nΩt,where An and Bn are the modal amplitude vectors corresponding to the nth harmonic and N is the number of terms considered in the approximation. For the linear damped case, it is only necessary to retain the second term in (2). Thus, the nodal displacement vector takes the form(3)ut=A1cos⁡Ωt+B1sin⁡Ωt.Introducing (3) into (1) results in(4)ΩCA1-Ω2MB1cos⁡Ωt-Ω2MA1+ΩCB1sin⁡Ωt=Ffcos⁡Ωt.TakingFf=F0 and Fi=Kut, where K is the linear stiffness matrix, and applying the HBM to equation (4), the following system is obtained:(5)C-ΩM+K-Ω2M+K-ΩCA1B1=F00.Equation (5) can be expressed in a more compact form as(6)K-ΩD--F-=0,where(7)K-Ω=ΩC-Ω2M+K-Ω2M+K-ΩC,D-=A1B1,F-=F00.As a result of the transformation of the system of equations to the frequency domain and the use of continuation techniques together with an arc-length constraint, the forcing frequencyΩ and the modal amplitude vector D- are the unknowns in (7) [29, 39, 40]. Even in the linear case, the transformation of (1) to the frequency domain leads to a nonlinear system of algebraic equations, with the use of nonlinear techniques such as the Newton-Raphson method being necessary.


## 2.3. HBM-Galerkin Methodology for Nonlinear Analysis

For the nonlinear analysis, considering that the steady-state response is periodic, it is convenient to use a periodic nondimensional time variable defined as(8)τ=Ωt.Substituting (8) into (1), the equation of motion can be rewritten as (9)Ω2Mu¨τ+ΩCu˙τ+Fiuτ=Fcos⁡τfand each component of the displacement vector can be approximated by(10)uiτ=∑n=0Nai,ncos⁡τ+bi,nsin⁡τ.For a damped system with quadratic nonlinearity, at least the first two terms of the Fourier series must be considered [39]:(11)uiτ=ai,0+bi,1cos⁡τ+ci,1sin⁡τ.For a damped system with only cubic nonlinearity, it is necessary to consider at least the following terms [39]:(12)uiτ=ai,1cos⁡τ+bi,1sin⁡τ.Plane frames exhibit quadratic and cubic nonlinearities.Define(13)ui=Cτd-iwith(14)Cτ=1,sinτ,cos⁡τ,d-i=ai,0,bi,1,ci,1T.Consider the following relation:(15)uτ=SτD-,where(16)Sτ=Cτ000⋯00Cτ00⋯000Cτ0⋯0000Cτ⋯⋮⋮⋮⋮⋮⋱0000⋯0CτN×3N,D-=d-1,d-2,…,d-N1×3NT.N is total number of degrees of freedom. Introducing (15) into (9), the following matrix equation is obtained:(17)Ω2MS¨τD-+ΩCS˙τD-+FiSτD-=F0cos⁡τ.Considering that the solution is periodic, multiplying both sides of (17) by S(τ)T and integrating the resulting expression over a period, one obtains [40](18)∫02πSτTΩ2MS¨τD-+ΩCS˙τD-+FiSτD-dτ=∫02πSτTF0cos⁡τdτ.In a compact form, (18) can be expressed as(19)Ω2M-D-+ΩC-D-+Fi-D-=F-,where(20)M-=∫02πSτTMS¨τ·dτ,C-=∫02πSτTCS˙τ·dτF-=∫02πSτTF0cos⁡τdτ,Fi-D-=∫02πSτTFiSτ·D-dτ.Equation (19) is a nonlinear system of algebraic equations, in which the unknowns are the forcing frequency Ω and the modal amplitude vector D-. The internal force vector in frequency domain Fi-D- is obtained from the transformation of Fiu(τ). In this context, an explicit expression for Fiu(τ) is necessary. The next step is the solution of the nonlinear system of algebraic equations.


## 2.4. Solution of the Nonlinear System of Equations

The transformation of the equations of motion from time to frequency domain results in a nonlinear system of algebraic equations both for the linear and nonlinear cases, as shown by (6) and (19), with the unknown variables being the forcing frequency and the modal amplitudes. In this section, the technique to solve the nonlinear system of equations must consider the possibility of frequency and amplitude limit points. In this context, (6) can be defined in a general form as(21)RΩ,D-=0.Through an incremental analysis, the dynamic equilibrium in theith step is given by(22)RΩi,D-i=0.The unknown variables in theith step are obtained by the incremental analysis as(23)Ωi=Ωi-1+ΔΩi,D-i=D-i-1+ΔD-i.Introducing (22) into (21) results in(24)RΩi-1+ΔΩi,D-i-1+ΔD-i=0.The frequency and amplitude incrementsΔΩi and ΔA-i can be calculated iteratively as(25)ΔΩik=ΔΩik-1+δΩik,ΔD-ik=ΔD-ik-1+δD-ik,where δΩik and δD-ik are the correctors which can be obtained through the first variation of (24) as(26)δRk=∂Rk-1∂D-iδD-ik+∂Rk-1∂ΩiδΩikor, in a more compact form, as(27)Km-k-1δD-ik+δΩikf-k-1=δRk.To solve (27), an additional constraint equation is necessary. In this work, a spherical arc-length constraint is used [41, 42]:(28)ΔD-iTΔD-i+ΔΩik2f-k-1Tf-k-1-Δli2=0.From (24), (25), (27), and (28), it is possible to obtain the iterative frequency corrector as(29)δΩik=±Δli2vTv+f-k-1Tf-k-1fork=1-rkTvkvkTvkfork>1,where(30)vk=Km-k-1-1f-k-1,rk=Km-k-1-1Rk-1.The signal of the first frequency corrector, also called frequency predictor, can be obtained using the positive work criterion as(31)sign⁡δΩi1=sign⁡vkT·f-k-1.For linear systems,Km- and f- can be defined by the following expressions:(32)Km-=ΩC-Ω2M+K-Ω2M+K-ΩC,f-=C-2ωM-2ωM-CD-.Analogously, for the nonlinear case,Km- and f- can be obtained from the following relations:(33)Km-=Ω2M-+ΩC-+∂Fi-D-∂D-,f-=2ΩM-+C-·D-.


## 2.5. Geometric Nonlinear Consideration

An explicit expression forFiu(τ) is necessary. Considering an incremental analysis, the nodal displacement vector can be expressed as(34)uti=uti-1+Δut.Analogously, the internal forces vector can be calculated as(35)Fiti=Fiti-1+ΔFit.F
i can be derived once the kinematic relations and constitutive law are defined. These relations depend on the adopted nonlinear strain measure formulation [7, 16]. In this paper, only the geometric nonlinearity is considered. In this context, the internal force increment vector can be approximated as(36)ΔFint=KL+KNLΔu,where KL is the linear elastic stiffness matrix and KNL is the nonlinear stiffness matrix. Considering a linearized formulation, the nonlinear stiffness matrix for an Euler-Bernoulli beam-column plane element in a local system of reference can be defined as [7, 16](37)KNL=PLM1+M2L20-PL-M1+M2L20M1+M2L26P5LP10-M1+M2L2-6P5LP100P10PL150-P10-PL30-PL-M1+M2L20PLM1+M2L20-M1+M2L2-6P5L-P10M1+M2L26P5L-P100P10-PL300-P10PL15,where P is the average internal axial force and M1 and M2 are the internal bending moments in the initial and the end nodes, respectively. For a prismatic bean-column element, the internal forces can be expressed as a function of the axial and transversal displacements, respectively, ux and vx, as(38)P=EA∂ux∂x+1L∫0L∂vx∂x2·dxM=EI∂2vx∂x2,where E is Young’s modulus, A is the cross-sectional area, and I is the moment of inertia. The displacement fields can be expressed by interpolation of nodal displacement components as(39)ux=∑i=12Hixui,vx=∑i=36Hixui,where Hi(x) is the shape function used for the finite element discretization. Analogous to (15), the increment of the nodal displacement vector can by calculated as(40)Δut=SτΔD-.Solving (15), (35), (36), and (40), the nodal internal force vector Fi at the ith increment can be obtained as(41)Fii=KLSτD-+Fii-1+KNLSτΔD-.Finally, the nodal internal force vector in frequency domain defined in (20) can be calculated as(42)Fi-iD-=K-LD-i+K-NLΔD-+Fi-i-1D-,where(43)K-L=∫02πSτTKLSτdτK-NL=∫02πSτTKNLSτdτFi-i-1D-=∫02πSτTFNLi-1dτ.In linear case,Fi-D-=K-LD-, and the resultant nonlinear system of equations defined by (19) is the same as the system obtained through the classical HBM (see (5)).


## 2.6. Equivalent Rotation Matrix

In order to assemble the global system of equations, the matricesM, C, KL, and KNL must be defined in global coordinates. The matrices M, C, and KL are the same for any reference system, since they are constant. On the other hand, KNL, defined by (37), is formulated in a local reference frame. Thus, its rotation to the global reference system is necessary. It is defined in global reference frame by(44)K-NL=∫02πSτTTTKNLTSτdτ,where T is a rotation matrix. The rotation matrix for a plane beam-column element can be found in [15]. In terms of the computational implementation, it is not convenient to rotate the stiffness matrix KNL before frequency domain transformation, since the elements of KNL have large expression which depend on nodal displacement, even in the local reference frame. So, it is convenient to rotate KNL after the transformation of the equations to the frequency domain, using an equivalent rotation matrix. In this case, it takes the form (45)K-NL=T′T∫02πSτTKNLSτdτT′,where T′ is the equivalent rotation matrix, which satisfies the following condition:(46)T′TSτT=SτTTT.Taking into account the block-diagonal characteristic of Sτ and since T is an orthogonal matrix, it is possible to find an explicit expression for T′ as a function of T through the following expression:(47)T′=t1,1′…t1,NGL′⋮⋱⋮tNGL,1′…tNGL,NGL′, where NGL is the number of degrees of freedom and ti,j′ is a submatrix given by(48)ti,j′=Ti,jI3,where I3 is an identity matrix of order three.


## 3. Numerical Examples

### 3.1. Linear Analysis Validation

As a first example, a simple tower model with a concentrated mass at the top and an elastic rotational support at the base is studied to validate the linear formulation. The finite element model is composed of 10 beam-column Euler-Bernoulli elements. The geometrical and material properties are shown in Figure1. The tower is under the action of a horizontal harmonic base displacement.Figure 1
Tower model with concentrated mass at top and elastoplastic support.Figure2 shows the variation of the maximum displacement at the top of the tower as a function of the forcing frequency Ω (resonance curve) for a forcing magnitude Ax=0.8g, where g is the acceleration of gravity, considering two values of the elastic support stiffness kr. The results obtained with the present formulation and continuation scheme agree with the results obtained from a time domain analysis. This shows the coherence of the result obtained with the proposed method.Figure 2
Variation of the maximum displacement at the top of the tower as a function of the forcing frequencyΩ obtained with the harmonic balance method (HBM) and with time domain integration. Ax=0.8g.


### 3.2. Nonlinear Analysis Validation

As a second example, a simple mass-spring system with cubic nonlinearity under harmonic excitation is studied to validate the nonlinear formulation. The schematic representation of the system and relevant properties are shown in Figure3. In this example, damping is not considered.Figure 3
1-DOF model with cubic nonlinear stiffness under harmonic excitation.Figure4 shows the nonlinear resonance curves obtained through the present HBM-Galerkin methodology and time domain simulations for positive (hardening) and negative (softening) cubic nonlinearity (knl = ±10 N/m3). The results obtained with the present HBM-Galerkin scheme and with time domain simulations are very close, validating the nonlinear formulation. It can be also observed that the algorithm is able to bypass the limit point associated with the saddle-node bifurcation. As expected, the time domain analysis is not able to trace the unstable branches of resonance curves.Figure 4
Nonlinear resonance curves. Maximum amplitude as a function of the forcing frequency.


### 3.3. Frequency Domain Analysis of Slender Frames

#### 3.3.1. Pitched-Roof Frame

A pitched-roof frame fixed at the base and with a constant cross section is studied to assess the behavior of this type of commonly used structure under horizontal and vertical base excitation. The geometry and material properties of the frame are shown in Figure5, as well as the configuration of the first four natural vibration modes, with the first and third modes being antisymmetric and the second and fourth modes being symmetric. The influence of the roof height on the response is also investigated. Table 1 shows the first four natural frequencies for the three values of H considered in the numerical analysis. The vibration frequencies decrease with H, as the frame slenderness increases. The structure is modeled with twenty bean-column elements: 8 elements of the same size for the columns and 12 elements for the roof.Table 1
Four first natural frequencies of the pitched-roof frame.
H (m)
Natural vibration frequency (rad/s)
1st mode
2nd mode
3rd mode
4th mode
0.0
21.82
48.97
130.54
161
3.0
19.61
47.14
108.23
159.80
6.0
15.40
36.42
83.21
92.08Figure 5
Pitched-roof frame model and four first natural vibration modes.Figure6 shows the variation of the norm of the horizontal and vertical components of the displacement at top of the frame, Dx and Dy, respectively, as a function of the forcing frequency, Ω, considering both linear (blue) and nonlinear (red) FE formulations for three values of the height of the roof pitches (H = 0 m, 3 m, and 6 m) and considering a horizontal excitation magnitude: Ax=0.96g. The horizontal base motion excites the first and the third modes only, which correspond to predominantly lateral deformation modes. The nonlinear effect is present in the first resonance region, where a hardening behavior is observed. The first peak and, consequently, the nonlinearity increase with increasing roof height, with the nonlinear maximum vibration amplitude being slightly lower than the linear one. After the first peak, the nonlinear and linear responses are practically coincident. The vertical displacement Dy is smaller than the vertical component and is overestimated in a linear analysis. Figure 7 displays the response for a vertical excitation magnitude: Ay=0.96g. The vertical base motion excites the second and fourth modes; these modes correspond to symmetric vibration modes with predominantly vertical motion. For vertical base excitation, the resonance peak of the vertical motion decreases when the roof height increases and the effect of nonlinearity is negligible. In this case, due to the symmetry of the displacement field, the lateral motion of the top node is zero.Figure 6
Variation of the norm of the horizontal and vertical vibration amplitudes at the roof top as a function of the excitation frequency for a horizontal base excitation and three values of the roof heightH. Ax=0.96g.Figure 7
Variation of the norm of the vertical vibration amplitude at the roof top as a function of the excitation frequency for a vertical base excitation and three values of the roof heightH. Ay=0.96g.Figure8 shows the norm of the horizontal and vertical displacement as a function of the forcing frequency for the frame excited in both the horizontal and vertical directions. This is a typical excitation in seismic analysis where both components are usually present. Following the suggestion of some seismic code [43], the intensity of the vertical base motion is adopted as 0.67 times the intensity of the horizontal base motion (Ay=0.66Ax). In this case, all vibration modes are excited. The influence of geometric nonlinearity is evident in the two first resonance regions, where a hardening behavior, leading to possible dynamic jumps, is observed. As the roof height increases, the magnitude of Dx corresponding to the first resonant peak increases, while the magnitude of second peak decreases. On the other hand, the magnitude of vertical displacement Dy is reduced in the two first resonance regions with the increase of the roof height. These nonlinear resonance curves exhibit amplitude and frequency limit points. This indicates the existence of fold (saddle-node) bifurcations, which can only be mapped by the use of continuation techniques. The solutions between the two limit points are unstable.Figure 8
Variation of the norm of the vertical and horizontal vibration amplitudes at the roof top as a function of the excitation frequency for simultaneous vertical and horizontal excitation and three values of the roof heightH. Ax=0.8g and  Ay=0.667Ax.


#### 3.3.2. Arched-Roof Frame

Now, an arched-roof frame fixed at the base and with a constant cross section is studied to assess the behavior of this typical structural form under horizontal and vertical base excitation. The structure has the same span, column height, and cross section as the pitched roof, as well as the same material properties, as shown in Figure9, to enable comparisons of the results. An incomplete circular arch with height H is considered. Figure 9 also shows the configuration of the first four vibration modes. Table 2 shows the first four natural frequencies for three values of the arch height, H. As in the previous example, the first and third modes are antisymmetric, while the second and fourth modes are symmetric.Table 2
Natural vibration frequencies of the arched-roof frame.
H (m)
Natural vibration frequency (rad/s)
1st mode
2nd mode
3rd mode
4th mode
1.0
21.02
51.73
122.60
167.05
3.0
17.81
48.51
97.77
151.60
5.0
13.20
40.79
74.02
113.41Figure 9
Arched-roof frame model and four first natural vibration modes.Figures10 and 11 show the variation of the norm of the vertical and horizontal displacements at top of the arch as a function of the forcing frequency for horizontal and vertical harmonic excitation, respectively, and considering both linear (blue) and nonlinear (red) formulations. For the horizontal base excitation, the first and the third modes are excited and the nonlinear effect is only noticed in the first peak leading to a hardening behavior and decrease of the maximum lateral displacement. As in the previous example, the first peak increases with the increase of the roof height. For vertical base excitation, only the second and fourth modes are excited and the nonlinear effect is slightly noted in the first peak only, where a softening behavior is observed. In this case, the first peak amplitude decreases when the roof height increases. Comparing the results of the pitched roof and arched roof for H = 3 m, one can observe that the arched roof is more flexible, leading to higher vibration amplitudes.Figure 10
Influence of the roof height on the nonlinear response of the arched-roof frame under horizontal excitation.Ax=0.96g.Figure 11
Influence of the roof height on the nonlinear response of the arched-roof frame under vertical excitation.Ay=0.96g.Figure12 shows the variation of the maximum horizontal (Dx) and vertical displacements (Dy) at the top of arch, considering both horizontal and vertical base excitations with Ay=0.667Ax and Ax=0.8g. In this case, all vibration modes are excited. The influence of the geometric nonlinearity is evident in the two first resonance regions, where a hardening behavior leading to possible dynamic jumps is observed. As the roof height increases, the magnitude of Dx corresponding to the first resonant peak increases, while the magnitude of second peak decreases. The magnitude of vertical displacement Dy decreases in the two first resonance regions with the increase of the roof height.Figure 12
Variation of the norm of the vertical and horizontal vibration amplitudes at the roof top as a function of the excitation frequency for simultaneous vertical and horizontal excitation and three values of the roof heightH. Ax=0.8g and  Ay=0.667Ax.


#### 3.3.3. Pitched-Roof Frame with a Long Span

Now, the pitched roof is again analyzed considering a span lengthL=16.0 m, increasing thus the frame’s slenderness. The other properties of the system are the same as those shown in Figure 5. Table 3 shows the first four natural frequencies for the three values of the roof height. Compared with the results in Table 1, a strong decrease in the natural frequencies is observed.Table 3
Natural vibration frequency of pitched-roof frame.
H (m)
Natural frequency of vibration (rad/s)
1st mode
2nd mode
3rd mode
4th mode
0.0
16.60
21.09
60.68
117.35
3.0
15.52
25.58
58.03
91.95
6.0
12.93
24.82
50.24
59.47Figure13 shows the resonance curves of the horizontal (Ax) and vertical (Ay) displacement at the top of the frame considering vertical and horizontal harmonic base excitation acting simultaneously with amplitude of the horizontal acceleration of Ax=0.8g and vertical acceleration of Ay=0.66Ax. The nonlinear effects are visible in the two first resonance regions, similar to the frame with L=10 m. However, the two first natural frequencies are close, especially for H=0 m, leading to modal interaction between the two first modes and a large excitation region, where large-amplitude vibrations occur. Compared with the results in Figure 6, for L=10 m, a marked increase in the lateral and vertical displacements is observed. The vertical displacements are particularly large around the second vibration frequency.Figure 13
Influence of the roof height on the nonlinear response of the frame under horizontal and vertical excitation.Ax=0.8g and  Ay=0.667Ax.


#### 3.3.4. Strong Base Motion

Now, the effect of strong base motion on the nonlinear dynamics of the pitched-roof frame (Section3.3.1), arched-roof (Section 3.3.2) frame, and pitched-roof frame with long span (Section 3.3.3) is assessed. For this, a parametric study considering the variation of the ground acceleration magnitude is presented. A maximum acceleration peak of 3.6g for a single direction excitation and of 3.0g for horizontal and vertical excitations acting simultaneously is adopted.Figure14 shows the variation of the maximum horizontal displacement at top of the pitched-roof frame (Figure 5) as a function of the forcing frequency for increasing base excitation magnitude. A sharp increase in the vibration amplitude and the hardening nonlinearity of the response in the first resonance region is observed. For strong base motions, the effect of the nonlinearity in the vicinity of the third vibration frequency begins to appear, being also of the hardening type. The second mode is not excited.Figure 14
Influence of increasing magnitudeAx on the nonlinear response of the frame under horizontal excitation. Norm of the horizontal displacement at top of pitched-roof as a function of the forcing frequency. L=10 and H=6 m.Figure15 shows the frame response due to vertical harmonic acceleration, considering Ax=1.8g and Ax=3.6g. The effect of the geometric nonlinearity is observed in the resonance regions associated with the second and fourth vibration modes, exhibiting softening behavior. The softening behavior is due to quadratic nonlinearities and is connected with the interaction between axial and transversal beam displacements. In such cases, the linear analysis underestimates the maximum vibration response. The softening effect due to negative quadratic nonlinearity can be observed in other slender structures under compressive loads like shallow arches and arises from the interaction between in-plane compressive forces and bending [25].Figure 15
Influence of increasing base motionAx on the nonlinear response of the frame under vertical excitation. Norm of the vertical displacement at top of pitched roof as a function of the forcing frequency. L=10 m and H=6 m.Figure16 shows the horizontal and vertical components of vibration as function of forcing frequency at top of pitched-roof frame (L=10 m and H=6 m), considering harmonic horizontal and vertical base excitations acting simultaneously and considering Ay=0.667Ax and two values of the acceleration amplitude: Ax=1.5g and Ax=3.0g. The nonlinear effect is present in all resonance regions for the two displacement components. A strong difference between the linear and nonlinear responses is particularly observed in the first resonance region.Figure 16
Horizontal and vertical displacements at top of the frame as a function of the forcing frequency due to simultaneous horizontal and vertical forcing.L=10 m and H=6 m.Finally, Figure17 shows the horizontal and vertical displacement components at top of pitched-roof frame as a function of forcing frequency for the frame with L=16 m and H=6 m, considering harmonic horizontal and vertical base excitations acting simultaneously with Ay=0.667Ax, for two values of Ax: 1.5g and 3.0g. Again, the strong base motion leads to a highly nonlinear response.Figure 17
Horizontal and vertical displacements at the top of the frame as a function of the forcing frequency due to simultaneous horizontal and vertical forcing.L=16 m and H=6 m.In all cases analyzed here, the resonance curves for the linear and nonlinear cases were easily obtained by the proposed numerical methodology, which can be easily applied to any structural system discretized by the FEM.


## 3.1. Linear Analysis Validation

As a first example, a simple tower model with a concentrated mass at the top and an elastic rotational support at the base is studied to validate the linear formulation. The finite element model is composed of 10 beam-column Euler-Bernoulli elements. The geometrical and material properties are shown in Figure1. The tower is under the action of a horizontal harmonic base displacement.Figure 1
Tower model with concentrated mass at top and elastoplastic support.Figure2 shows the variation of the maximum displacement at the top of the tower as a function of the forcing frequency Ω (resonance curve) for a forcing magnitude Ax=0.8g, where g is the acceleration of gravity, considering two values of the elastic support stiffness kr. The results obtained with the present formulation and continuation scheme agree with the results obtained from a time domain analysis. This shows the coherence of the result obtained with the proposed method.Figure 2
Variation of the maximum displacement at the top of the tower as a function of the forcing frequencyΩ obtained with the harmonic balance method (HBM) and with time domain integration. Ax=0.8g.


## 3.2. Nonlinear Analysis Validation

As a second example, a simple mass-spring system with cubic nonlinearity under harmonic excitation is studied to validate the nonlinear formulation. The schematic representation of the system and relevant properties are shown in Figure3. In this example, damping is not considered.Figure 3
1-DOF model with cubic nonlinear stiffness under harmonic excitation.Figure4 shows the nonlinear resonance curves obtained through the present HBM-Galerkin methodology and time domain simulations for positive (hardening) and negative (softening) cubic nonlinearity (knl = ±10 N/m3). The results obtained with the present HBM-Galerkin scheme and with time domain simulations are very close, validating the nonlinear formulation. It can be also observed that the algorithm is able to bypass the limit point associated with the saddle-node bifurcation. As expected, the time domain analysis is not able to trace the unstable branches of resonance curves.Figure 4
Nonlinear resonance curves. Maximum amplitude as a function of the forcing frequency.


## 3.3. Frequency Domain Analysis of Slender Frames

### 3.3.1. Pitched-Roof Frame

A pitched-roof frame fixed at the base and with a constant cross section is studied to assess the behavior of this type of commonly used structure under horizontal and vertical base excitation. The geometry and material properties of the frame are shown in Figure5, as well as the configuration of the first four natural vibration modes, with the first and third modes being antisymmetric and the second and fourth modes being symmetric. The influence of the roof height on the response is also investigated. Table 1 shows the first four natural frequencies for the three values of H considered in the numerical analysis. The vibration frequencies decrease with H, as the frame slenderness increases. The structure is modeled with twenty bean-column elements: 8 elements of the same size for the columns and 12 elements for the roof.Table 1
Four first natural frequencies of the pitched-roof frame.
H (m)
Natural vibration frequency (rad/s)
1st mode
2nd mode
3rd mode
4th mode
0.0
21.82
48.97
130.54
161
3.0
19.61
47.14
108.23
159.80
6.0
15.40
36.42
83.21
92.08Figure 5
Pitched-roof frame model and four first natural vibration modes.Figure6 shows the variation of the norm of the horizontal and vertical components of the displacement at top of the frame, Dx and Dy, respectively, as a function of the forcing frequency, Ω, considering both linear (blue) and nonlinear (red) FE formulations for three values of the height of the roof pitches (H = 0 m, 3 m, and 6 m) and considering a horizontal excitation magnitude: Ax=0.96g. The horizontal base motion excites the first and the third modes only, which correspond to predominantly lateral deformation modes. The nonlinear effect is present in the first resonance region, where a hardening behavior is observed. The first peak and, consequently, the nonlinearity increase with increasing roof height, with the nonlinear maximum vibration amplitude being slightly lower than the linear one. After the first peak, the nonlinear and linear responses are practically coincident. The vertical displacement Dy is smaller than the vertical component and is overestimated in a linear analysis. Figure 7 displays the response for a vertical excitation magnitude: Ay=0.96g. The vertical base motion excites the second and fourth modes; these modes correspond to symmetric vibration modes with predominantly vertical motion. For vertical base excitation, the resonance peak of the vertical motion decreases when the roof height increases and the effect of nonlinearity is negligible. In this case, due to the symmetry of the displacement field, the lateral motion of the top node is zero.Figure 6
Variation of the norm of the horizontal and vertical vibration amplitudes at the roof top as a function of the excitation frequency for a horizontal base excitation and three values of the roof heightH. Ax=0.96g.Figure 7
Variation of the norm of the vertical vibration amplitude at the roof top as a function of the excitation frequency for a vertical base excitation and three values of the roof heightH. Ay=0.96g.Figure8 shows the norm of the horizontal and vertical displacement as a function of the forcing frequency for the frame excited in both the horizontal and vertical directions. This is a typical excitation in seismic analysis where both components are usually present. Following the suggestion of some seismic code [43], the intensity of the vertical base motion is adopted as 0.67 times the intensity of the horizontal base motion (Ay=0.66Ax). In this case, all vibration modes are excited. The influence of geometric nonlinearity is evident in the two first resonance regions, where a hardening behavior, leading to possible dynamic jumps, is observed. As the roof height increases, the magnitude of Dx corresponding to the first resonant peak increases, while the magnitude of second peak decreases. On the other hand, the magnitude of vertical displacement Dy is reduced in the two first resonance regions with the increase of the roof height. These nonlinear resonance curves exhibit amplitude and frequency limit points. This indicates the existence of fold (saddle-node) bifurcations, which can only be mapped by the use of continuation techniques. The solutions between the two limit points are unstable.Figure 8
Variation of the norm of the vertical and horizontal vibration amplitudes at the roof top as a function of the excitation frequency for simultaneous vertical and horizontal excitation and three values of the roof heightH. Ax=0.8g and  Ay=0.667Ax.


### 3.3.2. Arched-Roof Frame

Now, an arched-roof frame fixed at the base and with a constant cross section is studied to assess the behavior of this typical structural form under horizontal and vertical base excitation. The structure has the same span, column height, and cross section as the pitched roof, as well as the same material properties, as shown in Figure9, to enable comparisons of the results. An incomplete circular arch with height H is considered. Figure 9 also shows the configuration of the first four vibration modes. Table 2 shows the first four natural frequencies for three values of the arch height, H. As in the previous example, the first and third modes are antisymmetric, while the second and fourth modes are symmetric.Table 2
Natural vibration frequencies of the arched-roof frame.
H (m)
Natural vibration frequency (rad/s)
1st mode
2nd mode
3rd mode
4th mode
1.0
21.02
51.73
122.60
167.05
3.0
17.81
48.51
97.77
151.60
5.0
13.20
40.79
74.02
113.41Figure 9
Arched-roof frame model and four first natural vibration modes.Figures10 and 11 show the variation of the norm of the vertical and horizontal displacements at top of the arch as a function of the forcing frequency for horizontal and vertical harmonic excitation, respectively, and considering both linear (blue) and nonlinear (red) formulations. For the horizontal base excitation, the first and the third modes are excited and the nonlinear effect is only noticed in the first peak leading to a hardening behavior and decrease of the maximum lateral displacement. As in the previous example, the first peak increases with the increase of the roof height. For vertical base excitation, only the second and fourth modes are excited and the nonlinear effect is slightly noted in the first peak only, where a softening behavior is observed. In this case, the first peak amplitude decreases when the roof height increases. Comparing the results of the pitched roof and arched roof for H = 3 m, one can observe that the arched roof is more flexible, leading to higher vibration amplitudes.Figure 10
Influence of the roof height on the nonlinear response of the arched-roof frame under horizontal excitation.Ax=0.96g.Figure 11
Influence of the roof height on the nonlinear response of the arched-roof frame under vertical excitation.Ay=0.96g.Figure12 shows the variation of the maximum horizontal (Dx) and vertical displacements (Dy) at the top of arch, considering both horizontal and vertical base excitations with Ay=0.667Ax and Ax=0.8g. In this case, all vibration modes are excited. The influence of the geometric nonlinearity is evident in the two first resonance regions, where a hardening behavior leading to possible dynamic jumps is observed. As the roof height increases, the magnitude of Dx corresponding to the first resonant peak increases, while the magnitude of second peak decreases. The magnitude of vertical displacement Dy decreases in the two first resonance regions with the increase of the roof height.Figure 12
Variation of the norm of the vertical and horizontal vibration amplitudes at the roof top as a function of the excitation frequency for simultaneous vertical and horizontal excitation and three values of the roof heightH. Ax=0.8g and  Ay=0.667Ax.


### 3.3.3. Pitched-Roof Frame with a Long Span

Now, the pitched roof is again analyzed considering a span lengthL=16.0 m, increasing thus the frame’s slenderness. The other properties of the system are the same as those shown in Figure 5. Table 3 shows the first four natural frequencies for the three values of the roof height. Compared with the results in Table 1, a strong decrease in the natural frequencies is observed.Table 3
Natural vibration frequency of pitched-roof frame.
H (m)
Natural frequency of vibration (rad/s)
1st mode
2nd mode
3rd mode
4th mode
0.0
16.60
21.09
60.68
117.35
3.0
15.52
25.58
58.03
91.95
6.0
12.93
24.82
50.24
59.47Figure13 shows the resonance curves of the horizontal (Ax) and vertical (Ay) displacement at the top of the frame considering vertical and horizontal harmonic base excitation acting simultaneously with amplitude of the horizontal acceleration of Ax=0.8g and vertical acceleration of Ay=0.66Ax. The nonlinear effects are visible in the two first resonance regions, similar to the frame with L=10 m. However, the two first natural frequencies are close, especially for H=0 m, leading to modal interaction between the two first modes and a large excitation region, where large-amplitude vibrations occur. Compared with the results in Figure 6, for L=10 m, a marked increase in the lateral and vertical displacements is observed. The vertical displacements are particularly large around the second vibration frequency.Figure 13
Influence of the roof height on the nonlinear response of the frame under horizontal and vertical excitation.Ax=0.8g and  Ay=0.667Ax.


### 3.3.4. Strong Base Motion

Now, the effect of strong base motion on the nonlinear dynamics of the pitched-roof frame (Section3.3.1), arched-roof (Section 3.3.2) frame, and pitched-roof frame with long span (Section 3.3.3) is assessed. For this, a parametric study considering the variation of the ground acceleration magnitude is presented. A maximum acceleration peak of 3.6g for a single direction excitation and of 3.0g for horizontal and vertical excitations acting simultaneously is adopted.Figure14 shows the variation of the maximum horizontal displacement at top of the pitched-roof frame (Figure 5) as a function of the forcing frequency for increasing base excitation magnitude. A sharp increase in the vibration amplitude and the hardening nonlinearity of the response in the first resonance region is observed. For strong base motions, the effect of the nonlinearity in the vicinity of the third vibration frequency begins to appear, being also of the hardening type. The second mode is not excited.Figure 14
Influence of increasing magnitudeAx on the nonlinear response of the frame under horizontal excitation. Norm of the horizontal displacement at top of pitched-roof as a function of the forcing frequency. L=10 and H=6 m.Figure15 shows the frame response due to vertical harmonic acceleration, considering Ax=1.8g and Ax=3.6g. The effect of the geometric nonlinearity is observed in the resonance regions associated with the second and fourth vibration modes, exhibiting softening behavior. The softening behavior is due to quadratic nonlinearities and is connected with the interaction between axial and transversal beam displacements. In such cases, the linear analysis underestimates the maximum vibration response. The softening effect due to negative quadratic nonlinearity can be observed in other slender structures under compressive loads like shallow arches and arises from the interaction between in-plane compressive forces and bending [25].Figure 15
Influence of increasing base motionAx on the nonlinear response of the frame under vertical excitation. Norm of the vertical displacement at top of pitched roof as a function of the forcing frequency. L=10 m and H=6 m.Figure16 shows the horizontal and vertical components of vibration as function of forcing frequency at top of pitched-roof frame (L=10 m and H=6 m), considering harmonic horizontal and vertical base excitations acting simultaneously and considering Ay=0.667Ax and two values of the acceleration amplitude: Ax=1.5g and Ax=3.0g. The nonlinear effect is present in all resonance regions for the two displacement components. A strong difference between the linear and nonlinear responses is particularly observed in the first resonance region.Figure 16
Horizontal and vertical displacements at top of the frame as a function of the forcing frequency due to simultaneous horizontal and vertical forcing.L=10 m and H=6 m.Finally, Figure17 shows the horizontal and vertical displacement components at top of pitched-roof frame as a function of forcing frequency for the frame with L=16 m and H=6 m, considering harmonic horizontal and vertical base excitations acting simultaneously with Ay=0.667Ax, for two values of Ax: 1.5g and 3.0g. Again, the strong base motion leads to a highly nonlinear response.Figure 17
Horizontal and vertical displacements at the top of the frame as a function of the forcing frequency due to simultaneous horizontal and vertical forcing.L=16 m and H=6 m.In all cases analyzed here, the resonance curves for the linear and nonlinear cases were easily obtained by the proposed numerical methodology, which can be easily applied to any structural system discretized by the FEM.


## 3.3.1. Pitched-Roof Frame

A pitched-roof frame fixed at the base and with a constant cross section is studied to assess the behavior of this type of commonly used structure under horizontal and vertical base excitation. The geometry and material properties of the frame are shown in Figure5, as well as the configuration of the first four natural vibration modes, with the first and third modes being antisymmetric and the second and fourth modes being symmetric. The influence of the roof height on the response is also investigated. Table 1 shows the first four natural frequencies for the three values of H considered in the numerical analysis. The vibration frequencies decrease with H, as the frame slenderness increases. The structure is modeled with twenty bean-column elements: 8 elements of the same size for the columns and 12 elements for the roof.Table 1
Four first natural frequencies of the pitched-roof frame.
H (m)
Natural vibration frequency (rad/s)
1st mode
2nd mode
3rd mode
4th mode
0.0
21.82
48.97
130.54
161
3.0
19.61
47.14
108.23
159.80
6.0
15.40
36.42
83.21
92.08Figure 5
Pitched-roof frame model and four first natural vibration modes.Figure6 shows the variation of the norm of the horizontal and vertical components of the displacement at top of the frame, Dx and Dy, respectively, as a function of the forcing frequency, Ω, considering both linear (blue) and nonlinear (red) FE formulations for three values of the height of the roof pitches (H = 0 m, 3 m, and 6 m) and considering a horizontal excitation magnitude: Ax=0.96g. The horizontal base motion excites the first and the third modes only, which correspond to predominantly lateral deformation modes. The nonlinear effect is present in the first resonance region, where a hardening behavior is observed. The first peak and, consequently, the nonlinearity increase with increasing roof height, with the nonlinear maximum vibration amplitude being slightly lower than the linear one. After the first peak, the nonlinear and linear responses are practically coincident. The vertical displacement Dy is smaller than the vertical component and is overestimated in a linear analysis. Figure 7 displays the response for a vertical excitation magnitude: Ay=0.96g. The vertical base motion excites the second and fourth modes; these modes correspond to symmetric vibration modes with predominantly vertical motion. For vertical base excitation, the resonance peak of the vertical motion decreases when the roof height increases and the effect of nonlinearity is negligible. In this case, due to the symmetry of the displacement field, the lateral motion of the top node is zero.Figure 6
Variation of the norm of the horizontal and vertical vibration amplitudes at the roof top as a function of the excitation frequency for a horizontal base excitation and three values of the roof heightH. Ax=0.96g.Figure 7
Variation of the norm of the vertical vibration amplitude at the roof top as a function of the excitation frequency for a vertical base excitation and three values of the roof heightH. Ay=0.96g.Figure8 shows the norm of the horizontal and vertical displacement as a function of the forcing frequency for the frame excited in both the horizontal and vertical directions. This is a typical excitation in seismic analysis where both components are usually present. Following the suggestion of some seismic code [43], the intensity of the vertical base motion is adopted as 0.67 times the intensity of the horizontal base motion (Ay=0.66Ax). In this case, all vibration modes are excited. The influence of geometric nonlinearity is evident in the two first resonance regions, where a hardening behavior, leading to possible dynamic jumps, is observed. As the roof height increases, the magnitude of Dx corresponding to the first resonant peak increases, while the magnitude of second peak decreases. On the other hand, the magnitude of vertical displacement Dy is reduced in the two first resonance regions with the increase of the roof height. These nonlinear resonance curves exhibit amplitude and frequency limit points. This indicates the existence of fold (saddle-node) bifurcations, which can only be mapped by the use of continuation techniques. The solutions between the two limit points are unstable.Figure 8
Variation of the norm of the vertical and horizontal vibration amplitudes at the roof top as a function of the excitation frequency for simultaneous vertical and horizontal excitation and three values of the roof heightH. Ax=0.8g and  Ay=0.667Ax.


## 3.3.2. Arched-Roof Frame

Now, an arched-roof frame fixed at the base and with a constant cross section is studied to assess the behavior of this typical structural form under horizontal and vertical base excitation. The structure has the same span, column height, and cross section as the pitched roof, as well as the same material properties, as shown in Figure9, to enable comparisons of the results. An incomplete circular arch with height H is considered. Figure 9 also shows the configuration of the first four vibration modes. Table 2 shows the first four natural frequencies for three values of the arch height, H. As in the previous example, the first and third modes are antisymmetric, while the second and fourth modes are symmetric.Table 2
Natural vibration frequencies of the arched-roof frame.
H (m)
Natural vibration frequency (rad/s)
1st mode
2nd mode
3rd mode
4th mode
1.0
21.02
51.73
122.60
167.05
3.0
17.81
48.51
97.77
151.60
5.0
13.20
40.79
74.02
113.41Figure 9
Arched-roof frame model and four first natural vibration modes.Figures10 and 11 show the variation of the norm of the vertical and horizontal displacements at top of the arch as a function of the forcing frequency for horizontal and vertical harmonic excitation, respectively, and considering both linear (blue) and nonlinear (red) formulations. For the horizontal base excitation, the first and the third modes are excited and the nonlinear effect is only noticed in the first peak leading to a hardening behavior and decrease of the maximum lateral displacement. As in the previous example, the first peak increases with the increase of the roof height. For vertical base excitation, only the second and fourth modes are excited and the nonlinear effect is slightly noted in the first peak only, where a softening behavior is observed. In this case, the first peak amplitude decreases when the roof height increases. Comparing the results of the pitched roof and arched roof for H = 3 m, one can observe that the arched roof is more flexible, leading to higher vibration amplitudes.Figure 10
Influence of the roof height on the nonlinear response of the arched-roof frame under horizontal excitation.Ax=0.96g.Figure 11
Influence of the roof height on the nonlinear response of the arched-roof frame under vertical excitation.Ay=0.96g.Figure12 shows the variation of the maximum horizontal (Dx) and vertical displacements (Dy) at the top of arch, considering both horizontal and vertical base excitations with Ay=0.667Ax and Ax=0.8g. In this case, all vibration modes are excited. The influence of the geometric nonlinearity is evident in the two first resonance regions, where a hardening behavior leading to possible dynamic jumps is observed. As the roof height increases, the magnitude of Dx corresponding to the first resonant peak increases, while the magnitude of second peak decreases. The magnitude of vertical displacement Dy decreases in the two first resonance regions with the increase of the roof height.Figure 12
Variation of the norm of the vertical and horizontal vibration amplitudes at the roof top as a function of the excitation frequency for simultaneous vertical and horizontal excitation and three values of the roof heightH. Ax=0.8g and  Ay=0.667Ax.


## 3.3.3. Pitched-Roof Frame with a Long Span

Now, the pitched roof is again analyzed considering a span lengthL=16.0 m, increasing thus the frame’s slenderness. The other properties of the system are the same as those shown in Figure 5. Table 3 shows the first four natural frequencies for the three values of the roof height. Compared with the results in Table 1, a strong decrease in the natural frequencies is observed.Table 3
Natural vibration frequency of pitched-roof frame.
H (m)
Natural frequency of vibration (rad/s)
1st mode
2nd mode
3rd mode
4th mode
0.0
16.60
21.09
60.68
117.35
3.0
15.52
25.58
58.03
91.95
6.0
12.93
24.82
50.24
59.47Figure13 shows the resonance curves of the horizontal (Ax) and vertical (Ay) displacement at the top of the frame considering vertical and horizontal harmonic base excitation acting simultaneously with amplitude of the horizontal acceleration of Ax=0.8g and vertical acceleration of Ay=0.66Ax. The nonlinear effects are visible in the two first resonance regions, similar to the frame with L=10 m. However, the two first natural frequencies are close, especially for H=0 m, leading to modal interaction between the two first modes and a large excitation region, where large-amplitude vibrations occur. Compared with the results in Figure 6, for L=10 m, a marked increase in the lateral and vertical displacements is observed. The vertical displacements are particularly large around the second vibration frequency.Figure 13
Influence of the roof height on the nonlinear response of the frame under horizontal and vertical excitation.Ax=0.8g and  Ay=0.667Ax.


## 3.3.4. Strong Base Motion

Now, the effect of strong base motion on the nonlinear dynamics of the pitched-roof frame (Section3.3.1), arched-roof (Section 3.3.2) frame, and pitched-roof frame with long span (Section 3.3.3) is assessed. For this, a parametric study considering the variation of the ground acceleration magnitude is presented. A maximum acceleration peak of 3.6g for a single direction excitation and of 3.0g for horizontal and vertical excitations acting simultaneously is adopted.Figure14 shows the variation of the maximum horizontal displacement at top of the pitched-roof frame (Figure 5) as a function of the forcing frequency for increasing base excitation magnitude. A sharp increase in the vibration amplitude and the hardening nonlinearity of the response in the first resonance region is observed. For strong base motions, the effect of the nonlinearity in the vicinity of the third vibration frequency begins to appear, being also of the hardening type. The second mode is not excited.Figure 14
Influence of increasing magnitudeAx on the nonlinear response of the frame under horizontal excitation. Norm of the horizontal displacement at top of pitched-roof as a function of the forcing frequency. L=10 and H=6 m.Figure15 shows the frame response due to vertical harmonic acceleration, considering Ax=1.8g and Ax=3.6g. The effect of the geometric nonlinearity is observed in the resonance regions associated with the second and fourth vibration modes, exhibiting softening behavior. The softening behavior is due to quadratic nonlinearities and is connected with the interaction between axial and transversal beam displacements. In such cases, the linear analysis underestimates the maximum vibration response. The softening effect due to negative quadratic nonlinearity can be observed in other slender structures under compressive loads like shallow arches and arises from the interaction between in-plane compressive forces and bending [25].Figure 15
Influence of increasing base motionAx on the nonlinear response of the frame under vertical excitation. Norm of the vertical displacement at top of pitched roof as a function of the forcing frequency. L=10 m and H=6 m.Figure16 shows the horizontal and vertical components of vibration as function of forcing frequency at top of pitched-roof frame (L=10 m and H=6 m), considering harmonic horizontal and vertical base excitations acting simultaneously and considering Ay=0.667Ax and two values of the acceleration amplitude: Ax=1.5g and Ax=3.0g. The nonlinear effect is present in all resonance regions for the two displacement components. A strong difference between the linear and nonlinear responses is particularly observed in the first resonance region.Figure 16
Horizontal and vertical displacements at top of the frame as a function of the forcing frequency due to simultaneous horizontal and vertical forcing.L=10 m and H=6 m.Finally, Figure17 shows the horizontal and vertical displacement components at top of pitched-roof frame as a function of forcing frequency for the frame with L=16 m and H=6 m, considering harmonic horizontal and vertical base excitations acting simultaneously with Ay=0.667Ax, for two values of Ax: 1.5g and 3.0g. Again, the strong base motion leads to a highly nonlinear response.Figure 17
Horizontal and vertical displacements at the top of the frame as a function of the forcing frequency due to simultaneous horizontal and vertical forcing.L=16 m and H=6 m.In all cases analyzed here, the resonance curves for the linear and nonlinear cases were easily obtained by the proposed numerical methodology, which can be easily applied to any structural system discretized by the FEM.


## 4. Conclusions

In this work, an incremental technique for the direct calculation of the nonlinear resonance curves of plane frames discretized by the finite element method and subjected to a base excitation is proposed. The transformation of discretized equations of motion, in the finite element method context, to the frequency domain is accomplished here through the classical harmonic balance method together with the Galerkin method. The resulting system of nonlinear equations in terms of the modal amplitudes and forcing frequency is solved by the Newton-Raphson method together with an arc-length procedure to obtain the nonlinear resonance curves. The formulation of the proposed method is validated comparing the present results with time domain simulations showing coherence and precision. The algorithm is able to map regions with coexisting stable and unstable solutions. A pitched-roof frame and an arched-roof frame under horizontal and vertical base excitation with varying roof height and spans are studied to illustrate the influence of the nonlinearity on the four first resonance regions. Horizontal base motion excites only the odd modes, while vertical base motion excites only the even modes. Under simultaneous vertical and horizontal excitation, all modes are excited. The roof height has a strong influence on the resonant peaks and degree of nonlinearity. The results also show the influence of increasing base excitation on the nonlinear behavior. Depending on the excitation magnitude, hardening or softening behavior can be observed. The proposed method is shown to be a useful method for the analysis of slender frame structures in frequency domain, since it is able to obtain the resonance curves with precision and small computational effort. This can be used as a preliminary design tool for frames under seismic and other types of base excitation, leading to a safer design.



---

*Source: 5281237-2017-03-30.xml*