The General Dispersion Relation for the Vibration Modes of Helical Springs

Abstract

A system of mathematical equations was developed for the calculation of the natural frequencies of helical springs, its predictions being compared with finite element simulation with ANSYS^®. Authors derive the general equations governing the helical spring vibration relative to the Frenet trihedral representing the normal, binormal and tangent unit vectors to the spring medium line. The dispersion relation $\omega =f\left(k\right)$ has been obtained to model a wave traveling along the axis of the wire.

1. Introduction

Helical springs are one of the most frequently used elastic elements in mechanical engineering. They are used in the most diverse structures as elastic energy accumulators supplementing damping devices, e.g., automobile or railroad suspensions, and in advance and return devices, e.g., camshafts and valves of internal combustion engines (Wahl [1], Shigley [2]; Kobelev [3]). One of the spring failure modes is caused by resonance vibration that occurs when the spring is excited with a periodic signal of frequency equal to its natural frequency.

Love [4] develops equations to study the static response of helical springs subjected to large deformations. Stokes [5] investigates vibrations in the case of a spring subjected to shock loads. Gironnet and Louradour [6] determine more precise expressions for the calculation of the natural frequencies of the helical spring. There are several numerical methods to determine the natural frequencies, the main ones being: the transfer matrix method (Yildirim [7]), the dynamic stiffness formulation (Lee [8]), and the pseudo-spectral method (Lee [9]). Becker [10] determines with the transfer matrix method the resonant frequencies of a spring subjected to an axial compression load. Jiang [11] studies forced vibrations and wave propagation in springs using the Laplace transform.

Wahl [1] already proposed in the middle of the 20th century the equations that are widely used today for the design of springs. Among these equations, Equation (1), which provides the stiffness constant of the spring subjected to an axial force, and Equation (2), which provides the natural frequency of the spring placed between two parallel flat plates, stand out.

$$k_{Wahl}=\frac{G{d}^4}{8D^{3}N}$$

(1)

where k_Wahl is the spring axial stiffness, G is the shear modulus of elasticity, d is the wire diameter, D is the helix diameter and N is the number of active coils.

Where f is the natural harmonic frequencies of a spring in Hz, and W is the mass of the spring in kg. The fundamental frequency is determined for $m=1$ and is usually the most important frequency in practice. Replacing in Equation (2) the value of the axial stiffness we obtain Equation (3) also called Harignx’s equation [12], who experimentally verified its validity.

$$f=\frac{m}{2}\sqrt{\frac{k_{Wahl}}{W}} m=1,2,3\dots .$$

(2)

$$f=\frac{d}{2\pi R^{2}N}\sqrt{\frac{G}{32\rho }}$$

(3)

R being the helix radius and ρ the density of the material.

Equation (3) indicates that for a given material, the fundamental frequency of a helical spring is proportional to the wire diameter and inversely proportional to the product of the helix diameter and the number of active coils, i.e., the resonance frequency of the spring depends on all the spring design parameters. If a higher frequency is desired for the same spring diameter, it is sufficient to increase the wire diameter (this increases the spring stiffness), and vice versa. If you want to vary the resonance frequency while keeping the wire diameter constant, just change the spring diameter or the number of coils of the spring. All the expressions shown above are raised with respect to an inertial reference system whose main X axis coincides with the symmetry axis of the helical spring.

By design of the machine itself, the spring is always inserted between two masses in relative oscillatory motion and its function is to keep them away from each other ensuring at all times that they do not come into contact. The design of the helical spring depends directly on the inertia between the masses at any given moment. In addition to the diameter of the wire and the spring itself, the number of coils and the pitch are responsible for achieving non-interaction between the mechanical elements in relative motion.

It is interesting to note that the behavior of the spring changes radically when its operation is associated with a mass much greater than that of the spring itself. For example, in the damping systems of a railway bogie such as the one in Figure 1. In that case, the mass of the spring is much lower than the mass of the railway car it is supporting and Equations (2) and (3) are no longer valid. Den Hartog [13] shows that in that configuration the natural frequency follows Equation (4).

$$f=\frac{1}{2\pi }\sqrt{\frac{g}{{\delta }_{est}}}$$

(4)

where g is the gravitational constant and δ_est is the static deflection of the spring under the weight of the suspended load.

Sauvage [15] and Campedelli [16] studied in detail the behavior of the springs of a bogie and simulated in finite elements the anomaly presented by the springs during axial compression called “effort de chasse”. This anomaly consists of a rotation and a lateral displacement appearing on the bearing surface of the spring. The “effort de chasse”, also called transverse flexibility of the spring, is caused by the forces and moments that appear at the ends of the spring during its axial compression. Campedelli [16] determines that the transverse flexibility is determined by both the geometry of the spring end and the way in which the gap between the coils closes during axial compression. The “effort de chasse” is an interesting anomaly in railway suspension springs since in their displacement the springs could interfere with other mechanical elements of the bogie generating interferences and breakages. To avoid this, the railway standard [17] requires a test for the springs of the highest category according to the scheme in Figure 2a. From the test, both the direction and the transverse bowing force (Φc) of the spring subjected to a defined axial load are obtained. The direction of the transverse direction is marked on the spring by a permanent system to take it into account in its mounting on the bogie. Figure 2b shows an example of a bogie spring assembly compensating for the spring bowing forces.

Kobelev [18] also studies another type of anomaly in the behavior of springs such as the transverse vibration of the spring once it enters into oscillations close to resonant frequencies. Kobelev [18] found that the fundamental natural frequency of transverse oscillations turns to be to zero when the lateral buckling of the spring occurs.

Yildirim [19] demonstrates that Equation (1), which is usually used to determine the spring stiffness, is only valid for spring angles α ≤ 10°. Figure 3a shows the spring geometry indicating the position of the angle α.

Yildirim [19] indicates that Equations (1) and (2) only take into account the effect of the torsion of the spring section when compressing it. However, as Yildirim [19] points out when the axial compression force acts, both stresses due to torsion and bending moment and forces normal and shear to the section appear. These last three effects are negligible for spring angles α ≤ 10° where the torsional effect dominates. Yildirim [19] proposes a global equation to determine the deformation of the spring subjected to an axial force. Additionally, Kato [20] proposes recently an equation outside the elastic range to determine large spring deformations.

In general, the previous authors focus the study of the spring behavior using a conventional Cartesian system in which the static or dynamic displacements are referred to as the first axis of symmetry of the helical spring, X. With this methodology the global behavior of a spring as a whole can be known with precision, i.e., constituted by a series of coils of a certain diameter and a certain thickness of the wire, contemplating in some of the occasions the ending effect.

This type of methodology, therefore, allows us to analyze the elastic effect of the spring as a whole, but it is not able to explain some mechanical behaviors intrinsic to its own constitution, i.e., it does not study the phenomenon occurring in the coil itself. In this way, anomalous spring behaviors such as those described above cannot be fully explained.

In order to achieve a reference system that meets these needs, the Frenet trihedral is used in this article, since it is intrinsic to the wire of the spring itself, configuring a reference system positioned according to the neutral line of the helix of the coil as shown in Figure 3b.

In the present study, the equations of motion of the coil of a helical spring are obtained in order to determine its modes of vibration. This is an original and novel approach since the usual practice is to obtain the modes of vibration of the complete spring referred to as a reference system according to its joint motion. Additionally, a sensitivity analysis has been performed for the main geometrical parameters of the coil: i.e., the helix diameter, the wire diameter and pitch. The obtained vibration modes have been validated satisfactorily with finite element calculations with ANSYS^® R17 (release 17).

2. Model of Coupled Bending and Torsional Oscillation in a Spring

2.1. Frenet’s Trihedral

For the approach of the oscillation model of the coil, the coordinate system formed by the Frenet trihedral, intrinsic to the neutral line of the spring, has been considered.

Figure 4 shows a differential element of the length of the curve defined by the midline of the coil of a helical spring. The Frenet trihedral or intrinsic trihedral is positioned on it. It is an orthonormal trihedron ($\overrightarrow{t},\overrightarrow{n,}\overrightarrow{b}$) whose components are named tangent, normal and binormal. The tangent vector is oriented to the neutral line of the loop. The normal is directed towards the bending center of curvature of the loop. The binormal vector follows the orientation defined by the other two. The vectors t and n form a tangent plane to the loop, containing the element ds, and perpendicular to the binormal vector which is directed towards the bending curvature of this element.

Considering as a scalar variable the arc length of the midline of the loop(s), we obtain the well-known Frenet–Serret equations, Equation (5), which give us the derivative with respect to s of the components of the trihedral as a function of the curvature of the midline of the coil.

$$\left(\begin{array}{c}\frac{d\overrightarrow{t}}{ds}\\ \frac{d\overrightarrow{n}}{ds}\\ \frac{d\overrightarrow{b}}{ds}\end{array}\right)=\left(\begin{array}{c}\overrightarrow{t\prime }\\ \overrightarrow{n\prime }\\ \overrightarrow{b\prime }\end{array}\right)=\left(\begin{array}{ccc}0& {\Omega }_3& -{\Omega }_2\\ -{\Omega }_3& 0& {\Omega }_1\\ {\Omega }_2& -{\Omega }_1& 0\end{array}\right)\left(\begin{array}{c}\overrightarrow{t}\\ \overrightarrow{n}\\ \overrightarrow{b}\end{array}\right)$$

(5)

where ${\Omega }_1={\tau }_0$, ${\Omega }_2=-k_0\sin {\varphi }_0$ and ${\Omega }_3=k_0\cos {\varphi }_0$. Where $k_0$ is the bending curvature, ${\tau }_0$ the torsional curvature and ${\varphi }_0$ the angle formed by the principal directions of inertia with the Frenet trihedron.

Equation (5) can be expressed as a function of the vector product of the Darboux vector $\overrightarrow{\Omega }\left(s\right)$ by the components of the intrinsic trihedron obtaining Equation (6).

$$\left(\begin{array}{c}\frac{d\overrightarrow{t}}{ds}\\ \frac{d\overrightarrow{n}}{ds}\\ \frac{d\overrightarrow{b}}{ds}\end{array}\right)=\overrightarrow{\Omega }\times \left(\begin{array}{c}\overrightarrow{t}\\ \overrightarrow{n}\\ \overrightarrow{b}\end{array}\right)$$

(6)

where the Darboux vector is defined according to Equation (7):

$$\overrightarrow{\Omega }\left(s\right)={\Omega }_1\overrightarrow{t}+{\Omega }_2\overrightarrow{n}+{\Omega }_3\overrightarrow{b}$$

(7)

The variation of any vector function, for example, the rotation and translation of the section, with respect to the length s of the arc follows the expression of Equation (8a) that we will simplify using the notation of Equation (8b) where we will call $\left[R\right]$ the tensor associated with the operation $\overrightarrow{\Omega }\times$.

$$\frac{d\overrightarrow{V}}{ds}=\frac{d\overrightarrow{V}}{ds}|{ }_S+\overrightarrow{\Omega }\times \overrightarrow{V}$$

(8a)

$$\left\{\overrightarrow{V}\right\}\prime =\left\{\overrightarrow{V}\prime \right\}+\left[R\right] \overrightarrow{V}$$

(8b)

$$\left[R\right]= \left(\begin{array}{ccc}0& {\Omega }_3& -{\Omega }_2\\ -{\Omega }_3& 0& {\Omega }_1\\ {\Omega }_2& -{\Omega }_1& 0\end{array}\right)$$

(8c)

where $\frac{d\overrightarrow{V}}{ds}|{ }_S$ is the derivative calculated as if the Frenet trihedron were fixed along the curve. That is, the variation of the function along the curve has two terms, the first represents the variation of the vector within the basis itself and the second represents the variation of the basis itself along the curve.

2.2. Rotation of the Spring Cross-Section

The spring in its natural geometrical configuration undergoes a deformation caused by external forces and moments. The plane cross-section remains plane and we neglect its small deformation. As shown in Figure 5, the cross-section of the deformed part is obtained by a translation $\overrightarrow{\eta }\left(s\right)$ of the center and a rotation $\overrightarrow{\psi }\left(s\right)$ of the section. Therefore, the point $P$ the intersection of the cross-section with the midline, will be transformed into the point $P_1$ (where $\overrightarrow{P{P}_1}=\overrightarrow{\eta }\left(s\right)$) and the perpendicular axes (j,k), coinciding with the normal and binormal of the intrinsic trihedron before the deformation, will become the axes $\overrightarrow{j_1} \mathrm{and} \overrightarrow{k_1}$ which we also assume to be perpendicular to each other. The axis $\overrightarrow{i_1}$ is obtained as the vector product of the previous ones, $\overrightarrow{i_1}=\overrightarrow{j_1}\times \overrightarrow{k_1}$.

Due to the shear produced by the shear forces, the i-axis is not transformed into the i-axis₁ by rotation and translation. Consequently, the tangent vector to a point on the midline of the part will no longer coincide with the first axis of the trihedron as it did in the initial equilibrium state but will form an angle $\overrightarrow{\beta }$, called the shear angle, proportional to the shear stress. If we do not consider the deformation due to the shear stress, the shear angle is canceled, so that the tangent to the midline of the part will be perpendicular to the section, or what is the same, will coincide with the first axis of the trihedron.

The derivative of $\overrightarrow{\psi }$ with respect to s is given to first order by Equation (9a,b).

$$\left\{\overrightarrow{\psi }\right\}\prime =\left[M_N\right] \overrightarrow{N}+\left[R\right] \overrightarrow{\psi }$$

(9a)

$$\left[M_N\right]=\left(\begin{array}{ccc}\frac{1}{k_{1}G{I}_1}& 0& 0\\ 0& \frac{1}{E{I}_2}& 0\\ 0& 0& \frac{1}{E{I}_3}\end{array}\right)=\left(\begin{array}{ccc}\frac{1}{I_1^{*}}& 0& 0\\ 0& \frac{1}{I_2^{*}}& 0\\ 0& 0& \frac{1}{I_3^{*}}\end{array}\right).$$

(9b)

where $M_N$ is the tensor associated with the torsional-flector torque increment, E is the modulus of elasticity, G is the shear modulus of elasticity, I_i is the geometric moment of inertia and k₁ is the torsion constant. $\overrightarrow{N}$ is the vector of internal torque increase. Its first component is the torsional moment, and the other two components are the bending moments in the normal and binormal directions. $I_1^{*}$ is the torsional stiffness, $I_2^{*}$ is the bending stiffness along the normal axis, and $I_3^{*}$ is the bending stiffness along the binormal axis.

2.3. Translation of the Spring Cross-Section

The derivative of the displacement vector $\overrightarrow{\eta }$ with respect to s is given to first order by Equation (10a–c).

$$\left\{\overrightarrow{\eta }\right\}\prime =\left[M_F\right] \overrightarrow{F}+\left[S\right] \overrightarrow{\psi }+\left[R\right] \overrightarrow{\eta }$$

(10a)

$$\left[M_F\right]=\left(\begin{array}{ccc}\frac{1}{EA}& 0& 0\\ 0& \frac{1}{k_{2}GA}& 0\\ 0& 0& \frac{1}{k_{3}GA}\end{array}\right)=\left(\begin{array}{ccc}{E}^{*}& 0& 0\\ 0& G_1^{*}& 0\\ 0& 0& G_2^{*}\end{array}\right)$$

(10b)

$$\left[S\right]=\left(\begin{array}{ccc}0& 0& 0\\ 0& 0& 1\\ 0& -1& 0\end{array}\right)$$

(10c)

where M_F is the tensor associated with the forces and k_i the shear constants, or the so-called Timoshenko’s k-factors [21]. Being for symmetric sections $k_2=k_3$. $\overrightarrow{F}$ is the vector of internal force increment. Its first component is the axial force, and the other two components are the shear forces according to the normal and binormal. $EA$ is the axial stiffness and $k_{2}GA$ is the shear stiffness.

Equations (9a,b) and (10a–c) constitute a system of six scalar equations for the twelve components of $\overrightarrow{\psi }, \overrightarrow{\eta }, \overrightarrow{F}, \overrightarrow{N}$. The remaining six equations are obtained from the dynamic study of the system for which the equations of the linear momentum and angular momentum of the part element comprised between s and s + ds are posed.

2.4. Linear Momentum

Considering the linear momentum, the derivative of $\overrightarrow{F}$ with respect to s is given by Equation (11).

$$\left\{\overrightarrow{F}\right\}\prime =\lambda \ddot{\overrightarrow{\eta }}+\left[R\right] \overrightarrow{F}$$

(11)

where $\lambda =\rho ·A$ is the linear density per unit length.

2.5. Angular Momentum

Considering the angular momentum, the derivative of $\overrightarrow{N}$ with respect to s is given by Equation (12a,b), which takes into account the effects of rotational inertia introduced by Rayleigh [22].

$$\left\{\overrightarrow{N}\right\}\prime =\frac{\lambda }{A}\left[I\right] \ddot{\overrightarrow{\psi }}+\left[S\right] \overrightarrow{F}+\left[R\right] \overrightarrow{N}$$

(12a)

$$\left[I\right]=\left(\begin{array}{ccc}{I}_1& 0& 0\\ 0& I_2& 0\\ 0& 0& I_3\end{array}\right)$$

(12b)

Table 1. System of Vectorial Equations to model spring behavior.

Vectorial Equations
$\left\{\overrightarrow{\psi }\right\}\prime =\left[M_N\right] \overrightarrow{N}+\left[R\right] \overrightarrow{\psi }$
$\left\{\overrightarrow{\eta }\right\}\prime =\left[M_F\right] \overrightarrow{F}+\left[S\right] \overrightarrow{\psi }+\left[R\right] \overrightarrow{\eta }$
$\left\{\overrightarrow{F}\right\}\prime =\lambda \ddot{\overrightarrow{\eta }}+\left[R\right] \overrightarrow{F}$
$\left\{\overrightarrow{N}\right\}\prime =\frac{\lambda }{A}\left[I\right] \ddot{\overrightarrow{\psi }}+\left[S\right] \overrightarrow{F}+\left[R\right] \overrightarrow{N}$

summarizes the complete system of equations obtained and which describes the behavior of the spring coil subjected to bending and torsion. It consists of four vector equations in linear partial derivatives, non-homogeneous and of second order.

2.6. Dispersion Relation

If in the Table 1 system of equations we change the variable using Equation (13a–d), we can simulate the movements in time of the spring and then perform a vibration analysis of its oscillatory behavior in the vector variables $\overrightarrow{\psi }, \overrightarrow{\eta }, \overrightarrow{F}, \overrightarrow{N}$. We will assume a wave traveling along the axis of the wire. Thus

$$\overrightarrow{\psi }\left(s,t\right)=Re \left(\overrightarrow{{\psi }_0}\left(s\right)·e^{i(ks-\omega t)}\right)$$

(13a)

$$\overrightarrow{\eta }\left(s,t\right)=Re \left(\overrightarrow{{\eta }_0}\left(s\right)·e^{i(ks-\omega t)}\right)$$

(13b)

$$\overrightarrow{F}\left(s,t\right)=Re \left(\overrightarrow{F_0}\left(s\right)·e^{i(ks-\omega t)}\right)$$

(13c)

$$\overrightarrow{N}\left(s,t\right)=Re \left(\overrightarrow{N_0}\left(s\right)·e^{i(ks-\omega t)}\right)$$

(13d)

The new parameters introduced are the wavenumber, k, and the angular velocity of the wave, ω (rad/s). Where the wavenumber k, is related to the length λ of the vibrating wave, with Equation (14).

$$k=\frac{2 \pi }{\lambda }$$

(14)

Substituting the system of equations in Table 1 we arrive at the following system in Equation (15a–d).

$$\left[\mathsf{\Lambda }\right] \overrightarrow{F}-\lambda {\omega }^2 \overrightarrow{\eta }=0$$

(15a)

$$\left[\mathsf{\Lambda }\right] \overrightarrow{\eta }+\left[M_F\right] \overrightarrow{F}+\left[S\right] \overrightarrow{\psi }=0$$

(15b)

$$\left[\mathsf{\Lambda }\right] \overrightarrow{\psi }+\left[M_N\right] \overrightarrow{N}=0$$

(15c)

$$\left[\mathsf{\Lambda }\right] \overrightarrow{N}-\frac{\lambda {\omega }^2}{A}\left[I\right] \overrightarrow{\psi }+\left[S\right] \overrightarrow{F}=0$$

(15d)

where the matrix $\left[\mathsf{\Lambda }\right]$ follows the expression of Equation (16).

$$\left[\mathsf{\Lambda }\right]=\left(\begin{array}{ccc}-ik& k_0& 0\\ -k_0& -ik& {\tau }_0\\ 0& -{\tau }_0& -ik\end{array}\right)$$

(16)

Solving for $\overrightarrow{\eta }$ in Equation (15a) and replacing it in Equation (15b), and similarly solving for $\overrightarrow{N}$ in Equation (15c) and replacing in Equation (15d) we obtain the following vector Equation (17a,b):

$$\left[A\right] \overrightarrow{F}+\left[S\right] \overrightarrow{\psi }=0$$

(17a)

$$\left[B\right] \overrightarrow{\psi }+\left[S\right] \overrightarrow{F}=0$$

(17b)

where the matrices $\left[A\right], \left[B\right]$ are calculated with Equations (18) and (19).

$$\left[A\right]=\frac{1}{\lambda {\omega }^2}{\left[\mathsf{\Lambda }\right]}^2+\left[M_F\right]$$

(18)

$$\left[B\right] \overrightarrow{\psi }+\left[S\right] \overrightarrow{F}=0$$

(19)

Equation (17a,b) constitutes a system of six equations with six unknowns ($F_1,F_2,F_3, {\psi }_1, {\psi }_2,{\psi }_3$) and with null independent terms. It is therefore a homogeneous system and for it to have a solution other than the trivial, null solution, it must be verified that its determinant is null, i.e., Equation (20).

$$\left|\begin{array}{cccccc}\lambda {\omega }^2{E}^{*}-\left(k^2+{\beta }_9\right)& {\beta }_4-2ik{\Omega }_3& {\beta }_5+2ik{\Omega }_2& 0& 0& 0\\ 2ik{\Omega }_3+{\beta }_4& \lambda {\omega }^2{G}_1^{*}-\left(k^2+{\beta }_8\right)& {\beta }_6-2ik{\Omega }_1& 0& 0& \lambda {\omega }^2\\ {\beta }_5-2ik{\Omega }_2& {\beta }_6+2ik{\Omega }_1& \lambda {\omega }^2{G}_2^{*}-\left(k^2+{\beta }_7\right)& 0& -\lambda {\omega }^2& 0\\ 0& 0& 0& k^2{I}_1^{*}+{\alpha }_4-{\omega }^2{\lambda }_1& -{\beta }_4{I}_3^{*}+ik{\alpha }_1& -{\beta }_5{I}_2^{*}-ik{\alpha }_2\\ 0& 0& 1& -{\beta }_4{I}_3^{*}-ik{\alpha }_1& k^2{I}_2^{*}+{\alpha }_5-{\lambda }_2{\omega }^2& -{\beta }_6{I}_1^{*}+ik{\alpha }_3\\ 0& 1& 0& {\beta }_5{I}_2^{*}-ik{\alpha }_2& {\beta }_6{I}_1^{*}+ik{\alpha }_3& -k^2{I}_3^{*}+{\omega }^2{\lambda }_3-{\alpha }_6\end{array}\right|=0$$

(20)

Operating the determinant and simplifying we obtain the alphanumerical Equation (21), which would be the dispersion relation, $\omega =f\left(k\right)$, which relates the length of a wave with its oscillation frequency in the spring coil.

$$m_1{\omega }^{12}+m_2{\omega }^{10}+m_3{\omega }^8+m_4{\omega }^6+m_5{\omega }^4+m_6{\omega }^2+m_7=0$$

(21)

where the coefficients of the polynomial are calculated according to the expressions given in Equation (A1) of Appendix A.

The real roots in ${\omega }^2$ of Equation (21) will give us the number of branches of the dispersion equation for each value of the wavenumber k.

If in the dispersion relation obtained for the spring we cancel the independent term, $m_7$ (which is equivalent to making $\omega \to 0$), we will have the solutions of the wave number shown in Equation (22a,b).

$$k^{*}=0$$

(22a)

$$k^{**}=\sqrt{{\Omega }_1^2+{\Omega }_3^2}$$

(22b)

That is to say, there exists a value of the wave number k different from zero (to which corresponds, therefore, a finite wavelength) for which one of the real and positive roots of the dispersion relation (degree 12) cancels, so that the branch corresponding to this solution, for any value of k, is tangent to the x-axis in $k=k^{**}$.

If we do $k\to 0$ we will obtain the number of branches of the dispersion relation (which will be equivalent to the number of real roots of the equation). For this case, the real and positive branches are determined by Equation (23a–c):

$$\omega =0$$

(23a)

$$\begin{gathered} {\rho }^2{I}_1{I}_2{\omega }^4-{\omega }^2\left[kG\left\{A\rho I_1+\left({\beta }_1+{\beta }_3\right)I_2\right\}+E{I}_2\rho \left(I_1{\beta }_1+I_2{\beta }_3\right)\right] \\ +kGE{I}_2\left[\left({\beta }_1+{\beta }_3\right)\left(I_1{\beta }_1+I_2{\beta }_3\right)+A{\beta }_3\right]=0 \end{gathered}$$

(23b)

$$\begin{gathered} {\rho }^2{I}_2{\omega }^4-{\omega }^2\rho \left[E{I}_2\left({\beta }_1+{\beta }_3\right)+kG\left[\left({\beta }_1{I}_2+A\right)+{\beta }_{3}kG{I}_1\right]\right] \\ +\left({\beta }_{3}kG{I}_1+{\beta }_{1}E{I}_1\right)\left({\beta }_{3}E+{\beta }_{1}kG\right)+{\beta }_{3}kGEA=0 \end{gathered}$$

(23c)

The wave numbers k obtained from Equation (23a–c) are indicators of the spring oscillation modes. Equation (23b) represents the bending oscillations of a differential element of spring length in plane 1, 2 determined by the principal directions of inertia of a cross-section (Figure 6). Equation (23c) represents the torsional oscillations of the differential element around the principal direction 1 of Figure 6. These oscillation modes will be explained also later in the finite elements section.

The result of the oscillation is bending accompanied by torsion which tends to make the ends of the coil open and close in time, oscillating according to the movement. This is composed of two modes: opening and closing of the ends, the coil remaining at all times inscribed in the original helix of the spring (of diameter D) and the other mode consists of the ends of the coil entering and leaving the original helix of diameter D, forming an ellipse. This last mode responds to the reason why the springs, under certain dynamic conditions, pull the coils out of their original helix, and may impact other elements of the machine, thus affecting its operation. Authors consider this can explain part of the “Effort de Chasse” anomaly described in Section 1.

In this article, we are in the case of a freely oscillating coil. In the case of several coils joined together, it would be necessary to consider the ligatures of their ends as shown in Figure 7 for the end of a squared and ground spring. That is to say, with a helix angle of zero degrees and a flat seat, which allows a better load transfer. The assembly thus defined becomes an oscillatory mode in which the spring as a whole adopts the barrel configuration. Pearson [23] studies the important contribution of the ends to the calculation of the spring natural frequencies.

The dispersion relation allows us to obtain the normal modes of oscillation by applying the appropriate boundary conditions. For example, in the case of a helical spring with one end clamped and the other free, the boundary conditions would be as follows: the rotation vector $\overrightarrow{\psi }$ and the translation vector $\overrightarrow{\eta }$ cancel out at the clamped end and the internal torque $\overrightarrow{N}$ and the internal stress $\overrightarrow{F}$ cancel at the free end. Assuming the clamped end for s = 0 and the free end for s = L, Equation (24) is satisfied which implies that the wave number k takes the values of Equation (25).

$$\overrightarrow{\psi }\left(0,t\right)=\overrightarrow{\psi }\prime \left(L,t\right)=0$$

(24)

$$k=\frac{n·\pi }{2L}=\frac{n·\pi }{2N\sqrt{h^2+{\pi }^2{D}^2}} n=1,3,\dots ..$$

(25)

Entering the value of k calculated according to Equation (25) in the dispersion Equation (21), the natural oscillation frequencies for the helical spring are obtained.

In order to validate the proposed theoretical model, a comparison with finite elements and with the experimental tests of a helical spring is made below.

3. Application of Finite Elements

The previous calculation performed to obtain the natural frequencies of the system has been compared with the results of the resonance frequencies obtained by numerical simulation. In addition, a complete modal analysis is performed. In order to compare the results, the oscillation modes of the one spring turn are calculated using standard computational methods. The simulations have been carried out by modeling the spring geometry by Finite Elements using ANSYS^® software.

A three-dimensional geometry of the helicoid representing the neutral fiber of a helical spring has been faithfully reproduced. The three-dimensional geometry has been generated from SOLID EDGE software, and it has been imported by the FEM software in the *.iges format. The SOLID EDGE software provides flexibility in the modelization and facilitates the actual geometrical definition compared to the CAD tool integrated in ANSYS. The spring is thus reproduced volumetrically.

We proceed to calculate in both cases, i.e., the theoretical model proposed in the previous section with one solved in finite elements, in order to compare the results.

In finite elements, a meshing has been performed using extruded hexahedral elements of 2 mm characteristic size, which represents the circumference of the wire cross-section (10 mm diameter) with 16 nodes in its contour. This results in 6000 nodes per turn, 54,000 nodes for the complete spring and 162,000 degrees of freedom.

The meshing strategy used consists of an extrusion method from the two-dimensional meshing (with quadrilaterals) of the cross-section, which guarantees a regular meshing that avoids geometric distortions of the elements and achieves a low bandwidth for the stiffness matrix, resulting in a better performance of the iterative methods of resolution.

The SOLID185 element has been selected in its homogeneous form, which is frequently used to model three-dimensional structures in ANSYS. The element consists of eight nodes with three degrees of freedom per node: translations in the three directions of space. SOLID185 uses the reduced selective integration method.

The material has been modeled linearly, with steel properties considering Young’s modulus and Poisson’s ratio.

The problem has been plausibly solved using the Lanczos method for the determination of eigenfrequencies and eigenmodes.

In free-free solid conditions, the different natural frequencies of one spring turn belonging to a helical spring, with different helix diameters and wire thickness have been calculated. In all cases, the pitch of the helix has always been the same.

In all cases, the corresponding oscillation modes are studied.

Table 2. FEM calculation of the first natural frequencies (Hz) of helical spring.

Natural Frequencies (Hz) Calculated by ANSYS
D (mm)	d (mm)	Pitch (mm)	Bending 1	Torsional 1	Bending 2	Torsional 2	Bending 3	Torsional 3
80	18	39.8	960	1218	2111	3184	4543	6322
80	22	39.8	1168	1468	2556	3733	5428	7414
80	28	39.8	1475	1823	3199	4435	6623	8254
108	18	39.8	539	664	1183	1815	2595	3721
108	22	25	667	792	1451	2155	3178	4339
108	22	39.8	657	806	1439	2162	3132	4349
108	22	50	648	817	1426	2169	3091	4355
108	28	39.8	832	1012	1815	2630	3899	5127
220	18	39.8	133	158	290	455	648	968
220	22	39.8	162	193	354	553	790	1169
220	28	39.8	207	245	450	695	1000	1454

N = 1 turn, E = 210 GPa, G = 84 GPa, $\rho$ = 7800 kg/m³, ν = 0.24407.

shows the results. Columns 4 to 9 show the frequency and the mode of vibration. Essentially, the first mode corresponds to an opening between ends, such that the original helix is transformed into an ellipsoidal helix (Figure 8a). In this case, the predominant deformation of the helix is bending. The second mode corresponds to the relative motion between these ends (Figure 8b) but preserves the circumference of the helix, which corresponds to a mode in which the torsional effect of the helix is more predominant than its bending. In this second case, the pitch of the coil would increase significantly.

From these two basic modes, the deformation profiles are reproduced with increasing inflection points. The third (Figure 8c) and fifth (Figure 8e) modes constitute double bending and triple bending, respectively. Always without leaving the imaginary plane that would contain the spiral if it were a ring in a plane and not helical, but distorting the circumference that defines it. The fourth (Figure 8d) and sixth (Figure 8f) modes constitute double torsion (with one inflection point) and triple torsion (with two inflection points), respectively.

If we continue to search for natural modes at increasing frequencies, increasingly complex deformation profiles with small displacements and numerous inflection points appear.

With respect to the results shown in Table 2, it can be deduced that, in the case of the first resonance mode (bending), keeping the helix diameter constant and varying the wire thickness, the resonance frequency increases with increasing thickness. This is logical because the stiffness of the spring increases. Figure 9 shows this trend. Moreover, as shown in Figure 9, as the diameter of the helix increases, the resonance frequency decreases. Despite maintaining the stiffness, the loop flexes with lower frequency values. A similar situation occurs with the first torsional mode.

Proceeding analogously for the first bending mode, the wire thickness is now kept constant, and the average diameter of the helix is varied (Figure 10). The frequencies decrease as the helix diameter increases. Moreover, as the wire diameter increases, the resonant frequency increases. The situation is repeated for the first torsional mode of vibration. It is a situation analogous to that of the previous paragraph.

Analyzing the variation of the frequency as a function of the variation of the average diameter of the helix, it is observed (Figure 10) that the first bending and torsional mode of the spring decrease as the mass increases due to a higher helix diameter. This has to do with the wire diameter/helix diameter ratio and it can be seen that there is a certain characteristic value of this ratio that linearizes the resonant behavior.

As shown in Figure 9 the relationship between frequencies and wire diameter (d) is linear. This result is consistent with Equation (3).

It has also been considered whether the pitch of the spring can be an influential design parameter when calculating resonances. Figure 11 shows that the pitch variation does not affect the frequencies. Figure 11 shows the different modes obtained and commented.

Table 3. Comparison of the numerical model with FEM for the first harmonic (Hz).

Fundamental Frecuency (Hz)
D (mm)	d (mm)	Pitch (mm)	Wavenumber k (m−1)	Bending 1 FEM	Bending 1 Numerical Model Equation (21)
80	18	39.8	24.69	960	989
80	22	39.8	24.69	1168	1190
80	28	39.8	24.69	1475	1495
108	18	39.8	18.39	539	555
108	22	25	18.47	667	680
108	22	39.8	18.39	657	670
108	22	50	18.32	648	665
108	28	39.8	18.39	832	840
220	18	39.8	9.08	133	145
220	22	39.8	9.08	162	170
220	28	39.8	9.08	207	211

N = 1 loop, E = 210 GPa, G = 84 GPa, $\rho$ = 7800 kg/m³, ν = 0.24407.

includes a comparison for the first harmonic of the results obtained with the finite element simulation and the numerical model presented in this paper. The results of the numerical model are obtained by solving $\omega$ in Equation (21) with the wavenumber k calculated with Equation (22b) or Equation (25). In general, it is observed that the numerical model overestimates the frequency values with an average difference of 3%.

4. Conclusions

The study of a single coil has been proposed since, in general, the authors who propose equations to model the vibration behavior of a spring consider the spring as a whole, including its ends, but do not analyze the behavior of the coil itself. From the analysis of a single coil, the singular behavior of the spring known as “effort de chasse” described in the introduction could be justified.

For the study of the vibratory behavior of a spring, an analytical model has been developed based on a reference system consisting of a Frenet trihedron whose tangent vector is the midline of the coil. One of the results obtained from the model is the dispersion relation that relates the wavelength to its oscillation frequency, $\omega =f\left(k\right)$. Applying the appropriate boundary conditions the dispersion relation allows us to calculate the normal modes of oscillation.

Another of the results obtained is that there is a value of the wavenumber $k=\sqrt{{\Omega }_1^2+{\Omega }_3^2}$ for which the roots of the dispersion equation cancel out. Being ${\Omega }_1$ and ${\Omega }_3$ the torsional curvature and the bending curvature, respectively.

Two clear modes of oscillation have been characterized, one in which the coil moves according to the torsion of the cross-section increasing the pitch of the coil, and another bending mode in which it moves by making an opening of its ends so that the original helix is transformed into an ellipsoidal helix.

A spring has been modeled by finite elements with ANSYS^® software in order to contrast results with the numerical model, obtaining an average difference of 3% in the fundamental oscillation frequency.