Rigid impacts of three-dimensional rocking structures

Studies of rocking motion aim to explain the remarkable earthquake resistance of rocking structures. State-of-the-art assessment methods are mostly based on planar models, despite ongoing efforts to understand the significance of three-dimensionality. Impacts are essential components of rocking motion. We present experimental measurements of free-rocking blocks on a rigid surface, focusing on extreme sensitivity of impacts to geometric imperfections, unpredictability, and the emergence of three-dimensional motion via spontaneous symmetry breaking. These results inspire the development of new impact models of three-dimensional facet and edge impacts of polyhedral objects. Our model is a natural generalization of existing planar models based on the seminal work of George W. Housner. Model parameters are estimated empirically for rectangular blocks. Finally, new perspectives in earthquake assessment of rocking structures are discussed.


Introduction
Rocking is the typical response of many structures to dynamic loads such as earthquakes. Rocking structures include masonry columns, arches, walls [7,15,16,23], innovative earthquake protection structures such as rocking shear-walls and frames [14,38], and some tanks, machines or vehicles [4,18,22]. Understanding the non-smooth dynamics of rocking motion is essential to properly design and to assess the safety of these structures.
Rocking motion of quasi-rigid objects is an example of contact-induced hybrid multibody dynamics, where episodes of continuous motion are interrupted by short impacts. Rocking is inherently nonlinear, which prevents the use of the theory of linear vibrations [21]. Moreover the time history of motion is notoriously unpredictable in many cases for two reasons [8,35]. The lack of reliable and accurate impact models is a central problem of research concerning rocking motion. Modeling friction and transitions between stick and slip is another crucial question, even though the motion of many rocking structures (such as slender blocks) appears to be free of slip, and in such cases impact models are the only source of difficulty.
Impacts between hard objects have a very short duration. Accordingly, impact models used in the context of rigid body dynamics often take the mathematical form of an instantaneous mapping assigning a post-impact velocity to given pre-impact velocity state. Such discrete models are phenomenological descriptions of a complex, multi-scale physical process, and thus they are often unable to provide accurate predictions. In the case of straight impacts at a point contact, the coefficient of restitution is the standard phenomenological parameter, describing the motion in the normal direction. In the case of oblique frictional impacts, additional parameters are required to describe how the tangential and rotational motion of the colliding objects change during the impact process. Among others, a Coulombtype coefficient of friction and the tangential coefficient of restitution are often used in this context. The parameters of impact maps are most often determined empirically and they are affected by various factors including material properties of the colliding objects (stiffness, strength, ductility, etc.), and other mechanical parameters (e.g., masses, velocity of collision, frictional characteristics of surfaces) [35].
The impacts of rocking objects occur in the presence of linear contact (Fig. 1a, b) or surface contact (Fig. 1c). Both contact setups are geometrically degenerate in the sense that points within a spatially extended region establish contact simultaneously. Degeneracy implies that the impact map is significantly affected by those geometric imperfections of the contact region, which are comparable in size to the deformations of the objects during an impact. In the case of light impacts of a quasi-rigid block on a hard surface (a common situation in the context of rocking motion), this mechanism manifests itself as extreme sensitivity to geometric imperfection, which is in sharp contrast with the more predictable motion of objects on a compliant surface [12]. For example, the slightest concavity of the base of a rigid, planar rocking block results in the concentration of contact forces at the vertices (Fig. 2a), whereas a slightly convex, polyhedral base may gives rise to a rapid sequence of very light impacts sweeping along the whole edge (Fig. 2b). Obviously, the two types of impact affect the angular velocity of the block differently.
Unlike other factors affecting the outcome of an impact, geometric imperfection varies among specimens in an unpredictable way. We believe that it has crucial role in the commonly observed unpredictability and irreproducibility of motion trajectories observed in physical experiments. The presence of such an effect limits the use of the traditional approach of impact modeling: empirically calibrated impact models cannot predict time histories. This feature of geometric imperfections is our motivation to develop an impact model that takes into account the effect of imperfections. We think that such a model is an important initial step toward the ability to assess the safety of rocking structure in the presence of unpredictability and irreproducibility.
Housner's classical impact model for a planar rocking block [23] assumes that the impulse transmitted during a rocking impact occurs at the corner of the block, which is consistent with the assumption of a slightly concave base surface (Fig. 2a). In contrast, physical experiments clearly indicate significant deviation from the predictions of Housner's model [3,5,20,28,31,32,36]. More recently, several authors proposed empirical model corrections [3,5,20,28,31], or a priori chosen geometric imperfections [36], to match the mean value of experimental measurements. Very recently, a few authors have proposed to consider the possibility of arbitrary geometric imperfections, in order to reproduce not only the mean, but also the observed variability of experimental results [10]. In the present paper, we revisit the problem of the planar rocking block, and present previous results using a unified approach in Sect. 2.
While the planar models uncover important aspects of rocking motion, real rocking motion is three dimensional, either because of out-of-plane excitation or  [23] (a). A slightly convex edge (b) induces no energy loss. The models proposed by Ther and Kollár [36] (c), and by Kallitzonis [24] (d) are consistent with other types of imperfection. For an arbitrary imperfection, the possible outcomes of the impact can be parameterized by a single scalar λ representing the location of the resultant ρ of the impulsive forces (e) as the result of spontaneously emerging out-of-plane motion even under planar initial motion and excitation.
Relatively, few works treat rocking as a 3D problem. The rocking motion of a rigid cylinder has been known for long time to be inherently three dimensional [34]. Rocking cylindrical columns have been studied numerically by using rigid models [25,37] as well as the discrete element method [1]. Some experimental investigations focusing on the spatial rocking behavior of ancient cylindrical columns have also been reported [19,30]. It is notable however that cylindrical blocks behave differently than polyhedral objects as they typically roll smoothly instead of undergoing impacts due to the lack of sharp vertices.
The first numerical model of 3D free rocking motion of a polyhedral block [26] focused on continuous motion, and did not propose an impact model. A similar analysis of an arbitrary rigid body with rectangular base by [40] proposed a natural 3D extension of Housner's impact model by assuming that the impact impulse is concentrated at a vertex of the object. In addition, the effect of slip was investigated in several papers [10,11,13].
Several authors including [11,17,40] pointed out that the overturning of a 3D block can occur under excitations which are lower than those which overturn a corresponding 2D block. Hence, 3D models are important for earthquake assessment. At the same time, none of these works attempted to systematically investigate the set of possible outcomes of three-dimensional impacts and the role of geometric imperfections. In order to fill this gap, we introduce three impact parameters in the case of facet impacts, and two parameters for edge impacts, which capture all physically relevant outcomes of an impact without slip. Thereby, we obtain for the first time a universal 3D rigid impact model. The impact parameters of some free rocking blocks are determined through fitting simulated trajectories to experimental data. We use these results as well as basic physical laws to estimate the relevant ranges of model parameters.
The rest of the paper is organized as follows. In Sect. 2, we review the most common rigid, planar models of rocking impacts in a unified framework focusing on the case of a single monolithic block, and on the role of geometric imperfection. Then, we present an experimental demonstration of transition from planarfree rocking motion into three-dimensional rocking by a spontaneous symmetry breaking. This experiment illustrates the crucial role of geometric imperfections. In Sect. 4, we develop the new, three-dimensional model of rocking impacts, which is highly analogous to planar models, but it can account for spontaneous symmetry-breaking. The experimental results are revisited in Sect. 5, and the parameters of the new impact model are fitted empirically. Section 6 concludes the work, and outlines future steps toward the successful application of spatial impact models in earthquake design.

A review of planar impact models
Consider a rigid, planar, rectangular object with mass m, and mass moment of inertia about the center of mass (COM) θ . Immediately before the impact, the block rocks around vertex V b , until the base V b V a hits the When the base of the block hits the ground, all points along V b V a come into contact simultaneously, and contact forces may emerge anywhere along V b V a . Housner's model [23] assumes that H1: immediately before the impact, the block rotates about V b H2: the rotational momentum of the block about the vertex V a colliding to the ground is preserved H3: the vertex V a stays in sticking contact with the ground after the impact These assumptions determine uniquely the post-impact angular velocity of the block in terms of the pre-impact value as follows. According to assumption H1, the velocity vector of yields where r = r Hou := θ + mh 2 − mb 2 θ + mh 2 + mb 2 (3) The scalar r will be referred to as angular velocity reduction factor (AVRF).
If we introduce the following notions -sticking impact: an impact where the impact point has 0 post-impact tangential velocity -inelastic impact: an impact where the impact point has 0 post-impact normal velocity then assumptions H2-H3 can be formulated in a more concise way as: H4: a single, perfectly inelastic sticking impact occurs at the vertex V a While Housner's model does not include explicit assumptions about geometric imperfections, assumption H4 is consistent with certain geometric imperfections ( Fig. 2a) but inconsistent with others. Many experimental works confirmed that the linear relationship (2) is a good approximation but the AVRF (3) derived using Housner's assumptions underestimates experimentally measured values [9,36]. The inaccuracy of the model has been attributed to the fact that H2 is not more than an a priori assumption [33]. To address this limitation, several works suggested improved models using either empirical modifications of the AVRF or modifications in Housner's initial assumptions [6,27,28]. For example, Kalliontzis et al. [24] proposed using a reduced effective width of the base, which was motivated by the fact that small deformations of the underlying surface give rise to a spatially extended contact region along the base during the impact. The effective width introduces a free parameter 0 ≤ ν ≤ 1 into the model and the corresponding value of r can be anywhere in the interval r ∈ (r Housner , 1). Ther and Kollár [36] pointed out the crucial role of geometric imperfections in the case of a hard underlying surface. In order to improve model accuracy, they proposed an a priori infinitesimally small geometric imperfection (Fig. 2c), which implies that a rocking impact consists of a sequence of two elementary impacts at points M and V a of Fig. 3. The smallness of the imperfection means that its effect on the physical properties (m, θ ) of the block is negligible. Furthermore, the points V b , M, V a are almost collinear, and the elementary impacts occur within very short time. This assumptions allowed them to estimate the effect of the two elementary impacts according to H5: The outcome of a sequence of elementary impacts can be calculated as if they occurred immediately after one another in the nominal impact configuration of the perfect block.
In addition, H4 was replaced by a more general assumption H4a: whenever a point of the imperfect surface hits the ground, a perfectly inelastic sticking impact occurs.
They showed that the corresponding AVRF r = r TK := θ + mh 2 θ + mh 2 + mb 2 (5) closely matched the mean value of many physical experiments. Wittich and Hutchinson [39] examined a model with arbitrary geometric imperfection and demonstrated how the corresponding value of r can be determined under assumptions H1, H4a, H5. Interestingly, the previously mentioned model of Kallitzonis [24] is also equivalent of assuming a hard surface with geometric imperfection in the form of a concave section surrounded by two identical convex segments (Fig.  2d). Chatzis et al. [10] seems to be the first one to notice that the linear momentum transferred along the base during the impact process can be replaced by an instantaneous resultant momentum [P x , P y ] T under assumptions H1, H3. All possible outcomes of the impact can be parameterized by a single scalar λ representing the position of the resultant impact force (Fig. 2e). Accordingly, they proposed to use assumptions H1, H3 together with H2a: the rotational momentum of the block about a specified internal impact point of the base given by the parameter λ is preserved The AVRF now becomes Importantly, the model of [10] is universal in the following sense: it captures all possible post-impact states through a single scalar parameter λ. On the other hand, the approach of [10] does not give a hint on how to choose λ. Clearly, many different factors may affect λ (e.g., material properties, impact velocity, etc.). Among these, the key role of geometric imperfections was first pointed out by [36,39].
One can also set up theoretical bounds on the parameter λ. The straightforward assumptions of nonnegative energy absorption implies 0 ≤ λ, whereas unilateral nature of contact forces implies −1 ≤ λ ≤ 1. Housner's assumptions correspond to λ = 1, and a slightly convex base (as in Fig. 2b) would induce an impact with λ ≈ 0 under assumptions H1, H4a, and H5. We notice without proof that all values in the interval λ ∈ (0, 1) can be realized by choosing an appropriate geometric imperfection. Thus, the resulting AVRF r Cha (λ) can take arbitrary values within the interval (r Hou , 1).
It is noticeable that some experimental works report AVRF values below r Hou , which would correspond to λ > 1 [5,20]. Impact parameter values beyond theoretical bounds are typically caused by unmodeled effects such as slip [5] or additional energy absorption of pin joints in the experimental device [20] .
In Sect. 4, we develop a similar model of threedimensional rocking impacts involving edge and surface contacts. Our contribution consists of a parameterization of all possible outcomes in order to ensure universality of the new model. Theoretical limits of the model parameters will also be identified.

The experimental setup
The free rocking motion of eight stone blocks ( Fig.  4) with identical nominal heights h and half width b 1 but different half depths b 2 has been investigated. The blocks were manufactured from the same solid granite slug. Examining the sliced surfaces, no cavities or cracks were found, so the blocks were considered to be homogeneous. The density of the granite blocks is ρ = 2692 ± 43 kg/m 3 . The motion of the blocks was recorded with the aid of an X-IMU inertial motion unit device produced by X-IO Technologies. We use an orthogonal coordinate frame fixed to the rocking block ( Fig. 4) such that the x-axis is parallel to the edge of length b 1 and the zaxis points vertically up. During the experiments, each block shown in Fig. 4 was tilted in the xz plane and released five times. Then, each block was turned upside down, and the experiments were repeated. The two settings will be referred to as positions A and B.
The bodies were released from a state of edge contact with initial angular velocity close to 0 and an initial inclination angle slightly below the neutral position of the block (ϕ n = 0.245 rad). The initial inclinations were for all experiments in the range of 0.17 rad < ϕ 0 < 0.24 rad. The support medium was a solid granite block of size 50 by 40 by 20 cm resting on a massive steel support.
Before the recorded experiments, some preliminary tests were run using 120 fps action cameras (GoPro Hero 8 Black) to investigate the effect of sliding and bouncing during impact. The front-and side-view of the block were recorded, and the frames of the movie were analyzed. No visible signs of sliding and or bouncing were observed. Figure 5 shows time history data of three experimental tests. The first two belong to granite block three, but the experiment was carried out in positions A and B, respectively. The last test belongs to block eight, which has a larger depth. Panels (a, d, g) show components of the angular velocity, where nonzero values of ω x or ω z are both signs of spatial motion. Panels (b, e, h) depict Euler angles ϕ x , ϕ y and ϕ z representing the attitude of the block using the zyx convention. (That is, a general rotation is the composition of three consecutive rotations about the local z, y, and x axes by angles ϕ z , ϕ y and ϕ x .) The X-IMU estimates time histories of the acceleration vector, angular velocity vector, attitude, and local magnetic field vector using on-board sensors including accelerometers, a triple-axis gyroscope, and magnetometers. Two operation modes are offered by the device. The attitude heading reference system (AHRS) modes combine all sensor data including the magnetometer. We observed that the big amount of steel equipment in the laboratory caused large systematic errors, and thus this mode could not be used. The inertial measurement unit (IMU) mode excludes magnetometer data. This mode provides sufficiently accurate angular velocity estimations, but large drift was observed in some of the attitude estimations (see Fig.  5e). Throughout the experiment, we used IMU mode with 256 Hz rate, and the drift of attitude data was compensated as described below (Fig. 5c, f, i). The modified Euler angles were determined as follows: -The local maximum points of |ω y | in the angular velocity diagram were located numerically. -These points mark the moment of an impact in which an edge of length b 1 is in contact with the ground. Thus, the ϕ y Euler angle should be zero -A piecewise linear error function was added to the Euler angles to enforce 0 value of the appropriate components at the impact times.
The same correction steps were also repeated for the ϕ x Euler angles.

Results of the experiment
Despite the xz plane being a plane of reflection symmetry, significant out-of-plane motion was observed in almost all cases (Fig. 5). In addition, large differences were observed between the A and B positions, whereas multiple trials in the same position yielded more similar results. Both phenomena are the results of manufacturing imprecisions of the object, and our observations suggest that minor imperfections of the block strongly affect rocking motion. For blocks with large depth, we recorded lateral motion emerging from time to time, but decaying very rapidly, so motion appeared to be planar most of the time (Fig. 5g).
In order to quantify these observations, two quantities were determined for every experimental trial: Then, an AVRF value corresponding to two adjacent maximum values was calculated as r exp,i = √ u i+1 /u i . The final estimated AVRF r EXP was the median of the nine values corresponding to one single experimental trial.
The measured AVRF and OPMF values for different values of block depth (b 2 ) are shown in Fig. 6. The out-of-plane motion factor was found to be as high as 0.1 for low values of b 2 , which is the result of spontaneous symmetry breaking induced by geometric imperfections. High dispersion of the AVRF data is likely also caused by imperfections, which can be explained by planar models. However the AVRF data also suggests that average values depend on b 2 , which cannot be explained by planar models. The amount of experimental data is insufficient to draw detailed conclusions about the role of b 2 ; however, the presented results clearly indicate that such an effect exists. This is our motivation to present a new three-dimensional impact model.

A new, three-dimensional impact model
This section is devoted to development of our new impact model. First, a general description of hybrid rocking motion is given to clarify the role and the types of impacts. Then, the basic idea of the new impact model is presented in Sect. 4.2. This is followed by a general formula applicable to all types of impact maps (Sect. 4.3). The initial formula is further developed into a fully explicit expression in terms of impact parameters, which is applicable to edge impacts (Sect. 4.4) and to two different types of facet impacts (Sects. 4.5-4.6). A unified map applicable to any facet impact is given (Sect. 4.7). Finally, we establish some theoretical bounds of impact parameters (Sects. 4.8-4.9) .

The role of impacts during rocking motion
Three-dimensional rocking motion of quasi-rigid, convex, polyhedral objects is a combination of continuous and discrete components. Contact-free motion is barely observed except for perhaps very small amplitude bouncing motion. Accordingly, models of rocking motion often ignore the possibility of free flight. Slip may or may not be significant depending on contacting materials and geometry of the objects. Here, we keep focusing on slip-free motion, which is the most common behavior of slender blocks unless the underlying surface is slippery. With these restrictions, the possible contact modes of motion are surface contact (i.e., rest), edge contact (rocking around the edge) and vertex contact (rocking possibly with spin component). Transitions between these modes occur through impacts. More specifically, rocking under vertex contact typically ends when all points along an adjacent edge of the block reach the ground simultaneously giving rise to an edge impact. The typical post-impact motion of slender rocking object is rocking about the other endpoint of that edge. In a similar fashion, rocking about an edge terminates when all points of the base facet reach the ground simultaneously, giving rise to a facet impact. Then, the typical post-impact mode is rocking about a vertex or another edge of that facet.
In addition, vertex contact terminating in a facet impact is also possible but non-generic, so it is not discussed here.
It is worth noting that sustained edge or facet contact (rest) can not be created by any of the single impact Fig. 7 Transitions during two-and three-dimensional rocking motion transitions listed above. One of the special features of rocking motion is the accumulation of infinite sequence of impacts in finite time intervals, which is often called Zeno behavior [2,8]. In planar models, a Zeno point occurs when motion stops. For three-dimensional models, a Zeno sequence of edge impacts may terminate in a state of rest or in a state of sustained edge contact (i.e., rocking). A Zeno sequence of facet impacts may only terminate in a state of rest. Figure 7 provides a schematic overview of the possible contact modes and generic transitions between them, including Zeno behavior.

Overview of the impact model
Similarly to planar impacts, three-dimensional rocking blocks have sustained contact at an edge or at a vertex immediately before an impact. Hence, it is reasonable to use an extension of assumption H1 (see Sect. 2). In addition, it is also reasonable to ignore the possibility of free flight and slip, similarly to assumption H3 in Sect. 2.
Analogously to planar rocking, we will approximate the impact by an instantaneous mapping of velocities. Then, the outcome of the impact is uniquely determined by the spatial distribution of impulsive forces transferred to the object along the edge or facet involved in the impact, and the exact time history of force during the impact process does not matter. Moreover, we can further simplify the description by considering the resultant of the impulsive forces in the spirit of Chatzis et al. [10].
There are however some differences between planar and spatial impacts. We have seen that an arbitrary distributed planar impulse is equivalent of a single resultant impulse ρ, which can be represented by its components ρ x , ρ y and its location x R (Fig. 8a). Two parameters are determined by the kinematic constraints associated with no slip and no liftoff, and thus the possible outcomes of a planar impact can be param- In contrast, the resultant of a spatially distributed three-dimensional impulse has six independent parameters, and representation of the distributed impulse by an equivalent single resultant impulse (which has only five parameters) is in general not possible. In the field of robotics, a distributed force is often represented by a 'wrench': a resultant force vector ρ = [ρ x , ρ y , ρ z ] along an appropriately chosen line of action (specified by two parameters) and an additional torque vector of magnitude τ parallel to ρ. This representation is unique. Impulses can be represented in a completely analogous way by a resultant impulse and an angular impulse, which can be parameterized by three components ρ x , ρ y , ρ z , two coordinates x R , y R of the resultant ρ as well as by the magnitude of the torque τ (Fig.  8b).
For convenience, we will use a slightly different representation of the distributed impulsive impact forces by using a 'modified wrench', which consists of the resultant impulse ρ = (ρ x , ρ y , ρ z ) and an additional vertical angular impulse vector of magnitude τ . The coordinates x R , y R are used to specify the line of action of ρ (Fig. 8c, d).
The modified wrench representation is unique provided that ρ z = 0. Its parameters can be calculated in the case of facet impact as follows. Assume that the distributed impulse acts over a facet F lying in the x − y coordinate plane and it is given by the function p(x, y) with components p x (x, y), p y (x, y), p z (x, y). Then, the parameters of the modified wrench are calculated by surface integrals: For edge impacts, p is distributed along that edge, and analogous linear integrals are used to determine the resultant. If the distribution of the impulse is discrete, then summation can be used instead of integration.
The main advantage of the modified wrench representation is the fact that for arbitrary p z (x, y) ≥ 0, the possible locations of the point R : (x R , y R , 0) are exactly the points of the facet F according to (8)- (9). Similarly, for an edge impact, the possible locations of R : (x R , y R , 0) are exactly the points of the edge V b V a involved in the impact. The latter will allow us later to use a single scalar parameter λ instead of x R and y R .

General formula of the impact map
Consider a three-dimensional impact of a rocking block with mass m and moment of inertia tensor Θ. Let u x , u y , u z denote unit vectors along the coordinate axes. The impact map ω − → ω + can be expressed as follows. The kinematic constraints of rocking motion immediately before and after the impact yield The angular impulse momentum theorem applied to point R yields One can replace the cross products in (11)-(13) by matrix multiplication, using the matrix and two other matrices R R , R a composed in the same way: Then, ω + can be expressed as

Explicit formula for edge impacts
In this case, V b , and V a are the two endpoints of the edge involved in the impact, and R is a point along the edge. Hence, one can introduce the scalar parameter λ such that which also implies that R R = (R a + R b )/2 + λ(R a − R b )/2. Thereby, we obtain from (18) a closed-form expression of the impact map in terms of known quantities and two impact parameters λ and τ :  block is not 0, then τ also affects the third (in this case: y) component.
Later, Sect. 5 shows that one can fix τ = 0, and the remaining impact parameter λ still allows a reasonable fitting of simulated trajectories to experimental measurements.

Explicit formula of facet impact followed by vertex contact
Let V i , i = 1, 2, . . . , n denote the vertices of the facet involved in the impact and r i the corresponding position vectors. Before a facet impact, the object rocks about an edge adjacent to that facet in general, which means two immobile vertices. {1, 2, . . . , n}) be the pre-impact axis of rotation.
(The index b − 1 should be understood modulo n). Immediately after the impact, the object rocks either about a vertex or about an edge of the base facet (see Fig. 7). Here, we discuss the first scenario.
So now, V a , a ∈ {1, 2, 3, . . . , n} is the unique vertex contacting the ground after the impact. As before, the impact is given by the general formula (18).
In this case, r R cannot be expressed in terms of a single scalar parameter λ. Instead, one can use two parameters x R and y R . Alternatively, if the facet is rectangular, then it is convenient to introduce two dimensionless parameters λ lon , λ lat according to Fig. 10b such that Then, the explicit expression of the impact map in terms of the impact parameters is given by (18), and (21) using the notation (14) for fixed values of a.
Importantly, the index a of the post-impact point of contact V a cannot be chosen a priori in the case of a facet impact. The feasibility of all values of a can be tested one by one by evaluating the impact map, and by checking the requirements associated with unilateral contact: for i = 1, 2, . . . , n. The first formula reflects that unilateral contact forces should increase u T z v, and the second one means that none of the vertices penetrates into the ground after the impact.
Ideally, this process should always provide a unique feasible solution; however, this is not true in general. As illustration, we show numerical results of this process for a rectangular block of mass m = 1 kg, and sizes b 1 = 5 cm; b 2 = 3 cm; h = 20 cm in Fig.  11. We assumed b = 2 (i.e., rocking about V 1 V 2 before the impact) and pre-impact angular velocity ω − = [0, −1, 0] T . Each panel corresponds to a fixed value of τ . The figure shows numerically identified regions of the impact parameters where a = 3 and those where a = 4, i.e., the impact is followed by rocking about vertex V 3 and vertex V 4 , respectively. In addition, there are regions of the impact parameters where no feasible solution exists, and others where there are two solutions.
The observed non-existence of the solution means that there are non-trivial ranges of the impact parameters, which are physically impossible. In addition, nonuniqueness means that there are ranges of the parameter values where the chosen parameterization of impact is ambiguous. Fig. 11 Illustration of a facet impact of a block. Background colors show the index a of the vertex in contact with the ground after the impact. Dash-dotted lines depict the singularity of (18). Dashed rectangles depict the region (27); the energy bound (26) is satisfied in the region above the continuous curves In order to get rid of nonuniqueness and nonexistence, we will postulate that i.e., there is no 'spinning' motion immediately after an impact. This assumption can be used to express τ explicitly from (18) after both sides of (18) have been multiplied by vector u T z from left. In the case of a homogeneous cuboid block, we obtain the trivial result τ = 0. Equation (24) is merely an assumption which does not follow from physical laws; however, it is plausible, as dry friction is likely to eliminate the spinning component of motion when a whole facet is in contact with the ground. As we will see, this assumption is also in agreement with experimental results.
Finally, (24) guarantees a unique feasible solution as shown in Fig. 11a. In particular, the impact is followed by rocking about V 3 if λ lat > 0 and rocking about V 4 if λ lat < 0. A sketch of the proof of uniqueness is presented in "Appendix."

Explicit formula of facet impact followed by edge contact
Similarly to the previous type of facet impact, the explicit expression of the impact map is given by (18), and (21) using the notation (14). In contrast to the previous case, there are two immobile vertices immediately after the impact (say V a and V a−1 ); furthermore, the kinematics of rotation implies that (24) must be satisfied. In the case of a slender cuboid block, it is easy to see that V a−1 V a must be the edge of the base facet opposite to V b−1 V b . Furtermore, these constraints imply that λ lat = τ = 0, whereas λ lon can be arbitrary. Hence, we only have one free parameter. This is not surprising given that facet impact followed by edge contact is essentially an example of two-dimensional hybrid rocking motion. Hence, in this case, the impact map is identical to the one-parameter map of planar rocking impacts of [10] revised by us in Sect. 2.

Unified impact map of facet impacts
Our findings in the previous two subsections imply that the two types of facet impacts can be given by one single impact map as shown in Fig. 11a. The impact map has two parameters: λ lat , and λ lon . Every parameter value delivers a unique solution, and the type of post-impact motion is determined by the sign of λ lat . For example, if pre-impact motion is rocking about V 1 V 2 , then postimpact motion of rocking about vertex V 3 if λ lat > 0; rocking about the edge V 3 V 4 if λ lat = 0, and rocking about V 4 if λ lat < 0. The dependence of the impact map on the parameters is shown in Fig. 12. The map is continuous everywhere but non-smooth at λ lat = 0. The dominant effect of λ lon is to set the amount of reduction in the component of ω + , which is parallel to ω − , whereas λ lat controls the emergence of an angular velocity component perpendicular to ω − . Fig. 12 Contour plot of the x and y components of the post-impact angular velocity ω + after a facet impact of a block with m = 1 kg, b 1 = 5 cm; b 2 = 3 cm; h = 20 cm. The pre-impact angular velocity is ω = [0, −1, 0], which corresponds to rocking around edge V 1 V 2 . Note that the third component of ω + is 0 by (24). The component u y ω + is non-smooth function of the impact parameters at λ lat = 0 because the post-impact mode of contact depends on the sign of λ lat If a rocking block has large width to depth ratio (like blocks 7 and 8 in Fig. 4), the experimentally observed motion is usually almost perfectly planar. In contrast, our new impact model introduces large lateral angular velocity component unless λ lat is very close to 0. This apparent contradiction can be resolved by noting the behavior of simulated blocks after such an impact. The facet impact is followed by spatial motion involving a sequence of edge impacts. Large depth value implies that the lateral component of the motion decays very rapidly (relatively to the time-scale of the longitudinal motion) and disappears through a Zeno point even if λ lat was not very close to 0. After this short transient, our model predicts planar motion in agreement with experiments.

Theoretical limits of parameters of edge impacts
Similarly to planar impacts, the parameters λ, τ of edge impacts capture all possible outcomes, but it remains to find theoretical limits of these parameters. As we have seen, λ determines the position of the resultant ρ, which must belong to the segment V b V a . From this observation, we obtain the limits Secondly, the impacts may not increase the kinetic energy of the rocking block. The kinetic energy of a rigid body rotating by angular velocity ω about a fixed point X can be calculated as 1 2 where R x is the matrix composed from r x according to (14). Then, we have The region of the τ − λ parameter plane consistent with the last two contraints is always bounded, as λ is bounded by (25), and for each value of λ, ΔE is a quadratic function of τ with positive leading coefficient. It should also be noted that for τ = 0, (26) implies λ ≥ 0 in analogy with previous findings about planar impacts. Figure 13 illustrates these bounds for various values of ω − . It should be noted that the contour lines are not affected by multiplication of ω − by a constant factor, nor by adding an arbitrary y component to ω − . Another theoretical bound is provided by a singularity of the inverted matrix Θ − mR R R a in (18). We note without proof that such singularites do exist, but for slender blocks, they appear for large values of λ, which are irrelevant due to (25) .
There is at least one more theoretical bound for edge impacts, which reflects the fact that the torque τ originates from frictional forces generated by tangential motion during the impact, thus it cannot be arbitrary.
Finding the induced constraints of τ is beyond the scope of this work.
The energy bound (26) (solid curves in Fig. 11) as well as the singularity of (18) (dash-dotted curves in Fig. 11) are also applicable to facet impacts. Again, we note without proof that the singularity appears for irrelevant ranges of (λ lat , λ lon ) provided that the block under investigation is slender. In addition, there are bounds of τ similarly to edge impacts that we do not investigate in this paper.

Methods
We now revisit the experimental tests introduced in Sect. 2. The motion always starts with the following phases of motion: 1. inital rocking about edge V 1 V 2 2. facet impact 3. motion with sustained contact at V 3 and/or V 4 . This phase includes one or several episodes of rocking about vertex V 3 , vertex V 4 or about edge V 3 V 4 . Each such episode is followed by an edge impact at edge V 3 V 4 . The number of impacts may also be infinite in the case that a Zeno point occurs during this phase of motion. 4. an edge impact or a facet impact when V 1 and/or V 2 reaches the ground again In order to reproduce the observed motion numerically, we developed a custom-made code in MATLAB environment, which simulates rocking motion on any vertex or edge as dictated by the Newton-Euler equations of rigid body motion. This was complemented with event detection algorithms to detect impact times and Zeno points, as well as numerical implementations of the impact maps. The code follows the logical structure of Fig. 7. This code is capable of simulating rocking motion of a block if the physical properties of the block (size, mass, moment of inertia), the initial time of release t 0 , initial tilt angle α, the parameters of edge impacts at edge V 3 V 4 during phase 3 (λ, τ ) and the parameters of the first facet impact in phase 2 (λ lat , λ lon , τ ) are specified.
The aim of this code is to find optimal values of the initial conditions and impact parameters, matching the experimentally measured trajectories. To that end, we introduce the error function (28) where t end is the time when the sequence of episodes listed above ends. ω m , and ω s are measured and simulated angular velocity functions, and ω s depends on the initial conditions and model parameters introduced above. In order to reduce the number of unknown parameters, we fix τ = 0 both for edge and for facet impacts. In addition, we assume that each of the λ, λ lat , λ lon impact parameters take constant values along an individual rocking test. The unknown initial conditions (t 0 , α) and the model parameters λ, λ lat , λ lon have been determined by semi-automated numerical mini-mization of the error function : a rough fitting was first obtained by trial-and-error, the results of which was used as initial guess in the fminsearch numerical optimization algorithm of MATLAB software [29].
An example of trajectory fitting is shown in Fig. 14, and the optimum values are summarized in Fig. 15 for both position (A or B). The fitting of λ lat was omitted for blocks 7 and 8, for which lateral motion decayed very rapidly. In the next two subsections, two important observations about these results are presented.

Variability of impact parameters
The first observation from Fig. 15 is the presence of large variability in the parameter values. This is partly caused by inaccurate measurement data. For example, impacts often cause rapid oscillations in measured angular velocity values (see Fig. 14), which is likely an artifact caused by mechanical vibration of the gyroscopes. Variation of initial tilt angles in the experimental tests may also contribute to variance of fitted parameter values because the impact maps often depend on impact velocities in a nontrivial way. Phenomena neglected by the model (bouncing, slip, elastic deformations, etc.), may also contribute to inaccuracy of fitting and noise in the fitted parameter values. All these factors call for a more extensive and controlled experimental campaign to properly calibrate the impact model in the future.
However, it is important to note that sensitivity of impacts to geometric imperfections could in itself cause large variability. Experiments using positions A and B of the same columns involve different geometric imperfections. Even subsequent tests of the same column in the same position involve different imperfections, due to moving over different areas of the underlying surface or due to the possible presence of small grains under the rocking blocks. The observed variability is thus more than just a weakness of our model, the experiments, or the analysis. Sensitivity implies that the anticipated result of model calibration is a realistic range of model parameters rather than specific values. Accordingly, impact parameters should be treated as uncertain parameters during the analysis of rocking motion. In particular, earthquake response analysis should verify the safety of a rocking structure for a whole range of impact parameter values rather than for single values of those parameters. This can be done via repeated sim-ulations to assess the safety of one single block under a given excitation.

Violation of theoretical bounds
A closer look at the fitted λ lon parameter values reveals that they are in the range (0, 1.5). This is consistent with the bound (26); however, the limit (25) is often violated by λ lon values exceeding 1. In general, the violation of a theoretical bound indicates limitations of a model such as unmodeled effects (e.g., bouncing, slip, pivoting friction, or finite duration of impacts) or some simplifying assumption (e.g., using constant values of impact parameters during each individual motion sequence). In this particular case, large λ means unexpectedly high reduction in angular velocity. We suspect that this may be caused by neglected energy loss during continuous rocking motion due to pivoting friction or micro-slip. The same explanation applies to the unexpectedly high λ parameter values of edge impacts.
The fitted λ lat values are in the interval (−5, 5), which is far beyond the limit (27). As in the previous case, we have no rigorous explanation for the violation of (27); however, we suspect that the results are caused by neglected small-scale bouncing motion as explained below.
One of the widely established assumptions in the literature of planar rocking motion is perfectly inelastic impacts (see assumption H4 in Sect. 2) despite the obvious fact that this assumption is not justified the materials used in experimental studies. Why is then inelasticity considered as acceptable? We believe that episodes of bouncing are indeed present in the experiments; however, they are of short duration and small amplitude. More specifically, rocking motion is a behavior of blocks with relatively small width-to-height (b/ h) ratio (Fig. 3). The kinematics of rotational motion implies that the base of the block hits the ground with low velocity (compared to the translational velocity of other points of the block). Consequently, bouncing motion after the impact will also have low initial lift-off velocity. The emerging bouncing motion is similar to the motion of a rigid bouncing ball, where the process ends with a Zeno point followed by sustained contact. The maximum height of bouncing is proportional to the square of the liftoff velocity and the total duration of the bouncing sequence is linearly proportional to the initial liftoff velocity. That is, both quantities are small.  From a practical point of view, these properties imply that bouncing is a transient behavior that has little effect on the dynamics, and it is also difficult to observe.
The observations presented above can be extended to spatial rocking motion, nevertheless even small-scale bouncing may have large effect on the emerging out-ofplane angular velocity. To illustrate this phenomenon, imagine a facet impact, in which the resultant of the impact momentum acts near vertex V 3 (see Fig. 11). Then, the model predicts rocking about vertex V 4 after the impact. If the impact is not ideally inelastic, then more vertical momentum is transferred to the block than the momentum predicted by the inelastic model. The surplus momentum causes the block to jump in the air, and to undergo bouncing at V 4 instead of smooth rocking. Moreover, the surplus momentum at V 3 during the impact is compensated by less momentum transferred to V 4 after the impact than the momentum generated by pure rocking motion during the same time interval. The surplus momentum emerging at V 3 instead of V 4 amplifies the out-of-plane angular velocity of the block after the facet impact. This is similar to the effect of large |λ lat | values.
Hence, our second important observation is that the optimal fit of impact parameters to experimental results may yield parameter values contradicting previously found theoretical limits. Such violations indicate limitation of the model. We were able to identify neglected phenomena, which have similar effects to parameter values exceeding theoretical bounds. This finding suggests that the practical applicability of the model might be improved if one uses empirically determined ranges of impact parameters exceeding theoretical limits. At the same time, further exploration of this idea is left for future work.

Conclusions
The remarkable earthquake resistance of rocking structures inspires ongoing efforts to understand the characteristic features of rocking. The challenges of the analysis include strong nonlinearity (even at small amplitudes), hybrid dynamics (sudden impacts and continuous motion), and in the case of quasi-rigid blocks and supports, geometric sensitivity associated with edge and facet impacts. Models used for the analysis of rocking motion often use various simplifying assumptions, such as rigidity, inelastic impacts, slip-free motion, and planar motion.
Our work was inspired by the well-known observation that free rocking blocks on hard surfaces tend to transition from planar motion to spatial rocking via spontaneous symmetry breaking. This is caused by microscopic geometric imperfections of the block, which have macroscopic effect on the motion due to the extreme geometric sensitivity of edge and facet impacts of quasi-rigid objects. This phenomenon cannot be explained by currently available models, and calls for a new modeling approach, which is three dimensional and non-deterministic. As a first step toward such a model, we proposed a new universal, threedimensional impact model, which can be used in the context of rigid body dynamics. Universality means that all reasonable outcomes of the impact are described by a small set of impact parameters. Our new model is a natural extension of the planar impact model [10], which uses one single impact parameter.
The parameters of the new impact model have been fitted to some free-rocking experiments. The results show that the parameters are quite unpredictable, which is not surprising as they are influenced by the effect of unknown geometric imperfections. We made initial steps toward identifying possible values of the impact parameters; however, more experimental data will be required to accurately identify their relevant ranges and probability distributions.
The long-term goal of our research effort is to improve current methods of assessing the earthquake resistance of rocking structures. We foresee that the inherently three-dimensional nature of rocking motion affects the ability of a structure to resist earthquakes without overturning. Moreover, the unpredictability of impact parameters suggests that parametric studies may be necessary to assess the safety of a single structure.
It is straightforward to show by direct calculation that v + + ω + × r a is a linear function of ρ for fixed values of a. As a consequence, the no slip conditions and the inequality (22) determine the direction of the vector ρ (i.e., ρ/|ρ|) uniquely. At the same time, inelasticity (31) can be used to determine the length |ρ|. Moreover, for slender blocks and realistic values of impact parameters, u T z (v + +ω + ×r a ) is a strictly increasing function of |ρ|.
If (24) holds, then the no-slip conditions (29)-(30) associated with different values of a become equivalent. Hence, for all values of a, the direction ρ/|ρ| of the impact impulse will be the same. The length |ρ| is however determined by (31) and it will be different for each value of a. Among these values, the largest one will provide a unique feasible solution as it implies that (31) is true, whereas for all values i = a, the inequality (23) will be satisfied instead.
In contrast, if (24) does not hold, then the noslip conditions (29)- (30) are not equivalent, hence the direction of ρ may be different for each value of a. Thus, the arguments outlined above are not applicable.