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 κi (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:(4)κi= kcSigxixs,ρ,where kc is the contact stiffness, and the Sigmoid function Sig is defined as follows:(5)Sigxixs,ρ=11+exp−ρxi/xs−1.

This function (5) together with the linear ramp xi−xs 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 shown in Figure 1(b). In this paper, we use ρ=104 which is large enough to approximate a nonsmooth system without incurring the significant computational cost of an event-driven, piecewise, formulation.

Let us introduce some nondimensional spatial and temporal variables:(6)xi=uixs,xg=ugxs,τ=ωt,and parameters(7)ω2=km,β=kck,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 kc 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:(8)u″+2γCu′+Ku+βS⋅u−1=−1ug″,where the derivative with respect to scale time τ is expressed using •′=d•/dτ ,•′′=d2•/dτ2. A viscous classical Caughey [38, 39] orthogonal damping matrix 2γωC is introduced, where all linear vibration modes have the same small damping ratio of γ=0.05. Thus, the diagonal sigmoid contact matrix S in equation (8) is defined as follows:(9)S=Sigu1,ρ0L000Sigu2,ρ…00M…0…M00…Sigun−1,ρ000L0Sigun,ρ.

In this paper, we consider the case of monochromatic base excitation:(10)ug″=A sinΩτ,Ω=ωfω,where A and Ω are the nondimensional forcing amplitude and forcing frequency ratio, respectively.

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:(11)u=Φq=∑i=1nϕiqi,where model coordinates are defined as qT=q1,…, qn and the Euclidean normalised matrix of eigenvectors is defined as Φ=ϕ1,…,ϕn. Hence, the uncoupled modal equations are defined as follows:(12)qi″+2γϖiqi′+ϖi2qi=−ϕiT1ug″,where the Euclidean norm is employed for the eigenvectors and ϖi2 are the eigenvalues of the dynamic matrix K. These modal coordinates allow a standard explicit solution as follows:(13)qiτ=hiτ+giτ,where hiτ and giτ are the homogeneous and inhomogeneous solutions, respectively, and are defined as follows:(14)hiτ=a1ih1iτ+a2ih2iτ,h1iτ=exp−γiϖiτsinϖi1−γi2τ,h2iτ=exp−γϖiτcosϖi1−γi2τ,giτ=ciA sinΩτ+diA cosΩτ,ci=ϖi2−Ω2ϖi2−Ω22+2γiϖiΩ2,di=−2γiϖiΩϖi2−Ω22+2γiϖiΩ2.

By premultiplying (11) by the jth modal eigenvector and employing the orthogonality of eigenvectors, we can obtain ϕjTu=ϕjTϕjzj; thus, by differentiation and algebra, we can solve of modal integration constants:(15)a1ja2j=h1jh2jh1j′h1j′−1ϕjTu0−gjϕjTu0′−gj′.

Therefore, the integration constants a1j and a2j of the jth modal homogeneous solution can be derived from a given set of boundary conditions u0=uτ0and u0′=u′τ0 at a given time τ0.

The first-order ODE (state-space) form of equation (8) can be expressed by introducing z2=u′ and z1=u:(19)z1′z2′=−0−IK2γCz1z2−0βS⋅z1−1+1ug″.

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

The first example cases are shown in Figure 2, which represent a phase-space projection onto ℝ2 of individual coordinate pairs ui,ui′. These 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₁) 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 u3 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 uiτ≥100T, where forcing period T is 2π/Ω. These steady-state solutions are stroboscopically sampled at τ=jT at j∈ℕ 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 u3 impacts the stop while DoFs u1 and u2 do not. Again, we have a P₁ solution. We introduce a notation C_ijk that implies contact of DoFs i, j, k with the stop, e.g., C₁₃ 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₈ solution with a C₂₃ (contact of DoFs u2 and u3 during steady-state responses). It is interesting to note that while DoF u3 impacts the stop at each and every cycle, the DoF u2 does not. This suggests that, at each change in the number of impacts of DoF u2, there is a grazing-type bifurcation. Thus, it appears the bifurcation structure of this system is extremely complex.

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 Ω=0.445 and contact-to-storey 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₃ contact, (iii) is C₂₃ contact, (iv) indicate C₁₃ and C₂₃ coexisting solutions, (v) indicate C₁₂₃ and C₂₃ coexisting solutions, (vi) indicate two coexisting C₁₂₃ solutions, (vii) indicate C₁₂₃ and C₂₃ coexisting solutions, and (viii) indicates C₁₂₃ contact. The contact states can be obtained by considering solutions whose DoFs are greater than unity, ui≥1.

The first grazing bifurcation between region (i) and (ii) displays the typical chaos before settling down to a P₁, C₃ solution as the amplitude increases.

The second significant grazing bifurcation occurs when floor DoF u2 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 P₉, P₈, P₇, P₆, and P₄ plotted in Figure 4, the C₂₃ impact contains DoF u3 impacting each and every cycle, while the DoF u2 impacts only twice per period of the solution. The P₅ periodicity was not captured on the sweep up, but only on the sweep down. It coexists with a chaotic attractor.

Region (iv) is the first one to exhibit detected C₁₃ and C₂₃ coexisting solutions. The presence of a P₂, C₁₃ solution shows impacting of floor DoFs u1 and u3, but floor DoF u2 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 P₁, C₂₃ solution with its coexisting P₁, C₁₂₃ solution. Finally, at very high amplitudes region (viii), we have not observed coexisting solutions, and in this region, impact is for all floors DoFs.

The 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. 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.

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

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 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. This is a sister plot to that shown in Figure 3. The initial grazing bifurcation chaos, region (ii), is adverse, but after this, the P₁ 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 P₂, C₁₃ solution produces a larger δ than that of the P₄, C₂₃ solution. This solution P₂, C₁₃ 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 P₁ impacting solutions. In this case, most of these P₁ solutions are shown to be beneficial in controlling the growth of interstorey drift.