Numerical High-Order Model for the Nonlinear Elastic Computation of Helical Structures

In this work, we propose a high-order algorithm based on the asymptotic numerical method (ANM) for the nonlinear elastic computation of helical structures without neglecting any nonlinear term. The nonlinearity considered in the following study will be a geometric type, and the kinematics adopted in this numerical modeling takes into account the hypotheses of Timoshenko and de Saint-Venant. The finite element used in the discretization of the middle line of this structure is curvilinear with twelve degrees of freedom. Using a simple example, we show the efficiency of the algorithm which was carried out in this context and which resides in the reduction of the number of inversions of the tangent matrix compared to the incremental iterative algorithm of Newton-Raphson.


Introduction
Helical structures, particularly springs, are considered to be the most used mechanical elements in several industries such as aeronautics, automotive, and civil engineering. Due to their reversible elastic deformability, springs can store or release energy when they are subjected to different loads such as normal force, torque, or both [1,2].
Numerical modeling of the linear elastic behavior of helical springs has been studied by many searchers. Taktak et al. [3] proposed a mixed hybrid formulation, including shear effects, based on the use of a curvilinear finite element with two nodes and six degrees of freedom per node. Dammak et al. [4] have shown that it is possible to obtain the distribution of the various generalized stresses along the middle line of the helical spring with great precision by using a single curvilinear finite element with twelve degrees of freedom. Taktak et al. [5,6] also used this same finite element for the calculation of the natural frequencies and the dynamic response of a simple or assembled helical spring. Other researchers [7][8][9][10][11] have used theoretical and numerical methods for the dynamic behavior study as well as the buck-ling of coil springs. However, most of these studies are limited to the linear elastic case.
The objective of this work is to propose a high-order algorithm based on the asymptotic numerical method (ANM) [12][13][14][15][16][17] to calculate the displacements and rotations at a point of the middle line of a helical structure. The originality of our numerical approach is to calculate these unknowns without neglecting any nonlinear term of the strain. The finite element used in our numerical approach is of curvilinear type, defined along the middle line of the structure considered [18,19]. The kinematics adopted in our theoretical formulation takes into account the hypotheses of Timoshenko and de Saint-Venant [3][4][5][6]. The efficiency of the proposed algorithm lies in the reduction of the number of inversions of the tangent matrix compared to the incremental iterative algorithm of Newton-Raphson [20,21].

Problem Modeling
We consider a helical structure of a circular section having a middle radius R, a constant pitch h, and a length L as shown in Figure 1, where L = ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi R 2 + ðh/2πÞ 2 q φ max and φ max is an arbitrary value. The position vector of each point of the middle line is written, in the Cartesian coordinate system base, as follows: where φ is the polar angle.
The direct orthonormal basis f τ ! , n ! , b ! g of the Frenet coordinate system R F , defined at point M, is expressed in the Cartesian basis as follows: where ω 2 = R 2 + ðh/2πÞ 2 . The pass matrix from the Cartesian coordinate system to the Frenet coordinate system is given by The derivatives of these vectors with respect to the curvilinear abscissa are given by where ρ 1 and ρ 2 designate, respectively, the radius of curvature and the radius of torsion, defined by dφ. Using the Timoshenko and de Saint-Venant hypotheses [3][4][5][6], the components of the displacement vector at a point P in the cross section are written in the Frenet coordinate system as follows: where X 2 and X 3 are the Cartesian coordinates in the cross section. uðsÞ, vðsÞ, wðsÞ, θ 1 ðsÞ, θ 2 ðsÞ, and θ 3 ðsÞ are the displacements and rotations at a point in the middle line as shown in Figure 2. The gradient of the displacements at each point of the cross section is given by  The Green-Lagrange strain vector is given by the following formula: where , and γ 13 = γ tb + X 2 χ t3 + X 3 χ t4 , with e t1 , e t2 , and e t3 designating the membrane strains; χ t1 , χ t2 , χ t3 , and χ t4 are the strains in the plane ðX 2 , X 3 Þ; χ n and χ b , χ bn1 , χ bn2 , χ bn3 , and χ bn4 are, respectively, the strains in the planes ðs, X 3 Þ andðs, X 2 Þ; and γ tn and γ tb are the shear strains (see Appendix A).
The stationarity of the potential energy of the studied structure results in where W d ðUÞ represents the elastic strain energy, P ext ðδUÞ is the virtual work of the external forces, and λ is the loading parameter.
The expression of the elastic strain energy is given by where S is the second Piola-Kirchhoff tensor.
The elementary volume at any point of the considered structure is written as The variation of the elastic strain energy is given by the following expression: where fSg and fγg are, respectively, the vectors of the generalized stresses and strains given by where N is the normal effort, T n and T b are the shearing forces, M n and M b are the bending moments along the transversal axes, M i=1,2,3,4 are quantities equivalent to the moment, and M i=5, 6,7,8,9,10 are additional quantities (see Appendix B). The gradient vector of generalized displacements is written in the following form: Taking into account the expression (14), the vector of generalized strains can also be written in the following form: where ½H is a constant matrix which depends on the radius of curvature and ½AðθÞ is a matrix which depends on the gradient of the unknown of the considered problem (see Appendix C). The generalized behavior law in nonlinear elasticity is given by where ½D is the elastic constant matrix which also depends on the radius of curvature of the studied structure (see Appendix D). By virtue of the previous relations, problem (9) is written:

Resolution Strategy
The solution of the nonlinear problem (17) is obtained using a high-order algorithm based on the discretization of the study domain into finite elements, the development of the different unknowns in the Taylor series, and the continuation technique.
3.1. Finite Element Discretization. The discrete form of the problem (17) is obtained by applying the finite element method [18,19]. The degrees of freedom of each curvilinear element as shown in Figure 3 are represented by the displacement vector as follows:

Modelling and Simulation in Engineering
The interpolation of the displacement vector in each element is given by where the interpolation functions are given by where L e = L/ðNe − 1Þ is the length of the finite element "e" and Ne is the number of nodes. The elementary displacement vector fUg e and elementary displacement gradient fθg e are given by where ½NN and ½G are, respectively, the matrix of interpolation functions and the matrix of their derivatives with respect to the curvilinear abscissa (see Appendix E). The variations of elementary vectors are written: where fδqg e is the variation of the nodal displacement vector.
Thus, the discretized form of the global problem (17) is written as follows: where ½BðθÞ = ½B l + ½B nl = ð½H + ½AðθÞÞ½G and ½KðqÞ is a matrix which depends on the unknown vector fqg.

Taylor Series Expansion.
In this stage, we write the unknown vectors, as well as the loading parameter λ in the form of Taylor series developments truncated to the order p, that is to say, where fq 0 g, fS 0 g, and λ 0 are known solutions at the initial point and "a" is the parameter of the asymptotic development defined by [15] q By injecting the developments (24) into (23) and (25) and by identifying the coefficients according to the different Figure 3: Discretization of the middle line into curvilinear finite elements. 4 Modelling and Simulation in Engineering powers of the parameter "a," we obtain a succession of linear problems given as follows: Problem at order k = 1: Problem at order 2 ≤ k ≤ p: where the tangent stiffness matrix ½K t of the discretized structure, evaluated at the initial point, is given by fF nl k g is a vector which depends on solutions to the previous orders, given by and ½Ŝ 0 is the matrix of initial stresses (see Appendix F). Thus, the computation of each asymptotic branch of the nonlinear problem (23) requires a single inversion of the tangent matrix. The validity range of the series (24) is defined by the following criterion [17]: where ε is the tolerance parameter.

Continuation Technique.
For the computation of the entire solution of the nonlinear problem (23), we apply the continuation technique [17] which consists in updating the tangent matrix (28) at the new initial point defined by fq 0 g = fqða max Þg, fS 0 g = fSða max Þg, and λ 0 = λða max Þ.

Results and Discussions
In this part, we consider a helical structure whose mechanical and geometric characteristics are mentioned in Table 1.
One of the ends of this structure is embedded, the other end is subjected to external actions represented by the torsor fF = 10 3 N ; M = 50 N · mg as shown in Figure 4.
Using this example, we solve the nonlinear problem (23) by the proposed high-order algorithm and the Newton-Raphson algorithm. The calculated displacements are expressed in the Cartesian coordinate system and the rotations in the Frenet coordinate system. In this application, we first study the influence of the mesh and the truncation order of the Taylor series on the solution quality of the problem (23) expressed here by the decimal logarithm of the norm of the residual vector (see Figures 5 and 6). For that, we have chosen in the first case a truncation order p = 10 and in the second case a mesh with 59 elements and this by using the same tolerance parameter ε = 10 −10 .
As shown in Figure 5, the optimal mesh corresponds to N e = 60 which ensures a good quality of the solution of the considered problem, expressed here by ðlogÞ 10 kR ! esk ≤ −6. Figure 6 shows that if we use a mesh with 59 elements in the proposed algorithm, then the optimal truncation order corresponds well to the value p = 10 which ensures the same quality of the solution represented in Figure 5.
In the following, we are interested in the comparison of computation times of CPU used by the high-order algorithm and the iterative incremental algorithm of Newton-Raphson and this for various meshes as shown in Table 2. In this numerical experiment, we have fixed the number of ANM steps, that is to say, N step = 10.
As shown in this table, the use of the Newton-Raphson algorithm requires for each mesh a high number of inversions of the tangent matrix in comparison with the proposed algorithm.

Modelling and Simulation in Engineering
In this application, we also studied the variation of the degrees of freedom U x , U y , U z , θ 1 , θ 2 , and θ 3 of the loaded node according to the parameter of loading λ by using the following data: Ne = 60, p = 10, N step = 10, and ε = 10 −10 for the proposed algorithm and ε = 10 −4 for the Newton-Raphson algorithm.
From Figures 7 and 8, we note the good concordance between the numerical solutions obtained by our highorder modeling and those calculated by the iterative Newton-Raphson algorithm.
In Figure 9, we represent the initial state and the middle line deformation of the considered structure, and in Figure 10, we give the evolution of the solution quality of the solved problem according to the loading parameter by using both algorithms.
From Figure 10, we note that our numerical modeling based on the ANM allows us to calculate the displacements and rotations of the middle line of the considered structure with a better quality compared to the results obtained by the Newton-Raphson algorithm.

Conclusion
In this work, we propose a high-order numerical modeling for the nonlinear elastic computation of a thin structure of helical type without neglecting any nonlinear term for the first time.  Modelling and Simulation in Engineering

Modelling and Simulation in Engineering
The algorithm realized in this context is based on the asymptotic numerical method steps. The finite element used in the discretization of the middle line of the structure considered is of curvilinear type with twelve degrees of freedom. Using a simple example, we have shown the efficiency of our numerical approach in comparison with the Newton-Raphson incremental iterative algorithm. First, we discussed the optimum choice of the mesh and the truncation order of the series allowing having the solution of the chosen example with a better quality. In the second stage, we presented a numerical analysis showing that the proposed high-order algorithm significantly reduces the number of inversions of the tangent matrix, and thus the CPU computation time, to have the entire solution of the chosen problem and this in comparison with the reference algorithm mentioned previously.  with