Discrete model for 3D motion of helical springs

Large displacements of coiled springs are affected by both contact between coils and geometric non-linearities. The two phenomena also interact with each other and the only technology that efficiently accounts for both is Finite Elements Method (FEM). Such analysis can become very expensive, especially when a great number of elements and contacts are involved, as it is often the case. A discrete 3D model for fast and efficient simulation of non-linear springs is presented. The discrete model is based on equivalent beam theory and each node is associated with a contact region placed on the helix. Thus contact is accounted on its real occurrence and fast non-linear solution procedure is implemented. A study case is presented and comparison with FEM is discussed.


Introduction
Contact between coils is often used as a mean to increase spring stiffness and provide additional damping. One of the most common applications of such strategy are valve springs in internal combustion engines and dual-rate or progressive springs in vehicle suspensions [1]. When coils come in contact, they almost stop deforming elastically, causing a reduction in the number of active coils and an increase in overall stiffness. Simulation of systems where non-linear springs are involved requires an appropriate modelling strategy. A full simulation of the coiled wire, e.g. using Finite Elements Method (FEM), would be too expensive to run and the correct managing of contact algorithm is non-trivial. Therefore, most commercial packages provide several tools to account for non-linear springs at low computational cost. A common approach is the multi mass model, where the spring is represented by several masses connected with a sequence of stiffness and damping elements [2]. The parameters are computed from coils properties, and non-linearity arises as the gaps between consecutive masses close. It is also possible to directly provide the non-linear spring characteristics, obtained from previous computations (analytically or via finite elements) or direct experiments. However, these methods only account for vertical motion and it seems to the authors that there is a lack in literature about the influence of coil contacts on flexural and lateral behaviour, while clash between coils not only affect the stiffness, but also changes the direction of spring motion and the volume occupied during operation. Starting from Haringx work on instability [3], several researchers have studied the lateral motion of helical springs. Analytical solutions have been obtained using equivalent beams under large displacements, [4] [5], but these models cannot manage coil contacts. On the other end, the dynamics of the wire has been addressed by many authors using numerical methods: finite element [6], transfer matrix method [7], finite difference [8], dynamic stiffness matrix method [9]. However, also in those cases it is still hard, and in some  cases unpracticable, to account for coil contacts and large deflections. A modelling strategy for two dimensional motion of non-linear springs with coil contact has been presented in [10]. The proposed model takes advantage of the equivalent beam formulation and represents the spring using several non-linear beam-like elements placed along the axis of the helix. Meanwhile, contact surfaces are related to the nodes and it is possible to introduce contact relations between adjacent coils. Such a model proved to be a good tool for quick and accurate prediction of the component behaviour and it is now extended to the general 3D case, where orientation of coils is defined by three consecutive finite rotations. The model is described in in section 2, focusing on rotation parametrization, kinematics of the elements and computation of nonlinear restoring forces. Contact is described in section 3 and the solution procedure is explained in 4. A numeric example is provided, considering FEM as benchmark in section 5 and finally conclusions are drawn.

3D equivalent beam 2.1. Outline of the model
The model is based on Haringx equivalent beam [3], [11], [12], a beam-column structure with elastic properties derived from the characteristics of the spring. Nevertheless, it is worth to mention two main assumptions in Haring's work that differ from classical theory of beam-column instability: (i) beam stiffness is kept constant, thus cross-section properties update accordingly with deformed length (ii) only bending rotations are considered for orientation of axial and shearing forces Thus it is possible to discretize the equivalent beam into a series of beam-like elements, that reflect previous assumptions, each node has six degrees of freedom (three translations and three rotations). The full helix in 3D can be reconstructed from nodal positions and orientations, in this manner it is possible to associate a contact zone on the helix to each node and verify the volume occupied by the spring in working conditions. A picture of the modelling strategy is given in figure 1.

Rotation parametrization
In 3D motion it is important to adopt a consistent set of rotation parameters that define the orientation of nodal frames of reference. For the present work it is found convenient to express the orientation through three consecutive rotations, as shown in figure 2: where e 1,J , e 2,J , e 3,J are the unit vectors of the moving frame of reference on node J when node I is assumed to be fixed. In the undeformed configuration all moving frames coincide with the fixed global one thus: Unit vectors in the deformed configuration are change as: where R(η, ζ, ξ) is the rotation matrix, given by [13]:

Coil kinematics and internal forces
The aim of this section is to show the procedure that allows to express generalized nodal forces as a function of generalized nodal displacements. According to Haringx assumptions, the deformation of the coils is expressed in terms of bending, twisting, extension/compression and shearing. Since it is assumed that the elastic properties of the equivalent beam are axisymmetric, bending can be expressed as a finite rotation about a generic axis R b (u, φ), that orient unit vector e 1,I towards e 1,J . Calling e 1,b the unit vector after a bending rotation is applied, it turns that: where Bending angle φ is computed as: φ = arccos( e 1,I · e 1,J ||e 1,I || ||e 1,J || ) (8) and bending axis u is: Bending rotation matrix is then [13]: where V φ = 1 − cos(φ). Usually (eq. 6) will not be satisfied for unit vectors e 2 and e 3 , thus a further rotation about e 1,J is required. The required angle ψ is the twisting angle, computed as either: The applied twisting rotation R t (e 1,J , ψ) is given again by (eq. 10) with appropriate substitutions. Of course, it is verified that: According to Haringx assumptions, the elongation and shearing directions are defined by bending rotations. For every element the elongation direction e b is computed considering the half-way bending rotation: The length l of the element is given by the component of nodal distance a along e b , and shearing v is the other component.
Obviously elongation ∆l of the element is given by: Having completely defined coil kinematics, it is possible to compute the elastic forces arising as a consequence of deformations. There are an axial force, f ax directed along e b , a shearing force f sh parallel to v, a bending moment m b applied on u and a twisting moment m t applied on e b : (20) Note that the magnitude of each force is directly proportional to the corresponding deformation through a stiffness constant that only depends on spring geometry and material properties. Their values are computed taking into account the most important component of stress in the wire [12] when an ideal flat coil is loaded: where cos α is a correction factor that accounts for the effect of the helix angle and n is the number of elements per each coil. For non-uniform coils average values are used. Finally, nodal forces on nodes I, and J are obtained by equilibrium:

Definition of contact
The algorithm implemented is based on penalty function. In fact, especially when several coils are involved in contact, which is often the case, penalty based methods are much easier to converge than Lagrangian multipliers. Moreover, local small compenetrations allow to represent distributed contacts. Contact regions are spheres of the nominal wire diameter, placed at the end of each coil element, see figure 1. When deformations occur, the position of the spheres is updated according to nodal displacements and orientations. Thus contact regions are placed along the helix, while their positions are computed from the kinematics of the equivalent beams. If p I and p K are the position vectors of nodes I and K, which define a contact pair, the position vectors of the corresponding contact spheres are given by: where R is the rotation matrix defined by rotation parameters as given in (eq. 5), and q 0 is the relative position of the sphere with respect to its node.If, for example, four elements per coil are being used, q 0 writes: However, spherical contact surfaces may eventually lead to unrealistic events, with the spheres slipping on each other along the tangent to the helix. For this reason the distance vector between spheres centres is projected on a plane perpendicular to tangent unit vector t, which is defined as the average between the tangent (helix angle effect is neglected) at both locations: where R AB is the orientation of the contact pair, taken as: and t 0 is the tangent direction in the undeformed configuration, perpendicular to q 0 : Contact force f C is directed as the distance vector h and its magnitude is a power function of compenetration: where K c is an appropriate stiffness constant, k is an exponent slightly greater than 1 and d * is the equivalent wire diameter. The application point of the force is the midpoint between sphere centres: Knowing contact force and its application point, it is trivial to compute corresponding forces and moments on nodes I and K.

Nonlinear solution
When accounting for the contribution of all elements, a system of equation is obtained, which gives internal forces as function of nodal coordinates and rotation parameters: When external loads and boundary conditions are applied, the generic static problem writes: where Σ represents the set of degrees of freedom where boundary conditions s are applied. To find the solution Newton-Raphson procedure is implemented, where Jacobians of elastic and contact forces are efficiently computed applying derivation rules to the equations described in the previous sections. Stiffness matrix K is obtained from Jacobians, then incremental displacements are computed from residual loads df and assumed displacements ds given by boundary contitions: where subscript f f refers to rows and columns corresponding to free dofs and subscript f c refers to free rows and constrained columns.

Numerical example
A numerical example is given and the results of the proposed discrete model are compared to those obtained by finite elements (FE). The study case is a conical spring, under compression, bending and a lateral force perpendicular to the bending plane. Spring ends are close coiled, so that one inactive coil is considered at each end. The FE model is made of 3D beam elements placed along the helix. Top and bottom coils are rigidly connected to a central node, located on spring axis at the beginning and the end of active coils. Loads and boundary conditions are applied to these central nodes.
The deformed configuration is shown in figure 3, where the 3D rendering of the discrete model solution with four elements per coil (yellow coils with blue and red contact spheres), is compared to FEM results (light blue coils).  The absolute displacements of the top centre node, where loads are applied, are plotted in figure ??. FEM and discrete models agree very well in the prediction of axial displacement. Some spread is observed when monitoring lateral displacements, which show a stronger dependence on the number of elements per coil. Solution given by the discrete model clearly changes when switching from three to four elements per coils and it seems not worthy to use more than six. A similar analysis on FE could not be run because of convergence problems. The best compromise sees four contact elements per coil. It would have been possible to use other contact types, but each one has its drawbacks, either in terms of accuracy or computational time.  Finite elements --3D Discrete with 3 elements per coil · · · · · · 3D Discrete with 4 elements per coil · · · · · · 3D Discrete with 6 elements per coil ♦· · · · · · 3D Discrete with 8 elements per coil • · · · · · · 3D Discrete with 12 elements per coil * · · · · · · IOP Publishing doi:10.1088/1757-899X/1214/1/012025 10

Conclusions
A 3D discrete model for non-linear springs is described and successfully applied to a case study. The model is based on Haringx equivalent beam theory, but take advantage of discretization in order to efficiently detect and simulate contact between spring coils. The peculiarity of the preposed modelling technique is the search and definition of contact on a projected geometry following deformation of the spring helix and the application of these effects on the equivalent beam nodes. The model is compared with FEM results; despite the challenges in contact managing it is shown that both modelling strategies, based on very different approaches, gave similar results, with an acceptable spread for practical purposes. The proposed 3D discrete model present several advantages in terms of computational time and model update procedures. The proposed modelling technique it is then considered a valuable basis for further developments aiming at the simulation of dynamic events.