Abstract
This paper explores the nonlinear dynamics of a multidegree of freedom (MDoF) structure impacting a rigid stop. The contact mechanics is simplified by continuous sigmoid function idealisation of a lossless spring. By introducing a smooth nonlinear formulation, we avoid the computational expense of event-driven, piecewise, nonsmooth dynamics. A large parametric study using high-performance computing is undertaken. The nondimensional equations of motion suggest one primary structural parameter, contact-to-storey stiffness ratio, and two excitation parameters, nondimensional ground amplitude and frequency. Bifurcation plots indicate an extremely rich and complex behaviour, particularly in the cases where at least two-floor degrees of freedom (DoFs) impact the stop and when the contact-to-storey stiffness ratio is large. When considering interstorey drift as a performance measure, period-1 impacting solutions are generally favourable when compared to an analogous nonimpacting case. This paper also discusses whether chaotic impacting can be favourable. Finally, we consider the question of whether higher modes are significantly excited, at a linear resonance, for impacting solutions to this system.
1. Introduction
Seismic pounding between closely spaced adjacent buildings results in an extremely complex nonlinear nonsmooth system dynamics. This problem is further exasperated by the uncertain time history of any future earthquake ground excitation, complicated 3D building geometries, and nonlinear material degrading behaviour. As an initial scoping mathematical study, we consider a reduced-order model that simplifies the system into one that is more amenable to parametric analysis.
Research into the pounding of buildings can be split into four main disciplines. The first is observational and began with the growth of close-packed and high-rise cities in earthquake-prone zones, most notably California. As computing capacity increased, the ability to solve pounding models with numerical methods became feasible. The third is research into mitigation methods, necessitated by inadequate or infeasible separation distances listed in building codes. The final is the extension of these models to other complex situations, such as bridge deck/abutment pounding. Observation of pounding as a phenomenon first received scientific attention in the 1970s and early 1980s [1–3]. Following a significant expansion of research in the area in the 1990s, Anagnostopoulos questioned the need for investigation into pounding at all [4]. He argued that the minute proportion of affected buildings damaged from pounding did not sufficiently justify the heavily conservative separation distances listed in design codes. In contrast, Maison performed an observational investigation into evidence of pounding damage from the 1989 Loma Prieta earthquake [5]. He stressed the need to address pounding hazards. Maison also found separation codes impractical, due to the need for efficient land use in cities, and failed to address pounding for existing buildings (those at most risk). Observation-based research is still required. Recent years have seen continued observational studies of pounding damage [6]. Jankowski [7] performed finite element analyses on pounding damaged and collapsed buildings. He found lighter/weaker buildings consistently at the most risk from pounding, especially when adjacent to a heavier or stronger one.
As observation-based research could explore individual pounding cases in little depth, theoretical models were required to better understand the problem. The first numerical models were developed in the late 1980s, encouraged by increased computing capabilities and a series of damaging earthquakes on the San Andreas Fault. Anagnostopoulos’s first publication on the subject [8] used a series of idealised SDOF systems, with pounding impact modelled by linear viscoelastic spring-dashpots (Kelvin–Voigt model). It was realised that the two main contributors to pounding damage were separation distances and differential structural properties (building heights and stiffnesses) between the structures considered. Many standard analytical methods were found to be inapplicable due to the nonlinearity of the problem. Maison and Kasai [9, 10] expanded this model to investigate the effect of building separation on factors such as displacement, shear, overturning moment, and interstorey drift, using ground motion based on past earthquakes. He found that although peak displacements were similar, all other stresses and moments were increased by pounding. It was suggested that due to the complexity of the problem, buildings should be evaluated on a case-by-case basis, rather than using separation distance codes. Anagnostopoulos and Spiliopoulos [11] and Papadrakakis et al. [12] confirmed the results of Maison’s model and stressed the importance of finding an alternative to separation codes.
Lin and Weng investigated [13] the probability of pounding based on separation distance using a lumped-mass cantilever beam model. These findings were applied to the dense urban area of the Taipei metropolitan district. He found that pounding probability was greatest when the difference in natural frequencies between adjacent buildings was highest. As natural frequency ratios were completely unaccounted for in the Taiwan building code (based on Californian codes), the suggested separation distance was likely to be ill-founded. This problem was explored in greater depth, numerically, by Pantelides and Ma [14, 15]. The Californian codes were found to be adequate but overly conservative. They suggested a relaxation of the codes if buildings increased active and passive building damping capacity. The first technical mitigation solution was proposed by Westermo, as a series of interstructural connections [16], to prevent pounding. Although eliminating pounding, it was found to increase stress/strain damage in the more common, nonpounding case. Anagnostopoulos proposed using shear walls extending across the building separation distance [17]. Although some degree of local damage in the collision area is unavoidable, the shear walls were found to prevent slabs pounding onto column midspans, the most damaging form of pounding. As a valid solution for buildings under construction and many existing buildings, shear collision walls appear to be the most feasible solution when compared to the “do nothing” case.
The most important extension of this work was the case of bridge decks pounding against abutments. Kawashima and Ruangrassamee used displacement response spectra to investigate the problem further [18] and found similar conclusions to building pounding: displacements were amplified, gap size between elements reduced effects, and differential masses caused significant damage in the smaller element. Vega et al. [19] confirmed these results and investigated them with respect to the Pacific Engineering Earthquake Response (PEER) formula framework and found it provided conservative forces for bridges. Komodromos investigated the effect that pounding has on the performance of seismically isolated buildings [20]. It was found that the increased accelerations and deflections may lead to higher modes of deformation rather than the rigid-body behaviour expected under seismic isolation. Papadrakakis and Mouzakis built on a lumped-mass model to three dimensions by treating each floor as a rigid diaphragm [21] to investigate the effects of friction in pounding. Utilising finite element analysis to model building pounding had not been considered until Jankowski [7] provided a modelling method more applicable to real-life design cases where building geometries are known.
A parallel research field to that of seismic pounding has been that of vibro-impact oscillators. In more forensic studies, chaos has been identified in these sinusoidally forced impacting systems. Foale and Bishop [22] investigated the bifurcation behaviour at the grazing point of vibro-impact oscillators with Hertzian impact. Ervin and Wickert [23] noted a cascade of period-doubling bifurcations in vibro-impact oscillator models, similar to the Poincaré sweeps discussed later. Modelling pounding as vibro-impact oscillators, Davis [24] and Chau and Wei [25] both noted chaotic effects at high differential stiffnesses when using a Hertz pounding model, and in 2004, Wagg and Bishop [26] mapped the regions of chaotic behaviour around the first modal frequency. All these studies have in common the use of Hertzian impact and sinusoidal forcing and are SDOF systems (thus chaos occurring at grazing). Andreaus has extended this to consider multisided impact with nonlinear impact [27, 28] with both experimental and simulational approaches [29, 30].
Work on MDOF systems has also been carried out by researchers such as John et al. [31] who found chaotic behaviour with a nonlinear impulse model for MDOF cantilevers. John et al. [31] modelled an MDOF system, with continuously distributed mass, that impacts only at a single point, namely, the top of the cantilever. In our paper, we explore an MDOF cantilever system, with discrete lumped storey masses, which has the potential for multiple points of impact. This is therefore a system with a richer range of coexisting nonsmooth bifurcation events.
In this paper, we consider a reduced-order model with a lossless spring contact approximated by a smooth sigmoid function. The general idea is to produce an analogous smooth simplified nonlinear system that shows almost all the esoteric features (chaos, bifurcations, and so on) of the nonsmooth piecewise system. By employing this reduced-order idealisation, we aim to parametrically explore a huge range of cases as the numerical solutions are more efficiently obtained. The aims of this paper are as follows:(i)What is the effect of stiffness differential, between structures, on the periodicity of pounding solutions for a range of forcing frequencies and amplitudes?(ii)Is pounding, at any periodicities including chaos, ever a desirable system feature?(iii)Does pounding amplify/attenuate higher-order structural modes?
2. Theoretical Modelling
Consider the case of a structure, shown in Figure 1(a), with lumped masses m at the n DoFs x1 to xn. These floor masses are supported laterally by “storey” stiffnesses k. The structure is subjected to horizontal ground motion . This structure is located at a separation distance from a stop. In this paper, the stop represents a massive “effectively rigid” secondary building whose response to the pounding of the small building is negligible. This would approximate the case of a small building neighbouring a very tall/large building without introducing the further nonlinear complexities of the stop dynamics as well. Thus, the system’s Lagrangian is as follows:where the first term represents the kinetic system energy, the second term represents the strain energy of the building, and the final term describes the work done by the contact forces , shown in Figure 1(b). This simplified formulation expresses a lossless Hertzian impact against a flexible medium.
(a)
(b)
The Euler–Lagrange equations of motions are as follows:where is a column vector of ones; the elastic stiffness matrix , nonlinear terms in the contact force vector , and DoF vector are defined as follows:
2.1. Sigmoid Contact Springs
Spring-dashpot impact models have been widely used since [32]. The linear spring and the stereomechanical [33, 34] are two common models. A linear spring increases linearly in force after the point of impact. A stereomechanical model is based not only on force but also on velocity. These have been used as the basis for the two major models of impact used in earthquake pounding: Kelvin–Voigt and Hertz. The Kelvin–Voigt model has been used since the late 1980s and is a combination of both these models resulting in a viscoelastic effect. Although first proposed by Davis [24], it was not until Chau and Wei [25] that the nonlinear Hertz model became widely used, possibly due to limitations on computational capacity. These models have been combined in numerous variants by Muthukumar and Jankowski [35]. As previously stated, the aim of introducing a nonlinear spring for the contact mechanics is to avoid the computational cost/complexity of nonsmooth [36], piecewise, dynamical simulations. Thus, we consider an approximation that results in a nonlinear smooth system. The stop stiffness (at the ith floor contact) is modelled using a sigmoid (smooth and continuous) function [37], equation (4), rather than a piecewise, nonsmooth, event-driven formulation:where is the contact stiffness, and the Sigmoid function is defined as follows:
This function (5) together with the linear ramp can produce a range of different types of behaviour. At , we have a simple linear, pretensioned spring, but as tends to infinity, equation (4) tends to a piecewise linear spring, as shown in Figure 1(b). In this paper, we use which is large enough to approximate a nonsmooth system without incurring the significant computational cost of an event-driven, piecewise, formulation.
2.2. Nondimensionalisation
Let us introduce some nondimensional spatial and temporal variables:and parameterswhere is a structural frequency parameter, the first modal frequency of the linear system (in the noncontact case) is simply , where is the smallest eigenvalue of dynamic matrix . Parameter represents the ratio of contact stiffness to structural storey stiffness .
In much of the literature on pounding, the seismic gap is highlighted as an important system parameter. In the lossless contact formulation presented in this paper, the seismic gap is subsumed into a nondimensional ground excitation amplitude, see equation (6). Hence, the nondimensional ground motion amplitude can be viewed as proportional to the ratio of the dimensional ground motion amplitude to the seismic gap.
Hence, the equation of motion (2) is re-expressed as follows:where the derivative with respect to scale time is expressed using . A viscous classical Caughey [38, 39] orthogonal damping matrix is introduced, where all linear vibration modes have the same small damping ratio of . Thus, the diagonal sigmoid contact matrix in equation (8) is defined as follows:
2.3. Defining Ground Excitation
In this paper, we consider the case of monochromatic base excitation:where and are the nondimensional forcing amplitude and forcing frequency ratio, respectively.
2.4. Modal Contributions in MDoF Responses
Modal contributions are heavily employed in earthquake engineering [40, 41]; therefore, it is useful to express the MDoF system in terms of a set of SDoF systems during the noncontact phases. The purpose here is to see whether there is any evidence suggestive of an increase/decrease in modal responses due to pounding. The solution of the system of equations (8) (in the case of no-contact) is expressed in terms of modal coordinates by the changes of coordinates:where model coordinates are defined as and the Euclidean normalised matrix of eigenvectors is defined as . Hence, the uncoupled modal equations are defined as follows:where the Euclidean norm is employed for the eigenvectors and are the eigenvalues of the dynamic matrix . These modal coordinates allow a standard explicit solution as follows:where and are the homogeneous and inhomogeneous solutions, respectively, and are defined as follows:
By premultiplying (11) by the jth modal eigenvector and employing the orthogonality of eigenvectors, we can obtain ; thus, by differentiation and algebra, we can solve of modal integration constants:
Therefore, the integration constants and of the jth modal homogeneous solution can be derived from a given set of boundary conditions at a given time .
2.5. Defining Response Measures
2.5.1. ith Modal Contribution to Overall Displacements
Computing the norm squared of the total displacement vector allows a contraction to a scalar function; thus,
The orthogonality of eigenvectors for is employed to simplify. The contribution of each mode can be assessed separately; thus,where is the ratio of the norm of the mode j displacement vector to the norm of the total structural displacement response vector (during the noncontact phase of the oscillations). These relative modal norms are time-varying functions; therefore, we will contact further by employing the mean of for the steady-state solutions.
2.5.2. Normalised Peak InterStorey Drift Ratio
The interstorey drift is defined as . It is a measure that is strongly correlated with the local maximum internal storey bending moments. This drift normally increases with increasing ground motion amplitude. Therefore, we propose, , the normalised peak interstorey drift to ground displacement amplitude ratio, defined as follows:where would be a constant for steady-state linear dynamical system responses, but for nonlinear systems, it can vary with forcing amplitude and frequency . Hence, the conservative/unconservative effects of nonlinearity can be described qualitatively.
2.6. State-Space Formulation and the Complete Set of System Parameter s
The first-order ODE (state-space) form of equation (8) can be expressed by introducing and :
The numerical solution of equation (19) is obtained using Matlab solver ode15s [42]. Therefore, in the final form of the equations of motion (19), we have four main system parameters:(i)Number of storeys, n(ii)Contact stiffness ratio, (iii)Nondimensional forcing amplitude, A(iv)Forcing frequency ratio,
3. Parametric Exploration for a Heuristic 3-Storey Structure
In this paper, we explore a heuristic case of a 3-storey structure (). Therefore, equation (8) represents a 3-DoF nonlinear system with a phase space of . This reduced-order model case approximates a small structure pounding a much larger more massive structural system.
3.1. Some Initial Examples of Phase-Space Responses
The first example cases are shown in Figure 2, which represent a phase-space projection onto of individual coordinate pairs . These figures are depicted for the case , which represents the linear first modal nondimensional, undamped, frequency of this 3-storey structure (in the case of no-contact), where is the square root of the eigenvalue, i.e., of the dynamic matrix .
(a)
(b)
(c)
Figure 2(a) depicts a period-1 (P1) steady-state solution for the case at very nearly linear resonance. This is at an amplitude A just lower than the first grazing bifurcation [43, 44] where the top floor DoF almost touches the stop. These steady-state solutions are obtained by discarding the transient parts of the solution, which we assume is approximated as the first 100 forcing cycles. Therefore, our steady-state solutions are , where forcing period is . These steady-state solutions are stroboscopically sampled at at and so Poincaré points are displayed as coloured dots on the solutions to highlight the periodicity of the solution. Figure 2(b) depicts the case where the top floor DoF impacts the stop while DoFs and do not. Again, we have a P1 solution. We introduce a notation Cijk that implies contact of DoFs i, j, k with the stop, e.g., C13 means floor DoF 1 and 3 which impact the stop during steady-state responses. Figure 2(c) depicts a much more complex case of a P8 solution with a C23 (contact of DoFs and during steady-state responses). It is interesting to note that while DoF impacts the stop at each and every cycle, the DoF does not. This suggests that, at each change in the number of impacts of DoF , there is a grazing-type bifurcation. Thus, it appears the bifurcation structure of this system is extremely complex.
3.2. Amplitude Sweep Bifurcation Diagrams at Linear Resonance
Given the potential complexity hinted at in the previous section, we perform a conventional amplitude sweep. Starting at a low amplitude A = 0.035 (a nonimpact case), we solve the system and determine the Poincaré points for the steady-state solution; then, we increment A. By using the previous solution’s Poincaré points as initial conditions for this new A, we path-follow (parametric continuation) in some primitive way the stable solution branch we are on. This is a pseudopath-following that would be attempted in an experimental setting. Because of the likely existence of fold bifurcations, we first sweep up (increasing) in amplitude and then sweep down (decreasing) in amplitude. This can identify some coexisting solutions but does not guarantee to identify all coexisting solutions.
Figure 3 displays a bifurcation plot (an amplitude sweep) at a constant forcing frequency ratio and contact-to-storey stiffness ratio . At the top of this figure, the Roman numerals (i) to (viii) indicate 8 distinct contact regions, where different combinations of DoFs impact the stop. (i) is a noncontact region, (ii) is C3 contact, (iii) is C23 contact, (iv) indicate C13 and C23 coexisting solutions, (v) indicate C123 and C23 coexisting solutions, (vi) indicate two coexisting C123 solutions, (vii) indicate C123 and C23 coexisting solutions, and (viii) indicates C123 contact. The contact states can be obtained by considering solutions whose DoFs are greater than unity, .
The first grazing bifurcation between region (i) and (ii) displays the typical chaos before settling down to a P1, C3 solution as the amplitude increases.
The second significant grazing bifurcation occurs when floor DoF starts to impact the stop at the beginning of region (iii). This spawns a very complex series of chaotic zones [45] interspersed with reducing periodicity “windows.” The phase space plots of these are described in Figure 4. For amplitude values between those depicted in Figure 4, there exist chaotic regions like that shown for A = 0.09877. It is interesting to note that, for attractors P9, P8, P7, P6, and P4 plotted in Figure 4, the C23 impact contains DoF impacting each and every cycle, while the DoF impacts only twice per period of the solution. The P5 periodicity was not captured on the sweep up, but only on the sweep down. It coexists with a chaotic attractor.
(a)
(b)
(c)
(d)
(e)
(f)
Region (iv) is the first one to exhibit detected C13 and C23 coexisting solutions. The presence of a P2, C13 solution shows impacting of floor DoFs and , but floor DoF does not impact the stop. This suggests that, during the noncontact phase, of each cycle, mode 3 is significant. This highlights the fact that impact significantly changes the ratio of the norms of system modal responses. At the end of region (vii), we have a fold bifurcation of the P1, C23 solution with its coexisting P1, C123 solution. Finally, at very high amplitudes region (viii), we have not observed coexisting solutions, and in this region, impact is for all floors DoFs.
3.3. Amplitude Sweep Bifurcation Diagrams away from Linear Resonance
The initial amplitude sweep at the undamped linear resonance shown in Figure 3 demonstrates the system complexity. In this section, we explore amplitude sweeps away from this value. Figure 5 shows frames from the video file that shows an animation of amplitude bifurcation plots as the frequency ratio is varied frame by frame. This video is one of the only ways to visually communicate the extremely rich bifurcation structure of this reduced-order system. The brute-force computational cost of so many multiple scans is very significant, of the order of 5 million individual time-history solutions to equation (19). This required the University of Bristol’s HPC BlueCrystal.
(a)
(b)
(c)
As an alternative, for each parameter in these sweep-up amplitude bifurcation plots, we identify the periodicity of the solution obtained (note that other coexisting solutions may also exist) and then plot a coloured 2D bifurcation plot displayed in Figure 6.
(a)
(b)
Figure 6(a) represents 2D-parameter space bifurcation plot which underlines the scope of the problem of predicting system performance, for a physically analogous system, in the case where two or more impacting points start to come into play. This plot was for a large contact-to-storey stiffness ratio, = 1000. Figure 6(b) repeats these analyses but for a smaller contact-to-storey stiffness ratio, = 100. This smaller value could be viewed as a reduction in contact stiffness or an increase in structure storey stiffness or some combination of the two. What is clear is that, at = 100, the complexity of the bifurcation parameter space is greatly reduced. Hence, if reduced complexity (and hence more predictability) is a design target, then lower values of should be one important criterion.
3.4. Beneficial/Adverse Effects of Impacting
A question considered here is whether we can conclude that impacting is generally deleterious to interstorey drift. As a general principle, it is important to repeat, as previously stated, that interstorey drift appears to always increase with increase in amplitude regardless of the impact regime. The normalised interstorey drift measure, equation (18), represents a quantity that should remain constant with increase in amplitude for the noncontact regimes. Therefore, a reduction in this normalised drift would indicate that contact is favourable, at particular parameter values, compared with the case where the stop is not present. Figure 7(a) displays the normalised drift during the amplitude sweep up and down for the case at the undamped linear resonance where . This is a sister plot to that shown in Figure 3. The initial grazing bifurcation chaos, region (ii), is adverse, but after this, the P1 impacting solution tends to reduce . The high-period “windows” and chaos of region (iii) sometimes produce a reduction and sometimes a small increase in . In region (iv), the P2, C13 solution produces a larger than that of the P4, C23 solution. This solution P2, C13 counterintuitively tends to excite the 2nd mode more during the noncontact phase of cycles. This change in modal norms during the noncontact part of cycles is highlighted in Figure 7(b). Exciting the system at the undamped resonance should produce almost entirely mode 1 response in region (i). As the amplitude increases and we move into contact regions (ii) to (viii), we observe an increase in mode 2 and 3 norms. This suggests that the typical seismic analysis using the participation factor concept for estimating modal contributions [38, 40, 41] is not suitable for an impacting system.
Figure 8 extends the analysis from Figure 7(b) to explore beneficial/adverse regions of parameter space for the sweep-up amplitude bifurcation plots (like those discussed previously in Figures 5 and 6). Figure 8(a) highlights the regions where chaos is marginally favourable in green and adverse in red. Both adverse and beneficial regions exist. Figure 8(b) repeats this analysis for P1 impacting solutions. In this case, most of these P1 solutions are shown to be beneficial in controlling the growth of interstorey drift.
(a)
(b)
4. Conclusions
The analyses in this paper indicate the extremely complex nature of the nonsmooth dynamics of this MDoF impacting oscillator. The requirement for extensive computational runs necessitated both the model simplification of the contact using a smooth sigmoid function and high-performance computing. This work has shown for the first time that smooth, continuous impact functions on relatively simple cantilever structures can exhibit the same effects as those seen in more computationally intensive impulse models.
The effect of the contact-to-structure stiffness ratio (the stiffness differential) is significant. As this contact-to-structure stiffness ratio reduces, the complexity of the bifurcation parameter space also reduces, i.e., the likelihood of encountering chaotic, high periodicity, and coexisting solutions is reduced. While at high contact-to-structure stiffness ratios, very complex and rich bifurcation structures are encountered. Hence, if predictability of response is a design objective, then lowering contact-to-structure stiffness ratio should be an aim.
A period-1 periodic pounding contact reduces the interstorey drift relative to the case where there is no stop. Although this is not completely universally true, there is a high probability that this is the case. At the initial grazing, bifurcation chaos appears to be adverse on interstorey drift. Apart from this case, there was no consistent evidence that chaotic impacting was either beneficial or adverse on interstorey drift. As the locations in the parameter space of any beneficial chaos are difficult to predict, it appears that designing for chaotic impacting is not a reasonable design objective.
One of the most important results is the effect of pounding/impact on higher modal amplitudes in the response. For the case of a linear system without any impact, modal amplitudes, at steady-state, are dependent on forcing frequency alone. At a linear resonance, there should be only ever a single mode in the response, at steady-state. For the impacting system, impacts generate forcing that excites the higher modes (during the noncontact phases on each forcing cycle) as demonstrated in this paper. This significant amplification of higher modes has important implications for the case of designing seismically exciting systems as traditional modal combination rules are no longer valid.
Data Availability
Data related to the research within this research article are available from the authors on request.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Authors’ Contributions
Peter J. Christopher and Nicholas A. Alexander contributed equally.