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

Abstract

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.

Keywords:
Conjugate Modal Force; Instability Analysis; Geometrically Exact FEM; Positional FEM

Graphical Abstract

1 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 non-conservative 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.

2 POSITIONAL MAPPINGS - KINEMATICS

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.

Observing Figure 1 it is easy to identify that the positional FEM is a Total Lagrangian formulation, as, by means of the initial mapping ${\zeta }^0$ (Eq. (1)) and the current mapping ${\zeta }^1$ (Eq. (2)) one writes the change of configuration function $\zeta$ 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:

${\zeta }^0(\xi )={\phi }^k({\xi }_1,{\xi }_2)x^k+\frac{h_0}{2}{\xi }_3{\phi }^k({\xi }_1,{\xi }_2)n^k$ (1)

${\zeta }^1(\xi )={\phi }^k({\xi }_1,{\xi }_2)y^k+\frac{h_0}{2}\left[{\xi }_3+{\phi }^m({\xi }_1,{\xi }_2)a^m{\xi }_3^2\right]{\phi }^k({\xi }_1,{\xi }_2)v^k$ (2)

in which $h_0$ is the initial shell thickness, ${\phi }^k$ is the shape function of node $k$, $x^k$ and $y^k$ are, respectively, nodal positions at initial and current configurations of node $k$, $n^k$ is the normal unit vector of node $k$ at the initial configuration, $v^k$ is the generalized vector of node $k$ at current configuration, $a^m$ is the thickness strain variation of node $m$ and $\xi =({\xi }_1,{\xi }_2,{\xi }_3)$ is the dimensionless coordinate vector. Furthermore, repetition of indexes $k$ and $m$ indicates summation over element’s nodes, which are highlighted in Figure 1.

Figure 1
Mapping of initial and current configurations from a dimensionless space.

As mentioned, the change of configuration function is written as:

$\zeta ={\zeta }^1\circ {({\zeta }^0)}^{-1}$ (3)

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:

$A=A^1{(A^0)}^{-1}$ (4)

in which $A^0=\nabla {\zeta }^0=\partial {\zeta }^0/\partial \xi$ and $A^1=\nabla {\zeta }^1=\partial {\zeta }^1/\partial \xi$.

From Eq. (4) one calculates the Green-Lagrange strain tensor (Ogden, 1984) as:

$E=\frac{1}{2}\left[A^{T}A-I\right]=\frac{1}{2}\left[{(A^0)}^{-T}{(A^1)}^T{A}^1{(A^0)}^{-1}-I\right]$ (5)

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 ${\zeta }^1$ (Eq. 2).

3 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:

$u=\frac{1}{2}E:ℭ:E$ (6)

in which $ℭ=2G\left[\Im +\nu /(1-2\nu )I\otimes I\right]$ is the fourth order elastic constitutive tensor, $G$ is the shear elastic modulus, $\nu$ is the Poisson ratio and $\Im$ 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:

$S=\frac{\partial u}{\partial E}=ℭ:E$ (7)

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.,

$\frac{{\partial }^{2}u}{\partial E\otimes \partial E}=ℭ$ (8)

4 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:

$\Pi =U+P$ (9)

in which $U$ is the internal strain energy and $P$ is the potential of external conservative forces, defined as:

$U={\displaystyle \underset{{\Omega }_0}{\int }u(\zeta )d{\Omega }_0}$ (10)

$P=-{\displaystyle \underset{{\Gamma }_0}{\int }{p}_0\cdot \overline{\zeta }d{\Gamma }_0}-{\displaystyle \underset{{\Omega }_0}{\int }{b}_0\cdot \zeta d{\Omega }_0}$ (11)

in which ${\Omega }_0$ is the solid domain at initial configuration, ${\Gamma }_0$ is the boundary of ${\Omega }_0$, $p_0$ is the distributed applied load vector,$b_0$ is the vector of volume forces and $\overline{\zeta }$ 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 ${\varphi }^k$, associated to node $k$ of a discrete element on ${\Gamma }_0$, in contrast with the previous notation ${\phi }^k$ used in Eq. (1) and Eq. (2). Then, the function $\overline{\zeta }$can be written as:

$\overline{\zeta }={\varphi }^k{y}^k$ (12)

The principle of stationary mechanical energy states that at equilibrium the first variation of Eq. (9) is null, i.e.:

$\delta \Pi =\delta U+\delta P=0$ (13)

From now on $y$ will be a nodal parameter vector that contains $y^k$, $v^k$ and $a^k$ of each node $k$, thus, the equilibrium condition can be set defining in Eq. (13) a variation regarding $\delta y$ as:

$\delta \Pi =\left(\frac{\partial U}{\partial y}+\frac{\partial P}{\partial y}\right)\cdot \delta y=0$ (14)

with

$\frac{\partial U}{\partial y}={\displaystyle \underset{{\Omega }_0}{\int }\frac{\partial u(\zeta )}{\partial y}d{\Omega }_0}={\displaystyle \underset{{\Omega }_0}{\int }S:\frac{\partial E}{\partial y}d{\Omega }_0}=f^{\mathrm{int}}$ (15)

$\frac{\partial P}{\partial y}=-{\displaystyle \underset{{\Gamma }_0}{\int }{p}_0\cdot \frac{\partial \overline{\zeta }}{\partial y}d{\Gamma }_0}-{\displaystyle \underset{{\Omega }_0}{\int }{b}_0\cdot \frac{\partial \zeta }{\partial y}d{\Omega }_0}=-f^{ext}$ (16)

being $f^{\mathrm{int}}$ and $f^{ext}$ 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 $\delta y$ is arbitrary and, to make possible the arc-length implementation, the external force $f^{ext}$ is written in Eq. (16) as a control parameter $\lambda$, and the system of equilibrium equations is written as:

$f^{\mathrm{int}}-\lambda f^{ext}=0$ (17)

in which $f^{int}$ 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 $\Pi$ (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:

$H=\frac{{\partial }^{2}U}{\partial y\otimes \partial y}=\frac{\partial f^{\mathrm{int}}}{\partial y}={\displaystyle \underset{{\Omega }_0}{\int }\frac{\partial E}{\partial y}:ℭ:\frac{\partial E}{\partial y}d{\Omega }_0}+{\displaystyle \underset{{\Omega }_0}{\int }S:\frac{{\partial }^{2}E}{\partial y\otimes \partial y}d{\Omega }_0}$ (18)

in which $ℭ$ is the constitutive stress tensor (Eq. (8)).

When non-conservative forces are present, the Hessian matrix becomes:

$H=\frac{\partial f^{\mathrm{int}}}{\partial y}-\lambda \frac{\partial f^{ext}}{\partial y}$ (19)

in which the last term is described in the following section.

4.1 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 thin-walled 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:

$\delta P=-{\displaystyle \underset{S}{\int }p(y)\cdot \delta \overline{\zeta }d\text{}S}$ (20)

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 $\overline{\zeta }$ is defined in Eq. (12).

For one element $\mathcal{l}$, Eq. (20) becomes:

$\delta P=-{\displaystyle \underset{S^{\mathcal{l}}}{\int }\left(pn(y)\right)\cdot \left({\varphi }^k\delta y^k\right)d\text{}S}$ (21)

in which $n(y)$ is the normal unit vector of the shell element mid surface.

As $\delta y^k$ is arbitrary, Eq. (21) can be added to the external force $f^{ext}$, see Eq. (17), and is written for a node $k$ of a finite element $\mathcal{l}$ as:

${(f_{\mathcal{l}}^{ext})}^k={\displaystyle \underset{S^{\mathcal{l}}}{\int }p{\varphi }^{k}n(y)d\text{}S}$ (22)

The unit normal vector is calculated as:

$n(y)=\frac{u\times v}{\Vert u\times v\Vert }$ (23)

being $u=(\partial {\varphi }^k/\partial {\xi }_1)y^k$ and $v=(\partial {\varphi }^k/\partial {\xi }_2)y^k$ tangent (to the shell surface) vectors and the symbol $\times$ vector product. The integral of Eq. (22) is performed here from the parametric dimensionless space $\xi =({\xi }_1,{\xi }_2)$. Denoting by $d\text{}\Lambda$ the infinitesimal area of the parametric space, knowing that $dS=J_{1}d\Lambda =\Vert u\times v\Vert d\Lambda$ and using Eq. (23), the non-conservative force becomes:

${(f_{\mathcal{l}}^{ext})}^k={\displaystyle \underset{{\Lambda }^{\mathcal{l}}}{\int }p{\varphi }^k(u\times v)d\Lambda }$ (24)

which simplifies the numeric operation. Applying index notation for node $k$ and direction $i$ of element $\mathcal{l}$, Eq. (24) is rewritten as:

${(f_{\mathcal{l}}^{ext})}_i^k={\displaystyle \underset{{\Lambda }^{\mathcal{l}}}{\int }p{\varphi }^k({\epsilon }_{ipq}{u}_p{v}_q)d\Lambda }={\displaystyle \underset{{\Lambda }^{\mathcal{l}}}{\int }p{\varphi }^k\left[{\epsilon }_{ipq}\left(\frac{\partial {\varphi }^{\alpha }}{\partial {\xi }_1}{y}_p^{\alpha }\right)\left(\frac{\partial {\varphi }^{\beta }}{\partial {\xi }_2}{y}_q^{\beta }\right)\right]d\Lambda }$ (25)

in which ${\epsilon }_{ipq}$ is the Levi-Civita cyclic permutation tensor of order 3 with $p,q=\mathrm{1,2,3}$. In Eq. (25) indexes $\alpha$ and $\beta$ 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 $y^z$ of node $z$ of a finite element $\mathcal{l}$, that is:

$\frac{\partial {(f_{\mathcal{l}}^{ext})}^k}{\partial y^z}={\displaystyle \underset{{\Lambda }^{\mathcal{l}}}{\int }p{\varphi }^k\left(\frac{\partial u}{\partial y^z}\times v+u\times \frac{\partial v}{\partial y^z}\right)d\Lambda }$ (26)

Considering Eq. (25) and after some algebraic manipulations over Eq. (26) results:

$\frac{\partial {(f_{\mathcal{l}}^{ext})}_i^k}{\partial y_{\gamma }^z}=p{\epsilon }_{i\gamma p}{y}_p^{\alpha }\left[{\displaystyle \underset{{\Lambda }^{\mathcal{l}}}{\int }{\varphi }^k\left(\frac{\partial {\varphi }^z}{\partial {\xi }_1}\frac{\partial {\varphi }^{\alpha }}{\partial {\xi }_2}-\frac{\partial {\varphi }^{\alpha }}{\partial {\xi }_1}\frac{\partial {\varphi }^z}{\partial {\xi }_2}\right)d\Lambda }\right]$ (27)

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.

5 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 $\lambda$ and $y$ as:

$g(y,\lambda )=f^{\mathrm{int}}(y)-\lambda f^{ext}(y)$ (28)

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,

$\Delta s^2=\Delta y\cdot \Delta y+{\psi }^2\Delta {\lambda }^2$ (29)

in which $\Delta s$ is the arc radius and $\psi$ is a factor that adjusts dimensions as $\Delta \lambda$ is dimensionless.

Figure 2
Circumference arc representing the additional restriction.

A tip to choose a value for $\psi$ 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:

$\psi =\frac{\Vert \Delta y_0\Vert }{\left|\Delta {\lambda }_0\right|}$ (30)

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.

5.1 Prediction stage

Let be $y_{k-1}$ and ${\lambda }_{k-1}$ a pair that corresponds to equilibrated position at an increment step $k-1$. In order to find a first prediction solution for step $k$, one increments this solution by $\Delta y_k$ and $\Delta {\lambda }_k$ found by a linearization of the residuum, Eq. (28), written as follows:

$g(y_{k-1}+\Delta y_k,{\lambda }_{k-1}+\Delta {\lambda }_k)\approx g(y_{k-1},{\lambda }_{k-1})+\frac{\partial g}{\partial y}(y_{k-1},{\lambda }_{k-1})\cdot \Delta y_k+\frac{\partial g}{\partial \lambda }(y_{k-1},{\lambda }_{k-1})\Delta {\lambda }_k=0$ (31)

in which

$H(y_{k-1},{\lambda }_{k-1})=\frac{\partial g}{\partial y}(y_{k-1},{\lambda }_{k-1})=\frac{\partial f^{\mathrm{int}}}{\partial y}(y_{k-1})-{\lambda }_{k-1}\frac{\partial f^{ext}}{\partial y}(y_{k-1})$ (32)

$\frac{\partial g}{\partial \lambda }(y_{k-1},{\lambda }_{k-1})=-f^{ext}(y_{k-1})$ (33)

As $y_{k-1}$ and ${\lambda }_{k-1}$ satisfy the equilibrium, one has $g(y_{k-1},{\lambda }_{k-1})=0$. The position increment $\Delta y_k$can be achieved from Eq. (31) as:

$\Delta y_k=\Delta {\lambda }_k{H}^{-1}(y_{k-1},{\lambda }_{k-1})\cdot f^{ext}(y_{k-1})$ (34)

however $\Delta {\lambda }_k$ is unknown and needs the radius of the arc to be solved.

At the beginning of the first step ($k=0$), $\Delta {\lambda }_0$ is prescribed by the user and $\Delta y_0$ is easily calculated by Eq. (34). Thus, the initial radius is simply given by:

$\Delta s_0=\sqrt{\Delta y_0\cdot \Delta y_0+{\psi }^2\Delta {\lambda }_0^2}$ (35)

For subsequent steps ($k>0$) the radius can be updated using the number of iterations of the previous step as a controlling factor when compared to an acceptable one $n_d$ (prescribed by the user), resulting:

$\Delta s_k=\Delta s_{k-1}{\left(\frac{n_d}{n_{k-1}}\right)}^{\gamma }$ (36)

in which $n_{k-1}$ is the number of iterations performed in step $k-1$ and $\gamma$ is a parameter that adjusts the desired rate.

Knowing the radius $\Delta s_k$, Eqs. (35) and (36), the current increment $\Delta {\lambda }_k$ is achieved substituting Eq. (34) into Eq. (29), i.e.:

$\Delta {\lambda }_k=\pm \frac{\Delta s_k}{\sqrt{y_{k-1}^T\cdot y_{k-1}^T+{\psi }^2}}$ (37)

in which $y_{k-1}^T$ (Eq. (38)) is the position in step $k-1$, given by:

$y_{k-1}^T=H^{-1}(y_{k-1},{\lambda }_{k-1})\cdot f^{ext}(y_{k-1})$ (38)

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):

$sign\left(\Delta {\lambda }_k\right)=sign\left(\Delta y_{k-1}\cdot y_{k-1}^T+\Delta {\lambda }_{k-1}{\psi }^2\right)$ (39)

After achieving $\Delta {\lambda }_k$, Eq. (37), the increment $\Delta y_k$ is achieved by Eq. (34). With $\Delta {\lambda }_k$ and $\Delta y_k$ known one updates the load factor ${\lambda }_k={\lambda }_{k-1}+\Delta {\lambda }_k$ and the position $y_k=y_{k-1}+\Delta y_k$ and goes to the correction stage.

5.2 Correction stage

At the prediction stage one finds the first trial for position and load factor. Using these values a new residuum $g(y_k,{\lambda }_k)$ appears and its nullity must be imposed using

$g(y_k+\delta y_k,{\lambda }_k+\delta {\lambda }_k)\approx g(y_k,{\lambda }_k)+\frac{\partial g}{\partial y}(y_k,{\lambda }_k)\cdot \delta y_k+\frac{\partial g}{\partial \lambda }(y_k,{\lambda }_k)\delta {\lambda }_k=0$ (40)

in which, recalling Eqs. (32) and (33), the following terms are detached:

$H(y_k,{\lambda }_k)=\frac{\partial g}{\partial y}(y_k,{\lambda }_k)=\frac{\partial f^{\mathrm{int}}}{\partial y}(y_k)-{\lambda }_k\frac{\partial f^{ext}}{\partial y}(y_k)$ (41)

$\frac{\partial g}{\partial \lambda }(y_k,{\lambda }_k)=-f^{ext}(y_k)$ (42)

For this correction stage, one keeps the radius $\Delta s_k$ 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:

$\Delta s_k^2=(\Delta y_k+\delta y_k)\cdot (\Delta y_k+\delta y_k)+{\psi }^2{(\Delta {\lambda }_k+\delta {\lambda }_k)}^2$ (43)

After some algebraic manipulations, Eq. (43) can be written as a quadratic equation regarding $\delta y_k$ and $\delta {\lambda }_k$, as:

$2(\Delta y_k\cdot \delta y_k+{\psi }^2\Delta {\lambda }_k\delta {\lambda }_k)+\delta y_k\cdot \delta y_k+{\psi }^2\delta {\lambda }_k^2=0$ (44)

Equation (44) can be solved neglecting high order terms of $\delta y_k$ and $\delta {\lambda }_k$, i.e.,

$\Delta y_k\cdot \delta y_k+{\psi }^2\Delta {\lambda }_k\delta {\lambda }_k=0$ (45)

From Eqs. (40), (41) and (42) one writes $\delta y_k$ as a function of $\delta {\lambda }_k$ as:

$\delta y_k=H^{-1}(y_k,{\lambda }_k)\cdot \left[\delta {\lambda }_k{f}^{ext}(y_k)-g(y_k,{\lambda }_k)\right]$ (46)

Substituting Eq. (46) into Eq. (45) and isolating $\delta {\lambda }_k$, one finds:

$\delta {\lambda }_k=\frac{\Delta y_k\cdot \left[H^{-1}(y_k,{\lambda }_k)\cdot g(y_k,{\lambda }_k)\right]}{\Delta y_k\cdot \left[H^{-1}(y_k,{\lambda }_k)\cdot f^{ext}(y_k)\right]+{\psi }^2\Delta {\lambda }_k}$ (47)

With the value of $\delta {\lambda }_k$, one achieves $\delta y_k$ by Eq. (46). With these vales $\Delta {\lambda }_k$ and $\Delta y_k$ are updated, as well as ${\lambda }_k$and $y_k$. The correction procedure continues until $\delta {\lambda }_k$ or $\delta y_k$ - defined by Eqs. (46) and (47) - become sufficiently small regarding an adopted tolerance. In this work the chosen tolerance test is done as follows:

$\frac{\Vert \delta y_k\Vert }{\Vert x\Vert }<tol$ (48)

in which the value adopted for the tolerance is usually $tol=1.0\cdot {10}^{-8}$.

6 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 $f^{per}$ 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:

$g(y_k,{\lambda }_k)=f^{\mathrm{int}}(y_k)-{\lambda }_k{f}^{ext}(y_k)-({\lambda }_k-{\lambda }_p)f^{per}$ (49)

where ${\lambda }_p$ is the value of $\lambda$ at step $p$ in which the perturbation force has been calculated. As ${\lambda }_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 $p<k\le p+n_{per}$, being $n_{per}$ the number of steps in which $f^{per}$ is allowed to act. For $k>p+n_{per}$ the perturbation force is removed and Eq. (28) is recovered. The value of $n_{per}$ can be given by the user or, alternatively, can be calculated by the same strategy defined to achieve the step $p$. Usually a small $n_{per}$ value is sufficient. After withdrawing the perturbation force the analysis continues on the corresponding post-critical branch.

As the new term ($f^{per}$) 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 $\lambda$ becomes:

$\frac{\partial g}{\partial \lambda }(y,\lambda )=-\left(f^{ext}(y)+f^{per}\right)$ (50)

In short, for step $k$ respecting the interval $p<k\le p+n_{per}$ the arc-length is applied considering Eq. (49) and all other equations presented in Sections 5.1 and 5.2 are used replacing $f^{ext}(y)$ by $f^{ext}(y)+f^{per}$.

7 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 $f^{per}$, see Eq. (50). In this section, it is presented our strategy to calculate $f^{per}$ 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 ${\delta }^2\Pi =0$ for some $\delta y\ne 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 $H^E$ called elastic part, the second is $H^G$ called geometric and the third part $H^L$ is related to the non-conservative load. Moreover, $H^G$ and $H^L$ are dependent of the load level. We introduce an auxiliary load parameter $\mu$ into the composition of the Hessian matrix only to calculate the following expression:

$\det \left[H^E+\mu \left(H^G+H^L\right)\right]=0$ (51)

In this expression, if $\mu =1$ the load level corresponds to the critical load, which implies $\det (H)=0$. If $\det (H)\ne 0$ the load level is subcritical and the imposition of Eq. (51) returns $\mu >1$, as the necessary load to enforce equality should be greater than the current one. In short, for Eq. (51) $\mu$ is an eigenvalue that assumes 1 for critical loads, i.e., corresponds to admit that the current load level $\lambda$ is the critical one. The calculation of $H^E$, $H^G$ and $H^L$ are made by the following formulae, already discussed:

$H^E={\displaystyle \underset{{\Omega }_0}{\int }\frac{\partial E}{\partial y}:ℭ:\frac{\partial E}{\partial y}d{\Omega }_0}$ (52)

$H^G={\displaystyle \underset{{\Omega }_0}{\int }S:\frac{{\partial }^{2}E}{\partial y\otimes \partial y}d{\Omega }_0}$ (53)

$H{\text{}}^L=-\lambda \frac{\partial f^{ext}}{\partial y}$ (54)

Equation (51) can be better represented as a generalized eigenvector problem, i.e.:

$\left[H^E+{\mu }_{(i)}\left(H^G+H^L\right)\right]u_{(i)}=0$ (55)

in which $({\mu }_{(i)},u_{(i)})$ corresponds to i-th eigenvalue/eigenvector pair associated to the step $k$ of the arc-length analysis, ${\mu }_{(0)}$ 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 ${\mu }_{(i)}$ at the end of each arc-length step $k$. The smaller eigenvalue ${\mu }_{(0)}$ is always greater than 1 and, as close it becomes to the unit, as near the critical load factor ${\lambda }_k$ becomes to the critical load. In this sense, the perturbation load $f^{per}$ is calculated here only once when ${\mu }_{(0)}<1.1$ and this load level is assumed to be ${\lambda }_p$ (used in Eq. (49)), or equivalently this step is the quasi-critical one, i.e., $p$.

Instead of calculating ${\mu }_{(0)}$ for all levels of $\lambda$, 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 ${\mu }_{(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 $u_{(i)}$ that corresponds to the chosen quasi-critical mode (usually related to the smaller eigenvalue ${\mu }_{(0)}$). To calculate this force it is necessary to create a fictitious configuration, as:

$y^{per}=y_p+\alpha \frac{u_{(i)}}{\Vert u_{(i)}\Vert }$ (56)

in which $i$ is the chosen mode (usually $i=0$), $\alpha$ is a scale factor that adjusts the intensity of the applied perturbation force and $y_p$ is the structure configuration at step $p$.

The perturbation force is the necessary force to reach configuration $y^{per}$ (Eq. (56)) from configuration $y_p$, i.e.:

$f^{per}={\displaystyle \underset{{\Omega }_0}{\int }{\left(S:\frac{\partial E}{\partial y}\right)|}_{y^{per}}d{\Omega }_0}-f^{\mathrm{int}}(y_p)$ (57)

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 non-conservative loads.

8 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 non-conservative 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 $tol=1.0\cdot {10}^{-8}$.

Meshes are generated by the software Gmsh (Geuzaine and Remacle, 2009) and visualizations are made by the software Paraview (Ahrens et al., 2005).

8.1 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: $E=2.1\cdot {10}^6$ and $\nu =0.3$.

Considering that $u$, $v$ and $w$ are displacements in $x$, $y$ and $z$ directions, respectively, the displacement boundary conditions are: $w(\mathrm{0,}y\mathrm{,0})=w(\mathcal{l},y\mathrm{,0})=v(\mathrm{0,0,0})=v(\mathcal{l}\mathrm{,0,0})=u(\mathcal{l}/\mathrm{2,0,0})=0$. In all cases the distributed loads are $q_x(x=0)=1.0$and $q_x(x=\mathcal{l})=-1.0$ while the load factor $\lambda$ varies automatically.

Figure 3
Static scheme and cross section details.

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$u_C$. 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.

Figure 4
This work discretization and imposed imperfections by Garcea (2001).

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.

Table 1
Eigenvalues for flexural imperfection.

Mode	Garcea (2001 )			Positional FEM	Dif (%)	Mode type
	Mixed	Frozen	Nastran
1	1291.6	1292.2	1298.2	1291.5	0.52	Flexural
2	1396.0	1192.5	1195.2	1188.9	0.53	Torsional
3	1994.4	1992.8	1919.4	1909.0	0.54	Local (3 half-waves)
4	2046.4	2044.7	2094.5	1955.1	6.66	Local (4 half-waves)

Table 2
Eigenvalues for local imperfection.

Mode	Garcea (2001 )			Positional FEM	Dif (%)	Mode type
	Mixed	Frozen	Nastran
1	1293.2	1289.9	1296.1	1289.4	0.52	Flexural
2	1409.5	1387.7	1391.2	1383.9	0.52	Torsional
3	3150.7	3134.3	3084.8	3083.1	0.06	Local (13 half-waves)
4	3150.9	3148.6	3097.4	3088.6	0.28	Local (14 half-waves)

Table 3
Eigenvalues for torsion imperfection.

Mode	Garcea (2001 )			Positional FEM	Dif (%)	Mode type
	Mixed	Frozen	Nastran
1	1264.9	1289.9	1296.0	1289.4	0.51	Flexural
2	1833.9	1711.8	1715.6	1706.3	0.54	Torsional
3	2808.5	2912.1	2842.4	2837.4	0.18	Local (7 half-waves)
4	2823.9	1930.0	2863.5	2855.5	0.28	Local (8 half-waves)

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., $n_d=4$, $\gamma =1.0$ and $\Delta {\lambda }_0=1.0$.

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 $\Delta s=0.005$.

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 ${\lambda }_p=\lambda$ and the perturbation position $y^{per}$ is achieved by Eq. (56), assuming $\alpha =0.001$. The automatic perturbation force is applied during 3 steps ($n_{per}=3$) 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.

Figure 5
Equilibrium paths for points A, B and C - flexural imperfections and without imperfection.

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.

Figure 6
Equilibrium paths for points A, B and C - local imperfections and without imperfection.

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.

Figure 7
Equilibrium paths for points A, B and C - torsion imperfections and without imperfection.

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.

Table 4
Critical loads achieved along the equilibrium path.

Imperfection	Garcea (2001 )			Positional FEM	Dif (%)
	Mixed	Frozen	Nastran
None	-	-	-	1289.2	-
Flexural	900.3	908.4	896.3	883.7	1.41
Local	1246.1	1284.8	1235.4	1223.3	0.98
Torsional	1128.7	1204.2	1117.0	1105.9	0.99

In Figure 8 the post-critical configurations for each imposed imperfection are presented. In all cases the current configurations present a composition of modes. Case (a) and case (b) present the composition of global flexural and local flanges modes. In case (c), there is a composition between global torsion and local flange modes.

Figure 8
Post-critical configurations for (a) flexural, (b) local and (c) torsional imperfections.

8.2 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 $r=6cm$ and thickness $t=0.12cm$. Three lengths are used in the analysis: $L_1=20cm$, $L_2=50cm$ (short pipes) and $L_3=400cm$ (long pipe). In all situations the transversal displacements at the supports are restricted (rigid diaphragm). The material presents the following properties: $E=210GPa$ and $\nu =0.3$. To calculate the stress distribution necessary to perform the linear analysis of stability it has been applied a radial external pressure of $1.0MPa$.

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:

$p_{cr\mathrm{,1}}=\frac{E}{1-{\nu }^2}\frac{t^3}{4r^3}\left(non-conservativepressure\right)$ (58)

$p_{cr\mathrm{,2}}=\frac{E}{1-{\nu }^2}\frac{t^3}{3r^3}\left(conservativeradialload\right)$ (59)

The adopted discretizations used here have $37296$, $40992$ and $180096$ degrees of freedom for lengths $20$, $50$ and $400cm$, 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.

Table 5
Critical loads ($MPa$) and corresponding instability modes.

$L(cm)$	Non-conservative pressure
	$p_{cr\mathrm{,1}}$	Basaglia et al. (2019)		Positional FEM
		GBT	ANSYS
20	-	3.4149	3.4232	3.3493
50	-	1.4097	1.4184	1.3995
400	0.4615	0.4624	0.4659	0.4626

Table 6
Critical loads ($MPa$) conservative pressure.

$L(cm)$	Conservative pressure
	$p_{cr\mathrm{,2}}$	Basaglia et al. (2019)		Positional FEM
		GBT	ANSYS
20	-	3.5273	3.5634	3.5570
50	-	1.5742	1.5712	1.5703
400	0.6154	0.6165	0.6176	0.6165

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 $\Delta {\lambda }_0=0.05$, kept constant until reaching the critical point vicinity. The intensity parameter of the perturbation force $\alpha =\pm 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.

Figure 10
Equilibrium paths of point A for the pressurized tube.

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.

8.3 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: $R=100cm$, $r=10cm$ and thickness $t=0.5cm$. The applied load $q=0.1\lambda kN/cm$ is distributed along the perimeter of the free end. The adopted material properties are: $E=210GPa$ and $\nu =0.3$.

Figure 11
Static scheme - torus subjected to a transverse load.

The presented discretization (Figure 11) corresponds to $1392$ curve cubic elements totalizing $44352$ degrees of freedom.

Adopting a small load factor, $\lambda =1.0$, the buckling loads (small displacements - linear analysis) are achieved and presented in Table 7. The corresponding modes are shown in Figure 12.

Table 7
Buckling loads for torus flexural / torsion.

Mode	Buckling Load ($\lambda$)
1	67.9609
2	68.4315
3	75.8721
4	76.1218

Figure 12
Buckling modes for torus flexural / torsion.

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.

Figure 14
Deformed configuration for the last analysis step ($\lambda =31.48$).

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.

9 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.

Acknowledgements

This study has three different funding: (i) Computational Resources 2020/05393-4 by São Paulo Research Foundation (FAPESP). (ii) PhD Grant 2018/19288-8 by São Paulo Research Foundation (FAPESP). (iii) Physical Space by Coordenação de Aperfeiçoamento de Pessoal de Nível Superior - Brasil (CAPES) - Finance Code 001.

References

• Ádány, S., Schafer, B.W. (2006a). Buckling mode decomposition of single-branched open cross-section members via finite strip method: Derivation. Thin-Walled Structures, 44:563-584.
• Ádány, S., Schafer, B.W. (2006b). Buckling mode decomposition of single-branched open cross-section members via finite strip method: Application and examples. Thin-Walled Structures, 44:585-600.
• Ádány, S., Schafer, B.W. (2008). A full modal decomposition of thin-walled, single-branched open cross-section members via the constrained finite strip method. Journal of Constructional Steel Research, 64:12-29.
• Ahrens, J., Geveci, B., Law, C. (2005). ParaView: An End-User Tool for Large Data Visualization, Visualization Handbook, Elsevier.
• Anbarasu, M., Sukumar, S. (2014). Local/Distortional/Global Buckling Mode Interaction on Thin Walled Lipped Channel Columns. Latin American Journal of Solids and Structures, 11(8):1363-1375.
• Basaglia, C., Camotim, D., Silvestre, N. (2019). GBT-based buckling analysis of steel cylindrical shells under combinations of compression and external pressure. Thin-Walled Structures, 144:106-274.
• Bazant, Z.P., Cedolin, L. (1991). Stability of Structures, Oxford University Press, New York, Oxford.
• Bebiano, R., Gonçalves, R., Camotim, D. (2015). A cross-section analysis procedure to rationalise and automate the performance of GBT-based structural analyses. Thin-Walled Structures, 92:29-47.
• Bleich, F. (1952). Buckling Strength of Metal Structures, McGraw-Hill, New York.
• Bonet, J., Wood, R.D., Mahaney, J., Heywood, P. (2000). Finite element analysis of air supported membrane structures. Computer Methods in Applied Mechanics and Engineering, 190(5-7):579-595.
• Camotim, D., Silvestre, N., Gonçalves, R., Dinis, P.B. (2006). GBT-based Structural Analysis of Thin-walled members: Overview, Recent Progress and Future Developments. Advances in Engineering Structures, Mechanics & Construction. Solid Mechanics and Its Applications, 140:187-204.
• Camotiom, D., Basaglia, C., Silvestre, N. (2010). GBT buckling analysis of thin-walled steel frames: A state-of-the-art report. Thin-Walled Structures, 48:726-743.
• Coda, H.B. (2003). An exact FEM geometric non-linear analysis of frames based on position description. In: 17th International Congress of Mechanical Engineering-COBEM, São Paulo.
• Coda, H.B. (2009). A solid-like FEM for geometrically non-linear 3D frames. Computer Methods in Applied Mechanics and Engineering, 198(47-48):3712-3722.
• Coda, H.B. (2015). Continuous inter-laminar stresses for regular and inverse geometrically non linear dynamic and static analyses of laminated plates and shells. Composite Structures, 132:406-422.
• Coda, H.B. and Paccola, R.R. (2008). A positional FEM Formulation for geometrical nonlinear analysis of shells. Latin American Journal of Solids and Structures, 5:205-223.
• Coda, H.B., Paccola, R.R., Carrazedo, R. (2017). Zig-zag effect without degrees of freedom in linear and nonlinear analysis of laminated plates and shells. Composite Structures, 161:32-50.
• Coda, H.B., Silva, A.P.O., Paccola, R.R. (2020). Alternative active nonlinear total Lagrangian truss finite element applied to the analysis of cable nets and long span suspension bridges. Latin American Journal of Solids and Structures,‏ 17(3):e268.
• Coda, H.B., Greco, M. (2004). A simple FEM formulation for large deflection 2D frame analysis based on position description. Computer Methods in Applied Mechanics and Engineering, 193(33-35):3541-3557.‏
• Coda, H.B., PACCOLA, R.R. (2007). An alternative positional FEM formulation for geometrically non-linear analysis of shells: Curved triangular isoparametric elements. Computational Mechanics, 40(1):185-200.
• Coda, H.B., PACCOLA, R.R. (2010). Improved finite element for 3D laminate frame analysis including warping for any cross-section. Applied Mathematical Modelling, 34(5):1107-1137.
• Crisfield, M.A. (1981). A fast incremental/iterative solution procedure that handles “snap-through”. Computers & Structures, 13(1):55-62.
• Dinis, P.B., Camotim, D., Silvestre, N. (2006). GBT formulation to analyse the buckling behaviour of thin-walled members with arbitrarily ‘branched’ open cross-sections. Thin-Walled Structures, 44:20-38.
• Feng, Y., Peric, D., Owen, D. (1996). A new criterion for determination of initial loading parameter in arc-length methods. Computers & Structures, 58(3):479-485.
• Fujji, F., Okazawa, S. (1997). Pinpointing Bifurcation Points and Branch-Switching. Journal of Engineering Mechanics, 123(3):179-189.
• Garcea, G. (2001). Mixed formulation in Koiter analysis of thin-walled beams. Computer Methods in Applied Mechanics and Engineering, 190(26):3369-3399.
• Geuzaine, C., Remacle, J.-F. (2009). Gmsh: a three-dimensional finite element mesh generator with built-in pre- and post-processing facilities. International Journal for Numerical Methods in Engineering, 79(11)1309-1331.
• Ghumare, S.M., Sayyad, A. S (2017). A New Fifth-Order Shear and Normal Deformation Theory for Static Bending and Elastic Buckling of P-FGM Beams. Latin American Journal of Solids and Structures, 14(11):1893-1911.
• Gonçalves, R., Dinis, P.B., Camotim, D. (2009). GBT formulation to analyse the first-order and buckling behaviour of thin-walled members with arbitrary cross-sections. Thin-Walled Structures, 47:583-600.
• Heydari, A., Jalali, A., Nemati, A. (2017). Buckling analysis of circular functionally graded plate under uniform radial compression including shear deformation with linear and quadratic thickness variation on the Pasternak elastic foundation. Applied Mathematical Modelling, 41:494-507.
• Iwata, K., Tsukimori, K., Kubo, F. (1991). A symmetric load-stiffness matrix for buckling analysis of shell structures under pressure loads. International Journal of Pressure Vessels and Piping, 45(1):101-120.
• Li, Z., Schafer, B.W. (2010). Application of the finite strip method in cold-formed steel member design. Journal of Constructional Steel Research, 66:971-980.
• Li, Z., Abreu, J.C.B., Leng, J., Ádány, S., Schafer, B.W. (2014). Review: Constrained finite strip method developments and applications in cold-formed steel design. Thin-Walled Structures, 81:2-18.
• Macneal, R.H. (1978). Simple Quadrilateral Shell Element. Computers & Structures, 8(2):175-183.
• Manta, D., Gonçalves, R., Camotim, D. (2020). Combining shell and GBT-based finite elements: Linear and bifurcation analysis. Thin-Walled Structures, 152:106665.
• Mirzaei, S., Azhari, M., Bondarabady, H.A.R. (2015). On the use of finite strip method for buckling analysis of moderately thick plate by refined plate theory and using new types of functions. Latin American Journal of Solids and Structures, 12(3):561-582.
• Mororo, L.A.T., Melo, A.M.C, Parente Junior, E. (2015). Geometrically nonlinear analysis of thin-walled laminated composite beams. Latin American Journal of Solids and Structures, ‏12(11):2094-2117.
• Murray, N.W. (1986). Introduction to the Theory of Thin-walled Structures, larendon Press, Oxford.
• Ogden, R.W. (1984) Nonlinear Elastic Deformation, Ellis Horwood, England.
• Paccola, R.R., Sampaio, M.S.M., Coda, H.B. (2016). Continuous stress distribution following transverse direction for FEM orthotropic laminated plates and shells. Applied Mathematical Modelling, 40:7382-7409.
• Pascon, J.P. and Coda, H.B. (2013). A shell finite element formulation to analyze highly deformable rubber-like materials. Latin American Journal of Solids and Structures, 10:1177-1209.
• Pastor, M.M., Roure, F. (2009). Open cross-section beams under pure flexural. II. Finite element simulation. Thin-Walled Structures, 47:514-521.
• Poorveis, D., Khajehdezfuly, A., Moradi, S., Shirshekana, E. (2019). A simple spline finite strip for buckling analysis of composite cylindrical panel with cutout. Latin American Journal of Solids and Structures, 16(8):e227.
• Ren, W.-X., Fang, S.-E., Young, B. (2006). Analysis and design of cold-formed steel channels subjected to combined flexural and web crippling. Thin-Walled Structures, 44:314-320.
• Riks, E. (1972). The Application of Newton’s Method to the Problem of Elastic Stability. Journal of Applied Mechanics, 39(4):1060-1065.
• Schwlizerhof, K., Ramm, E. (1984). Displacement dependent pressure loads in nonlinear finite element analyses. Computers & Structures, 18(6):1099-1114.
• Shahmohammadi, M.A., Azhari, M., Saadatpour, M.M., Sarrami-Foroushani, S. (2019). Stability of laminated composite and sandwich FGM shells using a novel isogeometric finite strip method. Engineering Computations, 37:1369-1395.
• Shi, J. (1996). Computing critical points and secondary paths in nonlinear structural stability analysis by the finite element method. Computers & Structures, 58:203-220.
• Silvestre, N., Duarte, A.P.C., Martins, J.P., Silva, L.S. da (2019). GBT Buckling Analysis of Cylindrical Panels Under Compression. Structures, 17:34-42.
• Siqueira, T.M., Coda, H.B. (2016). Development of Sliding Connections for Structural Analysis by a Total Lagrangian FEM Formulation. Latin American Journal of Solids and Structures, 13(11):2059-2087.
• Soares, H.B., Paccola, R.R., Coda, H.B. (2019). Unconstrained Vector Positional Shell FEM formulation applied to thin-walled members instability analysis. Thin-Walled Structures, 136:246-257.
• Sorensen, D.C., Lehoucq, R.B., Yang, C., Maschhoff, K. (2008). ARPACK, Rice University.
• Timoshenko, S.P., GERE, J.M. (1961). Theory of Elastic Stability, McGraw-Hill Book Co., New York.
• Wagner, W., Wriggers, P. (1988). A simple method for the calculation of postcritical branches. Engineering Computations, 5(2):103-109.
• Wempner, G.A. (1971). Discrete approximations related to nonlinear theories of solids. International Journal of Solids and Structures, 7(11):1581-1599.

Edited by

Editor:

Marco L. Bittencourt.

Publication Dates

• Publication in this collection
  15 Mar 2021
• Date of issue
  2021

History

• Received
  08 Sept 2020
• Reviewed
  18 Dec 2020
• Accepted
  22 Dec 2020
• Published
  08 Jan 2021