Abstract

This paper presents a modeling study of the dynamics of a helical spring element with variable pitch and radius considering both the static stiffness and dynamic response by using the geometrically exact beam theory. The geometrically exact beam theory based on the Euler–Bernoulli beam hypothesis is described, of which the shear deformations are ignored. Unlike the traditional spliced curved beam element method, the helical spring element is described with curvature vector and axial strain by establishing and spline-interpolating a function of the radius, the height, the polar angle, and the torsion angle of the whole spring. In addition, a model smoothing method is developed and applied in the numerical analysis to filter the high-frequency oscillation component of the flexible multibody systems, so as to correct the system dynamic equations and improve the calculation efficiency when solving the static equilibrium of the spring. This study also carries out five numerical trials to validate the above dynamic procedure of the helical spring element. The example of the spring static stiffness design shows that the proposed helical spring procedure enables one to deal with practical engineering applications.

1. Introduction

In mechanical systems, the helical spring plays a crucial role in the design of the entire system. As in the automotive engine valve train, the spring stiffness needs to meet the preload force requirements of the valve group subassembly. Considering that the helical spring contains many design parameters, the stiffness of the spring is difficult to calculate numerically due to the diversity of the irregular springs, so a more detailed model is needed for the spring analysis.

The study of helical spring modeling methods has gone through a long process. Love [1] developed the 12-order differential equation of motion for a spring with a very small helix angle and made the assumption that the axial force and the shear forces did not produce any deformation. Wahl [2] considered the spring wire as a round bar that withstands shear and torque. He used the correction factor to consider the curvature of the spring and ignored the coupling between axial and torsional deformation. Wittrick [3] modeled the helical spring by the Timoshenko beam theory incorporating the shear deformation and the rotary inertia and developed the governing partial differential equations of motion of the coil spring when the helix angle is not small. Mottershead [4] developed a special finite element method for the modeling of helical springs, which is adopted for dynamic analysis of helical rods. The element stiffness and mass matrix were based on an exact differential equation that governs the static behavior of infinitesimal elements. The element displacement function is obtained through the integral differential equation, which may be one or more turns of the spring or a part of a turn. Wittrick and Pearson [5] simplified the model by ignoring the rotational inertia and the shear deformation which has become a standard method converting a Timoshenko model to a classical Bernoulli–Euler model. The theory for calculating the dynamic stiffness matrix and the natural frequency was also proposed. Lee and Thompson [6] derived the free-wave motion equation in the coil spring as well as the dynamic stiffness matrix by means of the Timoshenko beam theory and the Frenet equation. The natural frequency was calculated from the dynamic stiffness matrix using the Wittrick–Williams method by applying appropriate boundary conditions. Taktak et al. [7, 8] proposed a two-node finite element spring model with six degrees of freedom at each node, which was able to model a three-dimensional isotropic helical spring with a constant curvature radius and torsional radius, taking account of the shear strain. In the above literature, the helical spring was set with constant pitch and radius. Chaudhury and Datta [9] proposed a method for designing circular sections of nonstandard coil springs with variable pitch and radius using analytical methods and had compared them with the finite element analysis. It shows that, using the straight beam element to simulate the spring and to analyze the deformation details of the helical spring with variable pitch and radius, the result is inefficient and inaccurate.

The helical spring is an elongated structure having an initial helix curvature. And it undergoes large displacements and finite rotations. The spring’s nonlinear dynamic behavior is usually simulated by using the geometric nonlinear spatial beam elements. As is well known, the geometrically exact beam theory allows us to formulate problems involving arbitrarily large displacements, rotations, and small strains. Simo and Vu-Quoc [10–12] introduced the global position and rotation vector as the nodal coordinates to independently interpolate the incremental displacement and the rotation fields. Zhang et al. [13] developed a two-node spatial beam element theory with the Euler–Bernoulli assumption for the nonlinear dynamic analysis of slender beams undergoing arbitrary rigid motions and large deformations. Le et al. [14] compared four nonlinear dynamic calculation formulas of three-dimensional beam elements, of which the inertia force vector and the tangential inertia matrix are derived from the total Lagrangian formulation and are calculated using two Gaussian points. In the Euler–Bernoulli beam theory, the normal of a cross-section is in parallel with the tangent of the central line. It means that the translational displacement and the rotation fields couple. Since the work of Simo [10], the geometrically exact initially curved beam theory has been developed. Adnan [15] discussed the finite element implementation of the three-dimensional curved beam elements of Reissner. The work of Christoph Meier et al. [16]was the development of a new finite element formulation for a curved beam according to the geometrically exact Kirchhoff theory of thin rods. Wenxiong Li et al. [17] developed mixed finite elements for geometrically nonlinear analysis of frame structures with curved members by the geometrically exact beam theory in order to construct the independent internal force field in the deformed configuration. Zupan et al. [18] presented the finite-element formulation of geometrically exact 3D beam theory based on the interpolation of strain measures. However, the adoption of a curved beam element to piece the helical spring based on the geometrically exact beam theory will increase the number of freedom degrees and reduce the calculation efficiency. At the same time, it is very difficult and inaccurate to use the existing curved beam element theory to establish the element function with the initial curvature of the helix. Therefore, it is necessary to find an analytic solution for the initial curvature of the helix when using the geometrically exact beam theory to construct the spring element.

Taguchi and Gen [19] transformed the optimal mass design problem of the helical spring with the allowable shear stress, effective winding number, and average spring wire section diameter as constraints into a nonlinear integer programming problem and solved the nonlinear constraints directly by using an improved genetic algorithm. Xiao et al. [20] introduced a spiral spring optimization algorithm based on particle swarm optimization algorithm. The optimization algorithm took the minimum mass of the coil spring as the objective function and the spring geometric parameters as the design variables, with shear stress, maximum shaft deflection, critical frequency, collision, fatigue strength, coil nontouch, and space limitation as constraints. Sun et al. [21] applied the static load by the ESL equivalent method, transformed the structural optimization nonlinear dynamic response into a static one, and solved the corresponding nonlinear static response structure optimization problem by nonlinear interior point method. Taktak et al. [22] proposed a helical spring optimization method based on the numerical formula proposed in the previous work, which combined dynamic parameters as constraints and objective functions to optimize the mechanical properties of the helical spring according to its geometric parameters, and calculated the dynamic response of the helical spring.

In this paper, we combine the geometrically exact beam elements with the Euler–Bernoulli beam hypothesis and the spring structural features to establish a model to simulate the dynamics of the helical spring with variable pitch and variable radius. Compared with other methods, the spring element not only describes the nonlinear dynamic behavior of the spring with a small number of description variables but also introduces the geometrically exact beam theory only with the assumption of rigid section in the modeling process to ensure the accuracy of the model. At the same time, the model smoothing method is introduced into the numerical calculation to improve the computational efficiency. All of these are beneficial to the design of the spring and the modeling and analysis of complex multi-flexible body systems that include spring components.

The outline is as follows: in Section 2, we describe the detail of the geometrically exact helical spring element using the Euler–Bernoulli beam assumption. A whole helical spring is provided based on the cubic spline interpolation, and the governing equations are derived based on the virtual power principle. In addition, several illustrative numerical examples are carried out in Section 3 to verify the validation of the developed helical spring element theory and to determine the designed static stiffness of the spring which meets the practical engineering application. Finally, some conclusions are presented in Section 4.

2. Development of Dynamics Helical Spring Element

During the deformation of the spring under compression or tension, it should be pointed out that the spring still maintains a spiral configuration. Considering the helical character, the spring can be modeled as a whole element, instead of being spliced by traditional curved beam elements. Based on the geometrically exact beam element theory described in the previous section, the curvature vector and axial strain are adopted to describe the flexible deformation within the helical spring element. It is also assumed that the spring is generated by a succession of rigid cross-sections orthogonal to the curvilinear axis.

2.1. Equations of the Helical Spring

The helical spring is three-dimensional and has an initial helix centroid line, defined in the global coordinate system . The equation of the curvilinear axis is introduced first and used to derive the equations of the movement of the helix spring. Figure 1 shows an initial configuration of the helix spring with variable pitch whose central axis forms a right-handed helix. In this paper, a superposed ‘bar’ denotes the parameters in the initial configuration. represent the helix radius, polar angle, and height of the line of centroids in the initial configuration.

Figure 1 

Initial geometric configuration of variable pitch helix.

The vectors of the helix with the global coordinate system can be defined bywhere indicate the helix radius, polar angle, and height, respectively, in the current configuration.

A family of orthonormal vector bases on the helix spring is introduced to describe the orientation of the cross-sections. and are oriented along the principal axes of inertia of the cross-section, respectively, and is normal to the cross-section: . For the Euler–Bernoulli beam assumption, the cross-section remains perpendicular to the tangent vector of the centerline during deformation. The relationship between the global base vector and cross-section base vector can be expressed aswhere R is the rotation matrix with CARDAN angles as follows:

Here, denote sine and cosine of Cardan angles, respectively.

2.2. Generalized Coordinates of the Spring Element

Helical spring comes in many types, such as variable pitch and variable radius. For any type, the configuration curve of the helix radius, height, polar angle, and torsion angle of the spring in the initial state is known, which are the initial polar coordinate’s functions. Figure 2 shows a schematic of the helix spring in the initial configuration and the current configuration with deformation. is the endpoint in the i-th circle. The parameters of the spring are changed by the stress, but the spring remains helical.

Figure 2 

The schematic of the helix spring in the initial and current configuration.

The position of the endpoints is shown in Figure 2. Taking the helix radius as an example, the configuration curve of the helix radius in the initial state and current state is as shown in Figure 3. In this paper, the symbol denotes the derivative with respect to the initial polar angle . is the helix radius of the endpoint . By employing the cubic spline interpolation method to discretize the configuration curve, the centroid vector at any point of the spiral can be obtained.

Figure 3 

Configuration curve of the variable helix radius in the initial state and current state.

The helix radius in the i-th circle, shown in Figure 3, is described aswhereare the Hermite shape functions with , .

For spline interpolation, the first derivative of the radius at the endpoint to the initial polar coordinate is expressed aswhere and and are the constant matrix. In this way, the first derivative of the intermediate node of the spring element can be expressed at both ends.

Thus, the helix radius corresponding to any point on the spiral iswhere is the shape function matrix as follows:

In equation (7), the vector of the endpoint coordinates is defined as

In the same way, the helix height, polar angle, and torsion angle of the wire section to any point, respectively, are written as

The generalized coordinates of the spring are given by

The discrete method enhances the continuity of the model while reducing the degree of freedom of the system, rather than splicing the spring with a section of the curved beam element. By substituting (7), (10), (11), and (12) into (1), the vector of the helix in the global coordinate system is obtained by the cubic spline interpolation method.

2.3. Normal Strain and Curvature of the Spring Element

Since the geometric nonlinear spring element takes into account the assumption of the Euler-Bernoulli beam, the strain only contains normal strain and curvature. If no shearing effect is considered, the normal vector of the cross-section is the tangent direction of the helix curvilinear axis in the current configuration and can be described as

Taking partial derivative of (1) with respect to the initial polar coordinate leads towhere and .

The arc length of an infinitesimal spatial element is defined by

According to expression (16), the derivative of the arc-length coordinate with respect to the initial polar coordinate can be expressed by

Expression (14) becomes with expressions (15) and (17)and we also havewhere , .

Combining the relationship between the CARDAN angles and the elements in the rotation matrix R, we get

Particularly, we note that is the helix angle of spring, so we also have and avoid the singular position of the CARDAN angles to describe the rigid-body rotation effectively. Since , the CARDAN angle can be derived aswhere and .

Substituting expression (18) into expressions (20) and (21), we can solve for and .

Taking the derivative of expressions (20) and (21), the derivative of the CARDAN angles with respect to the initial polar coordinate can be expressed bywhere and

Substituting expression (19) into expressions (22), we can solve for and .

The rotation of a rigid cross-section is described with the help of axes of rotation and corresponding CARDAN angles. The vectors of rotation axes are

The curvature equations provide the relationship of the curvature of the rigid cross-section fixed system relative to the initial system. The curvature vector can be expressed bywhere .

The components of in the rigid cross-section fixed system are extracted from this in the following:where and the CARDAN angle array can be expressed by .

Similar to the description of the current configuration, the components of in the initial configuration are extracted from this in the following:where and the CARDAN angle array in the initial configuration can be expressed by .

The curvature strain measures are then extracted from above, leading to bending strains and and torsion strain , respectively.

According to the definition of strain, the normal strain is given by

2.4. Kinematic Description and Virtual Variation of the Spring Element

In this paper, the symbol denotes the time derivative. According to the cubic spline interpolation of the global position vector of the helix, the velocity of the curvilinear axis vector and its virtual variation can be expressed bywhere .

The spring element generalized velocity and its virtual variation are given by

The accelerated velocity of the vectors of the helix can be obtained by the time derivation of

In order to facilitate the following derivation, the time derivative and virtual variation of arewhere , , , and .

The time derivatives of the normal strain and its virtual variation can be expressed by

The time derivatives of can be expressed bywhere , and , .

The time derivatives of the CARDAN angles and its virtual variation arewhere .

The second-order time derivatives of the CARDAN angles can be expressed bywhere , , , , , and .

The angular velocity of the cross-section is

The components of in rigid cross-section fixed configuration are extracted from this in the following:and its virtual variation can be written as

The angular acceleration iswhere , , and .

The time derivatives and virtual variation of can be expressed by

2.5. Inertial Force Virtual Power

The inertial force virtual power of the rigid section includes the translational inertial force virtual power and the rotational inertial force virtual power.

The translational inertial force virtual power is given bywhere is spring spiral density.

The rotational inertial force virtual power is given bywhere is rigid-section inertia moment tensor.

2.6. Internal Virtual Power and Average Stress

With the Euler–Bernoulli beam assumption, the internal virtual power of the helix spring, calculating the integral over the range , can be written aswhere n is spring coil number and is the internal virtual power of the infinitesimal element which is shown in Appendix.

The constitutive equation of the spring is defined aswhere is the stress array and is the strain array. The components of and are given relative to the cross-section basis .

The constitutive matrix iswhere is the axial stiffness, is the torsion stiffness, and , are the principal bending stiffnesses.

With the constitutive equation of the spring, the reduced expression of the internal power is given by

In terms of (32) and (40), the internal virtual power equation (46) can be written as

The spring system dynamics equation has strong rigidity, while the high-frequency elastic vibration mainly comes from the rapid change of strain. The model noise reduction method is to average the strain in a time interval . The is called the noise reduction factor [23]. In this way, the internal virtual power of the spring system is calculated.

A second-order Taylor series expansion is employed to represent the average strain. It can be given by

Modified internal virtual power is given by

Compared with the internal virtual power before modification, the inertia term and the damping term are added, so that the natural frequency of the spring system can be reduced and higher calculation efficiency can be obtained.

2.7. Governing Equation

The method of applying a test external force in the axial direction of the spring is adopted in this paper. The external force virtual power includes the gravitational virtual power and the test force’s virtual power.

The gravitational virtual power can be expressed as

The test external force’s virtual power can be written as

The axial test external load is shifted to the action node and a corresponding torque is generated aswhere is the mean coil radius at the end node.

The virtual power equation of the spring element system can be written as

The mass matrix of the system is defined by

The generalized force vector is written as

The governing equation of the spring system is calculated by Gaussian integration directly.

3. Numerical Results

In this section, several numerical examples are carried out to evaluate the static equilibrium calculation and dynamic response of the proposed dynamics helical spring element. The governing equation (53) is solved using an explicit Runge-Kutta method (ODE45) implemented in MATLAB. The relative error tolerance is , the absolute error tolerance is , and the maximum step size is .

3.1. A Numerical Example of the Variable Pitch Spring Model

In the example, the model of spring with variable pitches is investigated, shown in Figure 4. The geometrical properties and material properties of the spring are presented in Table 1.

Figure 4 

Model of a spring with variable pitches in the initial configuration.

The curvature of the spring with variable pitches in the initial configuration is shown in Figure 5, which is modeled by the proposed spring element.

(a)

(b)

(c)

(a)
(b)
(c)

Figure 5 

Curvature of the variable pitch spring in the initial configuration. (a) The torsion curvature of the spring in the initial state. (b) The bending curvature of the spring in the initial state. (c) The bending curvature of the spring in the initial state.

As can be seen from Figure 5, the curvature vector obtained by the spline interpolation maintains good continuity and smoothness inside the spring element.

A fixed constraint is applied to the bottom of the spring model, and the top plane is applied with a process pressure load with a maximum load . The dynamic load equation is

Figure 6 presents the height of each spring coil endpoint in the compression process. The location of the spring coil endpoints is shown in Figure 4, and the endpoint is fixed. The equilibrium state of the spring can be obtained by the spring model. This method can be used to calculate the static stiffness of the spring.

Figure 6 

Time evolution of the height of each spring coil endpoint in the compression process.

The shape of the graphs presented in Figure 7 is remarkable; it presents the torsion angle of the coil endpoint 4 in the end coil is the maximum, which is close to the fixed position.

Figure 7 

Time evolution of the torsion angle of each spring coil endpoint in the compression process.

Figure 8 shows the radius of the spring coil endpoints in the middle is greater than that of the endpoints at both ends. The analysis results are in line with the actual ones.

Figure 8 

Time evolution of the radius of each spring coil endpoint in the compression process.

3.2. Performance of the Noise Reduction Factor

The geometrical characteristics, material properties, and stress conditions of the variable pitch spring are consistent with the first example. The static equilibrium equation can be obtained from the system equation (53) by letting . It is calculated by using FSOLVE function in the MATLAB package.

The calculation results with the different noise reduction factor and the static equilibrium equation are shown in Table 2.

In Table 2, it is noticed that the influence of the noise reduction factor on the calculation result is very small, but the influence on the calculation time is obvious. The convergence speed of the model calculation using the static equilibrium equation is not as fast as that of the dynamic calculation model with the noise reduction factor. Setting the noise reduction factor can calculate the displacement of the spring at equilibrium. The possibility of spring design can be guaranteed by fast simulation.

3.3. Verification of the Helical Spring Element

This example is used to calculate quickly the displacement of the end coil at static equilibrium with the noise reduction factor and illustrate the correctness of the geometrically exact helical spring element. The studied spring is a helical spring with a circular cross-section, which has the mechanical and geometrical characteristics presented in Table 3.

The bottom of the spring is fixed, and an axial compression load with 1000 N is applied to the free extremity of the spring. The displacements of a spring under static compression load are given in Table 4, compared with ANSYS results.

In Table 4, it is found that when the pitch of the spring is 40 mm, the error is the smallest. With the increase of the pitch, the error value gradually increases. But the overall error is not very large. The correctness of the method is proved by the data in the table.

3.4. A Numerical Example of the Spring Static Stiffness

The static stiffness of the spring can be obtained by a static equilibrium to the spring element. Therefore, the dynamic geometrically exact helical spring model is used to solve the spring stiffness design by applying the dynamic load equation (56). The helical spring model is combined with an optimization program that defines the objective function and the constraint equation. These programs are based on three MATLAB optimization algorithms which are “trust-region-reflective” (TR), “Levenberg–Marquardt” (LM), and “sequential-quadratic-program” (SQP).

In this case, the valve spring installation requirements of an automobile engine are selected as the spring stiffness design objective. When the valve is closed, the compression distance of the spring is 10 mm and the preloading load is 120 N. When the valve is open to the maximum position, the compression distance of the spring is 17 mm and the compression load is 204 N. The objective is to design the helical spring with the linear stiffness, namely, least square method. After the setting of an initial point for each algorithm, the stiffness of spring is optimized. The active coil number is 4, which is supposed fixes. These parameters of the spring are the wire diameter , the middle helix diameter , and the helix pitches , , ,. The least square method is used to construct the objective function, namely, the square of the difference between the calculated height of the spring and the objective height under the six test loads which should be in the range . The designed parameters are shown in Figure 9.

Figure 9 

Drawings of the designed spring.

This leads to the generic form of the presented optimization problem, which can be formulated as follows.Findto minimizesubject to

The material properties of the spring designed are presented in Table 5.

The design results are shown in Table 6. The corresponding design variables, number of iterations, and calculation time for each algorithm are listed.

In Table 6, it is noticed that the TR algorithm quickly gives the spring design parameters in accordance with the target stiffness. The design parameters need to be brought into the system equation to verify that the stiffness characteristics are met. The designed stiffness of springs for each algorithm is presented in Figure 10.

Figure 10 

Designed stiffness of springs for each algorithm.

Figure 10 shows that the calculated results of the TR algorithm for this case are the most consistent with the target curve. Meanwhile, the calculated values of the other two algorithms are basically consistent with the target stiffness curve. Figure 11 shows the model of the designed spring. Figure 12 shows the height of each endpoint of the spring designed by the TR algorithm under the maximum compression load of 204 N.

Figure 11 

Model of the designed spring.

Figure 12 

Time evolution of each coil endpoint height of the spring designed by TR algorithm in the compression process.

According to the data in Figure 12, the minimum distance between the two adjacent coils is 5.39 mm, larger than the wire diameter of the spring. Meanwhile, no clash occurs between the two adjacent coils during the compression process of the spring. This is to meet the spring stiffness design requirements. So the optimal design is as follows: the wire diameter , the middle helix diameter , and the helix pitches , , , .

3.5. A Numerical Example of the Spring Dynamic Response

A spring with variable radii is subjected to a tensile load applied in the axial direction of spring. The geometrical properties and material properties of the spring are presented in Table 7. The model of the spring with variable radii is shown in Figure 13.

Figure 13 

Model of a spring with variable radii.

The helix radius function of the spring is written asin which and .