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

Authors: Fatima Boussaoui; Hassane Lahmam; Bouazza Braikat
Journal: Modelling and Simulation in Engineering
 (2021)

---

## Abstract

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.


---

## Body

## 1. 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–11] have used theoretical and numerical methods for the dynamic behavior study as well as the buckling 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–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–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].


## 2. Problem Modeling

We consider a helical structure of a circular section having a middle radiusR, a constant pitch h, and a length L as shown in Figure 1, where L=R2+h/2π2φmax and φmax is an arbitrary value.Figure 1
Geometric description of the helical structure.The position vector of each point of the middle line is written, in the Cartesian coordinate system base, as follows:(1)OM→φ=RcosφRsinφh2πφi→,j→,k→,where φ is the polar angle.The direct orthonormal basisτ→,n→,b→ of the Frenet coordinate system RF, defined at point M, is expressed in the Cartesian basis as follows:(2)τ→=1ω−RsinφRcosφh2πi→,j→,k→n→=cosφsinφ0i→,j→,k→b→=1ωh2πsinφ−h2πcosφRi→,j→,k→,where ω2=R2+h/2π2.The pass matrix from the Cartesian coordinate system to the Frenet coordinate system is given by(3)P=−RωsinφRωcosφh2πωcosφsinφ0h2πωsinφ−h2πωcosφRωThe derivatives of these vectors with respect to the curvilinear abscissa are given by(4)dτ→ds=ρ1n→,db→ds=−ρ2n→,dn→ds=ρ2b→−ρ1τ→,where ρ1 and ρ2 designate, respectively, the radius of curvature and the radius of torsion, defined by(5)ρ1=Rω2,ρ2=h2πω2,with ds=R2+h/2π2dφ.Using the Timoshenko and de Saint-Venant hypotheses [3–6], the components of the displacement vector at a point P in the cross section are written in the Frenet coordinate system as follows:(6)U1s,X2,X3=us−X2θ3s+X3θ2s,U2s,X2,X3=vs−X3θ1s,U3s,X2,X3=ws+X2θ1s,where X2 and X3 are the Cartesian coordinates in the cross section. us, vs, ws, θ1s, θ2s, and θ3s are the displacements and rotations at a point in the middle line as shown in Figure 2.Figure 2
Displacements and rotations of the middle line in the Frenet coordinate system.The gradient of the displacements at each point of the cross section is given by(7)∇¯U→PRF=11−ρ1X2∂U1∂s−ρ1U2+ρ1X3∂U1∂X2−ρ2X2∂U1∂X3∂U1∂X2∂U1∂X311−ρ1X2∂U2∂s+ρ1U1−ρ2U3+ρ2X3∂U2∂X2−ρ2X2∂U2∂X3∂U2∂X2∂U2∂X311−ρ1X2∂U3∂s+ρ2U2+ρ2X3∂U3∂X2−ρ2X2∂U3∂X3∂U3∂X2∂U3∂X3RF.The Green-Lagrange strain vector is given by the following formula:(8)γPRF=11−ρ1X2γ11γ12γ13γ1200γ1300RF,where γ11=e1t+X22et3+X32et2−X2χb+X3χn+X2X3χbn1+X23χbn2+X2X32χbn3+X22X3χbn4, γ12=γtn+X2χt1−X3χt2, and γ13=γtb+X2χt3+X3χt4, with et1, et2, and et3 designating the membrane strains; χt1, χt2, χt3, and χt4 are the strains in the plane X2,X3; χn and χb, χbn1, χbn2, χbn3, and χbn4 are, respectively, the strains in the planes s,X3 ands,X2; and γtn and γtb are the shear strains (see Appendix A).The stationarity of the potential energy of the studied structure results in(9)δWdU−λPextδU=0,where WdU represents the elastic strain energy, PextδU is the virtual work of the external forces, and λ is the loading parameter.The expression of the elastic strain energy is given by(10)WdU=12∫VS̿:γ̿dV,where S̿ is the second Piola-Kirchhoff tensor.The elementary volume at any point of the considered structure is written as(11)dV=1−ρ1X2dX2dX3ds.The variation of the elastic strain energy is given by the following expression:(12)δWdU=∫LδγSds,where S and γ are, respectively, the vectors of the generalized stresses and strains given by(13)γ=et1,γtn,γtb,χt1,χt2,χt3,χt4,χn,χb,et2,et3,χbn1,χbn2,χbn3,χbn4,S=N,Tn,Tb,M1,M2,M3,M4,Mn,Mb,M5,M6,M7,M8,M9,M10,where N is the normal effort, Tn and Tb are the shearing forces, Mn and Mb are the bending moments along the transversal axes, Mi=1,2,3,4 are quantities equivalent to the moment, and Mi=5,6,7,8,9,10 are additional quantities (see Appendix B).The gradient vector of generalized displacements is written in the following form:(14)θt=u,s,v,s,w,s,θ,s1,θ,s2,θ,s3,u,v,w,θ1,θ2,θ3.Taking into account the expression (14), the vector of generalized strains can also be written in the following form:(15)γ=H+12Aθθ,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(16)S=Dγ,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:(17)∫LδθHt+AθtSds−λPext=0,S=H+12Aθθ.


## 3. 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:(18)qe=ui,vi,wi,θi1,θi2,θi3i=1,2.Figure 3
Discretization of the middle line into curvilinear finite elements.The interpolation of the displacement vector in each element is given by(19)ues=N1sui+N2sui+1;θ1es=N1sθi1+N2sθi+11,ves=N1svi+N2svi+1;θ2es=N1sθi2+N2sθi+12,wes=N1swi+N2swi+1;θ3es=N1sθi3+N2sθi+13,where the interpolation functions are given by(20)N1s=si+1−sLe,N2s=s−siLe,where Le=L/Ne−1 is the length of the finite element “e” and Ne is the number of nodes.The elementary displacement vectorUe and elementary displacement gradient θe are given by(21)Ue=NNqe,θe=Gqe,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:(22)δUe=NNδqe,δθe=Gδqe,where δqe is the variation of the nodal displacement vector.Thus, the discretized form of the global problem (17) is written as follows:(23)∑e∫LeδqeBθTSds−λPext=0,S=DBθqe,orKqq−λFext=0,S=DBθqe,where Bθ=Bl+Bnl=H+AθG and Kq is a matrix which depends on the unknown vector q.


### 3.2. 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,(24)q=q0+∑k=1pakqk,λ=λ0+∑k=1pakλk,S=So+∑k=1pakSk,where q0, S0, and λ0 are known solutions at the initial point and “a” is the parameter of the asymptotic development defined by [15](25)q−q0q1+λ−λ0λ1=a.By injecting the developments (24) into (23) and (25) and by identifying the coefficients according to the different powers of the parameter “a,” we obtain a succession of linear problems given as follows:Problem at orderk=1:(26)Ktq1=λ1Fext,S1=DBθ0q1,q1q1+λ12=1.Problem at order2≤k≤p:(27)Ktqk=λkFext+Fknl,Sk=DBθ0qke+Sknl,q1qk+λ1λk=0,where the tangent stiffness matrix Kt of the discretized structure, evaluated at the initial point, is given by(28)Kt=∑e∫LeBθ0tDBθ0+GtS^0Gds.Fknl is a vector which depends on solutions to the previous orders, given by(29)Fknl=−∑e∫Le∑r=1k−1BrnltSk−r+12B0tD∑r=1k−1Aθrθk−rds,and S^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]:(30)amax=εU1Up1/p−1,where ε is the tolerance parameter.


### 3.3. 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 q0=qamax, S0=Samax, and λ0=λamax.


## 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:(18)qe=ui,vi,wi,θi1,θi2,θi3i=1,2.Figure 3
Discretization of the middle line into curvilinear finite elements.The interpolation of the displacement vector in each element is given by(19)ues=N1sui+N2sui+1;θ1es=N1sθi1+N2sθi+11,ves=N1svi+N2svi+1;θ2es=N1sθi2+N2sθi+12,wes=N1swi+N2swi+1;θ3es=N1sθi3+N2sθi+13,where the interpolation functions are given by(20)N1s=si+1−sLe,N2s=s−siLe,where Le=L/Ne−1 is the length of the finite element “e” and Ne is the number of nodes.The elementary displacement vectorUe and elementary displacement gradient θe are given by(21)Ue=NNqe,θe=Gqe,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:(22)δUe=NNδqe,δθe=Gδqe,where δqe is the variation of the nodal displacement vector.Thus, the discretized form of the global problem (17) is written as follows:(23)∑e∫LeδqeBθTSds−λPext=0,S=DBθqe,orKqq−λFext=0,S=DBθqe,where Bθ=Bl+Bnl=H+AθG and Kq is a matrix which depends on the unknown vector q.


## 3.2. 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,(24)q=q0+∑k=1pakqk,λ=λ0+∑k=1pakλk,S=So+∑k=1pakSk,where q0, S0, and λ0 are known solutions at the initial point and “a” is the parameter of the asymptotic development defined by [15](25)q−q0q1+λ−λ0λ1=a.By injecting the developments (24) into (23) and (25) and by identifying the coefficients according to the different powers of the parameter “a,” we obtain a succession of linear problems given as follows:Problem at orderk=1:(26)Ktq1=λ1Fext,S1=DBθ0q1,q1q1+λ12=1.Problem at order2≤k≤p:(27)Ktqk=λkFext+Fknl,Sk=DBθ0qke+Sknl,q1qk+λ1λk=0,where the tangent stiffness matrix Kt of the discretized structure, evaluated at the initial point, is given by(28)Kt=∑e∫LeBθ0tDBθ0+GtS^0Gds.Fknl is a vector which depends on solutions to the previous orders, given by(29)Fknl=−∑e∫Le∑r=1k−1BrnltSk−r+12B0tD∑r=1k−1Aθrθk−rds,and S^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]:(30)amax=εU1Up1/p−1,where ε is the tolerance parameter.


## 3.3. 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 q0=qamax, S0=Samax, and λ0=λamax.


## 4. Results and Discussions

In this part, we consider a helical structure whose mechanical and geometric characteristics are mentioned in Table1. One of the ends of this structure is embedded, the other end is subjected to external actions represented by the torsor F=103N;M=50N·m as shown in Figure 4.Table 1
Mechanical and geometric characteristics of the studied structure.
Young’s modulusPoisson’s ratioMiddle radiusRadius of circular sectionPitch of helixNumber of coilsLength of the structureE=2.11011Paν=0.28R=0.5mro=0.075mh=0.19mNc=2L=6.2947mFigure 4
Recessed-loaded helical structure.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.Figure 5
Influence of the mesh on the residue vector: evolution of the Euclidian norm of residue vector with respect to the loading parameter for different meshes.Figure 6
Influence of the truncation order of the Taylor series on the residue vector: evolution of the Euclidian norm of residue vector with respect to the loading parameter for different truncation orders.As shown in Figure5, the optimal mesh corresponds to Ne=60 which ensures a good quality of the solution of the considered problem, expressed here by log10R→es≤−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 Table2. In this numerical experiment, we have fixed the number of ANM steps, that is to say, Nstep=10.Table 2
Influence of the mesh on the CPU computation time used by high-order and Newton-Raphson algorithms for different meshes and for truncation orderp=10.
MeshMaximum loadANMNRCPUNRCPUMANTruncation orderpToleranceεNumber of inversionsCPUToleranceΕNumber of inversionsCPUNe=4016311010−10103.9810−43264160.0740.21Ne=5015051010−10104.7310−43012153.4832.44Ne=6014251010−10105.6910−42852152.4826.79Ne=7013421010−10106.4810−42686150.3423.20Ne=8012461010−10107.4910−42497149.7719.99Ne=9011761010−10108.3110−42354147.5017.74Ne=10011141010−10109.4010−42230145.3015.52Ne=11010561010−10109.3610−42114144.5815.44Ne=12010071010−101010.0610−42016143.3514.24Ne=130964.31010−101010.8810−41930142.9813.14As 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.In this application, we also studied the variation of the degrees of freedomUx, Uy, Uz, θ1, θ2, and θ3 of the loaded node according to the parameter of loading λ by using the following data: Ne=60, p=10, Nstep=10, and ε=10−10 for the proposed algorithm and ε=10−4 for the Newton-Raphson algorithm.From Figures7 and 8, we note the good concordance between the numerical solutions obtained by our high-order modeling and those calculated by the iterative Newton-Raphson algorithm.Figure 7
Load-displacement curves obtained by both algorithms.Figure 8
Load-rotation curves obtained by both algorithms.In Figure9, 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.Figure 9
Initial and deformed configuration of the middle line obtained using the ANM algorithm.Figure 10
Decimal logarithm of the norm of the residue vector according to the loading parameter obtained by both algorithms.From Figure10, 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.


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



---

*Source: 6655909-2021-03-18.xml*