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

Authors: Jian Zhang; Zhaohui Qi; Gang Wang; Shudong Guo
Journal: Shock and Vibration
 (2020)

---

## 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.


---

## Body

## 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 Section2, 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 systemg1,g2,g3. 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 g3 forms a right-handed helix. In this paper, a superposed ‘bar’ denotes the parameters in the initial configuration. ρ¯,θ¯,andz¯ 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 systemg1,g2,g3 can be defined by(1)r=ρcosθg1+sinθg2+zg3,where ρ,θ,z indicate the helix radius, polar angle, and height, respectively, in the current configuration.A family of orthonormal vector baseset,eb,es on the helix spring is introduced to describe the orientation of the cross-sections. et and eb are oriented along the principal axes of inertia of the cross-section, respectively, and es is normal to the cross-section: es=et×eb. 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 gi and cross-section base vector ei can be expressed as(2)etebes=Rg3g1g2,where R is the rotation matrix with CARDAN angles α,β,andγ as follows:(3)R=c2c3s1s2c3+c1s3−c1s2c3+s1s3−c2s3−s1s2s3+c1c3c1s2s3+s1c3s2−s1c2c1c2.Here,ciandsi 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. Figure2 shows a schematic of the helix spring in the initial configuration and the current configuration with deformation. Pi 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 Figure2. 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 Pi. 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 thei-th circle, shown in Figure 3, is described as(4)ρ=N1ρi+N2ρi+1+N3ρi′+N4ρi+1′,where(5)N1=1−ξi21+2ξi,N2=ξi23−2ξi,N3=2π1−ξi2ξi,N4=2πξi2ξi−1are the Hermite shape functions with ξi=θ¯−θ¯i/2π, ξi∈0,1.For spline interpolation, the first derivative of the radius at the endpoint to the initial polar coordinate is expressed as(6)ρi′=ρiS0+ρ0′ρn′S1,where ρi=ρ0ρ1...ρn and S0 and S1 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 is(7)ρ=Tiρ,where Ti is the shape function matrix as follows:(8)Ti=N1ei+N2ei+1+N3S0TS1Ti.+N4S0TS1Ti+1..In equation (7), the vector ρ of the endpoint coordinates is defined as(9)ρ=ρ0ρ1...ρnρ0′ρn′T.In the same way, the helix height, polar angle, and torsion angle of the wire section to any point, respectively, are written as(10)z=Tiz,(11)θ=Tiθ,(12)γ=Tiγ.The generalized coordinates of the spring are given by(13)q=ρ,z,θ,γT.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(14)es=∂r∂s=r′s′.Taking partial derivative of (1) with respect to the initial polar coordinate leads to(15)r′=ρ′e1+ρθ′e2+z′g3,where e1=cosθg1+sinθg2 and e2=−sinθg1+cosθg2.The arc length of an infinitesimal spatial element is defined by(16)ds=drdθdθ.According to expression (16), the derivative of the arc-length coordinate with respect to the initial polar coordinate can be expressed by(17)s′=r′=z′2+ρ′2+ρ2θ′2≜ψ.Expression (14) becomes with expressions (15) and (17)(18)es=ψ−1ρ′e1+ρθ′e2+z′g3,and we also have(19)es′=ψ−1b1e1+b2e2+z″g3−ψ′es,where b1=ρ″−ρθ′2, b2=ρ′θ′+ρ′θ′+ρθ″.Combining the relationship between the CARDAN angles and the elements in the rotation matrixR, we get(20)β=asing3Tes=asinz′ψ.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(21)α=atan2−g1Tes,g2Tes,where −g1Tes=−ψ−1ρ′cosθ−ρθ′sinθ and g2Tes=ψ−1ρ′sinθ+ρθ′cosθ.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(22)α′=λαTes′,β′=λβTes′,where λαT=−c2−1c1g1+s1g2T and λβT=c2−1g3TSubstituting 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(23)p1=s2es+c2c3et−c2s3eb,p2=s3et+c3eb,p3=es,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 by(24)κ=ψ¯−1α′p1+β′p2+γ′p3,where ψ¯≜z¯′2+ρ¯′2+ρ¯2.The components ofκ in the rigid cross-section fixed system et,eb,es are extracted from this in the following:(25)κsκtκb=ψ¯−1Aα′,where A=s201c2c3s30−c2s3c30 and the CARDAN angle array can be expressed by α=αβγ.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:(26)κ¯sκ¯tκ¯b=ψ¯−1A¯α¯′,where A¯=s¯201c¯2c¯3s¯30−c¯2s¯3c¯30and 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κt−κ¯t and κb−κ¯b and torsion strain κs−κ¯s, respectively.According to the definition of strain, the normal strain is given by(27)εs=dsds¯−1=ψ¯−1ψ−1.


### 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 by(28)r˙=Trq˙,δr˙=Trδq˙,where Tr=e1Ti,g3Ti,ρe2Ti,03×n.The spring element generalized velocity and its virtual variation are given by(29)q˙=ρ˙,z˙,θ˙,γ˙T,δq˙=δρ˙,δz˙,δθ˙,δγ˙T.The accelerated velocity of the vectors of the helix can be obtained by the time derivation of(30)r¨=Trq¨+αr,αr=2ρ˙e2−ρθ˙e1θ˙,In order to facilitate the following derivation, the time derivative and virtual variation ofψ are(31)ψ˙=Tψq˙,δψ˙=Tψδq˙,where Tψ=Fρ,Fz,Fθ′,01×n, Fz=z′/ψTi′, Fρ=ρ′/ψTi′+ρθ′2/ψTi, and Fθ′=ρ2θ′/ψTi′.The time derivatives of the normal strainεs and its virtual variation can be expressed by(32)ε˙s=ψ¯−1Tψq˙,δε˙s=ψ¯−1Tψδq˙.The time derivatives ofes can be expressed by(33)e˙s=Tesq˙,where Tes=Nρ,Nz,Nθ,03×n, Nρ=e1Ti′+e2θ′Ti−esFρ/ψ and Nz=g3Ti′−esFz/ψ, Nθ=−ρθ′e1Ti+e2ρ′Ti+ρTi′−eFθ′/ψ.The time derivatives of the CARDAN angles and its virtual variation are(34)α˙=Tαq˙,δα˙=Tαδq˙,where Tα=λαTTesλβTTesTiNγ.The second-order time derivatives of the CARDAN angles can be expressed by(35)α¨=Tαq¨+αα,where αα=λ˙αTe˙s+λαTαesλ˙βTe˙s+λβTαes01×n, λ˙αT=∂λαT/∂αα˙+∂λαT/∂ββ˙, ∂λαT/∂α=c2−1s1g1−c1g2T, ∂λαT/∂β=c2−1s2λαT, λ˙βT=∂λβT/∂ββ˙, and ∂λβT/∂β=s2c2−1λβT.The angular velocity of the cross-section is(36)ω=α˙p1+β˙p2+γ˙p3.The components ofω in rigid cross-section fixed configuration et,eb,es are extracted from this in the following:(37)ωsωtωb=Aα˙=ATαq˙=Tωq˙,and its virtual variation can be written as(38)δωsδωtδωb=Tωδq˙.The angular acceleration is(39)ω˙=Aα¨+A˙α˙=Tωq¨+αω,αω=Aαα+A˙α˙,where A˙=∂A/∂ββ˙+∂A/∂γγ˙, ∂A/∂β=c200−s2c300s2s300, and ∂A/∂γ=000−c2s3c30−c2c3−s30.The time derivatives and virtual variation ofκ can be expressed by(40)κ˙=ψ¯−1A˙α′+Aα′˙=ψ¯−1Tκq˙,δκ˙=ψ¯−1Tκδq˙.


### 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 by(41)δpv=ψ¯∫02nπδr˙Tr¨ρ^dθ¯≜δq˙TMvq¨+fv,Mv=ψ¯∫02nπδTrTTrρ^dθ¯,fv=ψ¯∫02nπδTrTαrρ^dθ¯,where ρ^ is spring spiral density.The rotational inertial force virtual power is given by(42)δpω=ψ¯∫02nπδωTJω˙+ω˜Jωρ^dθ¯≜δq˙TMωq¨+fω,Mω=ψ¯∫02nπδTωTJTωρ^dθ¯,fω=ψ¯∫02nπδTωTJαω+ω˜Jωρ^dθ¯,where J 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 range0,2nπ, can be written as(43)δpe=∫02nπδwedθ¯=ψ¯∫02nπδε˙sfs+δκ˙sms+δκ˙tmt+δκ˙bmbdθ¯,where n is spring coil number and δwe is the internal virtual power of the infinitesimal element which is shown in Appendix.The constitutive equation of the spring is defined as(44)σ=Dε,where σ=fs,ms,mt,mbT is the stress array and ε=εs,κs−κ¯s,κt−κ¯t,κb−κ¯bT is the strain array. The components of σ and ε are given relative to the cross-section basis et,eb,es.The constitutive matrix is(45)D=digEA,GJ,EIt,EIb,where EA is the axial stiffness, GJ is the torsion stiffness, and EIt, EIb are the principal bending stiffnesses.With the constitutive equation of the spring, the reduced expression of the internal power is given by(46)δpe=ψ¯∫02nπδε˙Tσdθ¯.In terms of (32) and (40), the internal virtual power equation (46) can be written as(47)δpe=ψ¯∫02nπδε˙Tσdθ¯≜δq˙Tfe,fe=∫02nπTψT+TκTσdθ¯.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 intervalt,t+h. The 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(48)σ¯≈σ+hσ˙2+h2σ¨6.Modified internal virtual power is given by(49)δpe≈ψ¯∫02nπδε˙Tσ¯dθ¯≜δq˙TMeq¨+Ceq˙+fe,Me=ψ¯h2D6∫02nπTψT+TκTTψ+Tκdθ¯,Ce=ψ¯D∫02nπTψT+TκTh2Tψ+Tκ+h26T˙ψ+T˙κdθ¯.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(50)δpg=−ψ¯∫02nπδr˙Tgρ^dθ¯≜−δq˙Tfg.The test external force’s virtual power can be written as(51)δpa=−δr˙nTfn−δωnTmn≜−δq˙Tfa.The axial test external loadfn is shifted to the action node and a corresponding torque is generated as(52)mn=ρng1cosθn+g2sinθn×fn,where ρn is the mean coil radius at the end node.The virtual power equation of the spring element system can be written as(53)δq˙TMq¨+F=0.The mass matrix of the system is defined by(54)M=Mv+Mω+Me.The generalized force vector is written as(55)F=fv+fω+fg+fa+fe+Ceq˙.The governing equation of the spring system is calculated by Gaussian integration directly.


## 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 systemg1,g2,g3. 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 g3 forms a right-handed helix. In this paper, a superposed ‘bar’ denotes the parameters in the initial configuration. ρ¯,θ¯,andz¯ 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 systemg1,g2,g3 can be defined by(1)r=ρcosθg1+sinθg2+zg3,where ρ,θ,z indicate the helix radius, polar angle, and height, respectively, in the current configuration.A family of orthonormal vector baseset,eb,es on the helix spring is introduced to describe the orientation of the cross-sections. et and eb are oriented along the principal axes of inertia of the cross-section, respectively, and es is normal to the cross-section: es=et×eb. 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 gi and cross-section base vector ei can be expressed as(2)etebes=Rg3g1g2,where R is the rotation matrix with CARDAN angles α,β,andγ as follows:(3)R=c2c3s1s2c3+c1s3−c1s2c3+s1s3−c2s3−s1s2s3+c1c3c1s2s3+s1c3s2−s1c2c1c2.Here,ciandsi 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. Figure2 shows a schematic of the helix spring in the initial configuration and the current configuration with deformation. Pi 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 Figure2. 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 Pi. 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 thei-th circle, shown in Figure 3, is described as(4)ρ=N1ρi+N2ρi+1+N3ρi′+N4ρi+1′,where(5)N1=1−ξi21+2ξi,N2=ξi23−2ξi,N3=2π1−ξi2ξi,N4=2πξi2ξi−1are the Hermite shape functions with ξi=θ¯−θ¯i/2π, ξi∈0,1.For spline interpolation, the first derivative of the radius at the endpoint to the initial polar coordinate is expressed as(6)ρi′=ρiS0+ρ0′ρn′S1,where ρi=ρ0ρ1...ρn and S0 and S1 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 is(7)ρ=Tiρ,where Ti is the shape function matrix as follows:(8)Ti=N1ei+N2ei+1+N3S0TS1Ti.+N4S0TS1Ti+1..In equation (7), the vector ρ of the endpoint coordinates is defined as(9)ρ=ρ0ρ1...ρnρ0′ρn′T.In the same way, the helix height, polar angle, and torsion angle of the wire section to any point, respectively, are written as(10)z=Tiz,(11)θ=Tiθ,(12)γ=Tiγ.The generalized coordinates of the spring are given by(13)q=ρ,z,θ,γT.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(14)es=∂r∂s=r′s′.Taking partial derivative of (1) with respect to the initial polar coordinate leads to(15)r′=ρ′e1+ρθ′e2+z′g3,where e1=cosθg1+sinθg2 and e2=−sinθg1+cosθg2.The arc length of an infinitesimal spatial element is defined by(16)ds=drdθdθ.According to expression (16), the derivative of the arc-length coordinate with respect to the initial polar coordinate can be expressed by(17)s′=r′=z′2+ρ′2+ρ2θ′2≜ψ.Expression (14) becomes with expressions (15) and (17)(18)es=ψ−1ρ′e1+ρθ′e2+z′g3,and we also have(19)es′=ψ−1b1e1+b2e2+z″g3−ψ′es,where b1=ρ″−ρθ′2, b2=ρ′θ′+ρ′θ′+ρθ″.Combining the relationship between the CARDAN angles and the elements in the rotation matrixR, we get(20)β=asing3Tes=asinz′ψ.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(21)α=atan2−g1Tes,g2Tes,where −g1Tes=−ψ−1ρ′cosθ−ρθ′sinθ and g2Tes=ψ−1ρ′sinθ+ρθ′cosθ.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(22)α′=λαTes′,β′=λβTes′,where λαT=−c2−1c1g1+s1g2T and λβT=c2−1g3TSubstituting 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(23)p1=s2es+c2c3et−c2s3eb,p2=s3et+c3eb,p3=es,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 by(24)κ=ψ¯−1α′p1+β′p2+γ′p3,where ψ¯≜z¯′2+ρ¯′2+ρ¯2.The components ofκ in the rigid cross-section fixed system et,eb,es are extracted from this in the following:(25)κsκtκb=ψ¯−1Aα′,where A=s201c2c3s30−c2s3c30 and the CARDAN angle array can be expressed by α=αβγ.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:(26)κ¯sκ¯tκ¯b=ψ¯−1A¯α¯′,where A¯=s¯201c¯2c¯3s¯30−c¯2s¯3c¯30and 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κt−κ¯t and κb−κ¯b and torsion strain κs−κ¯s, respectively.According to the definition of strain, the normal strain is given by(27)εs=dsds¯−1=ψ¯−1ψ−1.


## 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 by(28)r˙=Trq˙,δr˙=Trδq˙,where Tr=e1Ti,g3Ti,ρe2Ti,03×n.The spring element generalized velocity and its virtual variation are given by(29)q˙=ρ˙,z˙,θ˙,γ˙T,δq˙=δρ˙,δz˙,δθ˙,δγ˙T.The accelerated velocity of the vectors of the helix can be obtained by the time derivation of(30)r¨=Trq¨+αr,αr=2ρ˙e2−ρθ˙e1θ˙,In order to facilitate the following derivation, the time derivative and virtual variation ofψ are(31)ψ˙=Tψq˙,δψ˙=Tψδq˙,where Tψ=Fρ,Fz,Fθ′,01×n, Fz=z′/ψTi′, Fρ=ρ′/ψTi′+ρθ′2/ψTi, and Fθ′=ρ2θ′/ψTi′.The time derivatives of the normal strainεs and its virtual variation can be expressed by(32)ε˙s=ψ¯−1Tψq˙,δε˙s=ψ¯−1Tψδq˙.The time derivatives ofes can be expressed by(33)e˙s=Tesq˙,where Tes=Nρ,Nz,Nθ,03×n, Nρ=e1Ti′+e2θ′Ti−esFρ/ψ and Nz=g3Ti′−esFz/ψ, Nθ=−ρθ′e1Ti+e2ρ′Ti+ρTi′−eFθ′/ψ.The time derivatives of the CARDAN angles and its virtual variation are(34)α˙=Tαq˙,δα˙=Tαδq˙,where Tα=λαTTesλβTTesTiNγ.The second-order time derivatives of the CARDAN angles can be expressed by(35)α¨=Tαq¨+αα,where αα=λ˙αTe˙s+λαTαesλ˙βTe˙s+λβTαes01×n, λ˙αT=∂λαT/∂αα˙+∂λαT/∂ββ˙, ∂λαT/∂α=c2−1s1g1−c1g2T, ∂λαT/∂β=c2−1s2λαT, λ˙βT=∂λβT/∂ββ˙, and ∂λβT/∂β=s2c2−1λβT.The angular velocity of the cross-section is(36)ω=α˙p1+β˙p2+γ˙p3.The components ofω in rigid cross-section fixed configuration et,eb,es are extracted from this in the following:(37)ωsωtωb=Aα˙=ATαq˙=Tωq˙,and its virtual variation can be written as(38)δωsδωtδωb=Tωδq˙.The angular acceleration is(39)ω˙=Aα¨+A˙α˙=Tωq¨+αω,αω=Aαα+A˙α˙,where A˙=∂A/∂ββ˙+∂A/∂γγ˙, ∂A/∂β=c200−s2c300s2s300, and ∂A/∂γ=000−c2s3c30−c2c3−s30.The time derivatives and virtual variation ofκ can be expressed by(40)κ˙=ψ¯−1A˙α′+Aα′˙=ψ¯−1Tκq˙,δκ˙=ψ¯−1Tκδq˙.


## 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 by(41)δpv=ψ¯∫02nπδr˙Tr¨ρ^dθ¯≜δq˙TMvq¨+fv,Mv=ψ¯∫02nπδTrTTrρ^dθ¯,fv=ψ¯∫02nπδTrTαrρ^dθ¯,where ρ^ is spring spiral density.The rotational inertial force virtual power is given by(42)δpω=ψ¯∫02nπδωTJω˙+ω˜Jωρ^dθ¯≜δq˙TMωq¨+fω,Mω=ψ¯∫02nπδTωTJTωρ^dθ¯,fω=ψ¯∫02nπδTωTJαω+ω˜Jωρ^dθ¯,where J 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 range0,2nπ, can be written as(43)δpe=∫02nπδwedθ¯=ψ¯∫02nπδε˙sfs+δκ˙sms+δκ˙tmt+δκ˙bmbdθ¯,where n is spring coil number and δwe is the internal virtual power of the infinitesimal element which is shown in Appendix.The constitutive equation of the spring is defined as(44)σ=Dε,where σ=fs,ms,mt,mbT is the stress array and ε=εs,κs−κ¯s,κt−κ¯t,κb−κ¯bT is the strain array. The components of σ and ε are given relative to the cross-section basis et,eb,es.The constitutive matrix is(45)D=digEA,GJ,EIt,EIb,where EA is the axial stiffness, GJ is the torsion stiffness, and EIt, EIb are the principal bending stiffnesses.With the constitutive equation of the spring, the reduced expression of the internal power is given by(46)δpe=ψ¯∫02nπδε˙Tσdθ¯.In terms of (32) and (40), the internal virtual power equation (46) can be written as(47)δpe=ψ¯∫02nπδε˙Tσdθ¯≜δq˙Tfe,fe=∫02nπTψT+TκTσdθ¯.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 intervalt,t+h. The 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(48)σ¯≈σ+hσ˙2+h2σ¨6.Modified internal virtual power is given by(49)δpe≈ψ¯∫02nπδε˙Tσ¯dθ¯≜δq˙TMeq¨+Ceq˙+fe,Me=ψ¯h2D6∫02nπTψT+TκTTψ+Tκdθ¯,Ce=ψ¯D∫02nπTψT+TκTh2Tψ+Tκ+h26T˙ψ+T˙κdθ¯.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(50)δpg=−ψ¯∫02nπδr˙Tgρ^dθ¯≜−δq˙Tfg.The test external force’s virtual power can be written as(51)δpa=−δr˙nTfn−δωnTmn≜−δq˙Tfa.The axial test external loadfn is shifted to the action node and a corresponding torque is generated as(52)mn=ρng1cosθn+g2sinθn×fn,where ρn is the mean coil radius at the end node.The virtual power equation of the spring element system can be written as(53)δq˙TMq¨+F=0.The mass matrix of the system is defined by(54)M=Mv+Mω+Me.The generalized force vector is written as(55)F=fv+fω+fg+fa+fe+Ceq˙.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 εr=1×10−3, the absolute error tolerance is εa=1×10−6, and the maximum step size is Δt=1s.


### 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 Figure4. 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.Table 1
Geometrical properties and material properties of the variable pitch spring.
PropertiesYoung’s modulusPoisson coefficientDensityMiddle helix diameterWire diameterActive coil numberPitches (m)2.10 × 1011 Pa0.37800 kg/m30.040.0044[0.01, 0.02, 0.03, 0.04]The curvature of the spring with variable pitches in the initial configuration is shown in Figure5, which is modeled by the proposed spring element.Figure 5
Curvature of the variable pitch spring in the initial configuration. (a) The torsion curvatureκs of the spring in the initial state. (b) The bending curvature κt of the spring in the initial state. (c) The bending curvature κb of the spring in the initial state.
(a)(b)(c)As can be seen from Figure5, 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 loadfmax=200N . The dynamic load equation is(56)fa=fmaxt310−15t+6t2,0≤t≤1.5,fmax,t>1.5.Figure6 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 P0 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 Figure7 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.Figure8 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 q˙=0,q¨=0. 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 Table2.Table 2
Displacement of each coil with different noise reduction factors.
Factorh = 0.01 sh = 0.1 sh = 0.5 sStatics (mm)Coil 1 (mm)6.226.226.126.23Coil 2 (mm)20.9820.9820.7420.98Coil 3 (mm)45.3945.3945.0145.43Coil 4 (mm)81.6181.6181.1281.68CPU time (s)52.02949.1493.192153.53In Table2, 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 h=0.1s 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 factorh=0.1s 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.Table 3
Geometrical properties and material properties of the equal pitch spring.
PropertiesYoung’s modulusPoisson coefficientDensityMiddle helix diameter (m)Wire diameterActive coil number2.10 × 1011 Pa0.37800 kg/m30.1570.0136The 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 Table4, compared with ANSYS results.Table 4
Comparison of spring displacement results with ANSYS.
Pitch (mm)Solid95PresentDisplacement (mm)Displacement (mm)Error (%)4080.7180.450.326080.9079.681.548081.2378.673.2510081.4377.764.7212081.7077.105.9614082.3076.807.16In Table4, 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 diameterd, the middle helix diameter D, and the helix pitches p1, p2, p3,p4. The least square method is used to construct the objective function, namely, the square of the difference between the calculated height hi of the spring and the objective height h¯i under the six test loads which should be in the range 120N,204N. 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.Find(57)x=x1,x2,x3,x4,x5,x6=D,d,p1,p2,p3,p4,to minimize(58)y=∑16hi−h¯i2,subject to(59)40<x1<90,3<x2<5,10<x3,x4,x5,x6<14.The material properties of the spring designed are presented in Table5.Table 5
Material properties of the spring.
PropertiesYoung’s modulusPoisson coefficientDensity2.10 × 1011 Pa0.37800 kg/m3The design results are shown in Table6. The corresponding design variables, number of iterations, and calculation time for each algorithm are listed.Table 6
Corresponding design variables for each algorithm.
AlgorithmTRLMSQPInitial pointx1 (mm)808080x2 (mm)555x3–x6 (mm)[10, 10, 10, 10][10, 10, 10, 10][10, 10, 10, 10]Resultsx1 (mm)50.7780.3750.71x2 (mm)57.055x3–x6 (mm)[10.99, 10.99, 11.01, 10.91][11.38, 11.66, 11.72, 9.23][10.95, 11.15, 11.15, 10.67]Iteration101813CPU time3661 s6135 s5255 sIn Table6, 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.Figure10 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 Figure12, 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 d=5mm, the middle helix diameter D=50.77mm, and the helix pitches p1=10.99mm, p2=10.99mm, p3=11.01mm, p4=10.91mm.


### 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 Table7. The model of the spring with variable radii is shown in Figure 13.Table 7
Geometrical properties and material properties of the variable pitch spring.
PropertiesYoung’s modulusPoisson coefficientDensityWire diameter (m)Active coil numberPitch (m)2.10 × 1011 Pa0.37800 kg/m30.00460.0115Figure 13
Model of a spring with variable radii.The helix radius function of the spring is written as(60)ρi=ρ6−ρ6−ρ0i6,in which ρ0=0.03mm and ρ6=0.02mm.The endpoint of the spring element is fixed, and another endpoint is applied with a tensile load with a maximum loadfmax=100N. The load function is(61)fa=2fmaxt,0≤t≤0.5,0,t>0.5.The computations are carried out by using the ODE45 with the noise reduction factorh=0.001. The 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. This is consistent with the facts. It shows that the spring model is also suitable for the simulation of spring dynamic response.Figure 14
Time evolution of the displacement at each spring coil endpoint.


## 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 Figure4. 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.Table 1
Geometrical properties and material properties of the variable pitch spring.
PropertiesYoung’s modulusPoisson coefficientDensityMiddle helix diameterWire diameterActive coil numberPitches (m)2.10 × 1011 Pa0.37800 kg/m30.040.0044[0.01, 0.02, 0.03, 0.04]The curvature of the spring with variable pitches in the initial configuration is shown in Figure5, which is modeled by the proposed spring element.Figure 5
Curvature of the variable pitch spring in the initial configuration. (a) The torsion curvatureκs of the spring in the initial state. (b) The bending curvature κt of the spring in the initial state. (c) The bending curvature κb of the spring in the initial state.
(a)(b)(c)As can be seen from Figure5, 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 loadfmax=200N . The dynamic load equation is(56)fa=fmaxt310−15t+6t2,0≤t≤1.5,fmax,t>1.5.Figure6 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 P0 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 Figure7 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.Figure8 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 q˙=0,q¨=0. 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 Table2.Table 2
Displacement of each coil with different noise reduction factors.
Factorh = 0.01 sh = 0.1 sh = 0.5 sStatics (mm)Coil 1 (mm)6.226.226.126.23Coil 2 (mm)20.9820.9820.7420.98Coil 3 (mm)45.3945.3945.0145.43Coil 4 (mm)81.6181.6181.1281.68CPU time (s)52.02949.1493.192153.53In Table2, 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 h=0.1s 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 factorh=0.1s 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.Table 3
Geometrical properties and material properties of the equal pitch spring.
PropertiesYoung’s modulusPoisson coefficientDensityMiddle helix diameter (m)Wire diameterActive coil number2.10 × 1011 Pa0.37800 kg/m30.1570.0136The 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 Table4, compared with ANSYS results.Table 4
Comparison of spring displacement results with ANSYS.
Pitch (mm)Solid95PresentDisplacement (mm)Displacement (mm)Error (%)4080.7180.450.326080.9079.681.548081.2378.673.2510081.4377.764.7212081.7077.105.9614082.3076.807.16In Table4, 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 diameterd, the middle helix diameter D, and the helix pitches p1, p2, p3,p4. The least square method is used to construct the objective function, namely, the square of the difference between the calculated height hi of the spring and the objective height h¯i under the six test loads which should be in the range 120N,204N. 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.Find(57)x=x1,x2,x3,x4,x5,x6=D,d,p1,p2,p3,p4,to minimize(58)y=∑16hi−h¯i2,subject to(59)40<x1<90,3<x2<5,10<x3,x4,x5,x6<14.The material properties of the spring designed are presented in Table5.Table 5
Material properties of the spring.
PropertiesYoung’s modulusPoisson coefficientDensity2.10 × 1011 Pa0.37800 kg/m3The design results are shown in Table6. The corresponding design variables, number of iterations, and calculation time for each algorithm are listed.Table 6
Corresponding design variables for each algorithm.
AlgorithmTRLMSQPInitial pointx1 (mm)808080x2 (mm)555x3–x6 (mm)[10, 10, 10, 10][10, 10, 10, 10][10, 10, 10, 10]Resultsx1 (mm)50.7780.3750.71x2 (mm)57.055x3–x6 (mm)[10.99, 10.99, 11.01, 10.91][11.38, 11.66, 11.72, 9.23][10.95, 11.15, 11.15, 10.67]Iteration101813CPU time3661 s6135 s5255 sIn Table6, 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.Figure10 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 Figure12, 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 d=5mm, the middle helix diameter D=50.77mm, and the helix pitches p1=10.99mm, p2=10.99mm, p3=11.01mm, p4=10.91mm.


## 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 Table7. The model of the spring with variable radii is shown in Figure 13.Table 7
Geometrical properties and material properties of the variable pitch spring.
PropertiesYoung’s modulusPoisson coefficientDensityWire diameter (m)Active coil numberPitch (m)2.10 × 1011 Pa0.37800 kg/m30.00460.0115Figure 13
Model of a spring with variable radii.The helix radius function of the spring is written as(60)ρi=ρ6−ρ6−ρ0i6,in which ρ0=0.03mm and ρ6=0.02mm.The endpoint of the spring element is fixed, and another endpoint is applied with a tensile load with a maximum loadfmax=100N. The load function is(61)fa=2fmaxt,0≤t≤0.5,0,t>0.5.The computations are carried out by using the ODE45 with the noise reduction factorh=0.001. The 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. This is consistent with the facts. It shows that the spring model is also suitable for the simulation of spring dynamic response.Figure 14
Time evolution of the displacement at each spring coil endpoint.


## 4. Conclusions

This 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. The first example shows that the structural characteristic of dynamics helical spring element with geometrically exact beam theory can be simulated successfully. The effects of the noise reduction factor on the calculation result and time are also studied in the second example. The comparison results of the third example verify the correctness of the proposed helical spring element. Then, the method is applied to the spring static stiffness design, and the design parameters of spring are proposed to deal with practical engineering applications. The 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.



---

*Source: 8254606-2020-06-19.xml*