Nonlinear Dynamics of Swinging Clapper Bells under Arbitrary or Resonant Forcing Functions

Study of swinging clapper bells involves aspects encompassing sound and acoustic engineering, mechanical engineering, and structural engineering. From the musical point of view, clapper bells are directly played idiophone instruments, where the playing device, the clapper, although directly excited, is not explicitly controlled by the bell ringer. The achievement of a clear and optimal sound mainly depends on the acoustic characteristics of the bell and on the regularity of the clapper strokes, which is not only governed by the ringing style and the relevant parameters of clapper and bell but also by the real time corrections to the excitation introduced by trained bell ringers. In fact, despite centuries of experience allowed to optimize the bell performances, standardizing proportions and mounting arrangements, effective sound control requires some fine tuning of the forcing function. Another crucial topic, especially in view of assessing existing structures, regards the evaluation of time histories of the actions transmitted by the bell to the pivots and the study of the interactions between the bell and the supporting structures, belfries, and bell-towers. “Ringability” of swinging bells and bell-structure interactions are usually tackled in the framework of rigid body dynamics, so arriving at an initial value problem, governed by a system of two second order nonlinear ordinary differential equations (ODEs), whose solutions are piecewise-defined functions. In the relevant literature, numerical solutions of the system are commonly sought using built-in algorithms provided in advanced software packages; since the use of such general algorithms is subject to some restrictions, especially regarding the forcing functions, validity of the results is often limited. The present study focuses on an innovative procedure to solve the equations of motion. The method, extremely fast and effective, is based on original numerical explicit-implicit predictor-corrector integration algorithms with constant time step, duly validated reproducing the outcomes of relevant reference case studies. Each time the clapper strikes the bell a new “piece” of the solution is initialized, so avoiding user interventions in the elaboration phase. Independently on the oscillation amplitude and on the duration of the considered time interval, the algorithms can successfully manage undamped oscillations; friction and viscosity damped oscillations; free oscillations in transient and stationary phases; and can be applied also to solve stiff equations. Furthermore, the capability of the proposed methods to deal with arbitrary forcing functions is particularly innovative. The outcomes of relevant case studies, regarding the oscillations of the old tenor bell of the Great St. Mary church in Cambridge, confirm the potentialities of the method, also highlighting some topical issues, involving, for example, the assessment of damping equivalence. Finally, a pioneering feature of the algorithms is their ability to handle and to define “resonant” forcing functions, continuously tuning the frequency of the excitation to the natural frequency of the oscillation, according to the oscillation amplitude.


Introduction
Since their origin, dating in China back to around the 3rd millennium BC [1], bells have been widely used as musical percussion instruments for religious rites and other relevant events.
Starting from the early Middle Ages, bells knew a large diffusion also in Europe and in the Christian world, so that tower bells became a normal addition of churches and other important municipal or civil buildings.
From the musical point of view, according to the Musical Instrument Museums Online (MIMO) revision [2,3] of the Hornbostel-Sachs classification [4] the bells are directly struck idiophone instruments, that vibrate "without requiring stretched membranes or strings". The sophisticated categorization of bells given in [2] is a function of the mounting, of the position and of the way they are struck, but for the aim of the study a more rudimentary distinction is enough.
Depending on the position, bells can be classified as: a. resting bells, whose "cup is placed on the palm of the hand or on a cushion; its mouth faces upwards", not relevant for the aim of the study, or b.
Suspended bells are grouped in three main subcategories, depending on the way they are struck: b1. bells, often without internal clapper, struck from outside with a separate hammer, b2. clapper bells, characterized by an internal clapper, b3. bells with attached external clapper.
Oriental religions generally utilize bells belonging to type b1 or b3; on the contrary, in the Christian context, church bells generally fall in type b2, even if, in particular cases, or events, they can also be struck from outside.
Depending on the excitation, clapper bells can be rung in two different ways. In fact, the vibrations can be excited by clapper strokes indirectly caused by the swinging bell, or by chimes, i.e., leaving the bell quite motionless and directly striking the clapper against it. The latter technique is typical of the Orthodox Church tradition, while in the rest of Christianity swinging bells are the rule, even if in most cases, for conservation purposes, they are chimed via motorized devices.
Swinging clapper bells raise intricate mathematical and physical problems, covering sound and acoustic engineering, mechanical engineering, and structural engineering as well.
As already noted, the swinging bell is a peculiar musical instrument played with a device, the clapper, whose movements, although directly excited, are not explicitly controlled by the bell ringer. The achievement of a clear and optimal sound thus relies not only on the acoustic characteristics of the bell but also on the regularity of the clapper strokes, which depends on the parameters of clapper and bell governing the equations of motion. Such parameters have been empirically tuned over centuries of practical experience leading to standardized proportions, nevertheless, in some cases, the sound control could result extremely difficult even when the bells are driven by well experienced bell ringers. In fact, if on the one hand, in case the bells are manually rung, some odd behavior can be compensated by skilled bell ringers duly modifying the excitation, on the other, in case the bells are swung by motors, fine real time adjustments are practically precluded.
It must also be underlined that different layouts are adopted to suspend the bells, mainly depending on: the ringing technique; the weight of the joke, eventually contributing to counterweight the bell; the position of the bell pivot with respect to the center of gravity of the bell; the distance between the clapper pivot and the bell pivot. The discussion of the rationale and consequences of different arrangements is out of the scope of the present paper, but it must be stressed that distinctive features of bells are also influenced by the technique intended for ringing them.
As described in [5], generally, they can be distinguished: • the "flying clapper" ringing technique, mainly used in continental Europe, except the Iberic peninsula and north Italy, • the full circle techniques, like the English full circle ringing or change ringing technique [6][7][8][9], and the Italian "Bolognese", "Veronese", and "Ambrosiano" styles, • the Spanish system, typically used in Iberic peninsula.
In the "flying clapper" technique, the maximum oscillation of the bell, which is not counterweight, is around 90 • , although in some cases it can attain 135 • , and the clapper strikes the bell on the upper part of the mouth. Occasionally, in France, Belgium, and the Netherlands also the opposite "falling clapper" scheme is adopted, where the clapper strikes the lower part of the mouth.
The full circle techniques are based on alternate full circle oscillations, which facilitate the sound control. The bells are firstly driven in the upward position, where they are briefly stopped with the help of suitable devices or even manually, and then put in oscillations, which are maintained accurately adjusting the excitation. The English full circle and the Bolognese style use unbalanced bells, on the contrary "Veronese" and "Ambrosiano" styles adopt balanced bells.
In the Spanish system, balanced bells continuously rotate full circle in one direction [10][11][12][13]. From the above considerations, it clearly emerges that study of swinging bells can involve: • the assessment of the acoustic characteristics of the bell itself, intended as a musical instrument, • the quality and the precision of the sound, depending on the ease of controlling the clapper impacts during the oscillations, and • the interactions between the swinging bell and the supporting structures: belfry and bell tower.
The acoustics of the bells is a very specialized topic, outside the scope of the present study, therefore it will not be further discussed. The interested reader can find some basic information about theoretical, experimental, or numerical acoustic studies in [14][15][16][17].
Truly, the ease of controlling the impacts significantly influences the acoustic performance of the bell; in fact, as remarked in [18], "the "ringability" of bells depends not only on the linear acoustics of the modes and natural frequencies of the bells but also on the rotational dynamics of the bell and clapper, interacting through impact each time the clapper strikes the bell surface." The inertia forces due to the swinging of the bell play a crucial role in dynamic identification [19] or assessment of existing bell towers [8][9][10][11][12][13][20][21][22][23][24][25], especially when the frequency of the bell oscillations is in the neighborhood of the relevant natural frequencies of the tower. Some analytical solutions for bell-tower interaction are given in [13], considering free oscillations and different mounting systems. However, some studies focus on the seismic response of the bell itself [26].
The time histories of bell and clapper rotations and support reactions are commonly studied in the framework of the Lagrangian mechanics, following the approach proposed by Lagrange [27,28], which is described in detail in every modern book of theoretical mechanics [29][30][31]. The analytical model of the bell, based on rigid body dynamics, originally suggested by Veltmann [32] in the 19th century, was recovered in the late 20th century by Threlfall and Heyman [33], who also suggested a simple experimental method to derive the inertia forces transmitted by full circle ringing bells. Since then, several studies have been carried out, shortly discussed in the following.
A recurrent simplification disregards the presence of the clapper; the bell motion is thus governed by the ordinary differential equation (ODE) of the physical pendulum, depending on two independent initial conditions. But, as better illustrated in Section 2, more refined approaches lead to a system of two second-order nonlinear ODEs. The solutions of that system are the time histories of bell and clapper rotations, provided that four independent initial conditions are given.
It must be underlined that, except in the case of small amplitude oscillations, when the physical pendulum equation can be linearized, the ODEs are always nonlinear, even if springs and dampers are not present. That peculiarity significantly distinguishes the swinging bell model from the simple harmonic oscillator, whose behavior is nonlinear only if nonlinear terms are explicitly introduced, like, for example, nonlinear restoring forces or springs in the Duffing oscillator [34]. Sometimes Duffing oscillator has been also used to approximate the nonlinear pendulum equation, based on the McLaurin expansion of the sinus function, but these aspects are outside the scope of the study.
To solve the swinging clapper bell problem, many different approaches have been proposed, alternatively concentrating the attention on the bell-clapper interactions or on the bell-tower excitation. Frequently, aiming to obtain convenient expressions of the equations, drastic simplifications are introduced, so arriving at suitable analytical or numerical solutions.
Noteworthy simplifications, often combined, may include: • neglecting the clapper, so that the bell reduces to the already cited physical pendulum, suitably modifying the forcing functions to obtain "workable" equations from the analytical [22], or the numerical point of view [35][36][37][38], • adopting simplified models to describe the impact of the clapper on the bell.
However, if analytical solutions are available only for special cases, also usual numerical approaches can result unsatisfactory. In fact, numerical solutions reported in the relevant literature are commonly obtained adopting built-in algorithms for the solutions of ODEs, available in advanced software packages, like Matlab ® . These algorithms are frequently based on Runge-Kutta methods. The availability of such sound and validated algorithms, characterized by clearly defined limitations, is certainly beneficial, but resorting to them involves rather evident disadvantages, because these general methods are scarcely flexible and not easily adaptable to handle specific problems.
Moreover, hypothesizing viscous damping only and assuming special forcing functions, the field of application of these "traditional" approaches sensibly reduces.
With regards to the previous remarks, it must be stressed, on one hand, that other kinds of damping, like friction damping, can play a key role and on the other, that different forcing functions can be envisaged in stationary or transient phases, also depending on the type of excitation−manual or automatic, and on the type of the motorization: linear, or rotating, or even contactless, if any.
Since the pivots are generally equipped with bushes or roller bearings, in which dissipative forces depend on the radial reaction forces, friction damping appears as the most appropriate model. Nevertheless, it is common practice to adopt an "equivalent" viscous damping model, so drastically simplifying the equations. Unfortunately, as confirmed by the outcomes of the study, the equivalence holds on a limited range of oscillation amplitudes, so that the parameters of the equivalent model should be continuously adjusted and adapted, depending on the actual situation.
Regarding the forcing functions, it is noteworthy that a typical transient is the starting phase, in which the bell is set in motion and suitably excited until a pre-established configuration is achieved. Especially in the case that the bells are manually driven or contactless excitations are adopted, the forcing functions are often optimized to be "resonant", i.e., continuously and suitably tuned by an external device or by a trained bell ringer to match the instantaneous natural frequency of the system, which is a function of the amplitude of the oscillation.
For these reasons, an ideal and fully featured method to solve the ODEs of motion should allow: a. to take into account both viscous and friction damping, b.
to consider "arbitrary" forcing function, c.
to determine the solution in both stationary and transient phases, d.
to introduce or to deduce "resonant" excitations, automatically tuned on the instantaneous natural frequency of the system.
In addition, the method should work also with "stiff equations". The aim of the study is the setup of an original all-purposes general procedure to solve the system of ODEs describing the motion of the bell and of the clapper.
The proposed algorithms discussed in Section 6 are based on explicit-implicit predictor-corrector integration methods with constant time step, directly operating on the ODEs system. They do not require any preliminary manipulation of the equations themselves, with evident advantages. To highlight their flexibility, two slightly different versions of the algorithms are presented: • a moderately simplified version, in which the prediction-correction phase works only on the rotations of bell and clapper, • the complete, more refined, version, in which the convergence is improved thanks to an additional prediction-correction phase regarding the angular velocities and the angular accelerations too.
The paper is organized as follows: • Section 2 illustrates the deduction of the ODEs governing the problem. • Section 3 refers to the model of the clapper to bell impact and how it is considered in the equations. • Section 4 illustrates the friction and viscous damping models, while • Section 5, considering the dissipative moments corresponding to damping phenomena, describes the general form of the nonlinear system of ODEs. • Section 6 illustrates in detail the implementation of the proposed algorithms. • Section 7 deals with the validation of the proposed method, checking its capability to reproduce some relevant theoretical or experimental results. To enrich the comparison, the convergence of the proposed algorithms is assessed, and its dependence on the time step and other influencing parameters is critically discussed, also considering "stiff" equations. Of course, that section is not essential for the understanding of the remaining parts of the paper. • Section 8 considers a significant case study, consisting in the old tenor bell of the Great St. Mary's church in Cambridge, whose characteristics are precisely reported in [17]. Several relevant situations are analyzed, considering damped or undamped conditions; free oscillations; forced oscillations; stationary or transient phases; different initial conditions; demonstrating the main features and the potentialities of the proposed algorithms. • Section 9 summarizes the main conclusions, envisaging further developments.
The novel approach and the algorithms proposed in the paper allow one to overcome the limitations affecting the "traditional" approaches proposed until now; in fact, besides satisfying all the previously recalled requisites (a−c), they are also able to derive or to define resonant excitations and to solve stiff equations, as well. Independently on the oscillation amplitude and on the duration of the considered time interval, the algorithms can successfully manage undamped oscillations; friction and viscosity damped oscillations; free oscillations in transient and stationary phases; and can be applied also to solve stiff equations. Furthermore, the capability of the proposed methods to deal with arbitrary forcing functions is particularly innovative.
Typical outcomes of the method are time histories and phase plots in transient and stationary phases of relevant parameters under "arbitrary" forcing functions, including the excitations transmitted to the belfry and the bell tower through the supports. Even in the paper the attention is mainly concentrated on "unbalanced" swinging systems; the proposed methods apply also to the other system previously cited.

The Governing Equations
In the classical Lagrange-Euler approach, the equations of the bell and clapper dynamics can be easily derived starting from the Lagrangian function L of the two degrees of freedom system, where, • θ = θ(t) is the angle between the axis of the bell and the vertical axis, • ϕ = ϕ(t) the relative angle between the clapper axis and the bell axis (Figure 1), θ(t) and ϕ(t) their first derivatives with respect to the time t, i.e., their angular velocities, respectively, • T the kinetic energy of the system, • and U its potential energy.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 6 of 42 • = is the angle between the axis of the bell and the vertical axis, • = the relative angle between the clapper axis and the bell axis ( Figure 1), • = and = their first derivatives with respect to the time , i.e., their angular velocities, respectively, • the kinetic energy of the system, • and its potential energy.
Referring again to Figure 1, let be the pivot of the bell, the pivot of the clapper, and their distance, = − , being and the ordinates of and , respectively. An error affecting the position of the pivot of the clapper is hypothesized. Be such an error, which is the distance ′ between the actual position ′ of the clapper pivot and its theoretical position. Obviously, the positioning error can be expressed in other terms by means of the angle , = atan .
Although the positioning error, typically ranging in the interval = 10 − 20 mm, is significant only for small bells; it is considered here for the sake of generality.
where Referring again to Figure 1, let A be the pivot of the bell, B the pivot of the clapper, and r their distance, r = y B − y A , being y A and y B the ordinates of A and B, respectively.
An error affecting the position of the pivot of the clapper is hypothesized. Be δ such an error, which is the distance BB between the actual position B of the clapper pivot and its theoretical position. Obviously, the positioning error can be expressed in other terms by means of the angle α, Although the positioning error, typically ranging in the interval δ = 10 − 20 mm, is significant only for small bells; it is considered here for the sake of generality. Said: Finally, considering generalized forces, the system of nonlinear ordinary differential equations is obtained, where M e (θ, t) and m e (ϕ, t) are the external torque moments, if any, acting on the bell and on the clapper, respectively, M dis θ, . θ, t the dissipative moment at the bell pivot and m dis ϕ, . ϕ, t the dissipative moments at the clapper pivot, as illustrated in Figure 2.
It is quite obvious to remark that M dis θ, ϕ, t contributes to the motion of the bell. The application of Lagrangian mechanics is rigorously limited to nondissipative systems, even if, under certain conditions, it can be extended also to dissipative systems [39,40]. Nevertheless, since Newtonian mechanics can be applied also in case of dissipative systems, Equation (8) can be considered as dynamic equilibrium equations and the dissipative forces in the right sides as external forces.
Typical dissipative moments are associated to viscosity or friction. In the classical oscillation theory, dissipation is associated to linear viscous dampers, where resisting forces are proportional to velocities through equivalent damping coefficients. In the study, it has been considered also friction damping, characterized by resisting forces which are nonlinear functions of the displacements or of their time derivatives, as better illustrated in the following.
is obtained, where , and , are the external torque moments, if any, acting on the bell and on the clapper, respectively, , , the dissipative moment at the bell pivot and , , the dissipative moments at the clapper pivot, as illustrated in Figure 2. It is quite obvious to remark that , , and , , are oriented to resist the motion of the bell and of the clapper, respectively, while , , resists the motion of the bell only if < 0; on the contrary, when > 0, , , contributes to the motion of the bell. After some manipulations, omitting variables in the moment expressions, Equation (8) becomes, which is a system of two nonlinear second order ordinary differential equations (ODEs), whose solution depends on the given initial conditions. It must be stressed that the solution is represented by piecewise-defined functions, each one continuously differentiable in appropriate subspace domains corresponding to subintervals of the closed interval where the clapper can freely swing ( Figure 3): When positioning errors are considered and odd clapper impact occurs, the interval becomes As illustrated in Figure 3, in Equation (11) ϕ * and ϕ * * are given by: where When the swing angle of the clapper attains a limit value, the clapper impacts the bell, the angular velocities of the clapper and of the bell suddenly change, and a new "piece" of the solution, duly connected with the preceding one, is initialized. Obviously, if two consecutive clapper impacts occur on opposite sides of the bell, the corresponding "piece" of the solution is defined on the closed interval given by Equation (10) or by Equation (11); otherwise, it is defined on a left closed or a right closed subinterval of it, depending on the evolution over time of the relative clapper-bell configurations. The modeling of the impact is discussed in the next section, but it can be anticipate that the use of coefficient of restitution (COR) in terms of relative velocity, like in the classical Newton approach, seems to be the most appropriate for the purposes of the present study.
When positioning errors are considered and odd clapper impact occurs, the interval becomes − * * ≤ ≤ * .

Modelling the Clapper Impact
Over the years, several approaches have been proposed to model the clapper to bell impact, also depending on the aim of the specific study. In classical approaches, the impact is modelled from macroscopic point of view, based on the coefficient of restitution [41]. The coefficient of restitution is a simple and synthetic macroscopic parameter, which summarizes not only the intricate interactions occurring during the impact, but also the influence of the clapper and bell geometry and of the ringing style on the impact itself. The coefficient of restitution can be alternatively defined: • in kinematic terms, as ratio between the relative orthogonal velocities after and before the impact, according to the Newton definition, • in kinetic terms, as ratio between the normal impulses during the restitution phase and the compression phase, according to the Poisson definition, or • in energetic terms, as recently suggested in [42], as ratio between the works done during the restitution phase and the compression phase.
In case of the bell, the COV "measures" also the efficiency of the energy transfer between clapper and bell, in exciting bell vibrations. Manifestly, different definitions of COR are equivalent only if the impact is frictionless; on the contrary, as remarked in [42], if friction between the impacting bodies is relevant, hypothesizing the equivalence leads sometimes to illogical results.
Recently, more sophisticated models have been proposed, particularly focusing on the influence of the impact on the integrity [43] or the acoustic performance of the bell, but a careful examination of these sophisticate models is outside the scope of the present work.
In the present case, the kinematic approach in terms of relative velocities appears particularly comprehensible. The impact phase is studied applying the law of conservation of angular momentum of the whole system with respect to the bell pivot A in the configuration assumed by the system when the impact occurs ( Figure 4). According recurrent assumptions and considering that the impact is quite instantaneous, the impulse of ordinary forces, if present, the energy losses at the pivot, and the displacements of body particles, have been disregarded during the impact phase. In Equation (16) ′ is the distance between the bell pivot and the clapper center of gravity in the impact position ( Figure 4), and + cos * − is the measure of its projection on the clapper axis. Said the moment of inertia of the clapper around in that position, the moment of the inertia around of the whole system, and * the term pertaining to in (16) the law of conservation of the angular momentum for an isolated system gives being and the angular velocities of the bell before and after the impact, respectively. Recalling Equation (14), is easily obtained from Equation (21) as Finally, the variations of the angular velocities during the impact result: where . ϕ i and . ϕ f are the relative angular velocities of the clapper at the time of impact and at the end of the impact, respectively. As known, K = −1 corresponds to an impact against massless bell, K = 0 to the perfectly plastic case, and K = 1 to an ideal perfectly elastic impact against a fixed surface, in the latter case, the kinetic energy of the system is obviously conserved.
The kinematic COR can be easily measured in practice by means of high-speed video cameras, using super slow-motion video recording capabilities of such kind of devices.
As already said, the clapper strikes the bell each time it holds one of the conditions holds. For example, referring to the case ϕ = ϕ * , the angular momentum L of the system around A is: In Equation (16) b is the distance AG c between the bell pivot and the clapper center of gravity in the impact position ( Figure 4), and b + r cos(ϕ * − α) is the measure of its projection on the clapper axis. Said I cA the moment of inertia of the clapper around A in that position, I bt the moment of the inertia around A of the whole system, and I c * the term pertaining to . ϕ in (16) the law of conservation of the angular momentum for an isolated system gives being . θ i and . θ f the angular velocities of the bell before and after the impact, respectively.
Recalling Equation (14), . θ f is easily obtained from Equation (21) as Finally, the variations of the angular velocities during the impact result: Equations (14) and (22) allow to consider the impact effects in a very effective way, since they constitute, together with the pertinent condition on ϕ (Equation (15)) and the corresponding value of θ f at the end of the impact, the initial conditions for the subsequent "piece" of the piecewise-defined solution. From another point of view, these variations can be simulated setting the angular accelerations, .. ϕ j , at the j − th step in the form in which . ϕ j is the relative angular velocity of the impacting clapper and ∆t the time step adopted in numerical integration. Similar conclusions can be derived from the impulse-momentum theorem.
As said before, refined evaluation of the coefficient of restitution would require much more sophisticated modelling of the impact phase, nevertheless the simplified approach proposed above can be easily improved and corrected; on one hand, to avoid physical paradoxes, on the other, to properly consider additional phenomena associated with the impact.
First, it must be remarked that the bounds on K given in Equation (14) can be finer tuned, since, on the light of what previously discussed, in a swinging bell K = 1 infringes the law of conservation of energy. For physical coherence, the effective relative angular velocity between the clapper and the bell after the impact, ∆ω, involving not only that, recalling Equations (14) and (22), finally provides in which the upper bound corresponds to a perfectly elastic impact. It must be stressed that, while in the first clapper stroke(s) a considerable part of energy is spent to induce vibrations in the bell, in the subsequent ones the transfer is less relevant, being only responsible for maintaining the vibrations themselves. The first impacts are thus characterized by lower coefficients of restitution. Furthermore, in some cases bell ringers cover a side of the clapper with felt pads to muffle the sound as a function of the side of impact. In any case, these aspects can be simply managed, varying the COR as a function of the strike sequence or of the impact side.
At this point, it could be argued that the bell is not an isolated system, because part of energy, associated with the acoustic waves, leaves the system. This facet can be tackled suitably modifying the second Equation (23), and expressing the variation of the bell angular velocity ∆ . θ as in which ξ s is the quota of the angular momentum, ∆L cb , transferred by the clapper to the bell during the impact, which is spent to excite bell vibrations.

Modelling the Dissipation
In the analysis of simple oscillators, dissipation is usually modelled by viscous dampers, so that dissipative forces and moments are assumed proportional to the velocity .
x or to the angular velocity .

θ.
Consequently, dissipative effects are described in translational mechanics by the force and in rotational mechanics by the moment being c the damping coefficient, and C the angular damping coefficient. Although very simple and effective, these models are not satisfactory when more complex dissipation phenomena are involved, like, for example, in rotational mechanics. For that reason, in the present study, a more convincing friction dissipative model is also introduced. In that model, friction forces at the pivot depend on the radial reactive forces at the pivot itself.
Referring to a pivot of radius r p , the dissipative moment due to friction forces can be written as being F r the radial force and µ f the friction coefficient. In the present case, F r are the forces resulting from inertia forces and self-weight at the bell pivot and at the clapper pivot. Adopting the previously mentioned notation (Figure 1), after some trivial mathematical passages, the centrifugal reactive force at the clapper pivot, B, results Furthermore, the centrifugal reactive force at the bell pivot, A, results Appl. Sci. 2020, 10, 5528 13 of 44 Finally, also considering the self-weights, Equations (31) and (32) become

Dissipative Moments at Clapper Pivot and at the Bell Pivot
As already noted in Section 2, dissipative moments at the pivots resist the motion of the element directly connected to the pivot, but, in some cases, depending on their orientation, they contribute to the motion of other elements. Let C c and C b be the angular damping coefficient at the clapper and bell pivot, respectively; considering the already recalled viscous model (Equation (29)), the resisting moments at the clapper pivot B , M vB , and at the bell pivot A, M vA , result On the contrary, considering friction according to Equation (30), the friction moments are then at the clapper pivot, B and at the bell pivot, A, being r c the radius of the clapper hinge, r b the radius of the bell hinge, and µ f c and µ f b the respective friction coefficients. In Equation (37) → F rB is the modulus of the radial force acting in B (Equation (33)); in Equation (38) → F rA is the modulus of the radial force acting in A (Equation (34)).
Finally, in the most general case, considering both dissipative sources, the governing system of ODEs can be written as

Numerical Integration of the Nonlinear ODEs System
Considering the peculiar structure of the nonlinear system of ODEs (Equation (39)), a very effective numerical explicit-implicit predictor-corrector integration method with constant time step can be envisaged. Since numerical methods to solve ODEs are the subject of huge specialized literature, they are no further discussed; the interested reader can refer, for example, to [44].
In the proposed method, at a given time step, the angular velocities are calculated first, using the forward Euler method, starting from the angular accelerations at the previous step; subsequently, from the angular velocities at the current and at the previous step, the rotations are obtained, via the Heun's explicit trapezoidal rule; finally, the angular accelerations at the current step are predicted solving the resulting system of two linear algebraic equations. The solution is then corrected: the angular velocities and the rotations are recalculated using the Heun's explicit trapezoidal rule, while the angular accelerations are recalculated by solving again the linear system of algebraic equations.
The main steps of the practical implementation of the method are: 1.
Define the relevant parameters of the problem: masses, moments of inertia, geometry of bell and clappers, limit values of clapper oscillation, −ϕ * * and ϕ * , angular damping coefficients, C c and C b , friction coefficients, µ f c and µ f b , coefficient of restitutions K m .

2.
Define typology and expressions of the external excitations, M e and m e , if any.

4.
Define the time interval [0, t max ] to be investigated.
Considering M e0 and m e0 , substitute the initial values in Equations (39)

7.
Start the iteration at the considered time step i, t i = t 0 + i ∆t.

8.
Calculate the predicted values of the angular velocities . θ ip and . ϕ ip at the time t i by means of the Euler forward formula: 9.
Calculate the predicted values of the rotations θ ip and ϕ ip at the time t i by means of the Heun's explicit trapezoidal formula: 12. Calculate the corrected values of θ i and ϕ i at the time t i using the Heun's explicit formula: 13. Check the value of ϕ i : if −ϕ * * < ϕ i < ϕ * go to step 16. 14. If ϕ i ≤ −ϕ * * , set ϕ = −ϕ * * , else, if ϕ i ≥ ϕ * , set ϕ = ϕ * .
• Register the clapper strike number, n, and the step number, j, • Take note of the time, t nj = t 0 + j ∆t, • Select the pertinent coefficient of restitution, K n , depending on n and on the side of impact, Initialize the new clapper oscillation, t 0 = t nj , and assume as initial values (Equation (22)): • Considering M e0 and m e0 , substitute the initial values in Equations (39), and evaluate the unknowns ..  17. If t i = t max then end the elaboration. 18. Go to step 7.
Obviously, to optimize the computer memory usage, the procedure can be repeated as many times as necessary, duly modifying the initial time and starting from the previously obtained initial conditions. The proposed procedure can also be slightly simplified, eliminating the correction steps 10, 11 and 12, so jumping indeed directly from step 9 to step 13.
The validation of the method is discussed in the next section, referring to some relevant cases.

Validation of the Methods
For the sake of validation, the methods have been applied to reproduce significant experimental and theoretical results with increasing degree of complexity, also discussing their convergence.
The considered cases are: 1.
Frequencies of small free oscillations of bell and clapper, considering the old tenor bell of the Great St. Mary's, Cambridge, and a laboratory bell, indicated in the following as bell 1, both described in detail in [18], whose parameters are summarized in Table 1. In the table L b and c are the equivalent lengths of the bell and of the clapper, and T a and T b the periods of the free small oscillations of the bell and of the clapper, respectively: 2.
Bell oscillations and clapper strikes in one cycle for the bell 1, for which also convergence is discussed, and for other two bells (bell 2 and bell 3) studied in [18] as well (see Table 1), 3.
Damped oscillations of the of the Great St. Mary's bell.
The above mentioned cases have been selected in [18] suitably modifying the most relevant parameters characterizing the bell behavior, namely the ratios r/L c and L b /L c . By means of the algorithm proposed before, the small oscillations frequencies of clapper and bell of the of Great St. Mary's tenor bell and of the laboratory bell 1 have been numerically evaluated assuming, for the clapper, the initial conditions: and, for the bell, the initial conditions: The diagrams, obtained considering a massless clapper in case of bell excitation, and r = 0 in case of clapper excitation, are reported in Figure 5a,b for bell and clapper excitation of St. Mary's tenor bell and in Figure 6a,b for bell and clapper excitation of bell 1. The natural periods so obtained fit very well the theoretical ones, as evident from the comparison (see Table 2).

Bell Oscillations and Clapper Strikes of Three Reference Bells
As recalled before, the clapper strikes on one cycle of a free oscillating bell have been extensively studied in [18], where results are diagrammatically shown for the three reference bells.
Assuming the properties already recalled in Table 1, the cycles, starting at an angle of 171°, have been reproduced with the proposed algorithm, considering two cases, according as the bell is starting ringing right or wrong. The bell is said to be ringing right if the clapper impacts are "in phase" with the bell oscillations, i.e., the clapper impacts the bell on the upper part of the mouth, or it is said ringing wrong if the clapper impacts are "in opposition of phase" with the bell oscillations, i.e., the clapper impacts the bell on the lower part of the mouth. Analytically, the bell starts ringing right if the first clapper impact satisfies the condition: it starts ringing wrong, if the first clapper impact satisfies instead. In the considered cases, the corresponding initial conditions are thus

Bell Oscillations and Clapper Strikes of Three Reference Bells
As recalled before, the clapper strikes on one cycle of a free oscillating bell have been extensively studied in [18], where results are diagrammatically shown for the three reference bells.
Assuming the properties already recalled in Table 1, the cycles, starting at an angle of 171 • , have been reproduced with the proposed algorithm, considering two cases, according as the bell is starting ringing right or wrong. The bell is said to be ringing right if the clapper impacts are "in phase" with the bell oscillations, i.e., the clapper impacts the bell on the upper part of the mouth, or it is said ringing wrong if the clapper impacts are "in opposition of phase" with the bell oscillations, i.e., the clapper impacts the bell on the lower part of the mouth. Analytically, the bell starts ringing right if the first clapper impact satisfies the condition: it starts ringing wrong, if the first clapper impact satisfies instead. In the considered cases, the corresponding initial conditions are thus for the bell ringing right and for the bell ringing wrong. According to the assumption made in [18], in both cases the kinetic coefficient of restitution has been hypothesized K = √ 0.2, independently on the impact sequence.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 18 of 42 and = obtained in the present study are reported in Figure 8c,e, regarding bell 2, characterized by = 4.137 s and in Figure 8d,f, regarding bell 3, characterized by = 4.235 s.  It must be considered that in [18] specific information about the period of the cycles as well as on the magnitude of the angular velocity of the clapper is not given, therefore they have been considered in relative sense. In Figure 7c,d they are reported the θ = θ(t) and ϕ = ϕ(t) diagrams obtained with the proposed algorithm in the two considered cases, while the corresponding . ϕ = . ϕ(t) clapper angular velocity diagrams are shown in Figure 7e,f.
Adopting a time step of ∆t = 0.025 ms, the calculated periods are T = 2.883 s for the bell ringing right and T = 2.706 s for the bell ringing wrong.
To test the convergence, the case of bell 1 ringing right has been analyzed with the complete version of the algorithm, as well as with the simplified version, where the second correction phase, corresponding to steps 10, 11, and 12 in Section 6, is lacking. In the tests, time steps varying in the range ∆t = 0.025 − 5 ms have been adopted. The results in terms of minimum and maximum angular velocities of bell and clapper are summarized in Table 3, while in Table 4 are reported the analogous results obtained with the simplified version. (1) the solution becomes numerically unstable (θ(t > 3 s) > θ 0 ).
The analysis of the results demonstrates that the complete version converges satisfactorily even adopting a rather coarse time step ∆t = 5 ms and that using ∆t ≤ 1 ms the errors are insignificant. Using the simplified method, satisfactory results are obtained, as expected, only adopting time steps ∆t ≤ 1 ms, while errors become negligible only when ∆t ≤ 0.1 ms. In addition, adopting a time step ∆t = 5 ms, the simplified method leads to a numerically unstable solution, since the bell exhibits a full rotation, in fact θ(t > 3 s) > θ 0 .
ϕ(t) obtained in the present study are reported in Figure 8c,e, regarding bell 2, characterized by T = 4.137 s and in Figure 8d Finally, to complete that analysis, the case of bell ringing wrong has been studied for bell 2, characterized by = 4.342 s, and for bell 3, characterized by = 4.514 s. Results in terms of =  Finally, to complete that analysis, the case of bell ringing wrong has been studied for bell 2, characterized by T = 4.342 s, and for bell 3, characterized by T = 4.514 s. Results in terms of θ = θ(t) and ϕ = ϕ(t) are reported in Figure 9a for bell 2 and in Figure 9b for bell 3 respectively, while the diagrams of the angular velocities of the clapper are illustrated in Figure 9c,d, respectively.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 20 of 42 and = are reported in Figure 9a for bell 2 and in Figure 9b for bell 3 respectively, while the diagrams of the angular velocities of the clapper are illustrated in Figure 9c,d, respectively.

Damped Oscillations of St. Mary Old Tenor Bell
As previously remarked, one of most important features to be assessed in studying numerical methods for solution of ODEs is their ability to manage "stiff" equations. Although there is no precise definition of equation stiffness, it is commonly agreed that differential equations manifest stiffness if their integration requires unnecessarily short time steps.
Since equation stiffness increases raising the damping coefficient, to assess the capability of the proposed methods to manage this kind of equations, they have been applied to St. Mary tenor bell to evaluate the value of the critical damping coefficient; to determine the structural response in critical or overdamped conditions as well as to estimate the relevant conditions for convergence. In the analysis, time steps have been varied in the range [0.025 − 10 ms].
As known, said its natural angular frequency, the critical value, , of the angular damping coefficient of a pendulum characterized by mass , equivalent length , and moment of inertia , is ϕ(t) curve obtained with the proposed algorithm for bell 3.

Damped Oscillations of St. Mary Old Tenor Bell
As previously remarked, one of most important features to be assessed in studying numerical methods for solution of ODEs is their ability to manage "stiff" equations. Although there is no precise definition of equation stiffness, it is commonly agreed that differential equations manifest stiffness if their integration requires unnecessarily short time steps.
Since equation stiffness increases raising the damping coefficient, to assess the capability of the proposed methods to manage this kind of equations, they have been applied to St. Mary tenor bell to evaluate the value of the critical damping coefficient; to determine the structural response in critical or overdamped conditions as well as to estimate the relevant conditions for convergence. In the analysis, time steps have been varied in the range [0.025 − 10 ms].
As known, said ω n its natural angular frequency, the critical value, C cr , of the angular damping coefficient of a pendulum characterized by mass M, equivalent length L, and moment of inertia I, is Being that a G is the distance between the center of gravity of the pendulum and the pivot. From practical point of view, the damping coefficient can be more conveniently represented by a relative damping factor, c, whose critical value, depends only on L, or by the dimensionless damping ratio ζ which is the ratio between the actual damping coefficient and the critical one. As known, underdamped oscillators, characterized by ζ < 1, tend to the equilibrium position by damped oscillations, overdamped systems, for which ζ > 1, smoothly return to the equilibrium position never overpassing it, critically damped systems, ζ = 1, require the minimum time to regain, without oscillations, a pre-established position.
Applying Equations (54) and (55) to the bell and the clapper of Sty. Mary's tenor bell (Table 1), it results C cr,b = 7069.44 kg m 2 s −1 , for the bell and C cr,c = 78.44 kg m 2 s −1 , c cr,c = 5.503 m s −1 , for the clapper. The critical damping obtained for the St. Mary's bell and clapper small oscillations with the proposed methods are reported in Table 5. In the elaborations, damping has been considered in turn on the bell and on the clapper. In addition, the influence of the clapper on the bell behavior has been neutralized assuming the clapper as massless and c c = 0; while the influence of the bell on the clapper behavior has been minimized suitably increasing c b : The first condition in Equation (60) is dictated by the convergence condition of both methods, which result in the adoption of a constant time step ∆t satisfying the inequality ∆t c b ≤ 2.6 m. (61) In the present exercise, initial conditions have been set to: for small amplitude bell oscillations and to θ(0) = 0; . .
for small amplitude clapper oscillations.
In Figure 10a the time−rotation curve of the critically damped bell in semilogarithmic scale is represented, while in Figure 10b the time−rotation curve of the critically damped clapper is represented. The first condition in Equation (60) is dictated by the convergence condition of both methods, which result in the adoption of a constant time step Δ satisfying the inequality Δ ̅ ≤ 2.6 m. (61) In the present exercise, initial conditions have been set to: for small amplitude bell oscillations and to 0 = 0; 0 = 0; 0 = 2°; 0 = 0 , for small amplitude clapper oscillations.
In Figure 10a the time−rotation curve of the critically damped bell in semilogarithmic scale is represented, while in Figure 10b the time−rotation curve of the critically damped clapper is represented. Results clearly demonstrate that, if time step complies with Equation (61), the proposed methods can manage stiff equations too.
To further stress that conclusion, also two overdamped cases have been studied, for the initial values 0 = 171°; 0 = 0; 0 = 32°; 0 = 0 , for bell and clapper, respectively; corresponding to damping ratios: Results clearly demonstrate that, if time step complies with Equation (61), the proposed methods can manage stiff equations too.
Damping coefficients have been set to: for bell and clapper, respectively; corresponding to damping ratios: To enrich the study, the solutions obtained with a rather coarse time step, ∆t 1 = 10 ms, nearly corresponding to the limit condition for convergence, are compared with the reference solution, obtained with a significantly smaller time step, ∆t 2 = 0.025 ms, in Figure 11a,b, which refers to initial values given by Equations (64)  To enrich the study, the solutions obtained with a rather coarse time step, Δ = 10 ms, nearly corresponding to the limit condition for convergence, are compared with the reference solution, obtained with a significantly smaller time step, Δ = 0.025 ms, in Figure 11a,b, which refers to initial values given by Equations (64)  In Figure 11, the dashed lines, referring to Δ , practically coincide with the solid lines, referring to Δ , thus confirming the soundness and the reliability of the algorithms in handling stiff equations.

Some Relevant Case Studies Referring to St. Mary's Tenor Bell
The nonlinear ODE system and the proposed numerical algorithms have been applied to solve some relevant problems, referring to the St. Mary's old tenor bell. In these problems, the coefficient of restitution has been always assumed = √0.2, unless otherwise specified in the following.
More precisely, the considered cases studies consist in: 1. Large undamped free oscillations of the bell considering four significant cases, characterized by different sets of initial conditions: • Case I.a: initial conditions as per Equation (63)

Large Undamped Free Oscillations
As already said, large undamped free oscillations have been considered in two relevant cases, characterized by different initial angles : = 171° in case I or = 40° in case II. In Figure 11, the dashed lines, referring to ∆t 1 , practically coincide with the solid lines, referring to ∆t 2 , thus confirming the soundness and the reliability of the algorithms in handling stiff equations.

Some Relevant Case Studies Referring to St. Mary's Tenor Bell
The nonlinear ODE system and the proposed numerical algorithms have been applied to solve some relevant problems, referring to the St. Mary's old tenor bell. In these problems, the coefficient of restitution has been always assumed K = √ 0.2, unless otherwise specified in the following. More precisely, the considered cases studies consist in: 1.

2.
Assuming nearly equivalent logarithmic decays for small oscillations, comparison of curves resulting from viscous and friction dissipative models, considering initial conditions as per Case I.

3.
Considering initial conditions as per Case I.a, assessment of time histories of horizontal and vertical excitations on the pivots in three subcases: undamped oscillations; friction damped or "equivalent" viscous damped oscillations 4.
Synchronized resonant forced oscillations of the bell in transient phase.

Large Undamped Free Oscillations
As already said, large undamped free oscillations have been considered in two relevant cases, characterized by different initial angles θ 0 : θ 0 = 171 • in case I or θ 0 = 40 • in case II. It must be stressed that, despite that damping is nil, energy is not conserved, since in both cases, as anticipated, the clapper impact is not perfectly elastic K = √ 0.2 . The analysis of the results shows that, in the given conditions, the bell tends to maintain the initial ringing condition. More precisely, if the bell starts ringing wrong, it continues to ring wrong. Nevertheless, as clearly results from Figure 12b, while in a first phase the bell ringing wrong exhibits two well distinct impacts per semicycle, over time the impacts become more regular, and only one instantaneous and clear contact occurs per semicycle. Moreover, the number of significant impacts (Figure 12a In addition, said θ * i the amplitude and T * i the "apparent period" of the i − th cycle, defined as the time interval between θ * i and θ * i+1 , apparent logarithmic decrements of amplitude, δ * θ , and period, δ * T , can be defined. In the following, these apparent logarithmic decrements are evaluated considering the first five apparent periods, according to: The outcomes of the analysis are summarized in Table 6, which confirms that case I.a is characterized by smaller apparent logarithmic decrements. Although marginal in the present paper, another important question regards the aptitude of the bell to maintain the initial ringing condition, especially when initial condition is wrong. As already remarked, the adoption of a coefficient of restitution K = √ 0.2 implies that in case I.b the bell tends to ring wrong, but this is not always the case. Truly, hypothesizing higher values of the coefficient of restitution the result can be the opposite. For example, assuming K = 0.8, the curves θ = θ(t) and ϕ = ϕ(t) resulting in cases I.a and I.b become those illustrated in Figure 13a,b, respectively. Curves in Figure 13a confirm that, like in Figure 12a, the bell lasts in ringing right state over time; on the contrary, curves in Figure 13b show that, after some cycles, three in that case, the ringing state switches from wrong to right. It must be stressed that, despite that damping is nil, energy is not conserved, since in both cases, as anticipated, the clapper impact is not perfectly elastic = √0.2 . The tendency to stay in "wrong" condition thus depends not only on the relevant ratios r/L c and L b /L c but also on the COR. Moreover, as shown in the following, that tendency depends not only on the quotas of energy dissipated or transferred outside the system during the impact but also on the extent and on the nature of other dissipative phenomena.
to ring wrong, but this is not always the case.
Truly, hypothesizing higher values of the coefficient of restitution the result can be the opposite. For example, assuming = 0.8, the curves = and = resulting in cases I.a and I.b become those illustrated in Figure 13a,b, respectively. Curves in Figure 13a confirm that, like in Figure 12a, the bell lasts in ringing right state over time; on the contrary, curves in Figure 13b show that, after some cycles, three in that case, the ringing state switches from wrong to right.  Results regarding case II, undamped oscillations of the bell starting from θ 0 = 40 • , are reported in Figure 14a,b in terms of θ = θ(t) and ϕ = ϕ(t) curves, in Figure 14c In the figures, impacts obviously correspond to sudden variations, whose extents, recalling Equation (23) or (24), are inverse function of the adopted time step.
These remarks are further corroborated by Table 7, evidencing that in case II.a part of energy is lost, while in case II.b, lacking any kind of dissipation, the energy is conserved. The latter observation confirms that the proposed methods are symplectic integrators and that the energy drift is practically nil.
Recalling Equations (33) and (34), the radial accelerations at the clapper pivot B and at the bell pivot A, including gravity acceleration, can be finally obtained as (70) on the extent and on the nature of other dissipative phenomena.

Case II Undamped Free Oscillations θ = 40°
Results regarding case II, undamped oscillations of the bell starting from = 40°, are reported in Figure 14a,b in terms of = and = curves, in Figure 14c,d in terms of − curves, and in Figure 14e,f in terms of − curves, in analogy with the criterion adopted for case I.  Referring again to case II.b, the resulting time histories of the dimensionless horizontal and vertical component of the radial accelerations in pivots A and B are illustrated in Figure 16. More precisely, Figure 16a,b refers to horizontal components a rA,x /g and a rB ,x /g, while Figure 16c,d refers to vertical components a rA,y /g and a rB ,y /g. occur about every two oscillations. That results even more clearly from Figure 15a,b, which illustrates the time−angular acceleration curves = and = obtained in case II.b adopting ∆ = 0.1 ms. In the figures, impacts obviously correspond to sudden variations, whose extents, recalling Equation (23) or (24), are inverse function of the adopted time step.
These remarks are further corroborated by Table 7, evidencing that in case II.a part of energy is lost, while in case II.b, lacking any kind of dissipation, the energy is conserved.  The latter observation confirms that the proposed methods are symplectic integrators and that the energy drift is practically nil.
Recalling Equations (33) and (34), the radial accelerations at the clapper pivot ′ and at the bell pivot , including gravity acceleration, can be finally obtained as Referring again to case II.b, the resulting time histories of the dimensionless horizontal and vertical component of the radial accelerations in pivots and ′ are illustrated in Figure 16. More precisely, Figure 16a,b refers to horizontal components , / and , / , while Figure 16c,d refers to vertical components , / and , / . As expected, that amplitude does not induce significant excitation on the supports.

Viscous and Friction Damped Free Oscillations
Aiming to compare "equivalent" damping phenomena, case I.a has been studied considering viscous or friction damping as well. Preliminarily, the equivalence between viscous and friction damping has been assessed.
As known, the period, , of viscous damped small free oscillations is constant over time; it depends only on the damping ratio and on the period, , of undamped oscillations, being independent on the exponentially decreasing amplitude of oscillations On the contrary, the apparent period of friction damped oscillations rigorously depends, besides on the friction coefficient, also on the amplitude of the oscillations; consequently, it varies over time, although slightly and the system generally does not retrieve the rest position. Given that the decay rate of friction damped oscillations is not exponential, their logarithmic decrement is variable, therefore the parameters of the equivalent viscous damping depend on the equivalence criterion.
Thanks to its simplicity, widely used procedures to determine the damping coefficient are based As expected, that amplitude does not induce significant excitation on the supports.

Viscous and Friction Damped Free Oscillations
Aiming to compare "equivalent" damping phenomena, case I.a has been studied considering viscous or friction damping as well. Preliminarily, the equivalence between viscous and friction damping has been assessed.
As known, the period, T d , of viscous damped small free oscillations is constant over time; it depends only on the damping ratio ζ and on the period, T, of undamped oscillations, being independent on the exponentially decreasing amplitude of oscillations On the contrary, the apparent period of friction damped oscillations rigorously depends, besides on the friction coefficient, also on the amplitude of the oscillations; consequently, it varies over time, although slightly and the system generally does not retrieve the rest position. Given that the decay rate of friction damped oscillations is not exponential, their logarithmic decrement is variable, therefore the parameters of the equivalent viscous damping depend on the equivalence criterion.
Thanks to its simplicity, widely used procedures to determine the damping coefficient are based on the measurement of the period of damped oscillations or on the evaluation of the logarithmic decrement δ * .
The damping ratio ζ can be thus calculated using Equation (71) or by the formula: Following that procedure, the equivalence between friction and viscous damping has been established equalizing the logarithmic decrement of the viscous case and the apparent logarithmic decrement δ * eq of the friction case determined over five oscillations, using the first Equation (68), accordingly deducing damping ratios and damping coefficients. The amplitude of the first oscillation has been assumed, in turn, θ 0 = 5 • or ϕ 0 = 5 • .
For the St. Mary's tenor bell, assuming µ b = 0.02 and r b = 50 mm, at the bell pivot, and µ c = 0.04 and r c = 10 mm, at the clapper pivot, the equivalent damping ratios and the relative damping coefficient indicated in the last two columns of Table 8 have been derived. As already said, the equivalence of the clapper has been derived assuming r = 0. The corresponding "equivalent" time histories of bell and clapper rotation are illustrated in Figure 17a,b, respectively: dashed lines refer to friction cases, and continuous lines to viscous cases.
Looking at the figures, it is evident that, although first oscillations of friction and viscous curves are nearly coincident, over time friction decay tends to prevail.
Comparison of Large Damped Free Oscillations in Case I (θ 0 = 171 • ) Cases I.a and I.b have been studied also in damped conditions, adopting the damping equivalent parameters previously determined (Table 8), assuming the initial values given by Equations (64) and (65).
The results are compared in Figure 18 in terms of time histories of rotations, in Figure 19 in terms of . θ − θ curves, and in Figure 20 in terms . ϕ − ϕ curves. Figure 17a,b, respectively: dashed lines refer to friction cases, and continuous lines to viscous cases.  Looking at the figures, it is evident that, although first oscillations of friction and viscous curves are nearly coincident, over time friction decay tends to prevail.

Comparison of Large Damped Free Oscillations in Case I θ = 171°
Cases I.a and I.b have been studied also in damped conditions, adopting the damping equivalent parameters previously determined (Table 8), assuming the initial values given by Equations (64) and (65).
The results are compared in Figure 18 in terms of time histories of rotations, in Figure 19 in terms of − curves, and in Figure 20 in terms − curves. Figure 18a,b shows the time histories of bell and clapper rotations derived in case I.a-bell ringing right-considering friction damping and viscous damping, respectively; Figure 18c,d refer to the corresponding situations in case I.b−bell ringing wrong. Figure 19a,b compares − curves in case I.a e I.b, respectively, while Figure 19c,d compares the corresponding curves − . In these diagrams, curves pertaining to friction damping are in blue and curves pertaining to "5° equivalent" viscous damping are in red. The diagrams demonstrate that equivalent viscous damping coefficients determined considering small oscillations cannot be adopted for large oscillations. In fact, the comparison of diagrams in Figures 18 and 19 is a crystal-clear evidence that amplitude and period decrease more quickly for "5° equivalent" viscous damped oscillations than for friction damped oscillations.
The reason for the latter remark is that, as soon as the initial amplitude of the oscillations increases, viscous dissipation increases much more rapidly than friction dissipation. The observation can be further corroborated comparing the energy dissipation caused by friction or by viscosity in a group of cycles or in a given time interval.
Since energy dissipation is nothing else than the work done by the corresponding dissipative force, its calculation is not particularly difficult.
The energy dissipated by friction damping and the energy dissipated by "5° equivalent" viscous damping are compared in the diagrams of Figures 20 and 21a.
In Figure 20 the dissipated energies are plotted as a function of the bell rotation, , cycle by cycle: Figure 20a refers to the small oscillation case considered to assess the equivalence, Figure  20b to the case I.a.     ϕ − ϕ. In these diagrams, curves pertaining to friction damping are in blue and curves pertaining to "5 • equivalent" viscous damping are in red.
The diagrams demonstrate that equivalent viscous damping coefficients determined considering small oscillations cannot be adopted for large oscillations. In fact, the comparison of diagrams in Figures 18 and 19 is a crystal-clear evidence that amplitude and period decrease more quickly for "5 • equivalent" viscous damped oscillations than for friction damped oscillations.
The reason for the latter remark is that, as soon as the initial amplitude of the oscillations increases, viscous dissipation increases much more rapidly than friction dissipation. The observation can be further corroborated comparing the energy dissipation caused by friction or by viscosity in a group of cycles or in a given time interval.
Since energy dissipation is nothing else than the work done by the corresponding dissipative force, its calculation is not particularly difficult.
The energy dissipated by friction damping and the energy dissipated by "5 • equivalent" viscous damping are compared in the diagrams of Figures 20 and 21a. that the equivalence parameters between friction and viscous dissipation is sensibly dependent on the initial amplitude. This aspect is not surprising since the energy dissipation due to friction depends, apart from the oscillation amplitude, on the increase of the friction forces due to the increase of radial forces, while the energy dissipation due to viscosity depends, besides the oscillation amplitude, on the angular velocity.
Owing the fact that the mean modulus of the radial acceleration at the pivot , which is | | ≈ for = 5° oscillations, becomes | | ≈ 2.2 for = 171° oscillations, and that the mean angular velocity, which is ̅ ≈ 0.15 s for = 5° oscillations, becomes ̅ ≈ 2.5 s for = 171° oscillations, the energy dissipations in = 171° oscillations can be roughly estimated, starting from the dissipation evaluated considering = 5° oscillations. As the mean value of the dissipated energy for = 5° oscillations is around 5° ≈ 4 J/cycle (Figure 20a), the dissipation caused by friction in = 171° oscillations, which is approximately , 171° ≈ 42.234 ≈ 300 J/cycle , results sensibly smaller than the energy dissipated by "5° equivalent" viscosity, which is around , 171° ≈ 41634 ≈ 2200 J/cycle. It must be remarked that these values, although roughly estimated, satisfactorily fit the calculated values for the first cycles (Figure 21a).
To conclude the analysis, the "171° equivalent" viscous damping parameters have also been determined, reported in Table 9. These values, which are one order of magnitude lower than the ones determined for small oscillations (Table 8), lead to the time−dissipated energy curve illustrated in Figure 21b (solid line), substantially coinciding with the dashed curve regarding the friction case.   In Figure 20 the dissipated energies E dis are plotted as a function of the bell rotation, θ, cycle by cycle: Figure 20a refers to the small oscillation case considered to assess the equivalence, Figure 20b to the case I.a.
To facilitate the reading, curves in Figure 20b are parameterized considering different time intervals: the dashed magenta curve is the energy dissipated by friction in the time interval 0 − 25 s, while the solid blue line is the energy dissipated by "5 • equivalent" viscosity in the time interval 0 − 6 s. Finally, the E dis − t curves obtained for case I.a are compared in Figure 21a, in which solid line refers to "5 • equivalent" viscosity and dashed line to friction. Figure 20a clearly shows that, as expected, for 5 • oscillations the energy dissipations due to friction or equivalent viscosity practically coincide, but Figures 20b and 21a underline that using "5 • equivalent" viscosity leads to huge overestimations of the dissipated energy. That confirms that the equivalence parameters between friction and viscous dissipation is sensibly dependent on the initial amplitude. This aspect is not surprising since the energy dissipation due to friction depends, apart from the oscillation amplitude, on the increase of the friction forces due to the increase of radial forces, while the energy dissipation due to viscosity depends, besides the oscillation amplitude, on the angular velocity.
Owing the fact that the mean modulus of the radial acceleration at the pivot A, which is |a rA | ≈ g for θ o = 5 • oscillations, becomes |a rA | ≈ 2.2 g for θ o = 171 • oscillations, and that the mean angular velocity, which is . θ ≈ 0.15 s −1 for θ o = 5 • oscillations, becomes . θ ≈ 2.5 s −1 for θ o = 171 • oscillations, the energy dissipations in θ o = 171 • oscillations can be roughly estimated, starting from the dissipation evaluated considering θ o = 5 • oscillations. As the mean value of the dissipated energy for θ o = 5 • oscillations is around E dis (5 • ) ≈ 4 J/cycle (Figure 20a), the dissipation caused by friction in θ o = 171 • oscillations, which is approximately E dis, f (171 • ) ≈ 42.234 ≈ 300 J/cycle, results sensibly smaller than the energy dissipated by "5 • equivalent" viscosity, which is around E dis,v (171 • ) ≈ 41634 ≈ 2200 J/cycle. It must be remarked that these values, although roughly estimated, satisfactorily fit the calculated values for the first cycles (Figure 21a).
To conclude the analysis, the "171 • equivalent" viscous damping parameters have also been determined, reported in Table 9. These values, which are one order of magnitude lower than the ones determined for small oscillations (Table 8), lead to the time−dissipated energy curve illustrated in Figure 21b (solid line), substantially coinciding with the dashed curve regarding the friction case. In addition, the entity of the dissipative forces can play an important role in discriminating the aptitude of the bell to ring right or wrong, as demonstrated by the St. Mary's tenor bell behavior when large free oscillations are considered. In fact: • on the one hand, as long as dissipative forces are low, over the time the bell tends to maintain the initial ringing state (right or wrong) (see Figure 18a,c), • on the other, if dissipative forces increase, the bell tends to switch, independently on the initial state, into a ringing wrong state, in which each cycle is characterized by a single instantaneous impact (see Figure 18b,d).

Time Histories of Excitations at the Pivots in Case of Large Oscillations
Adopting the approach already described in Section 8.1.2. for case II-undamped free oscillations with θ 0 = 40 • , the time histories of the excitations at the bell and clapper pivots in case of full circle oscillations (θ 0 = 171 • ) have also been determined, considering undamped and friction damped oscillations. Recalling that the "equivalent" viscous damping case, whose parameters accord with Table 9, leads to outcomes practically coinciding with the ones obtained for friction damping, this case has not been explicitly considered.
The time histories of the accelerations are compared in Figure 22, where solid lines refer to the undamped situation and dashed lines to the damped situation.
of full circle oscillations = 171° have also been determined, considering undamped and friction damped oscillations. Recalling that the "equivalent" viscous damping case, whose parameters accord with Table 9, leads to outcomes practically coinciding with the ones obtained for friction damping, this case has not been explicitly considered. The time histories of the accelerations are compared in Figure 22, where solid lines refer to the undamped situation and dashed lines to the damped situation.  In analogy with Figure 16, Figure 22a,b refers to horizontal components a rA,x /g and a rB ,x /g, while Figure 22c,d refers to vertical components a rA,y /g and a rB ,y /g. The diagrams show that the apparent period of damped large oscillations decreases over time, as a consequence of the reduction of the amplitude of the oscillations, which largely compensates the small increase due to the (low) damping considered in the present case. Obviously, this aspect is particularly relevant not only for the assessment of the bell behavior but also for the evaluation of the actions to be considered dealing with bell-structure interactions.

Synchronized Resonant Forced Oscillations of the Bell in Transient Phase
As known, manual excitation of a ringing bell in the starting phase is a typical example of resonance, where the forcing frequency is suitably tuned by the bell player in such a way that it is constantly synchronized with the natural frequency of the bell, which decreases as the oscillation amplitude increases. From the mathematical point of view, modelling this operation is not trivial, due to the complexity of the relationship between the natural frequency of the oscillations, their amplitude, the damping parameters, and their variations over time. For that reason, the study of this relevant transient phase is often disregarded or drastically simplified, for example considering easy to handle periodic forcing functions.
Dealing with forced oscillations, external forces are usually modeled by means of cosinusoidal functions of the time, so opening the door for easier analytical treatment. A similar approach is often adopted also in numerical integration of the equations of motion, especially when the solution is achieved via standard built-in numerical algorithms. Obviously, if the problem involves resonance aspects, the frequency of the forcing cosinusoidal function should be continuously modified to catch the actual natural frequency of the system, as a function of the oscillation amplitude. That approach is often worthless since it requires systematic human intervention to proceed with the solution. A common line of action is to consider forcing functions characterized by fixed frequency, suitably chosen in the range of natural frequencies associated with the expected range of oscillation amplitudes. Difficulties further increase considering external excitations described by arbitrary functions of the time, like, for example, the stepwise constant real excitations applied by bell ringers or by linear or rotating motors.
The proposed algorithms overcome all these issues. Really, they are able to successfully manage quite arbitrary excitations, M e (θ, t) and m e (ϕ, t), of bell and clapper, also in case their dependence on time is not explicit. That important original feature easily allows to automatically synchronize the excitation with the varying instantaneous natural frequency of the system, reproducing persistently resonant conditions and, consequently, minimizing the peaks of the excitation.
A possible workaround is considering torques explicitly depending only on the rotations, so that the dependence on t is implicitly expressed through θ = θ(t) and . θ = . θ(t) or ϕ = ϕ(t) and . ϕ = .

ϕ(t).
It is worthwhile to remark that, even if the focus of the study is bell excitation, the approach is much more general, and it is not subject to special limitations.
Moreover, in defining the forcing function it should considered that, depending on the arrangement, external torque can be unidirectional or bidirectional, like in the case handstroke and backstroke consent to swing the bell through a full circle.
To assess the potentiality of the method, it has been studied the transient situation corresponding to the time interval needed to achieve the bell working position, starting from the rest condition. In the simulation, once achieved the bell working position, θ max ≈ 171 • , the excitation has been switched off, so subsequently leaving the bell free to oscillate.
Four cases have been analyzed, characterized by two different forcing functions and by undamped or damped oscillations. In all cases, the following initial conditions have been assumed: which correspond to an initial small rotation of the bell, with the clapper in vertical position. Aiming to reproduce a manual excitation, it has been adopted a piecewise-constant bidirectional forcing function, which is the more realistic one: or a truncated sinusoidal bidirectional forcing function: It must be remarked that the peak modulus of M e1 and M e2 , 250 Nm, matching the torque easily exerted by a person, is significantly smaller than the peak of the torque needed to lift statically the St. Mary's tenor bell in the upward position, M max ≈ 10 kNm.
Since "equivalent" viscous parameters strongly depend on the amplitude of the oscillations, the damped case refers to friction damping. Friction coefficients and pivot radii are those already recalled in Tables 8 and 9, leading to very low level of damping, corresponding to a damping ratio ζ ≈ 0.1%.
Relevant curves obtained in the undamped and in the friction damped case are compared in Figures 23-25. In the figures, panel (a) refers to M e1 (Equation (74)), and panel (b) to M e2 (Equation (75)); curves referring to undamped case are generally in blue and curves referring to damped case are generally in red, except ϕ − t curves in Figure 23, whose colors are opposites.
It must be remarked that the peak modulus of and , 250 Nm, matching the torque easily exerted by a person, is significantly smaller than the peak of the torque needed to lift statically the St. Mary's tenor bell in the upward position, ≈ 10 kNm.
Since "equivalent" viscous parameters strongly depend on the amplitude of the oscillations, the damped case refers to friction damping. Friction coefficients and pivot radii are those already recalled in Tables 8 and 9, leading to very low level of damping, corresponding to a damping ratio ≈ 0.1%.
Relevant curves obtained in the undamped and in the friction damped case are compared in Figures 23-25. In the figures, panel (a) refers to (Equation (74)), and panel (b) to (Equation (75)); curves referring to undamped case are generally in blue and curves referring to damped case are generally in red, except − curves in Figure 23, whose colors are opposites.  Looking at the time histories of the rotations of bell and clapper (Figure 23), it results that: • in the undamped case, the time needed to reach the maximum amplitude, t d , i.e., the time during which the bell is driven, is t d ≈ 83 s, when the forcing function is M e1 , increasing to t d ≈ 93 s, when the forcing function is M e2 , • in the friction damped case, the time needed to reach the maximum amplitude is t d ≈ 110 s, when the forcing function is M e1 , increasing to t d ≈ 121 s, when the forcing function is M e2 , • also considering the transient phase, and without stopping it in the upward position before the free oscillation phase, the bell tends to ring wrong in all considered cases.
Considering forced oscillations, the comparison of the diagrams obviously confirms that the number of oscillations needed to reach a given state in the transient phase is bigger in damped conditions than in undamped conditions. Anyhow, it is interesting to remark that the amplitude and the period of these forced damped oscillations are generally smaller than those pertaining to forced undamped oscillations, even if, clearly, as already anticipated, the apparent period of a given large oscillation is higher in damped condition.  Considering forced oscillations, the comparison of the diagrams obviously confirms that the number of oscillations needed to reach a given state in the transient phase is bigger in damped conditions than in undamped conditions. Anyhow, it is interesting to remark that the amplitude and the period of these forced damped oscillations are generally smaller than those pertaining to forced undamped oscillations, even if, clearly, as already anticipated, the apparent period of a given large oscillation is higher in damped condition.
The latter observation is diagrammatically illustrated in Figure 26, where, referring to the present case study, the apparent periods T * of the forced resonant oscillations are plotted as a function of the maximum amplitude of the oscillations themselves θ * max . In the diagram, dashed line represents the damped case, solid line the undamped case. Evidently, due to the low damping ratio characterizing the case study, the increase of the period due to damping is significant only when the oscillation amplitude is large. The latter observation is diagrammatically illustrated in Figure 26, where, referring to the present case study, the apparent periods * of the forced resonant oscillations are plotted as a function of the maximum amplitude of the oscillations themselves * . In the diagram, dashed line represents the damped case, solid line the undamped case. Evidently, due to the low damping ratio characterizing the case study, the increase of the period due to damping is significant only when the oscillation amplitude is large.
Finally, adopting a representation in line with that adopted for Figures 23-25, in Figure 27 they are illustrated, as functions of the time, the resonant forcing functions resulting, in the considered cases, in damped and undamped conditions. Obviously, Figure 27a refers to (Equation (74)), Figure 27b to (Equation (75)). The graphs clearly demonstrate that the period of the forcing frequency is automatically tuned on the period of the system, so certifying the capability of the proposed approach to simulate also constantly resonant excitations.
To complete the case study, the energy dissipated in the transient phase of forced damped oscillations are represented as a function of the time in Figure 28. In the figure, the blue line refers to and the red line to . The total energy dissipated to attain the maximum amplitude results ≈ 5.15 kJ, if the forcing function is , increasing to ≈ 5.85 kJ, if the forcing function is , so indicating that, as expected, excitation is more efficient.    Obviously, Figure 27a refers to M e1 (Equation (74)), Figure 27b to M e2 (Equation (75)). The graphs clearly demonstrate that the period of the forcing frequency is automatically tuned on the period of the system, so certifying the capability of the proposed approach to simulate also constantly resonant excitations.
To complete the case study, the energy dissipated in the transient phase of forced damped oscillations are represented as a function of the time in Figure 28. In the figure, the blue line refers to M e1 and the red line to M e2 . The total energy dissipated to attain the maximum amplitude results E dis ≈ 5.15 kJ, if the forcing function is M e1 , increasing to E dis ≈ 5.85 kJ, if the forcing function is M e2 , so indicating that, as expected, M e1 excitation is more efficient. Evidently, the duration of the transient phase depends on shape and intensity of the forcing excitation, as well as on the initial conditions, but the proposed algorithms are not sensitive on it. In fact, as already remarked at the end of Section 6, the procedure can be repeated as many times as necessary, duly modifying the initial condition, until the wanted final condition is achieved.

Conclusions
Study of swinging bells involves topical arguments of structural and mechanical engineering as well as sound and acoustic engineering.
The research focuses, first, on the system of two nonlinear second order ODEs governing the problem. That system is derived under very general hypotheses, considering viscous damping, friction damping, and arbitrary excitations. Dissipative forces are treated like external forces interpreting the equation of motion of bell and clapper as dynamic equilibrium equations, in the spirit of Newtonian mechanics.
In the study, besides viscous forces proportional to the angular velocities, friction forces at the bell and clapper pivots are introduced as well. In that novel perspective, friction forces, commonly disregarded in the traditional approach, are realistically expressed as a function of the time dependent reactive radial forces at the pivots.
Clapper strokes are modelled as instantaneous impacts, resorting to the kinematic definition of the coefficient of restitution, which allows a synthetic, and experimentally measurable, quantification of the complex interdependency between the relevant parameters governing the impact. The suitable application of the law of conservation of the angular momentum provides the bell and clapper angular velocities at the end of the impact, which can be easily modified to consider also the acoustic energy leaving the system. Moreover, to preserve the physical coherence avoiding paradoxical results, the upper bound of the coefficient of restitution, corresponding to perfectly elastic impact, is suitably derived.
In the relevant literature about classical bell dynamics, the focus is alternatively concentrated on the bell-clapper interactions or on the bell-tower excitation. Whatever the problem, analytical solutions are available in a limited number of cases, therefore the governing ODEs are usually solved using numerical methods. The most frequently used numerical approaches rely on standard built-in algorithms, available in advanced software packages. Since adoption of these general algorithms is subject to some restrictions and often requires some additional manipulation of the equations, especially regarding the forcing functions, the significance of the outcomes derived in that way is limited. Moreover, the piecewise-defined structure of the solution makes the use of these classical methods even more complex. Evidently, the duration of the transient phase depends on shape and intensity of the forcing excitation, as well as on the initial conditions, but the proposed algorithms are not sensitive on it. In fact, as already remarked at the end of Section 6, the procedure can be repeated as many times as necessary, duly modifying the initial condition, until the wanted final condition is achieved.

Conclusions
Study of swinging bells involves topical arguments of structural and mechanical engineering as well as sound and acoustic engineering.
The research focuses, first, on the system of two nonlinear second order ODEs governing the problem. That system is derived under very general hypotheses, considering viscous damping, friction damping, and arbitrary excitations. Dissipative forces are treated like external forces interpreting the equation of motion of bell and clapper as dynamic equilibrium equations, in the spirit of Newtonian mechanics.
In the study, besides viscous forces proportional to the angular velocities, friction forces at the bell and clapper pivots are introduced as well. In that novel perspective, friction forces, commonly disregarded in the traditional approach, are realistically expressed as a function of the time dependent reactive radial forces at the pivots.
Clapper strokes are modelled as instantaneous impacts, resorting to the kinematic definition of the coefficient of restitution, which allows a synthetic, and experimentally measurable, quantification of the complex interdependency between the relevant parameters governing the impact. The suitable application of the law of conservation of the angular momentum provides the bell and clapper angular velocities at the end of the impact, which can be easily modified to consider also the acoustic energy leaving the system. Moreover, to preserve the physical coherence avoiding paradoxical results, the upper bound of the coefficient of restitution, corresponding to perfectly elastic impact, is suitably derived.
In the relevant literature about classical bell dynamics, the focus is alternatively concentrated on the bell-clapper interactions or on the bell-tower excitation. Whatever the problem, analytical solutions are available in a limited number of cases, therefore the governing ODEs are usually solved using numerical methods. The most frequently used numerical approaches rely on standard built-in algorithms, available in advanced software packages. Since adoption of these general algorithms is subject to some restrictions and often requires some additional manipulation of the equations, especially regarding the forcing functions, the significance of the outcomes derived in that way is limited. Moreover, the piecewise-defined structure of the solution makes the use of these classical methods even more complex.
In the paper an innovative procedure to solve the equations of motions of bell and clapper is discussed.
The novelty of the proposed method, which is described in Section 6, is an original numerical explicit-implicit predictor-corrector integration scheme with constant time step, implemented in two different slightly different versions, depending on the correction phase regards only the rotations or the angular velocities and the angular accelerations as well. The implementation is such that each time the clapper strikes the bell a new "piece" of the solution is automatically initialized, consequently user interventions in the elaboration phase are not necessary.
In addition, the structure of the algorithm allows to analyze arbitrarily long transient and stationary time intervals, until the required final situation is achieved, provided that the total time to be analyzed is divided in suitable subintervals. Evidently, the duration of each subinterval should be chosen according the available physical memory of the computer hardware, while the initial conditions of the system in each subinterval are obtained duly recalling the final conditions of the system at the end of the preceding subinterval.
These extremely fast and effective algorithms, which have been validated in Section 7 reproducing relevant theoretical, numerical, and experimental outcomes, allow one to successfully handle, independently on the oscillation amplitude: undamped oscillations; friction and viscosity damped oscillations; free oscillations; forced oscillations; transient and stationary phases, and are also suitable to solve stiff equations.
The method has been applied in Section 8 to study the oscillations of a relevant bell, the Great St. Mary's old tenor bell in Cambridge, hypothesizing several distinct operational conditions, also highlighting that the definition of damping equivalence requires nontrivial considerations.
The critical examination of the outcomes of the analyses demonstrates the potentialities of the method highlighting some original key issues, like its capability to handle arbitrary forcing functions.
A particularly significant and pioneering feature of the algorithms is their ability both to manage and to define "resonant" forcing functions, continuously tuned on the actual frequency of the oscillation, which depends on the amplitude of the oscillation itself. Outstanding examples of resonant forcing functions are, for example, the excitations that bell ringers apply on manually driven bells.
It must be stressed that the approach is absolutely general and can be applied independently on the adopted ringing style; as well as on dimensions, shape, and mechanical properties of bell and clapper.
Further refinements, especially regarding other sources of dissipation; linear and nonlinear external elastic restraints; and interaction with the supporting structures can be envisaged, which will be the subject of future studies.