Exploring the Dynamics of Base-Excited Structures Impacting a Rigid Stop

This paper explores the nonlinear dynamics of a multidegree of freedom (MDoF) structure impacting a rigid stop. The contact mechanics is simpliﬁed 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 stiﬀness 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-ﬂoor degrees of freedom (DoFs) impact the stop and when the contact-to-storey stiﬀness 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 signiﬁcantly excited, at a linear resonance, for impacting solutions to this system.


Introduction
Seismic pounding between closely spaced adjacent buildings results in an extremely complex nonlinear nonsmooth system dynamics. is 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. e 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. e third is research into mitigation methods, necessitated by inadequate or infeasible separation distances listed in building codes. e 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][2][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. e 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. ese 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. is problem was explored in greater depth, numerically, by Pantelides and Ma [14,15].
e Californian codes were found to be adequate but overly conservative. ey suggested a relaxation of the codes if buildings increased active and passive building damping capacity. e 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.
e 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 rigidbody 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 vibroimpact 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. is 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. e 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. e 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?

Theoretical Modelling
Consider the case of a structure, shown in Figure 1(a), with lumped masses m at the n DoFs x 1 to x n . ese floor masses are supported laterally by "storey" stiffnesses k. e structure is subjected to horizontal ground motion x g . is structure is located at a separation distance x s 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. is 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. us, 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 Figure 1(b). is simplified formulation expresses a lossless Hertzian impact against a flexible medium. e Euler-Lagrange equations of motions are as follows: where 1 is a column vector of ones; the elastic stiffness matrix K , nonlinear terms in the contact force vector F(x), and DoF vector x are defined as follows:

Sigmoid Contact Springs.
Spring-dashpot impact models have been widely used since [32]. e 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. ese have been used as the basis for the two major models of impact used in earthquake pounding: Kelvin-Voigt and Hertz. e 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. ese 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. us, we consider an approximation that results in a nonlinear smooth system. e stop stiffness κ i (at the ith floor contact) is modelled using a sigmoid (smooth and continuous) function [37], equation (4), rather than a piecewise, nonsmooth, eventdriven formulation: where k c is the contact stiffness, and the Sigmoid function Sig is defined as follows: Sig is function (5) together with the linear ramp (x i − x s ) can produce a range of different types of behaviour. At ρ � 0, we have a simple linear, pretensioned spring, but as ρ tends to infinity, equation (4) tends to a piecewise linear spring, as Mathematical Problems in Engineering shown in Figure 1(b). In this paper, we use ρ � 10 4 which is large enough to approximate a nonsmooth system without incurring the significant computational cost of an eventdriven, piecewise, formulation.

Nondimensionalisation.
Let us introduce some nondimensional spatial and temporal variables: and parameters where ω is a structural frequency parameter, the first modal frequency of the linear system (in the noncontact case) is simply ω 1 � ω �� λ 1 , where λ 1 is the smallest eigenvalue of dynamic matrix K. Parameter β represents the ratio of contact stiffness k c to structural storey stiffness k.
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 2cωC is introduced, where all linear vibration modes have the same small damping ratio of c � 0.05. us, the diagonal sigmoid contact matrix S in equation (8) is defined as follows: 2.3. Defining Ground Excitation. In this paper, we consider the case of monochromatic base excitation: where A and Ω are the nondimensional forcing amplitude and forcing frequency ratio, respectively.

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. e purpose here is to see whether there is any evidence suggestive of an increase/decrease in modal responses due to pounding. e 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 q T � [q 1 , . . . , q n ] and the Euclidean normalised matrix of eigenvectors is defined as Φ � [ϕ 1 , . . . , ϕ n ]. Hence, the uncoupled modal equations are defined as follows: where the Euclidean norm is employed for the eigenvectors and ϖ 2 i are the eigenvalues of the dynamic matrix K. ese modal coordinates allow a standard explicit solution as follows: where h i (τ) and g i (τ) 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 ϕ T j u � ϕ T j ϕ j z j ; thus, by differentiation and algebra, we can solve of modal integration constants: erefore, the integration constants a 1j and a 2j of the jth modal homogeneous solution can be derived from a given set of boundary conditions u 0 � u(τ 0 )and u 0 ′ � u ′ (τ 0 ) at a given time τ 0 .

i th Modal Contribution η j to Overall Displacements.
Computing the norm squared of the total displacement vector allows a contraction to a scalar function; thus, e orthogonality of eigenvectors ϕ T j ϕ i � 0 for ∀i ≠ j is employed to simplify. e contribution of each mode j can be assessed separately; thus, where η j 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). ese relative modal norms η j (τ) are time-varying functions; therefore, we will contact further by employing the mean of η j (τ) for the steady-state solutions.

Normalised Peak InterStorey Drift Ratio.
e interstorey drift is defined as u i − u i− 1 . It is a measure that is strongly correlated with the local maximum internal storey bending moments.
is drift normally increases with increasing ground motion amplitude. erefore, 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 A and frequency Ω. Hence, the conservative/unconservative effects of nonlinearity can be described qualitatively.

State-Space Formulation and the Complete Set of System
Parameter s. e first-order ODE (state-space) form of equation (8) can be expressed by introducing z 2 � u ′ and z 1 � u: e numerical solution of equation (19) is obtained using Matlab solver ode15s [42]. erefore, 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, Ω

Parametric Exploration for a Heuristic 3-Storey Structure
In this paper, we explore a heuristic case of a 3-storey structure (n � 3). erefore, equation (8) represents a 3-DoF nonlinear system with a phase space of R 2n+1 . is reducedorder model case approximates a small structure pounding a much larger more massive structural system.

Some Initial Examples of Phase-Space
Responses. e first example cases are shown in Figure 2, which represent a phase-space projection onto R 2 of individual coordinate Mathematical Problems in Engineering pairs (u i , u i ′ ). ese figures are depicted for the case Ω � 0.455, 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., �� λ 1 of the dynamic matrix K. Figure 2(a) depicts a period-1 (P 1 ) steady-state solution for the case at very nearly linear resonance.
is is at an amplitude A just lower than the first grazing bifurcation [43,44] where the top floor DoF u 3 almost touches the stop. ese steady-state solutions are obtained by discarding the transient parts of the solution, which we assume is approximated as the first 100 forcing cycles. erefore, our steady-state solutions are u i (τ ≥ 100T), where forcing period T is 2π/Ω. ese steady-state solutions are stroboscopically sampled at τ � jT at j ∈ N 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 u 3 impacts the stop while DoFs u 1 and u 2 do not. Again, we have a P 1 solution. We introduce a notation C ijk that implies contact of DoFs i, j, k with the stop, e.g., C 13 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 P 8 solution with a C 23 (contact of DoFs u 2 and u 3 during steady-state responses). It is interesting to note that while DoF u 3 impacts the stop at each and every cycle, the DoF u 2 does not. is suggests that, at each change in the number of impacts of DoF u 2 , there is a grazing-type bifurcation. us, it appears the bifurcation structure of this system is extremely complex.

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.
is 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. is 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 Ω � 0.445 and contact-tostorey stiffness ratio β � 1000. 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 C 3 contact, (iii) is C 23 contact, (iv) indicate C 13 and C 23 coexisting solutions, (v) indicate C 123 and C 23 coexisting solutions, (vi) indicate two coexisting C 123 solutions, (vii) indicate C 123 and C 23 coexisting solutions, and (viii) indicates C 123 contact. e contact states can be obtained by considering solutions whose DoFs are greater than unity, u i ≥ 1. e first grazing bifurcation between region (i) and (ii) displays the typical chaos before settling down to a P 1 , C 3 solution as the amplitude increases. e second significant grazing bifurcation occurs when floor DoF u 2 starts to impact the stop at the beginning of region (iii). is spawns a very complex series of chaotic zones [45] interspersed with reducing periodicity "windows." e 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 P 9 , P 8 , P 7 , P 6 , and P 4 plotted in Figure 4, the C 23 impact contains DoF u 3 impacting each and every cycle, while the DoF u 2 impacts only twice per period of the solution. e P 5 periodicity was not captured on the sweep up, but only on the sweep down. It coexists with a chaotic attractor.

Mathematical Problems in Engineering
Region (iv) is the first one to exhibit detected C 13 and C 23 coexisting solutions. e presence of a P 2 , C 13 solution shows impacting of floor DoFs u 1 and u 3 , but floor DoF u 2 does not impact the stop.
is suggests that, during the noncontact phase, of each cycle, mode 3 is significant. is 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 P 1 , C 23 solution with its coexisting P 1 , C 123 solution. Finally, at very high amplitudes region (viii), we have not observed coexisting solutions, and in this region, impact is for all floors DoFs.

Amplitude Sweep Bifurcation Diagrams away from Linear
Resonance.
e initial amplitude sweep at the undamped linear resonance Ω � 0.445 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. is video is one of the only ways to visually communicate the extremely rich bifurcation structure of this reduced-order system. e 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). is required the University of Bristol's HPC BlueCrystal.
As an alternative, for each parameter (A, Ω) 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. 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. is 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. is 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.

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. e normalised interstorey drift measure, equation (18), represents a quantity that should remain constant with increase in amplitude for the noncontact regimes. erefore, a reduction in this normalised drift would indicate that contact is favourable, at particular parameter (A, Ω) 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 Ω � 0.445.
is is a sister plot to that shown in Figure 3. e initial grazing bifurcation chaos, region (ii), is adverse, but after this, the P 1 impacting solution tends to reduce δ. e high-period "windows" and chaos of region (iii) sometimes produce a reduction and sometimes a small increase in δ. In region (iv), the P 2 , C 13 solution produces a larger δ than that of the P 4 , C 23 solution. is solution P 2 , C 13 counterintuitively tends to excite the 2nd mode more during the noncontact phase of cycles. is 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. is 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 sweepup amplitude bifurcation plots (like those discussed previously in Figures 5 and 6). where chaos is marginally favourable in green and adverse in red. Both adverse and beneficial regions exist. Figure 8(b) repeats this analysis for P 1 impacting solutions. In this case, most of these P 1 solutions are shown to be beneficial in controlling the growth of interstorey drift.   is 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. e effect of the contact-to-structure stiffness ratio β (the stiffness differential) is significant. As this contact-tostructure 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.

Conclusions
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. is 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.