Estimation of parameters of various damping models in planar motion of a pendulum

In this work, planar free vibrations of a single physical pendulum are investigated both experimentally and numerically. The laboratory experiments are performed with pendula of different lengths, for a wide range of initial configurations, beyond the small angle regime. In order to approximate the air resistance, three models of damping are considered—involving the three components of the resistive force: linear (proportional to velocity), quadratic (velocity-squared) and acceleration-dependent (proportional to acceleration). A series of numerical experiments is discussed, in which the damping coefficients are estimated by means of several computational methods. Based on the observed efficiency, a gradient method for optimization is treated as the main tool for determination of a single set of damping parameters, independent of the system’s initial position. In the model of resistive force, the term proportional to acceleration, associated with the empirical Morison equation, seems to be indispensable for the successful approximation of the real pendulum motion.


Introduction
The proper selection of a damping model and estimation of damping parameters is an essential problem in the area of dynamic simulation and analysis of real mechanical systems. Major difficulties result from a wide variety of energy dissipation mechanisms, their complexity and coexistence, as well as the approximate character of damping models. Typically, the inherent damping, i.e. naturally present within the system or its environment, is divided into three types: internal (material) damping, structural damping (at joints, supports or interfaces), and external damping (via fluid-structure interactions) [7].
There are several elementary models used to describe different forms of mechanical energy dissipation, such as the viscous and hysteretic damping, the Coulomb dry friction, or the velocity-squared drag model. Even if more sophisticated and physically detailed approaches to dissipation are available (e.g. within the theory of viscoelasticity or the theory of thermoelasticity for the case of internal damping), they may be too complex and therefore underused. In computational practice, relatively simple models are usually employed. The purpose is to approximately represent the overall damping in a given system, and capture only its general effects (e.g. the amplitude decay). In light of this, the basic models, especially the viscous damping or-a kind of its extension-the proportional damping, can serve as an equivalent or ''resultant'' representation of not only the internal, but also structural and external damping [2,3,7,27].
Such an approach with the limited interest in the actual damping sources and mechanisms can be mainly justified by the fact that the amount of energy dissipated within numerous real systems is very small (heavy machines and structures, conventional materials etc.). The magnitudes of forces associated with a certain type of damping may be regarded as negligible compared to the inertial and elastic forces. However, the effect of dissipative forces can be of high importance for relatively low-mass structures undergoing large motions. Good examples are flexible slender systems like chains, cables, ropes and highly flexible beams or blades [2,19]. Taking into account the nature of oscillations and the related geometric nonlinearity, it is quite natural to consider a pendulum to be a kind of archetype of such slender systems.
For many centuries, pendula of various types (normal/inverted, simple/physical, single/multiple, etc.) have been attractive objects of studies [1,15,22,30]. It is well known that they play a vital role in the field of theoretical and experimental mechanics, and are used to demonstrate essential concepts of, among others, Newtonian and analytical mechanics, small vibration theory, nonlinear dynamics and chaos, analysis of coupled oscillators, multibody dynamics (e.g. see [28,36,37,39]). Moreover, pendula became standard models and a source of test problems in other related areas such as automatic control and robotics, space technology, computational methods and algorithms, artificial intelligence, etc. When it comes to analysis and simulation of the slender systems, multiple pendula have been extensively used in investigations of chain dynamics [20,23,34,38,40], and sometimes as a basis of discrete (lumped-mass) models of elastic continua [9,10,32].
In the context of vibration damping, a single pendulum may also be seen as a very simple and convenient object of experimental-computational studies. Historical and mathematical background of the problem of a simple pendulum motion with air resistance was presented by Dahmen [6] and Gauld [12]. Nelson and Olsson [26] used a plane pendulum to measure the gravity acceleration, and discussed several corrections related to a rich variety of the physical features of a real pendulum, including the aerodynamic effect. The contribution from the drag on the string into the overall damping of a simple pendulum was analyzed by Mohazzabi and Shankar [24]. Ladino and Rondón [16] proposed a method to determine the damping coefficient of a simple pendulum with oscillating suspension point. The influence of both dry friction (at a support) and viscous damping on the pendulum amplitude was discussed by Zonetti et al. [42]. It should be noted that the mentioned studies are mainly focused on the case of a simple pendulum, and usually on a single damping model: the classical linear (velocity-proportional) or the quadratic (velocity-squared) damping.
In this article, planar free vibrations of a physical pendulum are investigated experimentally and numerically. The laboratory experiments are performed with pendula of different lengths, for a wide range of the initial swing angle-without the assumption of small oscillations. In a series of numerical simulations, three models of external damping with constant coefficients are used and the parameters are estimated by means of various computational methods. The aim of this study is to assess and compare the prediction accuracy of the three models, and to select the most appropriate one.
The paper is organized as follows. In order to present a gradual extension of the viscous damping model, the dynamic equations for the physical pendulum, based on the Lagrangian approach, are outlined in Sect. 2. The experimental setup and the details of the laboratory study are described in Sect. 3. Sections 4, 5 and 6 contain the results of parametric estimation for each of the three damping models, performed with the

Formulation of the problem
Let us consider a physical pendulum depicted schematically in Fig. 1. The system consists of a rod of length L and mass m, suspended from point O. The in-plane free oscillations of the pendulum can be described by angle uðtÞ, classically measured with respect to the vertical axis.
Neglecting any resistance forces, the equation of motion for the one-degree-of-freedom system can be easily derived using the Lagrange equation of the second kind: where u is the only generalized coordinate, _ u ¼ du dt is the angular velocity, and L ¼ T À V denotes the Lagrange function. The kinetic and potential energy of the system are given by where I O is the mass moment of inertia of the rod about z-axis (passing through point O). Now, applying Eq. (1) leads to Mathematically, this idealized model is a second-order ODE with constant coefficients. If energy dissipation in the mechanical system is taken into account, the Lagrange equation becomes d dt where Q d is the generalized damping force. Classically, the quantity is assumed to be derivable from the Rayleigh dissipation function which is a quadratic form in the generalized velocity: where c 1 is a constant coefficient. Consequently, the equation of motion includes a linear dissipative (viscous) term: The viscous damping model is the simplest one, readily used in vibration analysis, since it does not lead to any nonlinearity and computational difficulties. It is sometimes employed to approximate more complex types of internal and external damping in real systems (e.g. via the concept of equivalent viscous damping) [2,3,11]. However, such a linear model may be regarded as the lowest order approximation of the resistance coming from the medium surrounding a moving body. More generally, the damping force F d as a function of the oscillator velocity v can be expanded in power series in the velocity [8,17,24,35]: where a i are constants (i ¼ 0; 1; 2; . . .) and v ¼ jvj. Naturally, a 0 ¼ 0 because there is no resistance acting on the body at rest. In most textbooks on classical mechanics the topic of resisted motion is discussed using solely the first-or second-order term of the expansion, arguably for the sake of simplicity (e.g. see [14,36,37]). From the viewpoint of fluid mechanics, in turn, the well-known drag model is said to be pure quadratic. However, one should notice that the nondimensional drag coefficient involved is a function of, among others, the body Reynolds number: Since Re is proportional to the speed v, the final form of the resistive force may be more complicated (not necessarily purely quadratic in velocity). Particularly, at low Reynolds numbers, it may consist of the linear (or approximately linear) term [25,41]. All in all, the dynamic Eq. (6) is extended by the second-order term: For the purpose of the combined computationalexperimental studies, the linear-quadratic model included in Eq. (8) is enriched with an accelerationrelated part: Occurance of such a term in the expression for the resistance force should not be seen as unusual. In the case of oscillatory motion of elementary bodies (sphere, infinite cylinder etc.) in a viscous fluid, purely analytical solutions (for low Reynolds numbers) lead to the drag having two parts: proportional to velocity (the dissipative part) and proportional to acceleration (the inertial part) [18,33]. Also in the empirical Morison equation, commonly used in hydromechanics and dynamic analysis of offshore structures, the total wave force on a submerge body includes the drag and inertia forces, where the latter is related to the addedinertia concept [29,33]. Let us introduce the dimensionless time s ¼ x 0 t, where x 0 is the natural frequency of the undamped pendulum that undergoes small-amplitude vibrations (the linearized case): Consequently, equations of motion (6), (8) and (9) take the following non-dimensional forms: where now u ¼ uðsÞ and an overdot denotes differentiation with respect to s, while are the non-dimensional damping coefficients corresponding to the three considered components of the resistive force: linear (proportional to velocity), quadratic (velocity-squared) and acceleration-dependent (proportional to acceleration).

Experimental studies
The laboratory experiments were conducted with pendula of various length: L 2 f100; 150; 200; 300g mm. In each case the rod has a diameter d ¼ 4 mm and is made of an aluminium-magnesium alloy (AlMg3), which is often used for welding rods. The material density is q ¼ 2640 kg/m 3 . The ends of the rod are threaded in order to attach conical sleeves by connecting nuts (see Fig. 2a). The main purpose of the sleeves is to clamp a synthetic cord that plays a role of a connector between the rod and the support, or between two adjacent members of a multiple pendulum (in our future work, see Fig. 2b). Moreover, the nuts are a convenient location for tracking markers (see Fig. 2c). As shown in the photographs, the pendulum is mounted to the frame by the cord connection between the tip sleeve and the sleeve fixed to an auxiliary bracket. The relatively complex geometry of the pendulum link has been taken into account carefully to calculate mass and moment of inertia of the system. The inertia parameters for all the length variants are given in Table 1; additionally, the natural frequencies for the linear case are included.
Planar motion of the system was investigated by an optical measurement technique. A high speed camera was used (Photron 1024 PCI, max. resolution of 1024 Â 1024 pixels, max. frame rate of 10,000 fps). The experiments were performed with the frame rate of 250 fps, which allowed for recording motion lasting about 12 s. The control of the camera settings as well as the data acquisition process were done via a computer with a specialized software. The motion analysis, i.e. identification and tracking of the marker at the pendulum end, and determination of its position in particular frames, was conducted using the Photron Motion Tools application. The angle of pendulum swing was found based on the Cartesian coordinates ðx; yÞ of the marker and the pivot point. The experimental setup is presented schematically in Fig. 3.
The purpose of these experiments is to record behaviour of the pendula in a wide range of the swinging angle. Thus, we focused on four initial configurations of the system: Since the pendula were realeased from rest by hand, the starting angles were imposed on them only approximately. The values of u 0 calculated from the measured positions of the marker at the beginning of motion are given in Table 2. Here and in the further stages the initial configurations are numbered by k.

Application of finite difference schemes
Let us start with the linear damping model included in Eq. (10). The damping coefficient will be approximated by means of the following three methods: the finite difference schemes, the bisection method and the gradient method for optimization.
The first approach is focused on the approximation of angular velocity x ¼ _ u and angular acceleration e ¼ € u at a time point by the values of angle u measured experimentally. Since the marker coordinates as the experimental results contain errors, and the finite-difference representation of e may be especially sensitive to small disturbances, highly accurate (multi-point) difference schemes must be applied.   Consider a uniform grid, i.e. a sequence of equispaced time points (measurement instants). Central differece schemes are based on a stencil of n Ã ¼ ð2n þ 1Þ nodes: where f ðjÞ ðt i Þ is the jth time derivative of f, calculated at t i , whereas the remainder R p ðt; t i Þ satisfies the condition: Usually the method of undetermined coefficients is employed to derive an appropriate finite difference formula. With such an approach, expansion (14) up to term of order n Ã is used for the n Ã -point stencil. This leads to a n Ã Â n Ã Vandermonde system of linear equations that becomes poorly conditioned for large n Ã [21]. Here we apply another technique based on symbolic computation. We tend to use only three terms of series (14), i.e. just the derivatives that are to be determined. Thus, it is assumeed that f(t) can be represented by where play the role of basis functions, equivalent to the power functions from the Taylor series. Each coefficient d s , in turn, denotes the value of the derivative of order ðs À 1Þ of f at t i which can be approximated by the central difference scheme. Namely, where are the desired coefficients. Substitution of (16) and (18) into Eq. (15) produces the equation residual As in the Galerkin method, the quantities d s can be obtained by imposing the orthogonality conditions: where h; i is the inner product, which we define in a discrete sense, that is This gives a 3 Â 3 system of algebraic equations that can be written as where d ¼ ½d 1 ; d 2 ; d 3 T , and elements of A, b are given by Since the arguments and nodal values of f are kept in symbolic form ( , the symbolic solution vector d contains three linear combinations (18), and the desired finite difference coefficients b ðsÞ iþj can be easily extracted. Generally, for these studies, the central difference schemes based on n Ã ¼ 25 nodes (n ¼ 12) have been constructed, i.e. the first and second derivatives of uðtÞ are evaluated using the angle values measured at 25 time instants. In the case of first and last n time points, asymmetric schemes involving 49 nodes are used, because they are two times less accurate than the central-difference approximation.
Next, one of the simplest smoothing technique is applied to the calculated time series x i and e i , i.e. the moving average smoothing: where the parameter q is taken equal to 10. Consequently, at our disposal we have the observed time series u i and the smoothed onesx i ,ê i . Finally, the damping coefficient c 1 can be determined over time from Eq. (6): Additionally, the generalized viscous force is evaluated: The results for the case of L ¼ 100 mm are presented in Fig. 4. As can be seen, coefficient a 1 does not remain constant. Although very accurate finite-difference approximations have been used for the angular velocity and acceleration, the linear damping model exhibits some inadequacy. The results obtained for the other pendulum lengths and different initial angular position u 0 lead to scatter plots of the same character.

Bisection method
In order to determine certain values of the damping coefficient a 1 , some optimization methods can be used. The first approach utilized in this paper is the bisection method which basically belongs to the class of bracketing root-finding techniques [5,31]. We restrict our attention to the interval X ¼ ha L ; a R i, where a L , a R are two guesses that are supposed to bracket the optimal solution. In the first step, X is divided in half, that is the midpoint For both these values of the damping coefficient we integrate numerically the differential Eq. (10). The results are compared with the experimental time series, i.e. the total distance is calculated: where N is the number of frames recorded, while u exp i; k and u num i; k denote the experimental and numerical values of the pendulum angle for the kth initial configuration. The function f k expresses the distance in the sense of 2-norm. The minimum value of f k ða LM Þ, f k ða RM Þ indicates the subinterval within which the desired solution is located. Consequently, we redefine the search interval X as X L or X R , and simultaneusly a M is replaced with the initial estimate of the optimum, a LM or a RM . The process is repeated until a satisfactory result is achieved (the error falls below some threshold), or the difference between errors (24) in two subsequent iterations is small enough.
For the determination of a 1 it is reasonable to focus on the range h0; 1i, since the coefficient is a nonnegative and relatively small parameter. The results obtained for all the pendulum lengths are presented in Table 3. Instead of the total distance between the numerical solutions and the experimental data, f k , the root-mean-square (RMS) error and the average (AVG) error are given: Moreover, the total mean errors are calculated as follows: Again, at a constant length of the pendulum the damping coefficient takes different values as the initial position is varied. The errors, in turn, reach quite high values, especially in the case of short pendula (L ¼ 100 mm and L ¼ 150 mm). For example, for L ¼ 100 and k ¼ 3 the RMS error is e rms k % 8:2 , which gives a relative error of 6% when treating the initial angle u 0 as the reference one (see Table 2). Simultaneously, the average distance is e avg k % 0:2 , and its relative value is less than 0.2%.
A graphical comparison between the experimental and numerical time histories of the angular position a b Fig. 4 The finite-difference-based dependencies: a c 1 ðx i Þ, b Q 1 ðx i Þ. Results obtained for L ¼ 100 mm and k ¼ 1 uðtÞ is shown in Fig. 5. Slight differences in the oscillation amplitudes occur, particularly at the maxima, and a phase difference becomes more noticeable with time. However, the general character of the measured oscillations is well reflected by the numerical solution. To give a better insight into computational-experimental discrepancies, history plots focused on the first and the last 2 s of motion will be presented. As can be seen in Fig. 6, prediction quality of the model is not very poor, although the phase shift grows with time and the initial angle. Figure 7 presents a simple visualization of the pendulum motion for L ¼ 100 mm and k ¼ 3. The real motion was recorded with the frame time Dt ¼ 0:004 s, whereas the interval between the depicted configurations is usually equal to 10 Â Dt or 20 Â Dt. The increasing phase shift between the numerical and experimental results is well noticeable. The angular difference (in absolute value) between the configuration pairs shown in Fig. 7a ranges approximately from 0:1 to 5:8 whereas the difference in Fig. 7b ranges from 0:6 to 11:3 . Selected frames of this motion are presented in Fig. 8 in the form of the simulation-based configurations overlaid on snapshots of the real pendulum.

Gradient method for optimization
To determine a single value of the damping coefficient for a given pendulum length, that is independent of the initial configuration, the gradient method is used. Usually the gradient-based approach is applied to multidimensional unconstrained optimization. Here we perform a one-dimensional search with the first derivative evaluation [5,31].
We aim at minimization of the distance of type (24), but simultaneously for all the studied initial positions of the system. Thus, the objective function is constructed as follows: where f k has the form (24).
Firstly, we select an appropriate initial guess, a 1 . It can be done, for example, by taking the arithmetic mean of the results given in Table 2 for particular values of L. Additionally, the initial search direction d ð0Þ and the step size Da ð0Þ must be chosen. In the onedimensional case we move along a line (a 1 -axis) by the distance Da ðpÞ , and the direction d ðpÞ takes the value 1 or À 1.
In the optimization procedure, equations of motion (10) are solved at the starting value a 1 , and the error is evaluated, e ð0Þ . Next, the same computation is performed for the trial point a , and error f ð1Þ is obtained as a result. A new direction is determined by using the approximate gradient value: and the descent direction is Finally, the next step size can be expressed as where w 1 and w 2 are some control parameters from the interval ð0; 1Þ. Their role is to speed up the iteration process when the consecutive steps apparently tend to the minimum, or decrease the step size after a change of direction. In further stages, Eq. (30) becomes The process is stopped if a satisfactory error level is achieved, the difference between errors in two subsequent iterations is sufficiently small, or the maximum number of iterations is reached. Table 4 contains the obtained results for all the pendulum lengths. Apart from the optimal values of coefficient a 1 , the root-mean-square and average errors (25) for particular initial configurations are presented. As it is indicated by e rms k and e avg k values, if the time series for each initial angle u 0 is treated separately, the solutions are non-optimal, that is the bisection method leads to lower errors. The fourth initial configuration generates significantly greater errors than the others. Unlike in the latter case, the total errors are not considerably lower for the larger lengths of the pendulum. However, the main advantage of this approach is to find only one ''effective'' value of the damping coefficient (independent of u 0 ), and such attempts will be continued with the more advanced models.

Successive search and intersection technique
The results for the linear damping suggest that this model may be insufficient. Therefore, a series of computations have been conducted for the extended, linear-quadratic model included in Eq. (11). To find coefficients a 1 and a 2 , the optimization problem has been solved by means of two methods: the successive search in combination with the intersection technique, and the gradient method. The first approach is focused on a graphical representation of the objective function f ða 1 ; a 2 Þ, in the identical form as (27), over a two-dimensional domain. Firstly, the search is narrowed to the region X ¼ h0; ai Â h0; bi, where a; b [ 0. Next, the domain is discretized into a uniform grid of n 1 Â n 2 points (see Fig. 9), characterized by the spatial step sizes in the two directions Let a 1;i and a 2;j be the ith and jth values of the coordinates a 1 and a 2 on the plane, respectively (for i ¼ 0; 1; 2; . . .; n 1 and j ¼ 0; 1; 2; . . .; n 2 ). For every grid point fa 1;i ; a 2;j g equation of motion (11) has been solved, and the nodal value f i;j of the objective function has been evaluated. As an example, the results for L ¼ 100 mm are shown graphically in Fig. 10. In all the studied cases, a continuous representation of the the discretized surface f i;j ða 1;i ; a 2;j Þ exhibits an oblong valley, and one can suppose that there exist infinitely many minima. For a better insight, the straight line (e.g. described as a 2 ¼ Aa 1 þ B) in the longitudinal direction of such a ''trough'' (at its bottom) has been found numerically. Then, it can be observed that the valley direction varies with the pendulum initial angle u 0 . Finally, to determine a single representative pair of the damping coefficients for all the initial configurations, the intersection points of these lines have been analyzed within an unconstrained domain. From all the intersections, the one for which the function f takes the lowest value, is chosen as the optimal solution.
The results obtained with the intersection technique are given in Table 5. As can be seen, coefficient a 1 turns out to be negative in each case, but it does not have a physical meaning and indicates that the model is inappropriate.

Gradient method for optimization
The intersection technique has a relatively low accuracy, because a small disturbance of the slope coefficient of a line can lead to completely different results. Therefore, another attempt has been made to find a 1 and a 2 -by the use of the gradient method.
In the two-dimensional problem, the idea of the method is identical as in the one-dimensional case. However, the motion towards the minimum of the objective function takes place on a plane. Thus, the values of the decision variables at the pth step form a vector and, by analogy to Eq. (29), the new direction is given by where j j denotes the Euclidean norm. Equation (31), in turn, now takes the form where the right-hand side includes a dot product of the previous and new direction vectors. Table 6 contains the estimated values of the parameters and the related errors for every pendulum length. Again, the damping coefficient a 1 has negative values, which excludes the application of the linearquadratic model.
In order to make full use of the two-term model, the purely quadratic damping is also examined. Thus, the gradient-based optimization outlined above has been performed for equation of motion (11) when a 1 ¼ 0. As shown in Table 7, coefficient a 2 takes smaller values now. In the case of all the pendulum lengths and initial configurations the truncated damping model, in combination with the gradient method, provides better prediction of the real motion than the linear model (see a b Fig. 10 Graphical representation of the objective function: a surface plot, b contour plot. Results for L ¼ 100 mm and k ¼ 2   Table 4). Moreover, as our tests indicated, the same results can be obtained using relatively simple methods for the constrained optimization problem with two decision variables, a 1 and a 2 (e.g. the projected gradient method, the penalty function method or the Monte Carlo approach [4,13]). Figure 11 illustrates the responses of the pendulum of length L ¼ 200 mm as k ¼ 1 and k ¼ 2. Although the values of e rms k and e avg k belong to the lowest ones, significant differences in the angular displacements can be observed. Here, the phase shift is much less noticeable than the amplitude mismatch-the oscillation decay is stronger in the real mechanical system. 6 Numerical results for velocity and acceleration dependent damping

Gradient method for optimization
As could be seen in Sect. 5, the enrichment of a model by an additional term not necesserily leads to the desired results. Consequently, the three-component damping model included in Eq. (12) have been considered with the hope that it will work in the full form, i.e. simultaneously with all the (non-zero) coefficients. In other words, the third, acceleration dependent term may be the missing element, and its presence can make the other terms play their proper role in the dynamic simulation of the system under real conditions. In this section, the coefficients a 1 , a 2 and a a are determined by using the gradient method. Thus, the objective function (27) is now regarded as dependent on three decision variables: f ða 1 ; a 2 ; a a Þ. As before, the optimization procedure is simply extended to the three-dimensional case. Since the vector of decision variables (32) becomes composed of the triple coefficients, vectors (33)- (35) are extended to have a third element, corresponding to a a .
The results obtained for each length of the pendulum are given in Table 8. In all the cases, the following relation holds: a 1 \a 2 \a a , which shows that the acceleration-related effect is not negligible. The total RMS error reach values below 5 while the total AVG errors are less than 0:1 (compare Tables 3, 4 and 7). In particular, the errors are at least 5 times smaller than  Results obtained with the gradient method those for the linear model and the gradient-based optimization. Apparently, the velocity and acceleration dependent damping model is most appropriate. As can be seen in Fig. 12, there are relatively small phase shifts between the compared results. One can also notice just slight differences in the peak amplitudes. The numerical solutions of the equation of motion coincide quite accurately with the experimental observations.
The sequential visualization of the pendulum motion for L ¼ 150 mm and k ¼ 2 is shown in Fig. 13. The phase shift between the numerical and experimental results turns out to decrease with time. The angular difference (in absolute value) between the configuration pairs shown in Fig. 13a ranges approximately from 0 to 5:3 whereas the difference in Fig. 13b ranges from 0 to 1 . It is worth noting that in the initial phase of motion the real pendulum is potentially exposed to certain types of disturbances such as out-of-plane motion (due to the manual release of the pendulum end) or torsional vibration (due to the uneliminated twist of the cord). In a further stage, the system motion is observed to stabilize.
In order to examine the importance of the damping effect, the total generalized dissipative force versus time is depicted in Fig. 14 for L ¼ 150 mm and two initial angles. Similarly to uðtÞ, the function Q d ðtÞ decays, and their peak values are larger for greater u 0 . For this example, the maximum absolute values are Q max d % 4:8 Â 10 À6 Nm and Q max d % 10:4 Â 10 À6 Nm when k ¼ 1 and k ¼ 3, respectively. It is worth noting that the generalized gravity forces reach much higher levels-the extreme absolute values are about 1100 and 605 greater than Q max d . Additionally, the particular components of the damping force   At the final stage, our studies are focused on the most effective damping model, but the aim is to define such damping parameters that are independent of the pendulum length, and to find their optimal values. By the term ''overall optimization'', we mean that the distance between numerical solutions and the experimental data is minimized simultaneously for all the values of both u 0 and L. First, we try to express the original (dimensional) coefficients c 1 , c 2 and c a as a product of two factors: affected and non-affected by the leading dimension L of the pendulum. Let us consider the distributions of the linear velocity v(r) and acceleration a(r) along the rod (0 r L). The moments of the frictional forces exerted on the pendulum, related to the three damping terms, are given by ð38Þ where x ¼ _ u and e ¼ € u, whileĉ 1 ,ĉ 2 ,ĉ a are the coefficients (dimensional) independent of L, i.e. common for pendula of various length. It is obvious that M 1 ¼ ÀQ 1 , M 2 ¼ ÀQ 2 , M a ¼ ÀQ a , and all the moments are present in the equation of motion (9). Thus, we can clearly see that Now, we introduce the radius of gyration r O of the rod to write its moment of inertia as I O ¼ mr 2 O . Taking into account the non-dimensional Eq. (12) together with definitions (13), and using the new expressions (40) for c 1 , c 2 and c a , we obtain Mathematical model (41) of the system is just the one in which the parametersĉ 1 ,ĉ 2 andĉ a should not depend on the pendulum length.
Such a reformulation of the dynamic equation allows us to solve the overall optimization problem. More precisely, the gradient-based procedure can be applied simultaneously to the whole set of time series collected experimentally. Thus, Eq. (41) must be integrated numerically at each step, and the new coefficients play the role of the decision variables. The objective function has the form Table 9 The optimal values of the damping coefficients a 1 , a 2 , a a , and the corresponding errors for various pendulum's length. Results obtained with the gradient method for the overall optimization ; and the total mean errors-calculated for all k and L, by analogy to (43)-have the values:ê rms tot ¼ 0:165 rad, e avg tot ¼ 0:0034 rad. Finally, we can return to the dimensional damping coefficients by using relationships (40). The optimization results are presented in Table 9. The total errors for each pendulum length, e rms tot and e avg tot , are given for comparison purposes. As can be seen, these quantities reach higher values than before (compare to Table 8). The distances e rms k and e avg k are also larger almost for every k.
The prediction quality of this approach can be assessed from Fig. 15, where the numerical and experimental responses of the system are shown for L ¼ 100 mm and L ¼ 300 mm at k ¼ 1; 2; 4. This time, the plots are focused on the first 6 s of motion. When it comes to the shortest pendulum, slight phase and amplitude differences can be observed for every initial configuration. In the second case, the amplitude discrepancy is higher, and the phase shift between the angular displacements intensifies with time, a c b Fig. 15 Response of the pendula-numerical solution (solid line) and experimental data (markers) for L ¼ 100 mm (left column) and L ¼ 300 mm (right column): a k ¼ 1, b k ¼ 2, c k ¼ 4. Results obtained with the gradient method for the overall optimization particularly for larger values of u 0 . It seems that the overall optimization leads to a worse match of the model to the experimental data, which can be regarded as the cost of determining just a single set of parameters values (ĉ 1 ;ĉ 2 ;ĉ a ), uniform for every pendulum.
6.3 Some remarks on the models quality The described procedures for estimation of the damping parameters have been based on the Euclidean distance between data points: the experimental and numerical results. Also, the quality of the optimal solutions has been assessed by means of the related RMS and AVG errors. Such an approach seems to be natural and most widely used.
However, the classical distance means a point-topoint comparison of oscillations. It represents a very precise difference between two data sets, calculated sequentially along the time axis. Unfortunately, the general, shape-related similarity of the displacement time histories is not evaluated in this way. Consequently, even a slight phase shift can produce relatively high values of the RMS error.
In order to cast light on the quality of approximation of the real oscillations decay by numerical solutions, we will use the concept of the amplitude envelope [27,37]. More precisely, we assume that the decaying character of the pendulum motion can be approximated by the time-dependent exponential function: where A, b, B are constants. Thus, given the peaks of the system responses found experimentally and numerically, u exp and u num , their upper envelopes are fitted in the least-squares sense: Now, the difference between the curves can be measured in many ways. By analogy to the RMS and AVG errors involved up to now, the following integral forms of the envelope errors are introduced: Table 11 The envelope errors corresponding to the results in Table 9 k L¼ 100 mm L ¼ 150 mm   where t e denotes the end time of motion. It should be emphasized that, in contrast to e rms k and e avg k , the new quantities (47) are non-dimensional, since they represent the relative errors of U num .
The percentage values of e e rms k and e e avg k corresponding to the results discussed in Sects. 6.1 and 6.2 are presented in Tables 10 and 11, respectively. The integrals occurring in formula (47) are evaluated numerically. The total mean errors are specified in the same manner as given in Eq. (26). As can be seen, the RMS and AVG distances between the envelopes for the numerical and experimental results are much lower in the case of the optimization conducted separately for each pendulum length. The e e rms k values lie in the range between 0.4 and 9.1%, but most of them are less than 2%. The interenvelope distance per sample, in turn, does not exceed 1%. For the overall optimization, the ranges for the envelope errors are about two times wider. The e e rms k values lie between 4.7 and 18.6% (many of them exceed 10%), whereas e e avg k \2%.
Thus, the results generally confirm the previous observations. However, it is worth noting that the envelope-based errors do not coincide with the point-to-point ones. Usually the greatest values of e e rms k and e e avg k for some L are related to the smallest values of e rms k and e avg k , and vice versa (e.g. see L ¼ 100 mm, k ¼ 1; 4 in Tables 8 and 10, or L ¼ 300 mm, k ¼ 1; 4 in Tables 9 and 11). These two types of criteria for evaluation of similarity between u exp and u num seem to contradict each other, at least partially. As results from our trials to incorporate the interenvelope distance into the optimization procedures, this kind of objective function could slightly improve the amplitude match between the model and the experimental data, but at the cost of phase difference. The classical point-to-point distance, on the other hand, leads to a compromise between the amplitude and phase match.  Fig. 16 Amplitude envelopes based on the numerical solution (black solid line) and the experiment (markers) for L ¼ 100 mm and k ¼ 4: a optimization conducted separately for each pendulum length, b the overall optimization k ¼ 4 and L ¼ 300 mm, k ¼ 1. Apart from the exponential curves over the whole time range (0 t 12 s), their clipped views for the last 3 s of motion are presented. When dealing with the separate optimization, the difference between envelopes U num and U epx is practically negligible in both cases. The overall optimization, in turn, leads to a worse approximation of the oscillation decay: the interenvelope distance grows gradually with time. In the second example (see Fig. 17b), the difference becomes very large (e e rms k ¼ 18:6%).

Conclusions
In this work, we have focused our attention on planar free vibrations of a single physical pendulum. Such a system has been investigated both experimentally and numerically. The laboratory experiments were performed with pendula of different lengths, for various initial swing angles. It should be emphasized that our studies have not been restricted to the small-vibration case. The system motions were recorded by a high speed camera. In order to approximate the air resistance, three models of damping have been examinedinvolving the three components of the resistive force: linear (proportional to velocity), quadratic (velocitysquared) and acceleration-dependent (proportional to acceleration). A series of numerical simulations has been discussed, in which the corresponding damping coefficients have been estimated by means of several computational methods. The root-mean-square distance and the average distance per sample, evaluated between the numerical and experimental results, have been treated as the error measures.
The simplest, linear damping model is associated with the highest errors between the numerical solutions and the experimental data. The linear-quadratic model leads to non-physical results, i.e. negative values of one damping parameter; the only option is to use the velocity-squared part. The third variant of the resistance force, that is velocity and accelerationdependent, turns out to be the most appropriate. The term proportional to acceleration, that may be tied to a b Fig. 17 Amplitude envelopes based on the numerical solution (black solid line) and the experiment (markers) for L ¼ 300 mm and k ¼ 1: a optimization conducted separately for each pendulum length, b the overall optimization the empirical Morison equation, seems to be indispensable for the successful approximation of the real pendulum motion.
Obviously, there are certain discrepancies between the numerical and experimental results. In particular, the overall optimization leads to significant errors, higher than for the optimization performed separately for each length variant. In some degree, the differences might be attributed to friction at the pendulum support. However, due to a small diameter of the cord and its smooth surface, this effect is assumed to be much weaker than the air resistance. But the capability of the cord to twist and the out-of-plane vibration of the pendulum, especially in the initial phase of motion, may be of much greater importance. And last but not least, all the damping coefficients have been assumed constant, independent of the Reynolds number that is proportional to velocity. All in all, despite such strict assumptions, the considered models allow for quite accurate approximation of the pendulum motion.
When it comes to the computational approaches to the problem, the gradient method for optimization has been chosen as the major one on the path towards determination of a single set of damping parameters, independent of the system's initial configuration. The other techniques, such as the bisection method or the successive search and the intersection technique, are quite well-suited for the simpler (one-and twodimensional) estimation problems.
All the calculations have been conducted using the Wolfram Mathematica software. Many tasks, such as numerical integration of differential equations of motion, the least-squares curve fitting, or numerical evaluaton of integrals, have been accomplished by means of the built-in functions. However, the optimization procedures have been implemented by the authors.
The path of the computational studies reported in this paper, and the presented results indicate that both the modelling of damping and the parameters estimation are not a trivial task, even for relatively simple (but nonlinear) systems. Our future works will be devoted to the problem of this type for the case of three-dimensional motion, and for more complex systems such as multiple pendula.
Funding This study was partially funded by 0612/SBAD/ 3558 and 0612/SBAD/3566 grants from the Ministry of Science and Higher Education in Poland.

Compliance with ethical standards
Conflict of interest The authors declare that they have no Conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.