A conjugate modal force strategy for instability analysis of thin-walled structures: an unconstrained vector positional finite element approach

The importance of knowing critical loads and post-critical behavior of thin-walled structures motivates the development of several scientific and practical studies. Most references are concerned with stability analysis for small displacements (first order approach), or with second order stability analyzes, a less precise geometric nonlinear strategy. However, very flexible structures or ones that present small loss of stiffness after the first critical load need more careful analysis. Here we present a shell numerical formulation capable of carrying out stability analysis of thin-walled structures developing large displacements. This formulation uses generalized unconstrained vectors as nodal parameters instead of rotations. To make possible a complete stability analysis using unconstrained vectors, we present an original strategy that imposes a Conjugate Modal Force at the vicinity of critical points, allowing an accurate choice of post-critical paths. Non-conservative forces are also considered and results are compared with literature benchmarks, demonstrating the accuracy and capacity of the proposed formulation.


INTRODUCTION
The stability of thin-walled structures is studied by several important works of today and the last decades. These studies are justified by the search for light structures to support high loads, providing savings in natural resources and ensuring safety for users. In this sense, the knowledge of the sensitivity of the mechanical behavior of structural elements related to imperfections is a research area in constant improvement. Sources of imperfections may come from load application, geometry and material defects. Classical studies (Bleich, 1952;Timoshenko and Gere, 1961;Murray, 1986;Bazant and Cedolin, 1991) analyzing structures with simple geometry and boundary conditions have shown the considerable effect of the presence of imperfections in the achieved critical loads. Among analytical and semi-analytical approaches to solve robust plate buckling, Heydari et al. (2017) solved the buckling of circular functionally graded plates on elastic foundations, providing proper adhesive functions satisfying boundary conditions.
The use of numerical methods to analyze structural stability is indispensable when more general geometry and boundary conditions are considered. In this sense one may cite, among others, works related to the Finite Strip Methods (FSM) and Generalized Beam Theory (GBT). Regarding FSM, Ádány and Schafer (2006a) originally proposed the buckling mode decomposition for single-branched open cross-section members. In Ádány and Schafer (2006b) applications related to the proposed technique are performed. In Ádány and Schafer (2008) this strategy of modal decomposition was extended to general open cross sections. Finally, in Li and Schafer (2010) this formulation is applied in cold-formed steel members design. A good review of this subject can be seen in Li et al. (2014).
Regarding GBT, in Dinis et al. (2006) one can see the formal equations of GBT for general branched cross section, the separation of instability buckling modes and applications. In Camotim et al. (2006) an interesting overview of the possibilities of future applications of GBT is presented. Gonçalves et al. (2009) present the first order buckling behavior of thin-walled members of any cross section using GBT. In Camotiom et al. (2010) one can see a review of GBT, in which all aspects of this formulation are mentioned in detail. Bebiano et al. (2015) rationalize and automate the GBT strategy creating bases for software configurations. These methods, applied in instability analysis, received a lot of contributions in the last years, as for example, Silvestre et al. (2019), Basaglia et al. (2019) and Manta et al. (2020) for GBT and Mirzaei et al. (2015), Poorveis et al. (2019) and Shahmohammadi et al. (2019) for FSM. Some interesting works using classical Finite Element Method (FEM) and solving stability problems may be cited (Mororo et al., 2015;Pastor and Roure, 2009;Ren et al., 2006). Among them, applications related to the stability of thin-walled structures are present, for example, in Anbarasu and Sukumar (2014), Ghumare and Sayyad (2017) and Soares et al. (2019). These works verify the presence of limit points and bifurcations along the equilibrium path. Almost all computational codes use the Newton-Raphson and/or Arc-Length techniques to find the equilibrium path. The instability points along the equilibrium path are usually identified by some supporting strategy, as the verification of the tangent stiffness matrix determinant or eigenvalues.
In the present study, we are interested in developing a specific strategy to deal with equilibrium path and bifurcation points by the position based FEM. As far as the authors' knowledge goes, the first study to describe a position based FEM is the paper of Bonet et al. (2000) in which pressurized membranes are analyzed. In Coda (2003) this approach is applied to solve simple hyperelastic 2D solids and, after that, the last author developed, with co-workers, several variations of this formulation. Only mentioning the works that are related with the present study, in Coda and Greco (2004) a 2D frame dynamic formulation (still using angles as variables) was presented. When developing its extension to 3D applications the unconstrained vector idea arises as an alternative to Euler-Rodrigues or Quaternions. The first work using the unconstrained vector that presents (with details) the geometric description is Coda and Paccola (2007). Another work that clarifies the understanding of the unconstrained vector kinematics is Coda (2009). Three are the main advantages of using the positional unconstrained vector to model 3D bars and shells: (i) As the formulation is total Lagrangian, inverse problems are easily solved, see for instance Coda (2015). (ii) As unconstrained vectors are Lagrangian variables, there is no incremental updating (necessary when using co-rotational formulations, as an example). Thus, the equilibrium is achieved for the full variable, reducing accumulated errors, see example 11.2 at Coda and Paccola (2010), (iii) The mass matrix of unconstrained vector formulation is constant, therefore simple time integration algorithms can be used in transient analysis, see Coda (2009) or Coda (2015. Other advantages are the natural consideration of cross section shear strain; the presence of a numerical chain rule that facilitates the description of initially curved elements; and the simple teaching. More recently, the unconstrained vector kinematics has been successfully applied in various structural problems (Coda et al., 2020;Siqueira and Coda, 2016;Pascon and Coda, 2013). Other recent works using the positional FEM based on unconstrained vectors instead of finite rotations may also be mentioned (Coda, 2015;Paccola et al., 2016;Coda et al., 2017). The use of unconstrained vectors in the analysis of structures composed of non-smooth surfaces -as occurs in the usual analysis of metal profiles -exposes the problem of the non-uniqueness of the normal versor at corners and edges. This problem has been solved in Soares et al. (2019) enabling the technique to solve more complex structural problems.
In order to make possible the stability analysis using unconstrained vectors, in this study we present an original strategy that imposes a Conjugate Modal Force at the vicinity of structural critical points, allowing an accurate choice of post-critical paths by the arc-length method including non-conservative loads (Crisfield, 1981;Feng et al., 1996). The proposed strategy uses the primary idea of imposing perturbations at the vicinity of instability points (Wagner and Wriggers, 1988;Shi, 1996;Fujji and Okazawa, 1997) to find post-critical paths. In this work, the perturbation is originally expressed as forces conjugated to the instability modes at the critical point. The chosen conjugate force is appropriately applied and removed with the only objective of overcoming bifurcations.
The organization of the work begins with the kinematics description of the adopted finite element. Then the constitutive law is described and the equilibrium equations (in their weak form and considering the presence of nonconservative forces) are presented. The applied Arc-Length method is also described using a background of the positional formulation and the proposed perturbation methodology. The perturbation force calculated from eigenvectors is described at the end of the theoretical sections. Numerical examples explore equilibrium paths that have bifurcations, including pressurized tubes, proving the accuracy and applicability of the proposed technique. Coda and Paccola (2008) first presented the isoparametric positional shell element, adopted here. The nodal parameters are generalized vectors, positions and transverse strain variation. It is worth noting that the transverse strain variation is introduced to make the kinematics free from volumetric locking when using 3D complete constitutive models.

POSITIONAL MAPPINGS -KINEMATICS
Observing Figure 1 it is easy to identify that the positional FEM is a Total Lagrangian formulation, as, by means of the initial mapping 0  (Eq. (1)) and the current mapping 1  (Eq. (2)) one writes the change of configuration function  by the composition presented in Eq. (3), which states that the initial configuration is the reference one. The initial and current mappings are given by: in which 0 h is the initial shell thickness, k  is the shape function of node k , k x and k y are, respectively, nodal positions at initial and current configurations of node k , k n is the normal unit vector of node k at the initial configuration, k v is the generalized vector of node k at current configuration, m a is the thickness strain variation of node m and  As mentioned, the change of configuration function is written as: In order to achieve objective strains (Ogden, 1984) it is important to calculate the gradient of the change of configuration function. This gradient is denoted here by symbol A , and can be written as follows: in which (4) one calculates the Green-Lagrange strain tensor (Ogden, 1984) as: in which I is the second order identity tensor. The unknowns of the problem are nodal current positions, nodal unconstrained vectors and nodal strain variation values, clearly present in 1  (Eq. 2).

CONSTITUTIVE RELATIONS
Due to the thickness strain variation, the use of a complete 3D constitutive law does not imply locking. Thus, the Saint-Venant-Kirchhoff constitutive Law is adopted and the corresponding specific strain energy is written as: is the fourth order elastic constitutive tensor, G is the shear elastic modulus,  is the Poisson ratio and ℑ is the identity tensor of fourth order. The energy conjugate of the Green-Lagrange strain (Eq. (5)) is the second Piola-Kirchhoff stress S , defined from the derivative of Eq. (6) regarding the Green-Lagrange strain as: As usual, see (Ogden 1984), the constitutive elastic tensor can be recovered by the second derivative of the specific strain energy, Eq. (6), regarding Green-Lagrange strain, i.e.,

EQUILIBRIUM EQUATIONS
Before introducing non-conservative loads, the equilibrium equations can be obtained from the principle of stationary mechanical energy. The total mechanical energy is written as: in which U is the internal strain energy and P is the potential of external conservative forces, defined as: Latin American Journal of Solids and Structures, 2021, 18(2), e344 5/24 in which 0  is the solid domain at initial configuration, 0  is the boundary of 0  , 0 p is the distributed applied load vector, 0 b is the vector of volume forces and  is the change of configuration function at surfaces or lines. It is important to stress that non-conservative forces will be introduced only at Eq. (20).
The shape functions used to model surfaces or lines will be called k  , associated to node k of a discrete element on 0  , in contrast with the previous notation k  used in Eq.
(2). Then, the function  can be written as: The principle of stationary mechanical energy states that at equilibrium the first variation of Eq. (9) is null, i.e.: From now on y will be a nodal parameter vector that contains k y , k v and k a of each node k , thus, the equilibrium condition can be set defining in Eq. (13) a variation regarding y as: with being int f and ext f the internal and external forces, conjugate quantities of strain and external potential energies, see Eqs. (10) and (11) and S is Piola-Kirchhoff stress (Eq. (7)).
As y is arbitrary and, to make possible the arc-length implementation, the external force ext f is written in Eq. (16) as a control parameter  , and the system of equilibrium equations is written as: in which int f is defined in Eq. (15). Equation (17) is nonlinear regarding y and can be solved by iterative strategies, as Newton-Raphson or Arc-Length methods. In this study, the Arc-Length is employed in order to make possible the analysis of the complete equilibrium path, passing through critical and bifurcation points. The specific use of the arc-length method with the positional formulation is presented in Section 5.
It is important to show the Hessian matrix H of the mechanical energy  (Eq. (9)), necessary to develop the solution process, as well as, to calculate critical loads. When only conservative forces are present, the Hessian matrix is written as: vector positional finite element approach : in which C is the constitutive stress tensor (Eq. (8)). When non-conservative forces are present, the Hessian matrix becomes: in which the last term is described in the following section.

Non-conservative load -follower pressure
One of the objectives of this work is to study the influence of follower pressure on the instability behavior of thinwalled pipes. Therefore, it is necessary to extend the formulation presented in Coda and Paccola (2008) -briefly described above -to consider non-conservative forces. It is assumed that the pressure is always normal to the shell surface in the current configuration and that its intensity p is a constant value. Non-conservative forces can be introduced into equilibrium (Eq. (17)) via the energy variation, Eq. (14), as: in which ( ) p y is the following pressure that depends on the current position y . S is the area of the current configuration, over which the pressure is applied, and  is defined in Eq. (12). For one element  , Eq. (20) becomes: in which ( ) n y is the normal unit vector of the shell element mid surface.
As k y is arbitrary, Eq. (21) can be added to the external force ext f , see Eq. (17), and is written for a node k of a finite element  as: The unit normal vector is calculated as: y tangent (to the shell surface) vectors and the symbol × vector product.
The integral of Eq. (22) is performed here from the parametric dimensionless space . Denoting by d the infinitesimal area of the parametric space, knowing that and using Eq. (23), the nonconservative force becomes: which simplifies the numeric operation. Applying index notation for node k and direction i of element  , Eq. (24) is rewritten as: in which ipq  is the Levi-Civita cyclic permutation tensor of order 3 with , 1, 2, 3 p q  . In Eq. (25) indexes  and  indicate summation on element nodes.
The non-conservative contribution to be considered in the Hessian matrix is achieved from the derivative of Eq. (24) regarding positions z y of node z of a finite element  , that is: Considering Eq. (25) and after some algebraic manipulations over Eq. (26) results: Equation (27) contains both symmetric and antisymmetric parts to be added to matrix H . As mentioned by Schwlizerhof and Ramm (1984) and Iwata et al. (1991) the antisymmetric part of H has negligible effect in the solution of the linear system of equations present in the non-linear solution strategy. It is worth clarifying that the Hessian matrix only regulates the convergence rate of the iterative solver. Moreover, Schwlizerhof and Ramm (1984) and Iwata et al. (1991) also conclude that the antisymmetric part of Eq. (27) has a negligible effect on eigenvalues/eigenvectors determination. From this reasoning, in the present work, only the symmetric part of Eq. (27) is added to the total Hessian matrix and the conclusions of Schwlizerhof and Ramm (1984) and Iwata et al. (1991) are verified at the example section.

SOLUTION PROCESS -ARC-LENGTH
In specialized literature (Wempner, 1971;Riks, 1972;Crisfield, 1981), the arc-length strategy has been used and improved to solve equilibrium paths that present snap-troughs and snap-backs. In this section we show how to use it together with the positional FEM.
Considering the possibility of having non-conservative external forces and considering that vector y is a position trial, the residual vector (that is null at equilibrium, see Eq. (17)) is given as a function of  and y as: In the arc-length method, the following restriction equation -that represents a circumference arc ( Figure 2) -is provided to limit the range of position correction and load intensity, in which s  is the arc radius and  is a factor that adjusts dimensions as   is dimensionless.

Figure 2
Circumference arc representing the additional restriction.
A tip to choose a value for  in Eq. (29) is to use the ratio between the first position and force increments (Eq. (30)), because at this stage structural problems are well behaved, thus one writes: The iterative procedure consists of two stages: the prediction stage and the correction stage. The prediction stage is performed in the first iteration of each incremental step and determines the arc radius to be used in the iterative process. The correction stage starts from the known arc and performs the search for the solution of Eq. (28) using the Newton-Raphson method including the arc effect (Eq. (29)). In what follows these two stages are given in detail.

Prediction stage
Let be 1 k  y and 1 k   a pair that corresponds to equilibrated position at an increment step 1 k  . In order to find a first prediction solution for step k , one increments this solution by k y and k   found by a linearization of the residuum, Eq. (28), written as follows: Latin American Journal of Solids and Structures, 2021, 18(2), e344 9/24 however k   is unknown and needs the radius of the arc to be solved.
At the beginning of the first step ( 0 k  ), 0   is prescribed by the user and 0 y is easily calculated by Eq. (34). Thus, the initial radius is simply given by: For subsequent steps ( 0 k  ) the radius can be updated using the number of iterations of the previous step as a controlling factor when compared to an acceptable one d n (prescribed by the user), resulting: in which 1 k n  is the number of iterations performed in step 1 k  and  is a parameter that adjusts the desired rate.
Knowing the radius k s  , Eqs. (35) and (36) (38)) is the position in step 1 k  , given by: The signal choice in Eq. (37) should ensure an advance in the equilibrium path. A safe and automatic way to find this signal (Eq. (39)) is using the recent history of the equilibrium path evolutions, as proposed by Feng et al. (1996):

Correction stage
At the prediction stage one finds the first trial for position and load factor. Using these values a new residuum ( , ) k k  g y appears and its nullity must be imposed using y g y y y y 0 y (40) in which, recalling Eqs. (32) and (33), the following terms are detached: For this correction stage, one keeps the radius k s  fixed. Therefore, the additional constraint equation imposes that the correction to be determined remains on the circumference arc. This is done by rewriting Eq. (29) as: After some algebraic manipulations, Eq. (43) can be written as a quadratic equation regarding k y and k  , as: Equation (44) can be solved neglecting high order terms of k y and k  , i.e., From Eqs. (40), (41) and (42) one writes k y as a function of k  as: Substituting Eq. (46) into Eq. (45) and isolating k  , one finds: With the value of k  , one achieves k y by Eq. (46). With these vales k   and k y are updated, as well as k  and k y . The correction procedure continues until k  or k y -defined by Eqs. (46) and (47) -become sufficiently small regarding an adopted tolerance. In this work the chosen tolerance test is done as follows: in which the value adopted for the tolerance is usually 8 1.0 10 tol    .

SOLVING BIFURCATION POINTS
At bifurcation points of the equilibrium path, or even at simple critical points, the Hessian matrix, which largely controls the solution search, becomes singular. Thus, in order to obtain paths beyond a bifurcation points, the conjugate perturbation force strategy is applied here in the context of the positional FEM. This perturbation is imposed as external forces per f energetically conjugate to the current instability modes. In this section, the strategy to apply the chosen perturbation force in the Arc-Length solution process is presented and the achievement of the correct perturbation force is described in the next section.
In order to keep the good behavior of the solution process, this force is gradually introduced and the residuum of a step k (see Eq. (28)) is written as: where p  is the value of  at step p in which the perturbation force has been calculated. As k  grows the influence of the perturbation force increases. After Eq. (57) the criterion to determine step p is established.
The perturbation force is applied at the interval per p k p n    , being per n the number of steps in which per f is allowed to act. For per k p n   the perturbation force is removed and Eq. (28) is recovered. The value of per n can be given by the user or, alternatively, can be calculated by the same strategy defined to achieve the step p . Usually a small per n value is sufficient. After withdrawing the perturbation force the analysis continues on the corresponding postcritical branch.
As the new term ( per f ) does not depend upon y , the derivative of the residuum (Eq. (49)) regarding y does not change, see Eq. (19). However, the derivative regarding the load factor  becomes: In short, for step k respecting the interval

FINDING THE PERTURBATION FORCE
As mentioned in the previous section, the solution of the bifurcation points is accomplished by the temporary introduction of the perturbation force per f , see Eq. (50). In this section, it is presented our strategy to calculate per f as generalized forces conjugate to the instability modes near critical points over the equilibrium path. Thus, one starts describing how these modes are calculated using the positional FEM and then the calculation of perturbation forces are presented.
A bifurcation (or a critical point) occurs when 2 0    for some   y 0. This means that at critical points the Hessian matrix H , Eq. (41), must have at least one null eigenvalue, which implies that its determinant is necessarily zero. It is important to recall that the Hessian matrix H can be divided into three parts, see Eqs. (18) and (19). The first part is E H called elastic part, the second is G H called geometric and the third part L H is related to the non-conservative load.
Moreover, G H and L H are dependent of the load level. We introduce an auxiliary load parameter  into the composition of the Hessian matrix only to calculate the following expression: In this expression, if 1   the load level corresponds to the critical load, which implies det( 0 )  H . If det( 0 )  H the load level is subcritical and the imposition of Eq. (51) returns 1   , as the necessary load to enforce equality should be greater than the current one. In short, for Eq. (51)  is an eigenvalue that assumes 1 for critical loads, i.e., corresponds to admit that the current load level  is the critical one.
corresponds to i-th eigenvalue/eigenvector pair associated to the step k of the arc-length analysis,  is the smaller eigenvalue and internal Hessian tensors are given by Eqs. (52), (53) and (54).
The numerical solution of Eq. (55) is made by the solver ARPACK (Sorensen et al., 2008), that has efficient routines for sparse matrices.
The criterion used to define the quasi-critical level p in this work is calculating the eigenvalues ( ) i  at the end of each arc-length step k . The smaller eigenvalue (0)  is always greater than 1 and, as close it becomes to the unit, as near the critical load factor k  becomes to the critical load. In this sense, the perturbation load per f is calculated here only once when (0) 1.1   and this load level is assumed to be p  (used in Eq. (49)), or equivalently this step is the quasicritical one, i.e., p .
Instead of calculating (0)  for all levels of  , as it is done in this work, the value of the determinant of the total Hessian matrix (usually provided by solvers of linear systems) can be used as an indicative of the necessity of (0)  .
For k p  the described procedure in Section 6 allows the solution method to follow the best branch of the equilibrium path, and it is no longer necessary to solve the eigenvalue problem if no other bifurcation is expected.
The perturbation force applied as described in Section 6, and summarized by Eq. (49), is achieved here from the eigenvector ( ) i u that corresponds to the chosen quasi-critical mode (usually related to the smaller eigenvalue (0)  ). To calculate this force it is necessary to create a fictitious configuration, as: in which i is the chosen mode (usually 0 i  ),  is a scale factor that adjusts the intensity of the applied perturbation force and p y is the structure configuration at step p .
The perturbation force is the necessary force to reach configuration per y (Eq. (56)) from configuration p y , i.e.: Differently from what is described by Eq. (57), in Soares et al. (2019) for a very simple example the authors calculated the perturbation force at the beginning of the analysis, i.e., for small displacements, and kept this force along the entire equilibrium path analysis. However, poor results were achieved for the study of bifurcation points, post critical behavior and equilibrium paths. Here the force f per is calculated and applied in Eqs. (49) and (50), so that this force is used only at the vicinity of the current critical configuration. Therefore, this strategy enables the Arc-Length method to accurately overcome bifurcations and find post-critical branches of the equilibrium path, including the presence of nonconservative loads.

NUMERICAL EXAMPLES AND DISCUSSIONS
In this section, numerical examples show the effectiveness of the proposed formulation to determine equilibrium paths beyond a bifurcation point. In the first example, the consulted reference Garcea (2001) imposes imperfections via concentrated forces. In this work (using the strategy presented in Sections 6 and 7) the same problem is solved with and without previous imposed imperfections, easily finding post-critical paths. The second example is presented to consider the presence of non-conservative loads. The difference of equilibrium paths regarding the presence (or not) of the nonconservative following pressure is analyzed. The third example solves a more complicated thin-walled structure stability analysis. In this example, it is shown the large difference of the achieved critical loads when considering small and large displacements. The adopted tolerance for all examples (Eq. (48)) is 8 1.0 10 tol    . Meshes are generated by the software Gmsh (Geuzaine and Remacle, 2009) and visualizations are made by the software Paraview (Ahrens et al., 2005).

Channel section
This dimensionless example deals with the compression analysis of a channel section profile. It is used to verify the determination of critical loads and to compare the resulting equilibrium paths. Garcea (2001) studied this example numerically and the adopted static scheme and dimensions are presented in Figure 3. Garcea (2001) also defined the following material properties: . Considering that u , v and w are displacements in x , y and z directions, respectively, the displacement boundary In all cases the distributed loads are x q x     while the load factor  varies automatically. In addition to the acting compression force, situations involving the presence of imperfections proposed by Garcea (2001) are analyzed. These imperfections are represented by two concentrated forces acting on points A and B of the central section of the structural element. Point C will be used only to measure the displacement C u . Three imperfections cases are considered and defined in Figure 4. In this figure it is also presented the adopted discretization with 5275 nodes and 1120 cubic triangular finite elements. Firstly, a small displacement stability analysis is performed. The achieved critical loads for the three imperfections are compared in Tables 1, 2 and 3. In order to calculate relative differences, the values achieved by Garcea (2001) using Nastran are taken as reference. Results present very good agreement if considering differences among kinematics. It is observed that, in general, the values obtained here are lower than the reference values, emphasizing the adequacy of the adopted discretization.   Using the arc-length strategy associated with the positional FEM, the equilibrium path is achieved for the three presented imperfections. In all cases the same arc-length parameters are adopted, i.e., For the case without imperfections, the procedure proposed in this work is used to go through the equilibrium path. It is important to mention that reference does not solve this case. For the stable initial branch, before the bifurcation point, a fixed radius is chosen based on the number of desired steps to overcome this branch. The adopted radius for this branch is 0.005 s   .
At each step of the initial branch, a stability analysis, based on the eigenvalue problem, is performed. When the smallest eigenvalue becomes smaller than 1.1 , the analysis is considered to be sufficiently close to the critical point, and an automatically calculated perturbation force is applied. In this step it is attributed p    and the perturbation position per y is achieved by Eq. (56), assuming 0.001   . The automatic perturbation force is applied during 3 steps ( 3 per n  ) overcoming bifurcation, as expected.
In this example, the result of the analysis without imposing the predefined imperfections is used for comparison with the other situations. In Figure 5 the equilibrium paths for flexural imperfection, achieved for points A , B and C (Figure 3) are presented. Note that, in this case, the value of the flexural imperfection force, imposed by Garcea (2001), is high, as the critical load achieved with imperfection is much smaller than without imperfection. It is also observed that the equilibrium paths obtained here are similar to the paths presented by the reference, validating the implemented formulation. In this example the importance of large displacement analysis is evident. As the two kinematic presented in Garcea (2001) -called Mixed or frozen -are based on selected instability modes, the reference formulation is not capable to compute a good post critical behavior. The QUAD4 element is a linear Kirchhoff theory element with several modifications to relax excessive constraints as reduced order integration for shear terms, enforcement of curvature compatibility, and the augmentation of transverse shear flexibility to account for a deficiency in the bending strain energy (MacNeal 1978). Therefore, the QUAD4 results serve as a reference value and the use of our formulation certainly gives results that are more accurate. To corroborate with accuracy, our element models initially curved surfaces, includes the influence of transverse shear and uses full integration. Moreover, as mentioned in the introduction, as unconstrained vectors are Lagrangian variables there is no incremental updating (necessary when using co-rotational formulations, as an example) and the equilibrium is achieved for the full variable, reducing accumulated errors.
Results for local imperfection, case 2, are shown in Figure 6. It is observed that the curves are closer to the obtained path without the presence of imperfections, indicating the choice of small imperfection values. It is noted again that the paths obtained here have good adherence with the reference results, especially when compared to Nastran, whose nonlinear procedure does not involve expansion of modes. Finally, graphics for torsion imperfections, case 3, are presented in Figure 7. The interpretation of results is similar to the other cases.
For all cases with imperfection the peaks of the equilibrium paths obtained here are located slightly below those obtained by the reference. This confirms the prediction observed from the analysis of eigenvalues. This effect is also seen in the post-critical paths, which are always below the paths of Nastran, which is the reference result that comes closest to the formulation proposed here. As mentioned before, this is due to the flexibility conferred to our model by the use of high-order elements and Mindlin-like kinematics.  Table 4 presents the values of the critical loads obtained from the stability analysis along the equilibrium path. The result achieved by Garcea (2001) using the software Nastran is used as a reference to calculate relative differences. In this table one observes that the proposed formulation presents more flexible results than all other analyses made by the reference, indicating the good adopted discretization. Thus, we conclude that the formulation is suitable for path following stability analysis. In Figure

Pipe under external radial pressure
This example is the analysis of a pipe subjected to a radial pressure and results are compared with the ones presented by Basaglia et al. (2019). The pipe has the following geometrical characteristics: Radius The material presents the following properties: . To calculate the stress distribution necessary to perform the linear analysis of stability it has been applied a radial external pressure of 1.0 MPa .
This example points out the importance of considering non-conservative loads when treating problems that have following pressure (or loads). In Iwata et al. (1991) analytical solutions for linear stability considering conservative and non-conservative loads are presented, and are reproduced here as follows: The adopted discretizations used here have 37296 , 40992 and 180096 degrees of freedom for lengths 20 , 50 and 400 cm , respectively, and use curved elements with cubic approximation. Results are compared with the ones presented by Basaglia et al. (2019) considering non-conservative pressure (Table 5) and conservative load ( Table 6). The instability modes for conservative and non-conservative forces are practically identical, thus only the no-conservative modes are depicted in Table 5. The achieved results are in accordance with the ones of Basaglia et al. (2019), presenting a more flexible behavior, due to the better discretization. Good results for non-conservative pressure confirm the affirmation of Schwlizerhof and Ramm (1984) and Iwata et al. (1991) that using the symmetric part of the non-conservative load contribution to the Hessian matrix implies sufficient accuracy. It is worth noting that length 400cm is considered a long cylinder in Table 5, and its critical load approaches the analytical solution.
In order to analyze the path following instability of a very long (infinity) tube, rigid body degrees of freedom are eliminated from the discretization presented in Figure 9. This discretization represents only 1cm of the cylinder with 5880 degrees of freedom. In order to reproduce the infinity tube behavior, the components of unconstrained vectors regarding X direction are restricted (no curvature due to Poisson effects allowed).

Figure 9
Tube discretization used to determine the equilibrium path.
To confirm the validity of the new model, a linear stability analysis is performed. The achieved critical loads, considering non-conservative and conservative loads are, respectively, 0.4615MPa and 0.6152MPa . These values are practically the same as the analytical solutions.
The analyses along the equilibrium paths follow the same procedure defined for the previous example. The initial load factor increment is 0 0.05    , kept constant until reaching the critical point vicinity. The intensity parameter of the perturbation force 0.01    is adopted and its sign defines the two possible post-critical paths after bifurcation. In Figure 10 the equilibrium paths (non-conservative and conservative) are shown. Point A (Figure 9) is chosen to depict results. As point A is on the symmetry of the deformed structure the displacement along direction Z corresponds to the radial displacement. From Figure 10 one observes that, as the problem does not have any imperfection and, before the loss of stability, only small displacements take place, the critical loads coincide with the ones achieved using linear analysis. By the difference of the critical loads and the post-critical behaviors of the analyzed structure, it is clear that it is important to consider the influence of non-conservative loads in the Hessian matrix (or tangent stiffness matrix) for stability analysis of thin-walled structures. It is important to mention that, GBT and the Mixed or Frozen methodologies of the previous example use selected instability modes, which may cause difficulty in solving general problems, like the one presented in the next section.

Torus under bending and torsion
This last example is a torus clamped in one end and subjected to a transverse force acting on the other end, as shown in Figure 11. The torus has the following dimensions:  The presented discretization ( Figure 11) corresponds to 1392 curve cubic elements totalizing 44352 degrees of freedom.
Adopting a small load factor, 1.0   , the buckling loads (small displacements -linear analysis) are achieved and presented in Table 7. The corresponding modes are shown in Figure 12. A path following analysis is carried out in order to obtain the equilibrium path and to have a more accurate determination of the real critical load. The equilibrium path is achieved by the arc-length strategy, and it is depicted in Figure 13 for the lower point of the free end.

Figure 13
Equilibrium path -displacement in z direction of the free end.
From the equilibrium path the achieved critical load factor is 37.7 . This value is 44.5% smaller than the first buckling load factor of Table 7. Therefore, for more general structures a simple linear analysis of buckling loads can lead to not safe results. As an illustration, in Figure 14 the last deformed configuration is presented. It is worth stressing that it is important to perform path-following large deformation analysis in design of structures because, in some cases, inconvenient reduction of stiffness may occur at smaller loads than the ones predicted by simple linear buckling analysis. Moreover, the design of bi-stable structures and the solution of inverse problems can be improved by the use of path-following strategies.

CONCLUSIONS
In this work, a total Lagrangian nonlinear geometric shell formulation with kinematics described by generalized vectors and positions (positional FEM) was proposed and applied in the analysis of equilibrium paths passing through bifurcation points and considering non-conservative loads. To this end, an original strategy based on Conjugate Modal Force at the vicinity of critical points of the structure's equilibrium path was developed. This original way of introducing perturbation forces in unconstrained vector Lagrangian FEM was implemented in association with the arc-length method. In addition, the necessary terms to consider non-conservative force contributions in the Hessian matrix of the structure were deduced for the positional FEM. Examples showed that the proposed formulation is accurate and robust, highlighting the use of high-order curved elements. From a practical point of view, two important aspects in the analysis of instability of thin-walled structural members have been confirmed, namely: (i) Non-conservative forces significantly influence structural critical loads and equilibrium paths. (ii) For flexible structures, linear instability analyses (buckling) are not sufficient to predict critical load levels. Future works involve the development of simple and suitable spatial connections to model thin-walled structural elements with shell finite elements, bypassing detailed connection discretization to model complex structures in geometrical nonlinear regime.