EFFICIENT HIGH-ORDER IMPLICIT SOLVERS FOR THE DYNAMIC OF THIN-WALLED BEAMS WITH OPEN CROSS SECTION UNDER EXTERNAL ARBITRARY LOADINGS

. This paper aims to investigate, in large displacement and torsion context, the nonlinear dynamic behavior of thin-walled beams with open cross section subjected to various loadings by high-order implicit solvers. These homotopy transformations consist to modify the nonlinear discretized dynamic problem by introducing an arbitrary invertible pre-conditioner [ K (cid:63) ] and an arbitrary path following parameter. The nonlinear strongly coupled equations of these structures are derived by using a 3 D nonlinear dynamic model which accounts for large displacements and large torsion without any assumption on torsion angle amplitude. Coupling complex structural phenomena such that warping, bending-bending, and ﬂexural-torsion are taken into account. Two examples of great practical interest of nonlinear dynamic problems of various thin-walled beams with open section are presented to validate the eﬃciency and accuracy of high-order implicit solvers. The obtained results show that the proposed homotopy transformations reveal a few number of matrix triangulations. A comparison with Abaqus code is presented.

1. Introduction. The demand for engineering thin-walled structures with open cross section is continuously increasing on account of their ability to increase stability, functional requirements and reducing weight and cost in many fields of civil construction to mechanical and aerospace structures. During these last decades, many mechanical aspects have been taken into consideration in the design of these types of structures in order to improve their performance and extend operating life. One of important aspect of the design process is the dynamic response of these structures subjected to external dynamic loadings. Therefore, an accurate prediction of dynamic behavior corresponding to a given strength applied to these types of structures with open cross section is of fundamental importance in the design A realistic prediction must be pass obligatorily by either linear or nonlinear analyses according to the displacement and torsion sizes. When the displacement and torsion components of the structure are small, a linear modal analysis is generally used to accede to the structure dynamical characteristics and responses. But, as the displacement components and torsion became large a nonlinear dynamic analysis is necessary because the nonlinear dynamic effects are important and the structure show very complex structural dynamic behavior on account of coupling complex phenomena such that warping, bending-bending, and flexural-torsion.
The nonlinear dynamic of thin-walled beams with open cross section undergoing large displacements and large rotations is an interesting attractive research topic. The analysis of the nonlinear dynamic response of these types of thin-walled structures is a task frequently encountered in the engineering practice. However, the variability of the cross section makes this analysis more complex from a mathematical point of view.
During the last two decades an extensive research has been carried out on linear and nonlinear dynamic of thin-walled beams with constant open cross section undergoing large displacement and torsion subjected to various loading and different boundary conditions based on the total Lagrangian finite element formulations [7,19] and the co-rotational formulations [6,16]. Many finite element formulations and modelisations have been also proposed to analyze the linear and nonlinear vibration of these thin-walled structures [14,20]. However, for thin-walled beams with open cross section, the research works which deal with nonlinear dynamics analysis has received less attention than thin-walled beams with constant open cross section. The research works are very limited concerning this topic [3,27].
We deal in this paper with the nonlinear dynamic analysis of thin-walled beams with open cross section. The shortening effect, pre-buckling deformation, large torsion and flexural-torsional coupling and effect of load eccentricities are included in this nonlinear study. The objective is to investigate the influence of homotopies. Two other new efficient homotopies are considered in this paper, which are different from that adopted in [8]. The choice of these proposed homotopy transformations leads to strong reduction of stiffness matrix triangulations. This high-order implicit algorithm was applied successfully to solving instationary nonlinear problems [15,28,18] and to structural nonlinear dynamic problems [8,9]. Two numerical examples of forced nonlinear dynamic problems of thin-walled beams with open cross section subjected to external arbitrary dynamic loads are analyzed to access the efficiency and the reliability of the developed high-order implicit algorithm by considering new homotopy transformations. It is proven that the used implicit highorder algorithm is more reliable and less time consuming than the classical iterative methods and than the developed algorithm in [8].
The outline of this paper is as follows. Section 2 presents the considered 3D nonlinear dynamic model of thin-walled beams with open cross section and the derivation of the governing dynamic equations of motion of these thin-walled structures. The section 3 is devoted to the proposed numerical approach based on a high-order implicit algorithm. Its performance and comparisons with the Abaqus code and literature results are performed on two typical examples will be given in section 4. Finally the efficiency and the concluding remarks are drawing in section 5.
2. Nonlinear dynamic model for thin-walled beams with open cross section.
2.1. Statement of the problem and its nonlinear kinematic relations. Let us consider a 3D thin-walled beam element of length L and open cross section of area A(x) which occupies a domain of volume Ω bounded by the boundary ∂Ω as depicted in figure 1a. The beam is made of an homogenous, isotropic and elastic material with Young modulus E, shear modulus µ and mass density ρ. The used rectangular Cartesian co-ordinates system is Gxyz of centre G such that G y and G z are the transversal axes and G x is the initial longitudinal axis. The co-ordinates of shear centre C located in Gyz plane are (y c and z c ), those of a point M on the variable section A(x) are denoted by y, z and ω; with ω is the sectorial co-ordinate which characterizes the warping of the section at point M for a nonuniform torsion (see figure 1b) [29]. In the framework of large displacements, large twist angle and small , v M (x, y, z), w M (x, y, z) of a point M on the section contour are expressed by the following nonlinear relations [8,19,21,12] : represents the axial displacement of G, the components v(x) and w(x) represent the displacements of shear point C in y and z directions, θ x is the torsion angle, ω represents the warping function and the two variables c and s are defined by the following trigonometric functions: The symbol () , in equation (1), denotes the derivation with respect to the coordinate x. Noting that the expressions (2) of trigonometric functions c and s are conserved in this model without any assumption on the torsion angle amplitude in both theoretical and numerical analysis. The equation (1) is strongly nonlinear incorporating the flexion and torsion coupling terms with the trigonometric functions c and s. The Vlasov's linear model [29] can be recovered from equation (1) by approximating the trigonometric functions c and s by 0 and θ x respectively and using linear assumptions. The considered nonlinear Green-Lagrange strain tensor, taking into account the large displacements, membrane effect, bending, nonlinear warping effect, has the following components: where is the membrane deformation component, k y and k z are beam curvatures about y and z axes and R is the distance between the point M and the shear centre C, expressed by: with ψ is the variable associated with membrane component given by:

Governing equation of motion.
To determine the governing nonlinear dynamic equation of motion, we use the principle of virtual work: δ(T +U −W ext ) = 0, applied in its variational form, where T is the work of acceleration forces, U is the strain energy and W ext is the work of the applied external forces. The variation of the work of acceleration forces is given by: where e y = y −y c and e z = z −z c are the eccentricities from the shear centre C. The components of the acceleration vectorü M ,v M andẅ M , obtained from equation (1), are: (7) and (8) in equation (6), the variation of the work of acceleration forces, after the integration on elementary cross section area dA writes as: I y (x) and I z (x) are the second moments of area about y and z axes, I ω (x) is the warping constant and I o (x) is polar inertia moment about shear centre. These last quantities are defined respectively by: The expression (9) can be written in terms of mass matrix [M (θ, α)] and gyroscopic matrix [C(θ, α)] in the following form: where {u} is defined by : The mass and gyroscopic matrices [M (θ, α)] and [C(θ, α)] can be decomposed in linear and nonlinear parts as follows: with The variation δU of the strain energy is given by the following expression [8,12,22,23]: (19) where N is the axial force, M y and M z are the bending moments about the principal axes y and z, B ω is the bimoment acting on the cross section, M sv is the St-Venant torsion moment (see figure 2). The higher order stress resultant M R is the Wagners moment, responsible of nonlinear warping and playing an important role in stability and dynamical nonlinear behavior of the structure often neglected in literature. The where β y , β z and β ω are Wagner's coefficients. The expression (19) can be written in the following compact form: where the stress vector {S} is connected to the elastic behavior matrix [D(x)] of the thin-walled beam with open section cross through the constitutive law: with {S}, [D(x)] and {γ} are defined by: with J(x) is the St-Venant torsion given by: Let us introduce the following vector: Due to large torsion context, the strain vector {γ} is highly nonlinear. It can be written in the form: where {γ l } is the linear part of {γ}, {γ nl ({θ}, {θ})} and {γ nlα ({θ}, {α})} are its nonlinear parts representing the quadratic and flexural-torsion coupling terms respectively and t {α} denotes the vector: which takes into account of large rotation and twist coupling flexion, where: Using (29), the strain energy variation (21) and the constitutive law (22) are given by the matrix systems: ]. All of these matrices are given in [22]. For arbitrary loadings, the external virtual work can be expressed, in cartesian co-ordinates, in terms of kinematic variables and resultant forces. Let us denote by t {f j } =< f xej , f yej , f zej > the concentrated vector force applied at point M (x j , y j , z j ) of the cross section A(x i ) (see figure 3a) and by t {F } =< F xe , F ye , F ze > the vector of distributed forces imposed on a part L f of the beam (see figure 3c).
The variation δW ext of the virtual work of external forces, proportional to the load parameter λ, is expressed in terms of displacement variations of by: where N c in (31) is the number of loaded cross sections of the beam and n is the number of concentrated forces applied on the cross section A(x i ), which can be written as the sum of the work of the distributed forces and the work of the concentrated forces as: Inserting the variation of kinematic relations (7) and taking into account for δc =  −sδθ x and δs = (c + 1)δθ x , one gets: Inserting (33) in (31) and replacing the 3D loadings by their statically equivalent forces on the mean fibre (see figures 3b and 3d), we obtain the following variations of works W d and W C of the distributed and concentrated forces respectively: and The efforts used in (34) and (35) are the forces (F ,f ), the bending and twisting moments (M ,m) and Vlasov's warping forces are denoted by (B,b). They are expressed by: To have a matrix form, let us introduce the following vectors: The virtual work δW ext (31) or (32) can be written in the following matrix form: Let us note that one can apply the external real forces away the middle line using the matrices [M x1 ], [M x2 ] and [M x3 ] which represent the flexion and twist from gravity center G(y G , z G ) and twist center C(y c , z c ) that contribute in the stiffness matrix. Finally, using the above definitions and notations, the stationary condition of the total potential energy gives the following system of the dynamic equations and material behavior: The highly nonlinear coupled system (41)  3. Space discretization and methodology of resolution.
3.1. Finite element discretization. The discrete form of the system (41) is obtained by recourse to the standard finite elements method approach [5,10,30,4]. The thin-walled beam of slenderness L is meshed in many 3D beam elements with two nodes and seven degrees of freedom per node. For the shape functions, we adopt linear functions for axial displacement u, and Hermite cubic functions for the variables v, w and θ x . Denoting by:   . Then the problem (46) reduces to: The nonlinear system (48) is usually solved numerically by incremental iterative methods based on the Newmark direct integration [26] and Newton-Raphson method. These last is very consuming in terms of computational time since they require updating the different matrices at each time step or at each iteration.
The aim of this work is to propose three homotopy transformations coupled with high-order implicit solver for analysing the forced nonlinear vibration of thin-walled beams with open cross section subjected to arbitrary external dynamic load taking into account the influence of load eccentricities.

Solution strategy.
The solution of the nonlinear discrete dynamic problem (48) is obtained by using efficient high-order implicit solvers derived by considering good homotopy transformations different from that used in [8]. These high-order implicit solvers are developed by combining the following five techniques: i) a time discretization by the Newmark scheme, ii) a change of variable, iii) appropriate homotopy transformations, iv) power series expansion technique and v) a continuation method. These techniques are detailed in next subsections.
3.2.1. Time discretization procedure. The time discretization of space discretized nonlinear dynamic problem (48) is performed by using unconditionally stable classical Newmark time scheme, widely used in dynamic problems [26]. Knowing the terms {r n g } of (48), the velocities {ṙ n g } and the acceleration {r n g } at time n∆t; with ∆t is the time step, the corresponding quantities {r n+1 g }, {ṙ n+1 g } and {r n+1 g } at time (n + 1)∆t are given by: where α and β are the Newmark constants such that 2α ≤ β ≤ 1 2 which determine both the accuracy and the stability of scheme. From equation (49), one can express the nodal acceleration {r n+1 g } at time (n + 1)∆t as follows: where the coefficients a 0 , a 1 and a 2 in (50) are defined by: The equation (53) can be written in a form of a linear part [K T ]{∆r g } and a nonlinear part {F nl g } given by: where [K T ] is the stiffness matrix including the stiffness, mass, gyroscopic terms and the Newmark coefficients, {F nl g } represents a nonlinear vector and {SM g } is right hand-side. To solve (54) with a reduced computational effort and to decrease the computational cost in matrix inversions, a robust and efficient numerical strategy is needed.
In the next, three high-order implicit solvers are proposed and compared. they are obtained by introducing some homotopy transformations.

The homotopy transformations.
The homotopy technique is one of the important techniques widely used to approximate, in mathematical physics, the solutions of nonlinear static or dynamic problems. This technique was originally proposed by [11] and used after by [2]. It consists to modify the time discretized problem (54) by introducing an adimensional artificial parameter "a" and a pre-conditioner arbitrary matrix [K * ] in this problem. This pre-conditioner arbitrary matrix is homogeneous to the tangent stiffness matrix. Let us note that there exist many manners to perform this modification according to choice of the homotopy transformation. Hence, an appropriate choice of the homotopy transformation is of a great importance because it leads to strong reduction of computational effort and cost. To solve the time discretized problem (54), we present, in this section, three different kinds of homotopy transformations which will be compared in order to select the better in terms of computational cost. a-First homotopy transformation. In order to overcome the inversion of the tangent stiffness matrix at each time step, we adopt an homotopy technique used in the previous work of the same authors [8]. Its consists to introduce an arbitrary invertible matrix and apply an homotopy transformation which acts both on quadratic term and on the right hand side as follows [8,15,28,9,13] : This homotopy transformation allows to transform continuously the system from the original, when a = 1, to the easier system to solve with a = 0. This algorithm is denoted by Alg 1 .
a.1 Perturbation technique. To solve the nonlinear problem (55), the perturbation method is used and "a" is considered as a perturbation parameter [25,24]. The unknowns ∆r g are sought in the form of an integro-power series with respect to homotopy parameter "a" starting from order 1 and truncated at order p as follows: Inserting the series representation (56) in the homotopy system (55) and after equaling like power of "a", we get at each order, the following set of sequence recurrent linear systems: Order p = 1: Order p ≥ 2: Each linear system ((57)-(58)) is solved numerically. Only the right-hand sides vectors are computed for each of these linear systems.
a.2 Validity range. The validity range of series representation (56) is [0, a max ], where the maximal homotopy parameter a max is estimated by the criterion which is expressed in terms of the truncation order of a defined tolerance parameter ε and right hand sides of linear equations verified by the terms of series (56), given by the following relationship: The parameter a max which depends of time allows us to determine the maximum time t max of validity range of series which is defined by the following inequality a max ≥ 1. This maximum time will be considered later as an initial condition for the next step of the continuation procedure. When this inequality becomes false, we increase the truncation order of development in Taylor series.
b-Second homotopy transformation. This second homotopy transformation acts only on the right hand side as follows: This algorithm is denoted by Alg 2 .
b.1 Perturbation technique. Introducing the polynomial representation (56) in the homotopy problem system (60), and equating like power of "a", a set of following linear problems is obtained at each order: Order p = 1: c-Third homotopy transformation. This last homotopy transformation acts only on quadratic term as follows: This algorithm is denoted by Alg 3 .
c.1 Perturbation technique. In contrary to previous cases, let us search the unknown ∆r g of the homotopy problem (63) in the from of the following series expansion with respect to homotopy parameter "a" truncated at order p: Substituting the series representation (64) in (63) and equating like power of "a", we get, at each order, the following linear problems: Order p = 0: We consider a clamped-free mono-symmetrical beam with U cross section of length L = 9m subjected to concentrated dynamical loading F z (t) applied at the section B (see figure 4a). The temporal evolution of dynamical loading F z (t) is given in the figure 4b. All numerical computations are carried out by the following material and geometrical data which are similar to those used in [8,17]. The material data are E = 210GP a, ν = 0.33, ρ = 7850kg/m 3   The nonlinear dynamical analysis is performed in time range [0, 4s] with a time step ∆t = 10 −3 s. To compare the results of high-order implicit solvers to those of the Abaqus code [1], we take as a truncation order p = 15 and a tolerance parameter = 10 −6 . The domain of the considered beam is discretized in 40 3D beam elements [12] for high-order implicit solvers and in 4800 C3D20R solid elements for the Abaqus code.
The time evolutions of vector components (u(L, t), v(L, t), w(L, t), θ x (L, t)) at node 41 (point B) compared to those given by the Abaqus code are depicted in figure 6a-d. These results show a good agreement between the results obtained by high-order implicit solvers and those of Abaqus code. Comparing the appearances of obtained curves with those given in the works of [8,17] concerning thin-walled beams with open constant cross section, we observe that there is a change in the vibrational behavior when the cross section is variable. The full response curves of figure 6 is obtained by high-order implicit solvers using three inversions of the matrix [K ]. Whereas, the obtaining of same response curves by the Abaqus code requires 8965 iterations. Subsequently, we will study the influence of some parameters on the effectiveness of high-order implicit solvers.  a-Effect of the time step. In this paragraph, we study the influence of time step on the parameters of high-order implicit solvers (optimal truncation order: optimal order, Residue: Log|Res|). In the table 1, we remark that the optimal order increases when the time step increases for a tolerance parameter fixed to = 10 −5 . From these results, we notice that the solver Alg 3 is better than the others because its optimal order is the smallest. In Table 2, we notice that the three solvers diverge for the truncation order p = 7, the solver Alg 3 starts to operate from the truncation order p = 8, the solver Alg 2 from p = 9 and the solver Alg 1 from p = 10. For the truncation order p = 10, for example, the solver Alg 3 requires less matrix inversions and less the time CP U than the others. We also notice that the number of matrix inversions IM , the number of right hand sides RHS and the CP U time increase with the truncation order p. In this example, we analyze the nonlinear dynamic behavior of a cantilever bi-symmetrical beam with steel I section subjected to an eccentric trapezoidal load F z (t) applied at end B. The considered data are the same used in [27]; the Young's modulus E = 210GP a, the shear modulus µ = 8.0769 10 7 kP a, the mass density ρ = 7850kg/m 3 , the length L = 8m, the flange thickness t f = 0.03m, the web thickness t w = 0.012m, the section varies according to the relationship h(x) = (2t f + 0.5(1.5 − x/L)) (see figure 8a) and a constant width b = 0.30m (see figure 7a). The beam clamped-free is subjected at free end to an eccentric trapezoidal load F z (t) plotted in the figure 7b. This beam  Figure 7. Cantilever bi-symmetrical beam with steel I cross section under eccentric loading and its time evolution is discretized in 40 3D elements for high-order implicit solvers and in 1020 3D solid elements C3D20R for the Abaqus code. The nonlinear dynamic analysis is made in the time range [0, 0.2s] with a time step ∆t = 0.001s. For the resolution, we take a truncation order p = 6 and a tolerance = 10 −3 for high-order implicit solvers. In this example, the load F z (t) is applied at point O with eccentricities e y = −0.03m and e z = 0.155m with respect to centroid of the cross section where y c = 0m and z c = 0m (see figure 8b). The time evolution curves of components (u(L, t), v(L, t), w(L, t), θ x (L, t)) are reported in figure 9, compared to result computed by Abaqus code and those obtained by Spountzakis et al. [27]. These response curves show a very good agreement with those obtained by the Abaqus code and by Spountzakis et al.. In table 3, we study the influence of time step on the truncation order and on the residue. We remark that the optimal order increases with the time step. The solver Alg 3 is always better than the others because its optimal order is less than the those of others. In the following  Table 3. Comparison between three solvers Alg 1 , Alg 2 and Alg 3 : Influence of time step on the number of matrix inversions IM , number of right hand sides RHS and the CP U time. From this table, we notice that the solver Alg 3 starts to operate for the truncation order p = 3 contrary to the solvers Alg 2 and Alg 1 begin to operate respectively for truncation orders p = 4 and p = 6. We remark also that the solver Alg 3 is always better than the others.