Axial Static Load Dependence Free Vibration Analysis of Helical Springs Based on the Theory of Spatially Curved Bars

This work addresses an accurate and detailed axial static load dependence linearly elastic free vibration analysis of cylindrical helical springs based on the theory of spatially curved bars and the transfer matrix method. For a continuous system, governing equations comprise coupled vibration modes namely transverse vibrations in two orthogonal planes, torsional and axial vibrations. The axial and shear deformation effects together with the rotatory inertia effects are all considered based on the first order shear deformation theory and their effects on the frequencies are investigated. The effects of the initial stress resultants on the frequencies are also studied. After buckling, forward-shifting phenomenon of higher frequencies is noticeably demonstrated. It is also revealed that a free/forced vibration analysis with an axial static load should not be performed individually without checking buckling loads.

consider the vibration of helical springs with a static axial load.Employing Wittrick's (1966) differential equations Mottershead (1982) extended the theory to large displacements through Lagrange-Green strain equations and derived consistent geometric stiffness matrices.Unlüsoy (1983) and Yalçın (1984) studied the free vibration of helical springs subjected to an axial static load by employing the finite element procedure.In their study the helix axis was modelled as a composition of a number of straight beam elements.Haktanır and Kıral (1992) worked on the free vibration of helical springs subjected to a static axial force with the help of the stiffness matrix formulation and straight beam elements.Tabarrok and Xiong (1992) derived the governing equations for the buckling of spatial rods by considering perturbations about the critical state.Xiong and Tabarrok (1992) summarized a comprehensive energy formulation for the vibration analysis of spatially curved and twisted rods under various applied loads.In addition to the initial axial force, this formulation took into account the initial moments and shear forces as well as initial deformations.Based on this formulation, a spatially curved and twisted rod finite element was developed and used for several examples.Becker and Cleghorn (1992) developed a set of dynamic, partial differential equations in canonical form for helical compression springs for just axial static force, and they solved them for fixed ends and circular cross-section by the transfer matrix method.Chassie et al. (1997) extended the work by Becker and Cleghorn (1992) for the buckling of compression springs to include the additional effect of torsion about the axis of the spring.Chassie et al. (1997) showed that both Mottershead's (1982) and Pearson's (1982) buckling equations do not reduce correctly to the well-known equations for simple rods.The results obtained in the Chassie et al.'s work (1997) also showed a good agreement with those from Tabarrok and Xiong (1992) and Haringx's (1949) theory.They stated that the governing equations derived were in fact a special case of the earlier and more comprehensive analysis of curved and twisted rods under load.Becker et al (2002) broadened their previous works on the buckling analysis for the free vibration analysis of a helical spring subjected to a static axial compressive load by incorporating inertial terms into previous buckling equations.They examined the free vibration analysis numerically based on the transfer matrix method.As they stated this study represents primarily the first detailed examination of the fundamental natural frequencies of helical springs, under static axial load with clamped ends, since the pioneering work by Haringx (1949).Frikha et al. (2011) proposed a physical analysis of the effect of axial load on the propagation of elastic waves in helical beams.The general equations proposed in References (Xiong and Tabarrok, 1992;Becker and Cleghorn, 1992;Chassie et al., 1997), which govern the small perturbations of a helical beam subjected to a static axial load, were employed.An eigen-system was obtained through a Fourier transform along the axis of spring with circular section.Frikha et al. (2011) presented dispersion curves of both loaded and unloaded springs for different parameters such as helix angle, helix index, Poisson coefficient, and axial strain.For high, medium and low-frequency ranges they highlighted the effect of load on wave propagation.A branch identification of dispersion curves was performed for both propagating and non-propagating modes.The effect of loading was found to be significant on the four propagating modes in a low-frequency range.The dispersion curve of the flexural mode oscillating in the normal direction shifted to higher frequencies under tensile loads and vice-versa for compressive loads.Yıldırım (2012) established governing equations of the helical springs with a static axial load based on the Özbek's work (1966).Using the linearized disturbance dynamic equations derived systematically in Yıldırım's study (2012), a few functional works were conducted (Yıldırım, 2009a;İbrikçi et al. 2010).In (Yıldırım, 2009a), the numerical buckling analysis of such springs was performed based on the dynamical approach.Ibrikçi et al. (2010) proposed the use of artificial neural networks (ANN) to predict perfectly the critical buckling loads of cylindrical isotropic and homogeneous helical spring with fixed ends and with circular sections, and also having large pitch angles.Then almost perfect weight values were obtained to predict the non-dimensional buckling loads.The theory in Yıldırım's study (2012) was extended to the free vibration of composite cylindrical helical springs under compression by Kacar andYıldırım (2011, 2016).In Kacar and Yıldırım's study (2016) the stiffness matrix method was used for the solution procedure.For a helical element with constant static axial load, the element stiffness matrix which has six degrees of freedom at each node was obtained with the help of the complementary functions method (Haktanır, 1995).Kacar and Yıldırım (2016) presented natural frequencies and buckling loads of barrel, hyperboloidal and conical unidirectional composite helical springs with circular and rectangular sections.Kobelev (2014) studied the load dependence of transverse vibrations for helical springs by developing equations, which were represented by two coupled second order differential equations.Kobelev's (2014) method was based on the concept of an equivalent column.He derived the explicit formulas for the fundamental natural frequency of the transverse vibrations of the spring depends on the variable length of the spring.The reduction of frequency with the load was demonstrated by Kobelev (2014).In this study, similar to the Reference (Yıldırım, 2009a), it was shown that when the frequency nullifies the spring buckling occurs.It was also shown that the first vibration modes of helical springs correspond to low-frequency motions and they are strongly affected by the presence of applied axial loads.
In the present work, an accurate free vibration frequencies of isotropic and homogeneous cylindrical helical springs with circular section are conducted numerically in a detailed manner with the help of the equations offered by Yıldırım (2012).The transfer matrix method is employed for the numerical analysis.The exact numerical overall transfer matrix is obtained by an effective numerical algorithm developed previously for the springs without initial loads (Yıldırım, 1996).Validity of this algorithm for even the springs with initial static loads is demonstrated.To determine the geometrical configuration of the spring subjected to a static axial force, an analytical expression taking into account for the whole effect of the stress resultants such as axial and shearing forces, bending and torsional moments is used (Yıldırım, 2012(Yıldırım, , 2016)).The present natural frequencies obtained by the method described above are then compared with the available literature.A good harmony is observed with related benchmark studies.The effects of shear/axial deformations and stress resultants on the natural frequencies are investigated.The present work also takes attention to the relationship between buckling loads and natural frequencies.Forward-shifting of frequencies in higher modes is made obvious.It is revealed that both free and forced vibration analyses with an axial static load should be performed after carrying out corresponding buckling analysis to achieve rational results.Yıldırım (2012) derived systematically and consistently the following governing equations in a vector form for a spatial bar to study the static, buckling and vibration problems in the framework of spatial Timoshenko beam theory.

Latin
where In the above, the undeformed cross-sectional area is denoted by A, the moments of inertia with respect to the normal and binormal axes are denoted by n I and b I , respectively.
Here,  is the mass per unit length of the rod, are the radii of gyration with respect to the (t, n, b) axes, respectively.
where  is the density of the rod material.Let  be the helix pitch angle,  be the angular coordinate, and R=(D/2) be the radius of the cylinder.Frenet-Serret relations for cylindrical helical springs are given by (Yıldırım, 1999a)

NON-DIMENSIONAL FREE VIBRATION EQUATIONS OF A CYLINDRICAL HELICAL SPRING SUB JECTED TO A STATIC AXIAL FORCE
where the curvature and tortuosity are as follows The free axial length of the helix is defined by In equation ( 8), the total number of active turns is denoted by n.Frenet components of the equivalent force-couple system at any section of the spring wire due to the initial static axial force and moment vectors will be in the following form of (Figure 2) In order to study the free vibration analysis, the following harmonic solution may be assumed for equations (1) for helical elements, and then using equation ( 9), a set of twelve linear differential scalar equations in Frenet trihedral, which govern the free vibration of a cylindrical helical springs subjected to a static axial load, are achieved as follows (Yıldırım, 2012).
These linear equations given above have constant coefficients and valid for isotropic and homogeneous cylindrical helical springs with uniform sections having double symmetry axes.The free vibration equations presented by Becker et al. (2002) coincide with the equations presented above in the absence of the first terms of Equations ( 11j), ( 11k) and (11l).The first terms of Equations ( 11j), ( 11k) and (11l) represent the shear deformation effects of the static axial load (Yıldırım, 2012).Equations (11) may be expressed in a matrix notation as where the state vector is given by is the dynamic differential matrix.In the present study Equations ( 11) are put in a non-dimensional form by using the followings and Latin American Journal of Solids and Structures 13 (2016) 2852-2875   n and using Equations ( 14) and ( 15), elements of the dimensionless differential matrix to be used in the present study for the free vibration analysis of the spring under a static axial force can be attained.

DETERMINATION OF THE GLOBAL TIP DEFLECTION OF THE SPRING UNDER A STATIC FORCE
In a true free vibration formulation of helical springs it is necessary to implement the true deflected configuration of the spring due to the applied static axial force.Otherwise the numerical free vibration studies cannot give satisfactory results according to the specifications of the helix.So a general formulation of the static deflection under an axial force is required for both static and dynamic analyses.
It may be noted that this is also achieved by a numerical analysis.
Even today the helical spring formulas derived in 1960s still continue to be used in the spring design.These formulas are used in special conditions/sections and do not consider the whole simultaneous effect of the stress resultants such as axial and shearing forces, bending and torsional moments on the deformations.So they cannot called global.The popular deflection formula for closed coiled springs,   10


, proposed by Wahl (1963) is given by Both Equations ( 16) and ( 17) account for just the effect of the torsional moment on the tip deflection.However, it is obvious from Equation ( 9) that when even a single static axial force acts on the center of the helix then the four different types stress resultants such as an axial force, a shearing force, a bending moment and a torsional moment appear at the center of the cross-section of the helix wire.One may readily reach the consequence that there must be contributions of all the remaining stress resultants omitted in the literature on the global tip deflection of the spring.
Under a static axial force, for the determination of the vertical tip deflection of cylindrical helical compression springs with any closed section having double symmetry axes the following equation is proposed by Yıldırım (2012Yıldırım ( , 2016) ) and Haktanır (1994) Definitions of the cross sectional rigidities in Equation ( 18) are given in Equation ( 2).In the above equation, first term represents the effect of the shearing force, the second term is the effect of the axial force, the third term corresponds to the effect of bending moment and finally the last term is the effect of the torsional moment.Equation ( 18) implies that even for springs having very small angles there is a contribution of each stress resultants on the total deflection.
While Equation ( 18), which was obtained based on the Castigliano's first theorem (Haktanır, 1994) , is valid for helical springs having doubly symmetric cross-sections such as solid circle, square, rectangle, ellipse, and hollow circle etc., both Equations ( 16) and ( 17) are valid for just solid circular cross sections having the following properties.

⁄ ;
2 ; A 4 From the last term of Equation ( 18), Equation ( 17) is effortlessly obtained as follows It is very clear that Equation ( 16) is a special case of Equation ( 17) with 1 cos 


. For closedcoiled springs, the effects of both the shearing force and torsional moment become leading whereas the effects of both the axial force and bending moment gain a gradually increasing impact on open coiled springs.Wahl (1963) also offered the following formula to consider both the shearing force and torsional moment effects on the axial tip deflection.
This equation is also achieved by Equation ( 18) , and , it may be revealed that the maximum effect of the axial force is around 1% and of shearing force is around 6%.For large helix angles, while the effect of the torsional moment on the tip deflection decreases to around 65%, the effect of the bending moment increases to around 35%. Consequently, the computation truly of the tip deflection of the spring under static axial load is one of the crucial step of the vibration analysis of such springs.It may be noted that those effects gain considerably importance for especially non-circular cross-sections.
In the following sections of the present study, equations given in (16-18) are all examined to get a strong idea about the free vibration response of cylindrical helical springs with solid circle sections.

NUMERICAL DETERMINATION OF THE OVERALL TRANSFER MATRIX CONSISTING OF A SINGLE HELICAL ELEMENT
As stated above, in the present study, the equations solved numerically based on the transfer matrix method with an effective numerical algorithm previously offered by the author (Yıldırım, 1996).This algorithm allows to model the whole continuous system by a single helical element without confronting any numerical difficulty and any numerical instability even if the effect of the static axial load is considered.Since it corresponds to the analytical series solution of the governing equation, the numerical results are assumed to be exact.In this section the transfer matrix method and the numerical algorithm for the determination of the element transfer matrix will be described briefly below.
In the transfer matrix method, the solution to Equation ( 12) is given by Pestel and Leckie (1963) as follows where S(0) is the state vector at section 0   and ) ( F is the overall transfer matrix.Solution of Equation ( 23) is applicable for the helix having either varying or constant cross-sections.The transfer matrix also satisfies the following differential equation in all cases.
with the initial conditions where I is the unit matrix.Numerical computation of the overall transfer matrix is a crucial step in the transfer matrix method.If it is obtained in an accurate manner, then the results should be acceptable as accurate.The analytical solution of the overall transfer matrix in (25) for constant crosssections is (Pestel and Leckie, 1963).
This solution is generally preferred in the literature due to its simplicity (Pearson, 1982;Becker et al., 2002).However, it is not possible numerically to take enough terms in the exact solution of Equation ( 26).So this leads to consider lots of helical elements to form the whole helix axis.In this case the element transfer matrix is computed for each division.A proper successive multiplication of the element transfer matrices renders the overall transfer matrix (Pestel and Leckie, 1963).
Usually, when applying this series solution it should be necessary to take some numerical preventative measures for especially springs having large helix angles and large active turns.One of this precaution is the application of the Cayley-Hamilton theorem to the series solution (26).Based on the theorem, Equation ( 26) may be put in the form of Latin American Journal of Solids and Structures 13 (2016) 2852-2875 where k( ) are functions of scalar infinite series in  .Utilization of each term of the series k( ) of Equation ( 27) corresponds to twelve terms in (26).The number of terms taken from the infinite series k( ) determines the accuracy of the solution.In the present study, 1000 terms are taken in each k( ) series of Equation ( 27) for each    to calculate the overall transfer matrix.This number of terms corresponds to the 12000 terms in Equation ( 27).The number of terms can be increased without any trouble to increase the accuracy of solution in the numerical algorithm previously developed by the author (Yıldırım, 1996).Letting m be the total number of terms to be taken in Equation ( 27), each  series may be put in the following form of where each element T in  series may be computed from the previous element and for 0  k for any i, When the coefficients qk of the characteristic polynomial given by the following equation are available, ) (i k T terms are calculated recursively and then the transfer matrix is obtained in an accurate manner.
As it is well known, every square matrix satisfies its characteristic equation as shown below The numerical algorithm in (Yıldırım, 1996) is valid for a square matrix D having the following property It may be proved simply that the differential transfer matrix formed by Equation ( 11) satisfies this condition for even cylindrical helical springs with a static axial load.Let us change the order of elements of the state vector as follows In this case the differential matrix may be subdivided into four parts as in the following form of where and In this situation, the coefficients of the odd powers of the differential matrix,  D , in its charac- teristic polynomial becomes zero.Due to the fundamental features of the transfer matrix, it is stated that the matrices D and  D are similar, and so they and their powers have the same trace and determinants.This reveals that the numerical algorithm devised for effectively computation of the helical element transfer matrix in the absence of the static axial force may continue to be used for successfully estimation of the helical element transfer matrix with a static axial load.
For a given axial static load and by using the corresponding helix geometry, after computation of F , the frequency equation can be obtained from the boundary conditions given at both ends ( =0 and  =2n) using Equation (23).For instance, the eigen value equation for fixed-free ends (=0, U=0) is reduced to the following As seen from Equation ( 39), the order of determinant is only six for the spring supported at both ends.As described, the transfer matrix method offers an exact solution with minimal computation memory requirement for dynamic problems if the overall transfer matrix is calculated numerically in an accurate manner.In this study, the free vibration frequencies are obtained by the method of searching determinant.Accordingly, numerical values are attributed to the natural frequency and the corresponding determinant is computed.The values making the determinant zero are assumed to be the natural frequencies of the helix.All numerical computations were performed using the doubleprecision arithmetic.

VERIFICATION OF THE RESULTS WITH/WITHOUT A STATIC AXIAL FORCE
To check the accuracy of the present results a benchmark spring which is not subjected to any static axial force is first considered.The material and geometrical properties together with boundary conditions of the spring are as follows: Lo/D=3.6, /Lo=0, n=7.6, C=D/d=10, d=1mm,  1).From Table (1) it is observed that the natural frequencies obtained by the transfer matrix method show a good harmony.The minor differences between the present study and Reference (Yıldırım, 1999a) probably steam from the number of terms considered in the numerical computation of the overall transfer matrix.The widely held formula for the fundamental frequency in axial mode, axial f (in Hz), in the literature (Haringx, 1949) is given by In the second example a cylindrical helical spring under an axial static load compressed to half of its free axial length is considered.Material and geometric properties of the spring having clamped supports at two ends are: Lo/D=5,  /Lo=0. 5,n=30, C=10, d=1mm In this example, since the helix pitch angle is very small (<<5 o ), the effects of the axial force and bending moment on the tip deflection may be negligible.However, the shearing effect should be considered.
The corresponding axial static load to the prescribed deflection is to be found as Po=8.2532Newton by using Equation ( 18).This value was found as 8.242 by Becker et al. (2002).Now, Let's check whether the applied axial force, 0 P , is equal to critical buckling load, cr P , by using Figure (3).
As shown from this figure N P cr 12  for this example.This verifies that the first frequency will be found is certainly to be the fundamental frequency of the spring loaded by a static axial load.
For the sake of the comparison, by using Equations ( 16) and ( 18), the first 16 natural dimensionless frequencies of the spring considered are presented in Table (2) in a comparative manner with benchmark results.Since the helix pitch angle is small, the additional terms considered in the present study will be gained an importance.As shown from Table (2) those terms affect just symmetric and asymmetric bending modes.As expected, for small helix angles, the analytical theory may give acceptable results.Again a good accordance between Becker et al's (2002) and present results is observed.It may be noted that Becker et al. (2002) computed those frequencies by employing Equation (17).By using Equation ( 18) which considers the whole effect of the resultant force on the tip deflection of the spring, the present first 16 natural frequencies for a cylindrical helical spring considered are tabulated in Table (3) based on the Timoshenko and Bernoulli-Euler hypotheses.As shown from the table, the differences between two theories become visible for 0 0  P .

Mode No
Mode Types (Becker et  ) for a cylindrical helical spring under a static axial load compressed to half of its free axial length.
For this example, effects of the axial and shear deformations on the first 16 natural dimensionless frequencies are also studied and presented in Table (4).In Bernoulli-Euler results of Becker et al (2002), rotational inertia effects are also neglected together with axial and shear deformation effects.As is well known, those effects will become more clear for small spring indices, C<10 and some types of cross sections such as rectangular shapes.

NUMERICAL EXAMPLE FOR BUCKLING AND FREE VIBRATION WITH A STATIC AXIAL FORCE
As it is well known, the critical buckling loads may also be determined in a dynamic manner (Yıldırım, 2012;Kobelev, 2014).In this method, for a given static axial load, from 0 P to cr P at which the fundamental frequency of the spring becomes zero, the fundamental frequencies of the spring are searched.
The effect of especially bending moment on the tip deflection increases with increasing helix pitch angle.Increasing the helix pitch angle results in an increase in Lo/D ratios while the number of active turns remain constant.To see the effect of bending moment on both the critical buckling load and the free vibration frequencies of a compressed spring, here a cylindrical helical spring having the same properties of the spring in the previous example, except n=5 and Lo/D=10, is studied., it is stated that while the percent contribution of the bending moment on the total tip deflection is about 24%, this value is about 76% for the torsional moment.0  cr  For this example present free vibration frequencies and critical buckling load which are both obtained by using Equation ( 18) are presented in From the Table (6) it is shown that after determining the critical buckling load, the fundamental frequency is lost and replaced by the lowest second frequency.This means that buckling occurred.For the spring with a static axial force, natural frequencies higher than the first are all shifted forward.The numerical free vibration analysis continues to produce natural frequencies for the post buckling case..43 1442.24 1442.23 1442.31 1442.34 1442.35 1442.35 1936.80 1936.82 1936.84 14 1886.91 1909.39 1921.25 1933.56 1936.08 1936.59 1936.79As stated above, Equation ( 17), consider the deflection due to just torsional moment for large helix angles.Using this equation, the present free vibration frequencies and critical buckling load are also computed and presented in Table (7).In this case natural frequencies are found smaller than frequencies given in Table (6).From table (7), it is again observed the forward-shifting phenomenon of the frequencies in higher modes.As it is well known, the common equation numbered by ( 16) consider the deflection due to the just torsional moment for small helix angles.Finally, to see the effect of Equation ( 16), Let's consider Table (8).In this situation, natural frequencies have the numerical values between Table (6) and Table (7).Figure (5) shows the critical buckling loads obtained in a static manner for this example.The numerical values of the critical buckling loads calculated by both static and dynamic manners are exhibited in Table ( 9) in a concise way.From this table, it may be concluded that critical buckling loads obtained by using Equation ( 18) gives the higher values.This means that critical buckling loads found by Equations ( 16) and ( 17) fell in a safe region in design for this example.(Yıldırım, 2009b

CONCLUSIONS AND DISCUSSIONS
Vibration analysis of the springs under axial static loads requires a precise study to understand the real behavior of the springs.To accomplish such an accurate analysis first the most comprehensive differential equations which comprise some additional terms for taking into account for the axial and shear deformation terms on the buckling and vibration were used, second the exact total tip deformation for a spring subjected to an axial force is calculated by using a global analytical formula considering the whole effect of the resultant forces on the tip deflection, third the overall exact transfer matrix is numerically calculated by using an effective numerical algorithm.As a result, those choices render an accurate numerical free vibration analysis of cylindrical helical springs under a static axial force based on the transfer matrix method.The vibration modes correspond to low-frequency motions are strongly affected by the presence of applied static axial loads as Kobelev stated (2014).In the present study the forward-shifting phenomenon of the natural frequencies in higher modes are observed after buckling.That is with a small increment of the value of the applied static load, one may get the next natural frequency after buckling occurs.Before determining the natural frequencies subjected to a certain static axial force, the axial force considered should be tested for the critical buckling load.The author offers that a free and/or forced vibration analysis with an axial static load should not be performed individually without checking buckling loads.
American Journal of Solids and Structures 13 (2016) 2852-2875 2 GOVERNING EQUATIONS OF A SPATIAL BAR SUBJECTED TO STATIC AXIAL LOADS Consider a curved spatial bar having the curvilinear position coordinate s, and let t be the time.Let vector in Frenet coordinates (t,n,b) (Figure 1).Initial internal static force and moment vectors are denoted by ) and moment vectors are illustrated by p(s,t) and m(s,t), respectively.

Figure 2 :
Figure 2: Geometry of a helix subjected to an axial compressive force and a torque.
1963) corrected formula concurs with the present closed coil formulas by pi.The first terms of Fixed-Fixed ends.The results are presented in Table (

Figure 4 :
Figure 4: Relation between frequencies and static axial load for the same spring.

Table 2 :
The first 16 natural dimensionless frequencies (

Table 3 :
The present first 16 natural frequencies for a cylindrical helical spring based on the Timoshenko and Bernoulli-Euler hypotheses by using Equation (18).

Table 4 :
Effects of the axial and shear deformations on the first 16 natural dimensionless frequencies (Po=8.242N).
This may be done by numerically as shown in Table(5).

Table 5 :
Percent effects of the each stress resultants on the relative compression for some applied axial forces.

Table 6 :
Free vibration frequencies and critical buckling load which are both obtained by using Equation (18).

Table 7 :
Free vibration frequencies and critical buckling load which are both obtained by using Equation (17).

Table 8 :
Free vibration frequencies and critical buckling load which are both obtained by using Equation (16).

Table 9 :
). Critical buckling loads and corresponding deformations of the spring.