Study on TVD parameters sensitivity of a crankshaft using multiple scale and state space method considering quadratic and cubic non-linearities

Talebitooti, R.; Morovati, M.

doi:10.1590/S1679-78252014001400007

Abstract

In this paper the effect of quadratic and cubic non-linearities of the system consisting of the crankshaft and torsional vibration damper (TVD) is taken into account. TVD consists of non-linear elastomer material used for controlling the torsional vibration of crankshaft. The method of multiple scales is used to solve the governing equations of the system. Meanwhile, the frequency response of the system for both harmonic and sub-harmonic resonances is extracted. In addition, the effects of detuning parameters and other dimensionless parameters for a case of harmonic resonance are investigated. Moreover, the external forces including both inertia and gas forces are simultaneously applied into the model. Finally, in order to study the effectiveness of the parameters, the dimensionless governing equations of the system are solved, considering the state space method. Then, the effects of the torsional damper as well as all corresponding parameters of the system are discussed.

Crankshaft; non-linear vibrations; torsional vibration damper (TVD); multiple scales method; state space method

Study on TVD parameters sensitivity of a crankshaft using multiple scale and state space method considering quadratic and cubic non-linearities

R. Talebitooti^I; M. Morovati^II

^IAsst. Prof., School of Mechanical Engineering, Iran University of Science and Technology, Tehran, Iran. e-mail: rtalebi@iust.ac.ir

^IIM.Sc. graduated, School of Automotive Engineering, Iran University of Science and Technology, Tehran, Iran. e-mail: mehdimorovati@auto.iust.ac.ir

ABSTRACT

In this paper the effect of quadratic and cubic non-linearities of the system consisting of the crankshaft and torsional vibration damper (TVD) is taken into account. TVD consists of non-linear elastomer material used for controlling the torsional vibration of crankshaft. The method of multiple scales is used to solve the governing equations of the system. Meanwhile, the frequency response of the system for both harmonic and sub-harmonic resonances is extracted.Â In addition, the effects of detuning parameters and other dimensionless parameters for a case of harmonic resonance are investigated. Moreover, the external forces including both inertia and gas forces are simultaneously applied into the model. Finally, in order to study the effectiveness of the parameters, the dimensionless governing equations of the system are solved, considering the state space method. Then, the effects of the torsional damper as well as all corresponding parameters of the system are discussed.

Keywords: Crankshaft; non-linear vibrations; torsional vibration damper (TVD); multiple scales method; state space method

Nomenclature

b Crankshaft linear damping coefficient

c TVD linear damping coefficient

Coefficient of excitation

I_c Crankshaft polar mass moment of inertia

I_d TVD polar mass moment of inertia

k_c Crankshaft linear spring sitffnes

k_d TVD linear spring sitffnes

k₂ Crankshaft non-linear quadratic non-linearities

k₃ TVD non-linear cubic non-linearities

R Dissipated energy

t Dimensionless time

T Kinetic energy

V Potential energy

General coordinate system

ω₁ System natural frequency

ε Small dimensionless parameter

Ω₁ External excitation frequency

Ω₂ External excitation frequency

σ₁ Detuning parameter

σ₂ Detuning parameter

σ₃ Detuning parameter

1 INTRODUCTION

The Vibration phenomenon is one of the most important issues should be treated to reduce the unpleasant shaking in various compartments. If the movement is left untreated, many consequent problems such as noise transmission into cavity of the car as well as fracture and failure of the compartments will be occurred. The crankshaft is one of those compartments should be attentioned. It is mostly manufactured with high mechanical strength cast iron. It should be strong enough to tolerate piston strokes without high torsion. In addition, it should be balanced very carefully to prevent vibrations generated from out of center weight of the crank.

Crankshaft has been investigated by many researchers. Espindola et al. (2010) indicated that a hysteretic model can be derived from a viscoelastic material based on four fractional parameters. Moreover, they derived generalized quantities of ordinary and pendulum type absorbers considering both viscoelastic and hysteretic materials. They have also compared the performance of the system with absorbers of viscoelastic and hysteretic nature. Mourelatos (2001) introduced a model to analyze the dynamic behavior of an internal combustion engine crankshaft. The model couples the crankshaft structure dynamics, the main bearing hydrodynamic lubrication and the engine block stiffness using a system approach. Smaili and Khetawat (1994) investigated the vibratory behavior of an automotive crankshaft. Their FEM model involved a new scheme for modeling the stiffness and damping properties of the journal bearings. Mourelatos (2000) developed a structural analysis using dynamic sub-structuring with Ritz vector to predict the dynamic response of an engine crankshaft, with the aid of FEM. Asfar (1992) investigated the effect of non-linearities in elastomeric material dampers to isolate torsional oscillations of internal combustion engines shafts. Boysal and Rahnejat (1997) have studied a detailed multi-body nonlinear dynamic model of a single cylinder internal combustion engine. Their model comprises all rigid body internal members, support bearings, joints, couplers and connections between the various engine components, as well as vibration dampers. Montazersadegh and Fatemi (2007) studied a dynamic simulation on a crankshaft from a single cylinder four stroke engine. They also performed finite element analysis to obtain the various stress magnitude at critical locations. Murawski (2004) investigated additional bending stresses in the crankshaft as well as possible vibrations of the ship's structure due to the reaction force in the thrust bearings. Giakoumis et al. (2008) evaluated the crankshaft angular deformations during turbocharged diesel engine operation owing to the difference between instantaneous engine and load torques by the aid of an experimentally validated diesel engine simulation code.

Furthermore, many authors is applied multiple scales method to solve various kinds of partial differential equations (PDEs). Multiple Scales method is one of perturbation method branches which is aimed at finding approximate analytic solutions to problems whose exact analytic solutions cannot be found. The setting where perturbation methods are applicable is where there is a family of equations P(ε), depending on a parameter ε « 1, and where P(0) has a known solution. Perturbation methods are designed to construct solutions to P(ε) by adding small corrections to known solutions of P(0). The singular aim of perturbation methods is to calculate corrections to solutions of P(0). Perturbation methods do not seek to prove that a solution of P(0), with corrections added, is close to a solution of P(ε) for e in some finite range with respect to some measure of error. Its sole aim is to compute corrections and to make sure that the first correction is small with respect to the chosen solution of P(0), that the second correction is small with respect to the first correction and so on, all in the limit when ε approaches zero. R. Ghaderi and Azin Nejat (2014) is used multiple scales method to analyze the frequency response of Nano-Mechanical Cantilever (NMC). Then, they applied the primary resonance excitations to show the softening phenomenon in frequency response. Moreover, Eissa and Bassiouny (2003) is applied the method of multiple scales to construct a second order uniform expansion of the non-linear rolling response of a ship in regular beam seas.

In most of analytical studies surveyed above, only one term is used to simulate the external force applied on crankshaft. The other external forces are neglected which may lead into improper results. In most of the researches, authors have not considered both inertia and gas forces in their model. These two terms are of high importance. In other word, if these are not modeled simultaneously could influence on crankshaft torsional vibration. In addition, the corresponding parameters of crankshaft and TVD are obtained with the aid of the FEM model developed in ABAQUS software and curve fitting method applied in a numerical procedure developed in MATLAB.

In this paper, nonlinear vibration of a crankshaft is studied using the multiple scales perturbation technique (MSPT) at harmonic and sub-harmonic resonance. The system consists of the main structure along with torsional vibration dampers. The external forces including both inertia and gas forces are simultaneously applied into the model. In order to study the effective parameters, the dimensionless governing equations of the system are solved, considering the method of state space in a steady state process. Then, the effects of the torsional damper as well as all corresponding parameters of the system are discussed.

2 GOVERNING EQUATIONS OF THE SYSTEM

The governing equations of the system including the crankshaft and elastomer have been considered. Therefore, in order to derive these equations a mathematical model is represented in this section, using rod and disk elements.

Lagrangian formula is used to calculate differential equations followed as in Rao (2010):

where L is representing the difference between kinetic and potential energy, R is dissipated energy, Q_i is external forces and q_i is general coordinate system denoted as and .With substituting Eqs. (2,3) into Eq. (1) the governing equations of the system are followed as:

where k_c and k_d are linear spring stiffness of the crankshaft and TVD respectively; I_c and I_d are polar mass moment of inertia of the crankshaft and TVD, respectively; k₂ and k₃ are non-linear quadratic and cubic Non-linearities of the crankshaft and TVD; b and c are linear damping coefficients of the crankshaft and TVD and and are the Amplitudes of the excitations. Also, in Eqs. (4,5) it is assumed that .

In order to get involve with standard equations, non-dimensional variables are introduced as follows:

where t indicates the dimensionless time and ω₁ is the natural frequency of the system including crankshaft and TVD.

Thus, Eqs. (4) and (5) are written as:

where ε is a small dimensionless parameters and (.) represents the derivative with respect to t^* and the other parameters are introduced as follows:

3 SOLVING THE GOVERNING EQUATIONS

In this section the method of multiple scales proposed by Nayfeh and Mook (1995) is used to solve Eqs. (7,8). The approximate solution is represented followed as:

where ε << 1 , T₀ and T₁ are the fast and slow time scales defined as:

In addition the time derivative becomes:

where D_n = ∂/∂T_n.

Substituting Eqs. (10), (11) and (13) into Eqs. (7) and (8) and reconstructing the equation with respect to the power of ε will result into two sets of equations:

The equations involving the zero order ε⁰ can be written as:

The equations involving the first order ε¹ can be written as:

It should be noted that the two sets of equations represented above are coupled. In other hand, solution of Eq. (17) depends on the solution of Eqs. (14-16). The similar procedures are applied in solving to Eqs. (15) and (16).

where A₁ is a function of T₁ at this level of approximation. Substituting Eq. (18) in Eq. (15) yields the solution of θ₂₀ as:

where A₃ = β²A₁/(β² - 1) and cc stands for the complex conjugate of the preceding terms and A₂ is an unknown function at this level of approximation. Substituting Eqs. (18) and (19) into Eq. (16) yields:

where NST stands for terms that do not produce secular terms. Therefore, any particular solution of Eq. (20) contains secular or small divisor terms depending on the resonance conditions. The detuning parameters σ₁, σ₂ and σ₃ are introduced as follows:

Substituting Eq. (21) in Eq. (20) and eliminating the terms that produce secular terms and small divisors in θ₁₁ yields the following expression:

where the prime denotes the derivative with respect to T₁ and also the overbar shows the complex conjugate. The uniform solution of Eq. (20) can now be written as:

In Eqs. (23), the terms proportional to exp(2iT₀), exp(2iβT₀), etc are ignored, as these terms are neutral in resonance cases.

Substituting Eq. (21) in Eq. (24) and eliminating the terms that produce secular terms and small divisors in θ₂₁ yields the following equation as:

Now, it is convenient to introduce A_n in polar form as:

where a_n and ψ_n are real components. Substituting Eq.(26) in Eqs. (22) and (25) and separating real and imaginary parts yields the governing following equations:

where:

Eqs. (27-30) are first order coupled differential equations which have to be solved simultaneously. Two cases of resonances including harmonic and sub-harmonic are studied in next sections.

4 HARMONIC RESONANCE

Considering the detuning parameters as β ≅ 2 , Ω₁≅ 1 and also taking the value of Ω₂ far from 2, Eqs. (27-30) are written as:

Steady state solutions of the system are obtained with considering in Eqs. (32-35) simultaneously. From Eq. (35) it can be obtained that ζ₂ = ζ₁Γ. Using steady state condition in

Eq. (31) results in:

Thus, the steady state equations in this case are followed as:

Using Eqs. (37-39), the frequency response equation of the system is constructed as:

5 SUB-HARMONIC RESONANCE

With assuming the detuning parameters as β ≅ 2 , Ω₂≅ 1 and also taking the value of Ω₁ far from 1, Eqs. (27-30) are reconstructed as:

Assuming the steady state conditions, Eqs. (41-43) are represented as:

Using the above equations (44-46), the frequency response equation of the system is followed as:

6 NUMERICAL METHODS

In order to numerically solve the system equations, they are written in non-dimensional state space forms. Therefore the four variables are defined as:

The above equation can be written as follows:

Substituting equation (48) in dimensionless equations of system (7,8) results in below equations:

6.1 Linear torsional spring stiffness (k_c) and non-linearity (k₂) of the crankshaft calculation

In order to estimate the linear and nonlinear quadratic terms of torsional spring stiffness, a torsion load is applied into the crankshaft. The, crankshaft specifications are shown in Table. 1. Firstly, the CAE model is built; then the FE model is developed with a high precision mesh in HyperMesh software. This is due to the fact that solving problems using finite element with hexahedral meshes are faster and more accurate than tetrahedral meshes.Â In addition, sometimes in finite element analysis, the usage of tetrahedral meshes for complex geometries will result in unreliable results. Therefore, the hexahedral meshes are devoted to this complex geometry of crankshaft. Then, the model is exported into ABAQUS software. The cast iron material with the properties as depicted in Fig. 2 is specified into the model. The crankshaft is cantilevered at the end and a 15 N.m torque is applied at the other side of crankshaft where the flywheel is placed. The main goal to apply the torque is to obtain the corresponding parameters of the crankshaft and TVD which is substantially independent to input applying torque. Therefore, the amount of input torque does not influence on the results of these corresponding parameters. Fig. 3 Shows isometric view of mesh elements of the crankshaft.

Thumbnail

Nonlinear analysis is chosen to model the stiffness factors. Therefore, the torque applied at the end of crankshaft; then it is plotted versus torsional angles as depicted in Fig. 4. It should be noted that the applications of numerical techniques in engineering often involve curve fitting of experimental data described by Mathews and Fink (1999). Therefore, a polynomials curve k_cx + k₂x² is fitted into the simulated data. Consequently, the coefficients of stiffness factors could be extracted.

The value of k_c, k₂ and their corresponding errors are listed in Table 2. The amounts of errors listed in the Table are negligible which approve the values of stiffness parameters.

Thumbnail

6.2 The torsional stiffness parameters of TVD model

Torsional vibration damper is applied to reduce the vibration in crankshaft. The TVD is placed at one side of crankshaft. It is known as a pulley placed on the crankshaft. It includes the different types of hub, rubber and ring components, as shown in Fig. (5).

Also, TVD components specifications are listed in Table 3.

Thumbnail

Then, considering TVD components into the model the analysis are performed similar as previous analysis discussed previously. There is a difference here and that the EPDM (Ethylene Propylene Diene Monomer) as depicted in Fig. 6. is specified for the rubber as well as cast iron which is specified for the hub and ring components. The values of the torque applied into the system are plotted versus angle of crankshaft and then the stiffness factors are extracted in the same way as described before (Fig. 7).

6.3 Other component specifications

With meshing crankshaft and TVD the components in HyperMesh software and then devoting the material properties for components in ABAQUS software and, some key component specifications could be obtained as listed in Table 5.

Thumbnail

The values of k_d, k₃ and their corresponding errors are shown in Table 4. The errors reported here are also negligible.

Thumbnail

It is obvious from Table 5 that the amount of TVD damping coefficient is higher than crankshaft damping coefficient.

6.4 Excitation force

Mainly, the crankshaft vibrations originate among two principal sources from internal combustion engine. These sources contain the cylinder pressure as a first source. These also include the mass inertia of components caused from reciprocating movement of the piston. These two factors make the excitation function be more complicated. Fig. 8 shows the pressure excitation versus crankshaft angle with different terms. In addition, Fig. 9 shows the mass excitation moment of the engine components which are significantly influenced by engine rotation.

7 NUMERICAL RESULTS AND DISCUSSION

In some cases for designers it is so important to know the effect of varying design parameters on behavior of the vibration system. Therefore, in this section the effect of changing these parameters on vibration behavior of the system is presented. It should be noted that, dimensionless parameters introduced depend on one or more physical parameters of the crankshaft vibration behavior. For instance, with improving the TVD linear damping coefficient, will lead into increasing the dimensionless parameter ζ₁. In addition, enhancing the crankshaft linear damping coefficient will result in increasing the dimensionless parameter ζ. Moreover, with increasing the TVD linear spring stiffness the dimensionless parameter γ₂ would be increased in this application. In a case of harmonic resonance from Eqs. (37-39) it should be considered the variation of n a versus the detuning parameter σ₂. Fig. (10) shows the variation of a_n versus σ₂ for real parameters value. The bold lines represent crankshaft vibration amplitude versus detuning parameters in both unstable (curve) and stable (straight line) form. Also, the dotted lines represent TVD vibration amplitude versus detuning parameter in both unstable (curve) and stable (straight line) form.

Figure. (11) Shows the effect of increasing dimensionless parameter ζ₁ on response of the main system when it is 8 times greater than its real value. It is clear that increasing this dimensionless parameter leads to decreasing the curve inclination of crankshaft vibration and then becomes multi-valued at zone of steady state. In other words, the system status goes from unstable to stable during the long time.

Figure (12) Shows the effect of increasing dimensionless parameter ζ on response of the main system when it is 15 times greater than its real value. It is clear that in this case by increasing this dimensionless parameter, the system tending to show more stable vibration than unstable vibration because of the value of detuning parameter. Also, by increasing this dimensionless parameter, the crankshaft becomes completely multi-valued at zone of steady state but TVD becomes single-valued at steady state zone.

Figure. (13) Shows the effect of increasing dimensionless parameter γ₂ on response of the main system when it is 6 times greater than its real value. It is clear that in this case the system becomes stable during the minimum time. Also, by increasing this dimensionless parameter, the crankshaft amplitude at zone of transient is decreased and both crankshaft and TVD have singlevalued at zone of steady state.

Figure (14) Shows the effect of increasing dimensionless parameter b on response of the main system when it is 6 times greater than its real value. It is clear that the amplitude of vibration at transient zone is decreased and crankshaft has single-valued at zone of steady state. Also, the TVD has tending to have single-valued at zone of steady state.

Figure (15) Shows the effect of decreasing dimensionless parameter Γ on response of the main system when Γ = -0.33 . In this case the behavior of the system is completely changed and both crankshaft and TVD shows divergence behavior. In addition, the system tends to be switched from steady (both dotted and bold straight) to unsteady (both dotted and straight curved) form.

The non-dimensional equations of the system have been solved numerically in section 6. Now, in order to investigate the sensitivity of designing parameters of the system, the influence of parameter variation on system behavior of crankshaft and TVD will be simulated graphically in a harmonic resonance case.

Firstly, the ratio of I₂/I₁ is assumed to be 0.1 of the case specified before. According to Fig.(¹⁶ ,¹⁷ ), the crankshaft and TVD systems show different vibration behavior. In other word, the vibration amplitude of crankshaft is smaller than that of TVD system. This manner is a desirable case for dissipating the vibration energy, as the purpose of using TVD is to reduce vibration amplitude of the system.

In Figs. (¹⁸ ,¹⁹ ), linear coefficients of crankshaft and TVD are 10 times of the reference case. It is well observed that, while the Non-dimensional time becomes 100, the amplitude of crankshaft and TVD with different vibration behavior approach into steady manner. Thus, two different zones in both diagrams called transient zone (Non-dimensional time is a value between 0 to 100) and steady state zone (Non-dimensional time is a value between 100 to 1000) are observed. In addition, the amplitude of TVD in steady state zone is a little lower than that of the crankshaft. However, with increasing linear coefficients of crankshaft and TVD to 10 times of their initial values, the amplitude of crankshaft are considerably increased which leads into undesirable vibration.

Figures (²⁰ ,²¹ ) shows this negative effect of linear crankshaft damping on system vibration behavior. It should be noted that, when linear crankshaft damping is assumed negative about 25% from its real value, both crankshaft and TVD are tending to be unstable. In other words, the system has no desirable vibration behavior at all.

Figures (²² ,²³ ) show the effect of neglecting damping ratio on vibration of the system. In this case, the crankshaft linear damping coefficient is reduced to 0.001. It is well depicted from these figures that both crankshaft and TVD behave in the same manner. In other word, their amplitude decreases during the non-dimensional time. Therefore, as a result of reduction of the system amplitude during the time and the tendency of the vibration amplitude of the system to zero, this reduction of crankshaft linear damping coefficient is desirable.

As depicted in Figs. (²⁴ ,²⁵ ), the effects of non-linear cubic TVD coefficients on vibration of the crankshaft and TVD are investigated .While non-linear cubic TVD coefficient is reduced to 0.001, crankshaft and TVD show the same vibration behavior. It is clearly observed that the system vibration amplitude is considerably high. In other word, the vibration of the system is not desirable although the system having the steady state vibration during the broad band time.

8 CONCLUSION

In this paper the vibration behavior of a system consisting of crankshaft and TVD has been considered. This system is described with second order non-linear differential equations. The method of multiple scales is applied to study the control of a combustion engine crankshaft vibration using a non-linear elastomeric material vibration damper under the interaction of external excitations originated from different sources. Also, the numerical solution is applied to solve the nondimensional equations of the system. Practically the following results are reported as:

The damping coefficient of the crankshaft could greatly influence on system behavior. The small or negative damping factor, leads into the worst behavior of the system, as it causes larger steady-states amplitudes or instability for both crankshaft and TVD.
Large magnitude of TVD non-linearities reduces TVD effectiveness.
While the ratio of (I₂/I₁) reduced to 0.1, the system behavior became steady. Therefore, this ratio is suggested in design process, if there are no other significant factors of limitation.
Linear Damping Coefficient Reduction leads into reduction of amplitude vibration in both crankshaft and TVD. This parameter is also very essential in designing the TVD for reducing the vibration of the whole system.
With increasing b1, the amplitude of vibration at transient zone becomes smaller. In addition, the crankshaft at steady state zone is single-valued.
As the magnitude of Γ has been taken the negative value, the system are going to be diverged.
With increasing the γ₂, the crankshaft amplitude is decreased in transient zone. However, both crankshaft and TVD have single-value in steady state zone.
With increasing the ζ, the crankshaft becomes multi-valued in steady state zone. However, the TVD becomes single-valued in steady state zone.
While Non-Linear Cubic TVD Coefficient is reduced, both TVD and crankshaft show the same steady state behavior. However, the vibration amplitude of the system is considerably undesirable.

Acknowledgments

The authors would like to express their gratitude to the respected referees for many valuable comments which improve the presentation of this paper.

Received 21.01.2014

In revised form 04.08.2014

Accepted 04.08.2014

Available online 17.08.2014

Asfar, K. R., (1992). Effect of Non-Linearities in Elastomeric Material Dampers on Torsional Vibration Control. International Journal of Non-Linear Mechanics 27: 947-954.
Bouchart, V., Bhatnagar, N., Brieu, M., Ghosh, A. K., Kondo, D., (2008). Study of EPDM/PP polymeric blends: mechanical behavior and effects of compatibilization. Comptes Rendus Mecanique 336:714-721.
Boysal, A., Rahnejat, H., (1997). Torsional Vibration Analysis of a Multi-Body Single Cylinder Internal Combustion Engine Model. Applied Mathematical Modeling 21: 481-493.
Eissa, M., El-Bassiouny, A. F., (2003). Analytical and numerical solutions of anon-linear ship rolling motion. Applied Mathematics and Computation 134: 243-270.
Esmaili, A. A., Khetawat, M. P., (1994). Dynamic Modeling of Automotive Engine Crankshaft. Mech. Mach. Theory 29: 995-1006.
Espindola, J., Bavastri, A., Lopes, M. O., (2010). On the Passive Control of Vibration with Viscoelastic Dynamic Absorbers of Ordinary and Pendulum Types. Journal of Franklin Institute 347: 102-115.
Ghaderi, R., Azin Nejat, A., (2014). Nonlinear mathematical modeling of vibrating motion of nanomechanical cantilever active probe. Latin American Journal of Solids and Structures 11: 369-385.
Giakoumis, E. G., Rakopoulos, C. D., Dimaratos, A. M., (2008). Study of Crankshaft Torsional Deformation under Steady-State and Transient Operation of Turbocharged Diesel Engines. International Mechanic Engineering 222: 17-30.
Griswold, F. D., Stephens, Jr. R. I., (1987). Comparison of fatigue properties of nodular cast iron production and Y-block casting. Int J Fatigue 413: 3-10.
Mathews, J.H., Fink, K,D., (1999). Numerical Methods Using MATLAB, Third Edition, Prentice Hall.
Montazersadegh, F. H., Fatemi, A., (2007). Dynamic Load and Stress Analysis of a Crankshaft. SAE International, 01-0258.
Mourelatos, Z. P., (2000). An Efficient Crankshaft Dynamic Analysis Using Substructuring with Ritz Vectors. Journal of Sound and Vibration 238: 495-527.
Mourelatos, Z. P., (2001). A Crankshaft System Model for Structural Dynamic Analysis of Internal Combustion Engine. Computers and Structures 79: 2009-2027.
Murawski, L., (2004). Axial Vibration of a Propulsion System Taking into Account the Couplings and the Boundary Conditions. Journal of Marine Science and Technology 9:171-181.
Nayfeh, A. H., Mook, D. T., (1995). Non-Linear Oscillations, John Wiley & Sons.
Rao, S., (2010). Mechanical Vibrations, Prentice Hall.