This paper investigates the complex dynamical behavior of a rigid block structure under harmonic ground excitation, thereby mimicking, for instance, the oscillation of the system under seismic excitation or containers placed on a ship under periodic acting of sea waves. The equations of motion are derived, assuming a large frictional coefficient at the interface between the block and the ground, in such a way that sliding cannot occur. In addition, the mathematical model assumes a loss of kinetic energy when an impact with the ground takes place. The resulting mathematical model is then formulated and studied in the framework of impulsive dynamical systems. Its complex dynamical response is studied in detail using two different approaches, based on direct numerical integration and path-following techniques, where the latter is implemented via the continuation platform COCO (Dankowicz and Schilder). Our study reveals the presence of various dynamical phenomena, such as branching points, fold and period-doubling bifurcation of limit cycles, symmetric and asymmetric periodic responses, and chaotic motions. By using the basin stability method, we also investigate the properties of solutions and their ranges of existence in phase and parameter spaces. Moreover, the study considers ground excitation conditions leading to the overturning of the block structure and shows parameter regions wherein such behavior can be avoided.

The dynamics of the rocking block has been studied for more than 60 years, and still, we can find new inspiring, complex phenomena in such a simple system. Numerous structures are modeled as rocking blocks, i.e., tall buildings, free-standing power transformers, cell towers, silos, and containers. The rocking block is used to mimic the slender structure, which is not, or badly, connected to the ground. Hence, the transversal movements of the foundation can lead to the appearance of the so-called rocking motion. The rocking motion will cause an overturn or structures damage for some sets of the system’s excitation parameters. We distinguish two main types of excitation. The first one is an earthquake which affects the system placed on the ground. In contrast, the second one is a sea wave acting on container ships. As is easy to see, both of them have significantly different properties. The earthquake has a relatively high and broad frequency spectrum that does not last long. Hence, it affects the structure in a short time with a very harsh forcing. Usually, it results in damage or overturns of the system. On the other hand, the sea wave has a lower frequency, lasts for a long time, and varies slowly. Therefore, in addition to the damage or overturn, we can observe a stable oscillatory motion. This paper focuses on the second type of excitation, and we model it as a periodic function acting on a rigid block. In extreme events, the sea wave excitation has similar or worse properties as the earthquake, and one should study such cases by assuming the forcing as a complex signal. The rocking motion is modeled as rotations with respect to the structure’s left and right bottom corners. The impacts occur during the motion, so the model is not only nonlinear but also piecewise. Due to discontinuity, we observe various dynamical phenomena, including bifurcations and the coexistence of solutions. We describe the bifurcation scenario leading from steady-state to chaos via grazing, fold, symmetry breaking, and period-doubling bifurcations. We show the bifurcation borders in two-parameter space, i.e., the amplitude and the frequency of the excitation. Such projection allows presenting ranges where the given solution is stable. In the considered system, we observe the coexistence of at least two solutions (overturn and equilibrium/periodic/chaotic). To reveal the structure of the multistable phase space, we investigate the system by using the sample-based methods. Presented results describe the properties of solutions and their ranges of existence in phase and parameter spaces. It shows the probability of the existence of solutions as a function of initial conditions and parameter values. Hence, we can evaluate the influence of those quantities on the appearance of solutions and define the ranges where the given state is the most or the less probable.

## I. INTRODUCTION

The rocking response of slender rigid structures due to the excitation at its base, e.g., during an earthquake, has been the subject of great interest for many researchers (see Refs. 1–3). Such a rocking motion in large engineering structures, such as buildings, free standing power transformers, and cell towers can cause catastrophic failures as a consequence of overturning, both from an economical and a human perspective. One of the first systematic studies in this area was pioneered by Housner,^{4} who has established a rocking block model that has been widely used by many researchers in the past. His study suggests that the survival of a tall slender structure during earthquakes depends on a scale effect making the structure more stable against overturning. The block equations of motion were described by piecewise-smooth equations depending on the sign of its rotation angle. Despite its apparent simplicity, the problem poses serious difficulties from an analytical point of view. Thus, the majority of previous works have employed numerical approaches to analyze the underlying structure (e.g., Refs. 5–13), most of which are based on studies relying on direct numerical simulations. In the present work, we will undertake a comprehensive parametric study of the block structure under ground excitation by using specialized numerical techniques, based on a numerical continuation approach.

Past numerical works in this research mainly focused on understanding the fundamental properties of the rocking block model. For example, Yim *et al.*^{5} studied the rocking response of rigid blocks subjected to earthquake ground motion and developed a numerical procedure to solve the nonlinear equations of motion governing the rocking motion. They also studied the problem from a probabilistic point of view with the ground motion modeled as a random process. Their numerical results suggest that the response of the block is sensitive to small changes in its size and slenderness ratio, as well as the features of the ground motion. Ishiyama^{6} developed a computer program to simulate the motion of a rigid body subjected to horizontal and vertical ground motions. It was found that the criteria for overturning of the body depended on the horizontal acceleration and the velocity of the ground. Hogan^{7} analyzed a rigid block undergoing harmonic forcing without damping. An in-depth analytical and numerical study has been performed for the slender structures because in this case, it is possible to linearize equations of the motion and find their exact solutions. However, such a specific ratio between the rocking block dimensions is not typical for real-world systems and do not give a comprehensive overview of most existing system’s dynamics. Hogan presented the stability of resonant periodic orbits in the amplitude–frequency of the excitation plane for different restitution coefficient values. The values of these parameters significantly influence ranges of the resonant periodic solutions’ existence. Hence, its proper selection is crucial for the correctness of modeling of the system’s behavior. Then, the author presented the numerical study of the nonlinear equations of motion for slender and non-slender rocking blocks. He showed the asymmetric solutions and period-doubling sequence leading to chaos. In such a way, we can find the basins of attraction with fractal properties, confirming that multiple solutions coexist in this system. Finally, the analysis of the system response to a short time excitation, which simulates a real earthquake, was conducted. Hogan assumed that an earthquake has a finite time of duration, so one can predict the conditions leading to non-toppling or toppling states. The following paper by Hogan^{14} extended the above-described results and focused on the coexistence of resonant periodic solutions in the system. He proved that for adjacent initial conditions, one can observe transitions to different solutions; hence, the system is multistable. At the end of the paper, the comparison of theoretical results with the experiments by Tso and Wong^{15,16} was presented. The matching of results is satisfactory, and simulations help us to explain experimental observations. Hogan also investigated the effect of damping on rigid block motion under harmonic forcing in Ref. 9. The study is analytical and bases on linearized equations of the motion by assuming a slender rigid block. The author showed several two-parameter plots showing the stability of resonant periodic solutions for four parameters of the system, i.e., the amplitude and the frequency of the forcing, restitution, and damping coefficients. Then, he supported his study with numerical computation of the nonlinear equations of the motion and showed the multiple solutions, asymmetric responses, period-, and impact-doubling cascades. This study shows that the rocking block under harmonic excitation has complex dynamics, and a plethora of different behavior can occur. It is a motivation to study the system with new numerical tools to have a more general overview of dynamics, the structure of the phase space, and bifurcation scenarios.

Lenci and Rega^{17} developed a method to eliminate the heteroclinic intersections embedded in the block’s dynamics by controlling the periodic excitation. In Ref. 10, a comprehensive experimental investigation was carried out to study the rocking response of four blue granite stones with different geometrical characteristics under free, harmonic, and random vibrations. In particular, two simulation approaches, one through Lagrangian formalism and the second one employing discrete element model, were used to identify the models’ parameters and study the rigid block structure’s dynamic response. In Ref. 11, Dimitrakopoulos and DeJong introduced new closed-form solutions and original similarity laws, which can shed light on the fundamental aspects of the rocking block. The focus was on the transient dynamics of the rocking block under finite-duration excitations. In Ref. 12, closed-form expressions were obtained for the dynamic response of a rocking block, whereas rigorous overturning criteria were established for conditions with and without impact. Later on, Brzeski *et al.*^{13} studied the influence of tuned mass absorbers on the dynamics of a rigid block in order to prevent overturning. The presented method was used to optimize tuned mass absorber’s parameters to achieve the highest chance that the block remains standing. The final results confirmed that it was possible to design and tune the stabilizing device in such a way that it was highly effective regardless of the rigid block asymmetry. The influence of the excitation’s initial phase on the response of the structure has been studied for four selected values of the phase by Lenci and Rega.^{18} The authors show how the dynamics of the system varies and which phase shift is the best for stabilization of the rocking block, while in Ref. 19, the response of the rigid block on sine or cosine impulse is presented. Moreover, the authors show the analysis of values of initial conditions on the response of the structure. However, the selection of possible values of initial conditions was limited to values close to zero.

This paper presents a detailed numerical study of dynamical response of the rocking block model with a special focus on its symmetric and asymmetric solutions. We are particularly interested in the periodic solutions of the system and the evolution of such responses under the changing parameters of base excitation. When non-smooth nonlinearity, e.g., impact in this system, is involved, the dynamical study of systems of this type (see Refs. 20–22) brings well-known challenges from a numerical point of view. In order to study the behavior of the rocking block, we will adopt two different types of numerical approaches, namely, direct numerical integration and path-following methods for non-smooth dynamical systems.^{23} For this kind of system, the state space can be divided into disjoint subregions, so the rocking block model in each region can be described by a smooth vector field. Hence, special care needs to be considered in order to get reliable numerical approximations of the system response in an efficient way.^{24}

Path-following (continuation) methods are well-established techniques in applied mathematics,^{25} which allows a systematic analysis of a system response under the variation of control parameters. For piecewise-smooth systems, specialized computational tools for continuation, such as SlideCont,^{26} TC-HAT,^{27} and COCO,^{28} have been developed, and the latter will be adopted in the current work for the numerical study of the rocking block model. The MATLAB-based numerical platform COCO supports the continuation and bifurcation detection of multi-segment periodic orbits, a feature that has been extensively used in the past for the numerical bifurcation study of piecewise-smooth dynamical systems. For example, it was used to study the robustness of a controller for controlling the multistability in a vibro-impact capsule system.^{29} In Ref. 21, it was adopted to determine the optimal control parameters in an impacting system in terms of its energy expenditure due to the control signal and transient behavior of the control error. Liu *et al.*^{30} employed the continuation platform COCO to investigate the dynamical response of a higher-order drifting oscillator. By using the platform, the authors were able to eliminate the bistability observed in the system and improve system’s performance with respect to its rate of penetration. In the current work, we will formulate the rocking block model in the framework of piecewise-smooth dynamical systems and employ the continuation platform COCO to study its complex dynamical response in detail. Our main concern is to understand the bifurcation structure of the model, which will reveal the presence of branching points, fold and period-doubling bifurcation of limit cycles, symmetric and asymmetric periodic responses, and chaotic motion. In addition, one of the main concerns is to identify ground oscillation modes leading to overturning and how to avoid them by considering suitable parameter regions where such events can be avoided. This work can then be used as a starting point for developing effective control strategies to prevent overturning during the rocking motion of the rigid block.

The second method applied to study the dynamics of the system considered in this paper is a basin stability method.^{31,32} It enables to quantify the stability of solutions based on the probability of reaching a given attractor from random initial conditions. To calculate the basin stability measure, one needs to perform a significant number of Bernoulli trials and classify the final solution reached in each trial. The efficiency of the basin stability method in comparison to the method of basins of attraction is obvious when the dimension of phase space is larger than two. In such a case, the classical approach requires much larger computational effort and interpretation of results is very complex (necessity of projections of results on two-dimensional cross sections of multi-dimensional phase space).^{33–35} The basin stability method gives statistical information on existing solutions (the volume of basins of attraction in relation to the total volume of considered phase space) and their location in the phase space. Moreover, drawing of initial conditions allows us to monitor the probability to reach a given solution. We are also confident that all existing attractors have been detected and no rare solutions have been overlooked.^{36–39} In our previous papers,^{40–42} we extended the basin stability approach by assuming that values of system parameters were random. It allows fast scanning of the initial conditions and the parameter spaces to find all meaningful solutions that exist in the system and the ranges of their existence. Thus, we can determine where the system is mono-stable/multi-stable or which perturbation or parameter mismatch can lead to destabilization of the system.^{43} The important data can be obtained from a two-dimensional diagram, which shows the probability of a given solution as a function of two variables. Such data are important in practical applications because in many cases, parameters of the system are varying or it is impossible to identify their precise values and we must be sure that there is no danger of a sudden jump to another solution or appearance of bifurcation.^{44} In this paper, we use this method to investigate the structure of the phase space and identified the regions of multi-stability and the dependence of a given type of solution on the system parameters.

The rest of the paper is organized as follows. In Sec. II, the physical model and equations of motion of the rocking block are derived. Section III introduces the mathematical formulation of the rocking block model as a piecewise-smooth dynamical system. Then, a detailed numerical investigation using the continuation methods via COCO, including one- and two-parameter studies, is presented in Sec. IV. In Sec. V, we show the solutions coexisting in the system and their probability of occurrence using the basin stability method. Finally, some concluding remarks are drawn in Sec. VI.

## II. PHYSICAL MODEL

### A. Equations of motion

In the current work, the system presented in Fig. 1 is considered, where a rigid block with the mass $M$ stands on the base subjected to periodic excitation with infinite rigidity. Here, we posit that only the horizontal motion of the base, $x\xa8=A\Omega 2sin\u2061(\Omega t)$, is considered, so $x=\u2212Asin\u2061(\Omega t)$, and the coefficient of friction between the block and the base is sufficiently large to prevent sliding. Therefore, the block stands still or performs oscillations around one of its corners (i.e., point $O+$ or $O\u2212$) that is currently in contact with the ground. The angular position of the block is given by the generalized coordinate $\phi $, which describes its angular deflection from the standing still position. We assume that the distances between the center of mass and its corners are equal, given by the length $R$. The slenderness of the block is described by $\theta b$, which is the angle between the lower edge of the block and the line that connects the block’s center of mass with its corner. If there is no external excitation, the block will overturn for angles larger than $\pi 2\u2212\theta b$, and we define it as the critical angle for the block.

When $\phi >0$, the block is rocking around the point $O+$ the equation of motion is given by

while for $\phi <0$,

where $cb$ is the viscous damping coefficient in the pivot points of the block and $g$ is the acceleration due to gravity.

The block will start rotation around one of the pivot points ($O+$ or $O\u2212$) depending on its acceleration $x\xa8$. This motion continues up to the moment when the block overturns or until falls back to the standing position $\phi =0$ knocking the ground. Here, we assume that there is no bouncing from the block, so after the collision, the rotation of the block continues around one of the pivot points, and the angular momentum of the system is conserved. To model the contact between the block and the ground, the impact law^{5} is adopted as

where $\phi \u02d9\u2212$ and $\phi \u02d9+$ are the angular velocities of the block just before and after the impact, respectively. The corresponding restitution coefficient $e$ is defined as

which is determined from the conservation of angular momentum^{5} and depends on the slenderness of the block $\theta b$ only. It is worth noting that the energy loss due to the contact depends on the material properties of the block and the base. However, for the idealized rigid block and base, formula (4) can give a good approximation of the restitution coefficient, which has been widely adopted. The mathematical model described above will be investigated for the parameter values given in Table I, which correspond to a realistic physical scenario.

Parameters . | Values/range . | Units . |
---|---|---|

M | 20 × 10^{3} | kg |

c_{b} | 101 | N m s |

g | 9.81 | m/s^{2} |

R | 3.764 | m |

θ_{b} | 0.8795 | rad |

A | (0, 4] | m |

Ω | $[\pi 4,2\pi ]$ | rad/s |

Parameters . | Values/range . | Units . |
---|---|---|

M | 20 × 10^{3} | kg |

c_{b} | 101 | N m s |

g | 9.81 | m/s^{2} |

R | 3.764 | m |

θ_{b} | 0.8795 | rad |

A | (0, 4] | m |

Ω | $[\pi 4,2\pi ]$ | rad/s |

For our numerical investigation, it is convenient to write the equations of motion in non-dimensional form according to the formulas

By using these transformations, we obtain a set of dimensionless equations for the motion of the rocking block as follows. For $\phi >0$, we have

while for $\phi <0$, it gives

where the prime symbol denotes differentiation with respect to the non-dimensional time $t~$. In what follows, tildes will be dropped for the sake of simplicity.

## III. DESCRIPTION OF THE MODEL AS A NON-SMOOTH DYNAMICAL SYSTEM

The rocking block model [(3), (6), and (7)] introduced in Sec. II belongs to the class of *non-smooth dynamical systems*.^{23} Such a class of system naturally arises when the underlying governing laws of a certain problem include non-smooth phenomena, for instance, impacts (hard or soft), friction, switches, etc. The equations of motion in these cases can be written in terms of a piecewise-smooth vector field in combination with a set of jump laws, hence producing discontinuities in the solution and in its first or higher-order derivatives. In the considered rocking block model, a discontinuity in the solution is produced by the restitution law assumed in (3), which accounts for the instantaneous velocity change when the block hits the ground ($\phi =0$).

In general, a non-smooth dynamical system can be characterized by three main components: a set of (smooth) vector fields, *event functions*, and *jump functions*. In this way, a system trajectory is defined by a collection of operation modes whose evolution is described by a specific smooth vector field. The termination condition of a given operation mode is defined in terms of the zeroes of a particular *event function* ($\phi =0$ in our case). Then, the initial point for the next operation mode is determined by a given *jump function*, which maps the terminal point of the current mode to the initial point of the next one. In our case, the jump function defines the angular velocity of the block just after an impact with the ground occurs, according to the restitution law (3). This formulation enables the application of specialized numerical techniques for the analysis of non-smooth dynamical systems, for instance, via path-following (continuation) methods.^{27,28}

For the study of the rocking block model via path-following methods, it is convenient to write the governing equations in an autonomous form. To this end, we will employ a self-excited, nonlinear oscillator given by Ref. 25

which possesses the asymptotically stable solution $r(t)=sin\u2061(\omega t)$, $s(t)=cos\u2061(\omega t)$, for all $t\u22650$. In this way, we can write the periodically forced rocking block model in an autonomous form, which then allows us to study the model via numerical continuation methods.

Let us define $\lambda :=(\omega ,\alpha ,\xi ,\theta b)\u2208(R+)4$ and $z:=(\phi ,\psi ,r,s)T\u2208[\u2212\pi 2+\theta b,\pi 2\u2212\theta b]\xd7R3$ as the parameters and state variables of the system, respectively, where $\psi $ will capture the angular velocity [see Eqs. (9)–(11)]. In this setting, the rocking block model [(3), (6), and (7)] introduced in Sec. II can be written as a non-smooth dynamical system as follows:

where the vector fields are defined as

In addition, for our numerical investigation, we will consider the following solution measure:

computed for a given $T0$-periodic solution. In this way, we will be able to detect critical points defined by the condition $\phi =\u2212\pi 2+\theta b$ or $\phi =\pi 2\u2212\theta b$. They define static boundaries where block overturning can occur. For the numerical implementation, the quantity defined in (12) is computed from the discretized solution, using the $max$ function.

## IV. NUMERICAL STUDY OF THE PERIODICALLY EXCITED ROCKING BLOCK

In this section, our main goal will be the numerical investigation of the rocking block motion subject to periodic base excitation, introduced in Sec. II. As mentioned earlier, the considered physical model may represent a number of practical scenarios, such as buildings, silos, containers, scaffoldings, free standing power transformers, etc., and one of the main questions in this matter is to determine the conditions under which overturning of the object occurs. For the considered physical configuration, we will pay close attention to solutions crossing the boundaries $\phi =\theta b\u2212\pi 2$ and $\phi =\pi 2\u2212\theta b$ (see Fig. 1), which can lead to overturning if, for instance, the base excitation suddenly stops. In our numerical investigation, the parameter values given in Table I will be employed, which correspond to a realistic physical scenario.

In Fig. 2, an initial reference solution of the rocking block model (9) is plotted and computed for the parameter values given in Table I. This corresponds to a periodic trajectory showing (odd) symmetry, as can be seen from the time plots given in panel (b). As expected, the resulting solution presents points of discontinuity, owing to the restitution law assumed in (3), thereby producing recurrent discontinuities in the angular velocity. At such points (occurring when the rocking block hits the ground $\phi =0$), the angular velocity is instantaneously reduced by a factor $e$, corresponding to the assumed restitution coefficient. This piecewise-smooth periodic solution will be used in our investigation as a starting point.

### A. One-parameter analysis of the system response

As mentioned earlier, the periodic solution depicted in Fig. 2 will be used as a starting point for our numerical study. It corresponds to a rocking block motion with odd symmetry, and we will investigate first how this solution is affected when the frequency of ground excitation $\omega $ changes. For this purpose, we will employ numerical continuation methods for non-smooth dynamical systems, implemented via the continuation platform COCO.^{28} The result of this process is presented in Fig. 3, showing the behavior of the solution measure Amp [see (12)] with respect to $\omega $. In panel (a), the resulting bifurcation diagram is depicted, where the black branch represents the continuation of symmetric periodic solutions, like the one shown in panel (b), computed at the test point P1 ($\omega =0.95$). The solid lines represent branches of asymptotically stable solutions, while the dashed lines stand for families of unstable orbits.

As can be seen in the bifurcation diagram [Fig. 3(a)], the stable branch containing the test point P1 is limited by the critical points BP1 and SN1, corresponding to a branching point and fold bifurcation of limit cycles, respectively. At the branching point, the original stable symmetric solution loses stability, while a solution branch of stable asymmetric periodic orbits is born (green curve), as the one shown in Fig. 3(c). This green branch persists for a wide window of $\omega $-values and terminates at another branching point BP2. Here, the asymmetric solutions disappear and the family of symmetric solutions (black branch) recovers stability. In addition, two critical points SB1 and SB2 are detected during the numerical continuation, corresponding to $\omega $-values for which the periodic solution makes tangential contact with the critical overturning boundary $\phi =\xb1(\theta b\u2212\pi 2)$. In this way, we can identify a branch in the bifurcation diagram (between SB1 and SB2) where overturning can occur, as discussed before. An example of such a critical solution is given in Fig. 3(c), computed at SB2 ($\omega \u22482.6635$).

Next, we will analyze the behavior of the rocking block when the amplitude of the base excitation $\alpha $ is varied while keeping the excitation frequency constant. As before, we will carry out the investigation via numerical continuation methods for non-smooth dynamical systems, using the solution shown in Fig. 2 as starting point. The result of this process is displayed in Fig. 4. In this case, we identify three different branches, plotted in black, green, and red. The black and green curves, as before (see Fig. 3), represent families of symmetric and asymmetric periodic solutions, respectively. These two curves intersect each other at a branching point BP3 found for $\alpha \u22480.6279$. On the green branch (asymmetric solutions) two additional critical points are found, labeled SB3 ($\alpha \u22480.6819$) and PD1 ($\alpha \u22480.7365$), where the periodic solution touches tangentially the overturning boundary $\phi =\xb1(\theta b\u2212\pi 2)$ and a period-doubling bifurcation of limit cycles occurs, respectively. At the PD1 point, the original period-1 solution loses stability, while a branch of stable period-2 orbits is born. The numerical continuation of these orbits is given by the red branch shown in Fig. 4. Further numerical investigations reveal that the period-2 solution bifurcates once more via another period-doubling point (PD$2$, $\alpha \u22480.7542$), hence producing periodic solutions with four times the original period.

In order to gain a deeper understanding of the dynamics of the rocking block model (9), we will carry out a parametric investigation of the system via direct numerical integration as follows. First, we fix an initial value for the amplitude ($\alpha =0.5$) and then integrate the system over 300 periods of excitation in order to eliminate transients. Then, we extend the numerical integration for a range of 100 periods and store samples of the obtained solution at every $t=0,2\pi ,4\pi ,\u2026$, after which $\alpha $ is increased by a small amount and the process is repeated, using the last sample as starting value for the next step. The result of this numerical process is presented in Figs. 5(a) and 5(b), showing the angular velocity $\psi $ on the vertical axis. The picture confirms the numerical predictions obtained via numerical continuation described above, in particular, the period-doubling bifurcations detected during continuation (labeled PD1 and PD2). As can be seen in Fig. 5(b), the numerical investigation reveals the presence of a period-doubling cascade producing periodic solutions of increasing period, leading to chaotic behavior, as shown in Figs. 5(c)–5(f).

### B. Two-parameter analysis of the system response

In Sec. IV A, our main concern was to investigate the behavior of the rocking block model (9) under one-parameter perturbations, in particular with respect to the frequency $\omega $ and amplitude $\alpha $ of the ground oscillation. This study revealed the presence of critical parameter values (corresponding, for instance, to period-doubling and fold bifurcations) upon which the system dynamics suffers significant changes. In the present section, our main goal will be to investigate how these critical points are distributed in the $\alpha $–$\omega $ plane so as to be able to classify the system dynamics with respect to the excitation parameters. For this purpose, we will employ COCO’s numerical routines for the two-parameter continuation of codimension-1 bifurcations. Specifically, we will carry out the two-parameter continuation of the critical points found in Figs. 3(a) and 4(a), with respect to $\omega $ and $\alpha $.

The outcome of the numerical procedure described above can be found in Fig. 6. Panel (a) presents the resulting bifurcation diagram showing the continuation of period-doubling (green curve), branching points (purple curve), and fold (orange curve) bifurcations. In addition, a red curve is included that corresponds to the two-parameter continuation of the points SB$i$ found before, which gives a parameter-dependent family of period-1 solutions making tangential contact with the overturning boundary $\phi =\xb1(\theta b\u2212\pi 2)$. Consequently, this red curve provides critical information regarding the rocking block dynamics, as in this way, we can determine combinations of frequency and amplitude of the ground oscillation that may lead to the overturning of the structure. Panels (c) and (d) in Fig. 6 illustrate the meaning of this curve, giving examples of period-1 solutions with and without overturning possibility, respectively. These solutions also demonstrate the meaning of the purple curve shown in the bifurcation diagram, which corresponds to the two-parameter continuation of the branching points (labeled BP) found in Figs. 3(a) and 4(a). As explained before, this bifurcation defines the connection between symmetric and asymmetric period-1 solutions, where the region producing such types of solutions is colored in yellow and gray in Fig. 6(a), respectively. In panel (b), a phase portrait of a period-2 solution is depicted, owing its presence to the period-doubling phenomena encountered in the system (green curve). Moreover, notice that the numerical continuation of period-doubling (green curve) and fold (orange curve) bifurcations define a region in the $\omega $–$\alpha $ space [see Fig. 6(a)] of stable period-1 responses without any crossing with the critical boundary $\phi =\xb1\pi 2$, where overturning takes place. Further numerical tests [shown in Fig. 6(e)] reveal that below the fold curve, the system presents aperiodic responses where overturning occurs, as confirmed at the test points OV1–OV6 shown in the figure. In this way, the computed fold curve can also serve as a critical boundary in the parameter space so as to determine under which oscillation conditions overturning phenomena can be avoided.

## V. BASIN STABILITY ANALYSIS

In this section, we analyze previously presented results using the basin stability approach. As aforementioned, in the case of a single degree-of-freedom system, considered in this paper, the basin stability approach can be used to calculate basins of attractions when we draw two initial conditions (for fix parameter values of the system). Before the analysis of obtained results, we list solutions types that exist in the system:

equilibrium—rocking block is not moving and $\phi =\psi =0.0$,

period $1$ symmetric solution—rocking block is oscillating periodically with symmetric orbit,

period $1$ asymmetric solution—period $1$ solution is shifted in a positive $\varphi $ direction,

period $1$ asymmetric solution—period $1$ solution is shifted in a negative $\varphi $ direction, and

period $2$ solution bifurcating from both asymmetric solutions (we do not distinguish the asymmetric solution that it bifurcates),

chaotic solutions or periodic with a period above $2$,

overturning solution where the block reaches $\varphi =\pi /2$, and

overturning solution where the block reaches $\varphi =\u2212\pi /2$.

Hence, we distinguish eight solutions. Due to the very narrow range of existence of solutions with a higher period than $2$, we decide to present them together with chaotic oscillations.

To have overview of the structure of the phase space, we perform analysis for six selected values of $\omega $ ($\omega =1.0$, $\omega =1.16$, $\omega =1.5$, $\omega =2.0$, $\omega =3.0$, and $\omega =4.0$) and a constant value of the amplitude $\alpha =0.66$ (see bifurcation diagram in Fig. 4 and horizontal line in Fig. 6). In Fig. 7, we show calculated basins of attractions. Initial conditions are drawn from ranges: $\phi \u2208(\u2212\pi /2+\theta b,\pi /2+\theta b)$ and $\psi \u2208(\u22121,1)$ and for each diagram we calculate $20000$ Bernoulli trials. For $\omega =1.0$ [panel (a)], we observe period $1$ symmetric solution shown in yellow color. Most likely, it occurs when the initial position is negative. For larger values of the initial velocity, we observe domination of the overturning solution ($\phi =\pi /2$, green area), while for negative initial position and velocity, we see the overturning of the block in the opposite direction ($\phi =\u2212\pi /2$) indicated by purple area. Moreover, we observe a large range where the block is stable in the equilibrium position and does not oscillate. In panel (b), the value of the frequency of excitation is equal to $\omega =1.16$, and it is above the BP bifurcation line. Hence, we observe two period $1$ asymmetric solutions (red and blue colors). Similarly, as for $\omega =1.0$, positive values of the initial velocity result in the overturning of the block. The range where the block reaches the equilibrium position shrinks significantly, and this solution is less probable. Next, the plot (c) is computed for $\omega =1.5$, and we see that with the increase of the excitation frequency, the probability of the overturning is higher, while the equilibrium solution is no longer observed. For negative values of the initial velocity, we see the domination of the period $1$ asymmetric solution. Next, for $\omega =2.0$ [panel (d)], we see that the range of the overturning in the positive direction ($\phi =\pi /2$) is still increasing, while the overturning via negative $\phi $ has a very low probability of occurrence. Other solutions are moved toward negative values of the initial velocity. The same tendency is observed with a further increase of $\omega $ presented in panels (e) and (f). The result presented in the last panel (f) is calculated for $\omega =4.0$, which is below the BP bifurcation line; hence, we observe period $1$ symmetric solution (yellow area). Summarizing, the basins of attraction are compact and only at their borders, we observe complex structures. We calculate the probability of reaching all considered solutions. For values of $\omega $ for which we observe asymmetric solutions, probabilities of their occurrence are equal; hence, neither of them is more likely to be observed assuming a random initial state. For the higher frequency of excitation, the probability of overturning the block increases, and this solution becomes dominant. For the analyzed set of the system parameters, many solutions coexist; hence, the system is multistable. Based on convergence analysis, we see that it is enough to calculate $5000$ Bernoulli trials to have precise information about the phase space structure (we show more to have a densely filled plot). Please note that an important factor that influences the final solution is the form of the excitation function. In this paper, the excitation amplitude is harmonically modulated by $sin\u2061(\omega t~)$ function, which in the first phase of motion moves block into the right direction; hence, when the initial velocity is positive, it can more easily cause the block to overturn.

In the next step, we analyze the response of the system in two-parameter space ($\alpha $, $\omega $) in the following ranges: $\alpha \u2208(0.1,0.9)$ and $\omega \u2208(0.5,4.0)$ for fixed initial conditions: $\phi =0.01$ and $\psi =0.0$. Hence, the block is slightly angled into the positive $\phi $ side. It allows us to study which solutions are reachable from values of initial conditions close to the equilibrium position, which has practical importance. In the description of the previous figure, we conclude that the initial phase of the excitation significantly influences the response of the system; thus, we also take it into account in our calculations. We introduce an initial phase of excitation $\beta $ to harmonic modulation of the excitation function: $sin\u2061(\omega t~+\beta )$ [see Eqs. (6) and (7)]. The diagram obtained for $\beta =0$ is shown in Fig. 8(a) and for $\beta =\pi /2$ in Fig. 8(b). To compare the ranges of solutions’ stability calculated based on the basin stability method and the path-following analysis, we copy bifurcations lines from Fig. 6—namely, the SN (orange), the BP (purple), and the PD (green) line. In both plots, we do not see the difference in the range where the equilibrium solution (black dots) exists. However, based on bifurcation lines, we see significant changes above the SN line. The basins of attractions of other solutions for the two considered values of $\beta $ differ significantly. In panel (a), above the equilibrium solution, we observe the range of period $1$ symmetric solution (yellow dots), which is bounded from the top by the BP bifurcation line and the overturning solution (green dots). Above the BP bifurcation line, we should observe the period $1$ asymmetric solutions (red and blue dots) but they are present only in a very narrow range. This confirms that the selected phase of excitation leads to the overturning of the block, especially for higher values of the excitation frequency. Above the PD bifurcation line, we see the appearance of period $2$ asymmetric solutions (we do not distinguish the) and then fast transition to chaos via period-doubling sequence. Above the chaotic range, we observe overturning. For the second value of the excitation phase $\beta =\pi /2$ [panel (b)] period $1$ symmetric solution is observed nearly in the whole range above the equilibrium solution and below the BP bifurcation line (only in a narrow range the overturning solution appears—purple dots). Between the BP and the PD lines, period $1$ asymmetric solutions (red and blue dots) exist above $\omega \u22481.7$, while below this value, we observe block overturning (purple dots). Crossing of the PD bifurcation line results in a transition to period $2$ solution and then via period-doubling sequence appearance of a chaotic solution. Further increase of $\alpha $ leads to overturning of the block (green and purple dots). Presented plots show that the initial phase of the excitation is crucial for the stability of the block. Hence, the initial phase of motion strongly affects the shape of basins of attractions and the final state of the rocking block. The exception is the equilibrium solution that is not sensitive to the initial phase of the excitation $\beta $.

The subsequent analysis is devoted to the study of probability to reach previously presented solutions. The obtained results are shown in Fig. 9. We draw initial conditions [$\phi \u2208(\u2212\pi /2+\theta b,\pi /2+\theta b)$ and $\psi \u2208(\u22121,1)$] and two parameter values [$\alpha \u2208(0.1,0.9)$, $\omega \u2208(0.5,4.0)$] for two values of the phase of the excitation [$\beta =0.0$ (right column) and $\beta =\pi /2$ (left column)]. For each value of $\beta $, we perform $1000000$ Bernoulli trials, and we show the probability plots in ($\omega \u2212\alpha $) plane for four solutions, namely, the period $1$ symmetric, the period $1$ asymmetric (two types together), the overturning (without distinguishing direction), and the equilibrium. We do not show periodic solutions with periods higher than $2$ as well as neither periodic solutions because their probabilities are low in comparison to the presented ones. To show the probability of occurrence of each solution, we divide the whole range into $900$ cells ($30\xd730$), and we calculate the probability to reach the given solution for each cell. We also add the SN, the BP and the PD bifurcations lines obtained via path-following. In panels (a) and (b), we present the period $1$ symmetric solution, for which the phase of the excitation $\beta $ has a significant influence. For $\beta =0$, it is less probable than for $\beta =\pi /2$. Moreover, for $\beta =0$, it exists closer to the SN bifurcation line than for $\beta =\pi /2$, which is equally probable in the whole existence range. Next, we consider the period $1$ asymmetric solutions [panels (c) and (d)] and we see that for $\beta =\pi /2$, they are more likely to occur and the probability is equally spaced between the bifurcations lines. Then, we investigate the occurrence of the overturning solution [panels (e) and (f)]. We see that for $\beta =0$ with an increase of $\alpha $, the overturning solution starts to dominate, and at the top of the plot, it is the only existing solution. In the case of $\beta =\pi /2$ the overturning is less probable and its higher probability is reached for $\alpha >0.8$ and between $\omega =1$ and $\omega =3$. Last two panels (g) and (h) show the equilibrium solution. We see that there is no visible difference between both plots; hence, we claim that stabilization of the block in the steady-state is independent of the phase of the excitation $\beta $. The main conclusion from this analysis is the importance of the initial phase of excitation that has a dominating role and significantly changes the response of the system.

In the last considered case, we assume the random initial phase of excitation. The results are shown in Fig. 10. We draw initial conditions [$\phi \u2208(\u2212\pi /2+\theta b,\pi /2+\theta b)$ and $\psi \u2208(\u22121,1)$] and three parameters [$\alpha \u2208(0.1,0.9)$, $\omega \u2208(0.5,4.0)$ and $\beta \u2208(0.0,2\pi )$] and we perform $2000000$ Bernoulli trials. We show the probability of five solutions (overturning, equilibrium, period $1$ symmetric, period $1$ asymmetric, and all other responses) as a function of the initial deflection [panel (a)], the initial velocity [panel (b)], and the phase of the excitation [panel (c)]. For that purpose, we divide the span of a given parameter into $50$ equal ranges, and in each range, we calculate probabilities to reach all considered solutions. The number of trials in each section is approximately the same ($40000$ with a standard deviation $242$). For both initial conditions, we observe an increase in the probability of the overturning solution when their values move away from zero. The opposite tendency is observed for period $1$ symmetric and equilibrium solutions. The last two solutions are nearly not affected by initial conditions imposed on the rigid block. In panel (c), we see the significant dependence of solutions’ probabilities on parameter $\beta $. It is visible that when the excitation starts from zero and increase in positive ($\beta \u22480$) or in negative ($\beta \u2248\pi $), the overturning occurs more likely than in case when the excitation starts from maximum ($\beta \u2248\pi /2$) or minimum ($\beta \u22483\pi /2$) values. A different response is observed for the period $1$ symmetric and asymmetric solutions. Equilibrium and other solutions are not affected by parameter $\beta $.

## VI. CONCLUDING REMARKS

This paper carried out a detailed numerical study of the dynamical response of a rocking block under ground motion. The considered physical model can represent several practical scenarios, such as buildings, silos, containers, scaffoldings, and free standing power transformers. One of the main concerns of this work was to determine the conditions under which overturning of the underlying structure may occur. Therefore, special attention was paid to the periodic solutions of the system and their evolution under the variation of base excitation. To study the behavior of the model, two different types of numerical approaches, namely, direct numerical integration and path-following techniques for non-smooth dynamical systems, were adopted. The path-following analysis via COCO revealed the presence of critical parameter values (corresponding, for instance, to period-doubling and fold bifurcations) upon which the system dynamics suffer significant changes. In addition, branching points leading to symmetry breaking were detected, where families of symmetric and asymmetric periodic solutions collide. During this study, two critical boundaries were defined [$\phi =\xb1(\theta b\u2212\pi 2)$, labeled SB$i$], which give the critical angles beyond which overturning can occur, for instance, when the base excitation suddenly stops. In this way, it is possible to identify ground oscillation regimes with a risk of overturning.

Our work significantly extended the results presented by Hogan.^{7,9,14} The majority of his investigations were performed for the slender structures. Such an assumption allows us to linearize the equations of the motion and applies analytical methods. In our work, we assume that the ratio between the breadth and the height of the structure is equal to $1.21$, which is a typical value for this class of systems. Comparing our work with Hogan’s results, we see that the structure’s slenderness significantly influences the system’s response, showing the dynamics of the real structure, thus one needs to study the full nonlinear model in this paper. Due to the limitation of the numerical methods adopted in Hogan’s works, he showed a few simulations of non-slender systems’ behavior. Hence, he did not present a complete view of the dynamics in the frequency and amplitude of the excitation space. Our numerical investigation revealed the complex dynamical behavior of the rocking block model, in particular, due to the presence of period-doubling cascades producing periodic solutions of increasing period, leading to chaotic motion. Further studies were carried out using the two-parameter continuation of the codimension-1 bifurcations found before, using the frequency and the amplitude of the base excitation as the main control parameters. A critical curve corresponding to the two-parameter continuation of the points SB$i$ found in the one-parameter analysis was identified. It gave a parameter-dependent family of period-1 solutions making tangential contact with the overturning boundary $\phi =\xb1(\theta b\u2212\pi 2)$. This curve provides critical information regarding the dynamics of the rocking block, since in this way, a combination of frequency and amplitude of the ground oscillation, which may lead to the overturning of the structure, can be determined. The results from the path-following method are supported with the basin stability analysis. We studied the coexistence of solutions and the influence of the system parameters and initial conditions on the probability of their occurrence. These are the further developments of this work based on Hogan’s study,^{7,14} where he reported the existence of such behavior in the rocking block. We found that in the considered range of the amplitude and the frequency of the excitation, the system is mostly multi-stable, and the final solution strongly depends on the initial conditions and the initial phase of excitation. It is especially essential from a practical point of view because we usually do not have control on the excitation signal. Thus, the parameters of the block should be carefully chosen based on the most probable excitation function.

This study provides crucial information regarding the periodic and non-periodic response of the rocking rigid block subjected to excitation function with varying parameters. The obtained results can be used for development and practical applications to obtain a given type of stable oscillation or prevent overturning under sea waves excitation.

## ACKNOWLEDGMENTS

This work was supported by the National Science Centre, Poland (Project No. 2018/31/D/ST8/02439) (P. Brzeski and P. Perlikowski). Dr Y. Liu was supported by EPSRC under Grant No. EP/P023983/1.

The authors declare that they have no conflict of interest concerning the publication of this manuscript.

## DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request.

## REFERENCES

*Piecewise-Smooth Dynamical Systems. Theory and Applications*, Applied Mathematical Sciences Vol. 163 (Springer-Verlag, New York, 2004).

*Numerical Continuation Methods for Dynamical Systems*, Understanding Complex Systems, edited by B. Krauskopf, H. Osinga, and J. Galán-Vioque (Springer-Verlag, 2007).

*Recipes for Continuation*, Computational Science and Engineering (SIAM, Philadelphia, 2013).