High-Efficiency Dynamic Modeling of a Helical Spring Element Based on the Geometrically Exact Beam Theory

)is 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.)e 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.)is study also carries out five numerical trials to validate the above dynamic procedure of the helical spring element. )e example of the spring static stiffness design shows that the proposed helical spring procedure enables one to deal with practical engineering applications.


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. e 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. e element stiffness and mass matrix were based on an exact differential equation that governs the static behavior of infinitesimal elements. e 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. e theory for calculating the dynamic stiffness matrix and the natural frequency was also proposed. Lee and ompson [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. e 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 twonode 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. e helical spring is an elongated structure having an initial helix curvature. And it undergoes large displacements and finite rotations. e 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][11][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. e 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. erefore, 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. e 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. e 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.

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.

Equations of the Helical Spring.
e helical spring is three-dimensional and has an initial helix centroid line, defined in the global coordinate system g 1 , g 2 , g 3 . e 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 g 3 forms a righthanded helix. In this paper, a superposed 'bar' denotes the parameters in the initial configuration. ρ, θ, and z represent the helix radius, polar angle, and height of the line of centroids in the initial configuration. e vectors of the helix with the global coordinate system g 1 , g 2 , g 3 can be defined by where ρ, θ, z indicate the helix radius, polar angle, and height, respectively, in the current configuration. A family of orthonormal vector bases e t , e b , e s on the helix spring is introduced to describe the orientation of the cross-sections. e t and e b are oriented along the principal axes of inertia of the cross-section, respectively, and e s is normal to the cross-section: e s � e t × e b . For the Euler-Bernoulli beam assumption, the cross-section remains perpendicular to the tangent vector of the centerline during deformation. e relationship between the global base vector g i and cross-section base vector e i can be expressed as where R is the rotation matrix with CARDAN angles α, β, and c as follows: Here, c i and s i denote sine and cosine of Cardan angles, respectively.

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. P i is the endpoint in the i-th circle. e parameters of the spring are changed by the stress, but the spring remains helical. e 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 θ. ρ i is the helix radius of the endpoint P i . By employing the cubic spline interpolation method to discretize the configuration curve, the centroid vector at any point of the spiral can be obtained. e helix radius in the i-th circle, shown in Figure 3, is described as For spline interpolation, the first derivative of the radius at the endpoint to the initial polar coordinate is expressed as where [ρ i ] � ρ 0 ρ 1 ... ρ n and S 0 and S 1 are the constant matrix. In this way, the first derivative of the intermediate node of the spring element can be expressed at both ends. us, the helix radius corresponding to any point on the spiral is where T i is the shape function matrix as follows: Shock and Vibration 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 e generalized coordinates of the spring are given by e 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.

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 crosssection 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 to where e 1 � cos θg 1 + sin θg 2 and e 2 � −sin θg 1 + cos θg 2 . e arc length of an infinitesimal spatial element is defined by According to expression (16), the derivative of the arclength coordinate with respect to the initial polar coordinate can be expressed by Expression (14) becomes with expressions (15) and (17) and we also have 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 0 < β < π/2 and avoid the singular position of the CARDAN angles to describe the rigid-body rotation effectively. Since cos(β) > 0, the CARDAN angle α can be derived as 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 by   Shock and Vibration Substituting expression (19) into expressions (22), we can solve for α ′ and β ′ . e rotation of a rigid cross-section is described with the help of axes of rotation and corresponding CARDAN angles. e vectors of rotation axes are e curvature equations provide the relationship of the curvature of the rigid cross-section fixed system relative to the initial system. e curvature vector can be expressed by where ψ ≜ �������������� (z ′ ) 2 + (ρ ′ ) 2 + ρ 2 . e components of κ in the rigid cross-section fixed system e t , e b , e s are extracted from this in the following: where A � Similar to the description of the current configuration, the components of κ in the initial configuration e t , e b , e s are extracted from this in the following: where A � e curvature strain measures are then extracted from above, leading to bending strains κ t − κ t and κ b − κ b and torsion strain κ s − κ s , respectively.
According to the definition of strain, the normal strain is given by

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 by where T r � [e 1 T i , g 3 T i , ρe 2 T i , 0 (3×n) ]. e spring element generalized velocity and its virtual variation are given by e 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 ψ are where e time derivatives of the normal strain ε s and its virtual variation can be expressed by e time derivatives of e s can be expressed by where where e second-order time derivatives of the CARDAN angles can be expressed by

angular velocity of the cross-section is
Shock and Vibration e components of ω in rigid cross-section fixed configuration e t , e b , e s are extracted from this in the following: and its virtual variation can be written as δω s e angular acceleration is where e time derivatives and virtual variation of κ can be expressed by (40)

Inertial Force Virtual
Power. e inertial force virtual power of the rigid section includes the translational inertial force virtual power and the rotational inertial force virtual power. e translational inertial force virtual power is given by where ρ is spring spiral density. e rotational inertial force virtual power is given by where J is rigid-section inertia moment tensor.

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 [0, 2nπ], can be written as where n is spring coil number and δw e is the internal virtual power of the infinitesimal element which is shown in Appendix. e constitutive equation of the spring is defined as where σ � [f s , m s , m t , m b ] T is the stress array and T is the strain array. e components of σ and ε are given relative to the cross-section basis e t , e b , e s . e constitutive matrix is where EA is the axial stiffness, GJ is the torsion stiffness, and EI t , EI b 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 (47) e spring system dynamics equation has strong rigidity, while the high-frequency elastic vibration mainly comes from the rapid change of strain. e model noise reduction method is to average the strain in a time interval (t, t + h).
e h 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 6 Shock and Vibration added, so that the natural frequency of the spring system can be reduced and higher calculation efficiency can be obtained.

Governing Equation.
e method of applying a test external force in the axial direction of the spring is adopted in this paper. e external force virtual power includes the gravitational virtual power and the test force's virtual power. e gravitational virtual power can be expressed as e test external force's virtual power can be written as e axial test external load f n is shifted to the action node and a corresponding torque is generated as m n � ρ n g 1 cos θ n + g 2 sin θ n × f n , where ρ n is the mean coil radius at the end node. e virtual power equation of the spring element system can be written as e mass matrix of the system is defined by e generalized force vector is written as e governing equation of the spring system is calculated by Gaussian integration directly.

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. e governing equation (53) is solved using an explicit Runge-Kutta method (ODE45) implemented in MATLAB. e relative error tolerance is ε r � 1 × 10 − 3 , the absolute error tolerance is ε a � 1 × 10 − 6 , and the maximum step size is Δt � 1 s.

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. e geometrical properties and material properties of the spring are presented in Table 1.
e curvature of the spring with variable pitches in the initial configuration is shown in Figure 5, which is modeled by the proposed spring element.
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 f max � 200 N . e dynamic load equation is (56) Figure 6 presents the height of each spring coil endpoint in the compression process. e location of the spring coil endpoints is shown in Figure 4, and the endpoint P 0 is fixed. e equilibrium state of the spring can be obtained by the spring model. is method can be used to calculate the static stiffness of the spring. e 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 8 shows the radius of the spring coil endpoints in the middle is greater than that of the endpoints at both ends. e analysis results are in line with the actual ones.

Performance of the Noise Reduction
Factor. e geometrical characteristics, material properties, and stress conditions of the variable pitch spring are consistent with the first example. e static equilibrium equation can be obtained from the system equation (53) by letting _ q � 0, € q � 0. It is calculated by using FSOLVE function in the MATLAB package. e 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. e 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 Shock and Vibration 7 the noise reduction factor h � 0.1s can calculate the displacement of the spring at equilibrium. e possibility of spring design can be guaranteed by fast simulation.

Verification of the Helical Spring Element.
is example is used to calculate quickly the displacement of the end coil at static equilibrium with the noise reduction factor h � 0.1s and illustrate the correctness of the geometrically exact helical spring element. e studied spring is a helical spring with a circular cross-section, which has the mechanical and geometrical characteristics presented in Table 3. e bottom of the spring is fixed, and an axial compression load with 1000 N is applied to the free extremity of the spring. e 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. e correctness of the method is proved by the data in the table.

A Numerical Example of the Spring Static Stiffness.
e static stiffness of the spring can be obtained by a static equilibrium to the spring element. erefore, the dynamic geometrically exact helical spring model is used to solve the spring stiffness design by applying the dynamic load equation (56). e helical spring model is combined with an optimization program that defines the objective function and the constraint equation. ese programs are based on three MATLAB optimization algorithms which are "trustregion-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. e 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. e active coil number is 4, which is supposed fixes. ese parameters of the spring are the wire diameter d, the middle helix diameter D, and the helix pitches p 1 , p 2 , p 3 ,p 4 . e least square method is used to construct the objective function, namely, the square of the difference between the calculated height h i of the spring and the objective height h i under the six test loads which should be in the range [120 N, 204 N]. e designed parameters are shown in Figure 9.
is leads to the generic form of the presented optimization problem, which can be formulated as follows.Find subject to e material properties of the spring designed are presented in Table 5.
e design results are shown in Table 6. e 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. e design parameters need to be brought into the system equation to verify that the stiffness characteristics are met. e designed stiffness of springs for each algorithm is presented in Figure 10. 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.
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. is is to meet the spring stiffness  x 1 x 6 x 5 x 4 x 3 Figure 9: Drawings of the designed spring.

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. e geometrical properties and material properties of the spring are  presented in Table 7. e model of the spring with variable radii is shown in Figure 13. e helix radius function of the spring is written as in which ρ 0 � 0.03 mm and ρ 6 � 0.02 mm. e endpoint of the spring element is fixed, and another endpoint is applied with a tensile load with a maximum load f max � 100 N. e load function is e computations are carried out by using the ODE45 with the noise reduction factor h � 0.001. e computed response of each spring coil endpoints is shown in Figure 14. It is found that the amplitude of the spring coil away from the fixed position is the maximum in the free vibration of the spring. is is consistent with the facts. It shows that the spring model is also suitable for the simulation of spring dynamic response.

Conclusions
is paper describes in detail the dynamics of a geometrically exact beam theory based helical spring element with variable pitch and variable radius using the Euler-Bernoulli beam assumption. Five numerical examples with different scales demonstrate that the proposed helical spring element is effective and efficient for both the static stiffness calculation and dynamic response analysis. e first example shows that the structural characteristic of dynamics helical spring element with geometrically exact beam theory can be simulated successfully. e effects of the noise reduction factor on the calculation result and time are also studied in the second example. e comparison results of the third example verify the correctness of the proposed helical spring element. en, the method is applied to the spring static stiffness design, and the design parameters of spring are proposed to deal with practical engineering applications. e final example shows the amplitude response of the spring with variable radii in the free vibration. In future work, the optimization design of the nonlinear dynamic response of the spring in a multibody system will be studied.