Large deflection and post-buckling of thin-walled structures by finite elements with node-dependent kinematics

Abstract

In the framework of finite elements (FEs) applications, this paper proposes the use of the node-dependent kinematics (NDK) concept to the large deflection and post-buckling analysis of thin-walled metallic one-dimensional (1D) structures. Thin-walled structures could easily exhibit local phenomena which would require refinement of the kinematics in parts of them. This fact is particularly true whenever these thin structures undergo large deflection and post-buckling. FEs with kinematics uniform in each node could prove inappropriate or computationally expensive to solve these locally dependent deformations. The concept of NDK allows kinematics to be independent in each element node; therefore, the theory of structures changes continuously over the structural domain. NDK has been successfully applied to solve linear problems by the authors in previous works. It is herein extended to analyze in a computationally efficient manner nonlinear problems of beam-like structures. The unified 1D FE model in the framework of the Carrera Unified Formulation (CUF) is referred to. CUF allows introducing, at the node level, any theory/kinematics for the evaluation of the cross-sectional deformations of the thin-walled beam. A total Lagrangian formulation along with full Green–Lagrange strains and 2nd Piola Kirchhoff stresses are used. The resulting geometrical nonlinear equations are solved with the Newton–Raphson linearization and the arc-length type constraint. Thin-walled metallic structures are analyzed, with symmetric and asymmetric C-sections, subjected to transverse and compression loadings. Results show how FE models with NDK behave as well as their convenience with respect to the classical FE analysis with the same kinematics for the whole nodes. In particular, zones which undergo remarkable deformations demand high-order theories of structures, whereas a lower-order theory can be employed if no local phenomena occur: this is easily accomplished by NDK analysis. Remarkable advantages are shown in the analysis of thin-walled structures with transverse stiffeners.

1 Introduction

Nowadays, in many engineering fields, increasingly sophisticated structures are employed to fulfill the tasks of demanding applications. Such structures require the adoption of enhanced calculation techniques to perform a prediction of their operative behavior, and then, for proper design. The need for a high level of accuracy and reliability from the structural simulation has led to the usage of high-performing three-dimensional (3D) models, with a drastic increase in the computational cost. In order to cut down this drawback, scientists have been pushed to the development of lighter one-dimensional (1D) models, with the goal of maintaining the same level of accuracy when compared to the heavier 3D tools.

The best-known 1D theories developed in history were made by Euler [1] and Timoshenko [2, 3]. These theories are also known as classical theories. The former ignores any cross-sectional deformation and neglects shear components. Thus, it fits only for 1D structures with a high slenderness ratio. The latter introduces uniform shear stress to overcome these limits, but the field of its application is locked in a narrow range. In fact, when dealing with structures whose cross-sectional deformation plays a crucial role, which is often the case for thin-walled sections, advanced structural theories must be considered for a 1D model, see the many chapters devoted to that in the classical book by Novizhilov [4]. Carrera et al. [5] and Kapania and Raciti [6, 7] presented a detailed discussion and review of refined models for the analysis of isotropic structures, whereas Reddy [8] reported an overview on the developed 1D Finite Elements (FEs) based on both classical and refined theories, including a discussion on the locking-free methods. As further examples of higher-order beam models proposed in the history, Vlasov [9] introduced warping functions to capture the deformations of beam cross sections. This approach found a great success between scientists, see the works by Ambrosini et al. [10], Mechab et al. [11] and Friberg [12], who make use of warping functions for thin-walled structures. A combination of the refined Vlasov model and the classical Euler-Bernoulli model was adopted by Kim and Lee [13] to analyze thin-walled beams made of functionally graded materials. The so-called generalized beam theory (GBT) was suggested by Schardt [14]. This theory allows the displacement field to be expressed as a linear combination of cross section deformation modes, which vary continuously along the member length. GBT found many applications in the literature, for example, by Peres et al. [15] for the analysis of curved thin-walled beams, and by Silvestre [16] for buckling problems. GBT was also adopted for the analysis of laminated materials, as presented by Silvestre and Camotim [17, 18].

In many real applications, local phenomena and large cross-sectional deformations occur in particular areas of the structure, for example, nearby of external loads or constraint conditions. In such cases, it would be convenient to develop a model with variable kinematics, namely, capable of refining only the areas where higher-order cross-sectional phenomena appear. In this way, the accuracy is still guaranteed, with a drastic decrease in the number of degrees of freedom and, subsequently, of the computational cost. When models with different kinematics have to be linked, the continuity of the displacements between the two regions has to be guaranteed. The modern literature about the issue of joining incompatible structural models is extremely vast. Wenzel [19] proposed an exhaustive state-of-the-art about this topic. For instance, the compatibility between different domains can be reached, making use of Lagrange multipliers. This approach was described and adopted by Prager [20] and Carrera et al. [21], among others. Another solution to this problem is the adoption of the global-local technique. Basically, this method consists of a multistep procedure, where, at first, a “global” analysis is carried out using a coarse mathematical model for the structure. Then, a refined FE model is applied separately in specific and more deformable subregions, and the compatibility is ensured by enforcing the continuity of the displacement in the interfacial or overlapping zones. Some authors adopted an iterative procedure between the global coarse mesh and the refined local one, see Whitcomb [22] who, with this method, retained the same accuracy than a refined global mesh. Noor [23] proposed a review on the global-local approaches for the nonlinear analysis of composite panels. The global-local approach found application in many engineering fields. For instance, Hanganu et al. [24] applied this method in the analysis of civil structures, for an accurate evaluation of the damage within the structure.

Particular attention must be given to local phenomena when structures are subjected to large deformation, e.g. large displacements and large rotations. In fact, for an accurate design of structures undergoing extreme loading conditions, a geometrically nonlinear analysis must be carried out. The contributions of scientists about the nonlinear analysis of 1D structures are uncountable. Most of the geometrically nonlinear models developed in the past are based on the Timoshenko beam theory, see, for example, Refs. [25,26,27]. However, by using variational asymptotic approaches, the 3D nature of such cases was divided into a two-dimensional (2D) analysis of the cross section and a 1D problem of the beam axis, see for instance the work by Yu et al. [28, 29]. Nevertheless, these solutions lack the ability to accurately catch higher-order phenomena, such as coupling or local effects, which may occur within a thin-walled structure. Thus, the adoption of refined model results is mandatory. For example, the post-buckling behavior of thin-walled beam structures was evaluated by an extension of GBT by Basaglia et al. [30], and using the Newton–Raphson linearization method along with the Ritz approach by Machado [31]. In both cases, the deformability of the structures in the large displacement field was caught by considering bending and warping effects over the cross section.

The present work intends to assess the benefits of adopting variable kinematics on a thin-walled beam structure in the geometrically nonlinear analysis. The solution proposed is the adoption of the node-dependent kinematics (NDK) in an FE framework based on the Carrera Unified Formulation (CUF) [32, 33]. Thanks to the scalable nature of CUF, any arbitrary expansion of the FE unknowns can be used to achieve the desired theory of structures. In other words, the primary unknowns of a given problem (which, for a 1D problem, are the displacements along the beam axis), are expanded using the expansion functions. They can change along the beam, and no problems about joining different expansion functions arise since the finite element method (FEM) is used between the two domains. NDK was used and validated in the past years by Carrera and Zappino [34] and applied to composite structures [35, 36], a two-dimensional (2D) plate [37, 38] and shell problems [39]. An application of NDK to the fluid-structure interaction can be found in [40]. The geometrically nonlinear 1D governing equations of the beam theory are obtained by means of the so-called fundamental nuclei, which allows the automatic employment of low- to higher-order theory, arbitrarily. This geometrically nonlinear solution was validated for isotropic and composite materials [41, 42], and, then, further extended to the dynamic [43, 44] and 2D plate [45] and shell cases [46, 47].

This paper is organized as follows: (i) NDK approach and the Green–Lagrange relations are presented in Sect. 2, along with the geometrically nonlinear FE equations ; (iii) then, numerical results are discussed for thin-walled beams in Sect. 3, with symmetric and asymmetric cross section and both transverse and compressive loadings; (iv) finally, the main conclusions are drawn.

2 Unified beam element with node-dependent kinematics

2.1 Preliminaries

Let us consider a generic 1D ‘beam-like’ structure as shown in Fig. 1. Clearly, the cross section $\varOmega$ lays on the xz-plane of a Cartesian reference system (x, y, z), and, consequently, the beam axis is placed along the y direction. Thus, the transposed displacement vector reads:

$$\begin{aligned} \varvec{\mathbf {u}}(x,y,z) = \left\{ \begin{array}{ccc} u_{x}&\quad u_{y}&\quad u_{z} \end{array}\right\} ^{T}. \end{aligned}$$

(1)

The 3D stress, $\varvec{\mathbf {\sigma }}$, and strain, $\varvec{\mathbf {\epsilon }}$, components are introduced in the following, with a vectorial notation:

$$\begin{aligned} \varvec{\sigma } = \left\{ \begin{array}{cccccc} \sigma _{xx}&\quad \sigma _{yy}&\quad \sigma _{zz}&\quad \sigma _{xz}&\quad \sigma _{yz}&\quad \sigma _{xy} \end{array}\right\} ^{T}, \quad \varvec{\epsilon } = \left\{ \begin{array}{cccccc} \epsilon _{xx}&\quad \epsilon _{yy}&\quad \epsilon _{zz}&\quad \epsilon _{xz}&\quad \epsilon _{yz}&\quad \epsilon _{xy} \end{array}\right\} ^{T}. \end{aligned}$$

(2)

These will hereafter coincide with the components of the Green–Lagrange and Second Piola–Kirchhoff stress tensors, respectively. In fact, the aim of evaluating higher-order coupling phenomena, such as torsion, bending, and extension, forced the adoption of the full Green–Lagrange strains. Thus, the geometrical relations take the following form:

$$\begin{aligned} \begin{array}{lllll} \varvec{\mathbf {\epsilon }}&{} = &{} \varvec{\mathbf {\epsilon }}_l + \varvec{\mathbf {\epsilon }}_{nl} &{} = &{} (\varvec{\mathbf {b}}_l + \varvec{\mathbf {b}}_{nl}) \varvec{\mathbf {u}} \\ \end{array} \end{aligned}$$

(3)

where $\varvec{\mathbf {b}}_l$ consists of the linear differential operators and $\varvec{\mathbf {b}}_{nl}$ is the nonlinear differential operator, namely the aforementioned Green–Lagrange strain tensor. Interested readers can find the expression of $\varvec{\mathbf {b}}_l$ and $\varvec{\mathbf {b}}_{nl}$ in [41].

A linear elastic isotropic metallic is analyzed in this work. Consequently, the constitutive relation reads as:

$$\begin{aligned} \varvec{\mathbf {\sigma }}= {\varvec{\mathbf {C}}}\varvec{\mathbf {\epsilon }}, \end{aligned}$$

(4)

where ${\varvec{\mathbf {C}}}$ is the material elastic matrix of homogenoeus and isotropic materials. The explict form can be found in any book of elasticity. An exhaustive finite element treatment of the topic can be found in [48, 49].

Fig. 1

One-dimensional structure defined over a Cartesian reference system

Within the framework of the Carrera Unified Formulation (CUF), the three-dimensional (3D) displacement field $\varvec{\mathbf {u}}(x,y,z)$ as well as its variation (denoted by $\delta$) of 1D models can be expressed as a general expansion of the primary unknowns, as follows:

$$\begin{aligned} \begin{array}{cccc} \varvec{\mathbf {u}}(x,y,z) &{} = &{} F_{\tau }(x,z) \; \varvec{\mathbf {u}}_{\tau }(y), \quad \ &{} \tau = 1,2,\ldots ,M, \\ \delta {\varvec{\mathbf {u}}}(x,y,z) &{} = &{}F_{s}(x,z) \; \delta \varvec{\mathbf {{u}}}_{s}(y), \quad &{} s = 1,2,\ldots ,M \end{array} \end{aligned}$$

(5)

where $F_{\tau }$ and $F_{s}$ (see Fig. 1 in red) represent the base 2D functions over the xz cross section plane of the 1D structure, $\varvec{\mathbf {u}}_{\tau }$ and $\varvec{\mathbf {\delta u}}_{s}$ are the vectors of the unknowns displacements along the beam axis, M denotes the order of expansion, and the summing convention with the repeated indicies $\tau$ and s is assumed. The 1D theory to be employed is strictly related to the $F_{\tau }$. In other words, the choice of the expansion functions determines the order of the adopted model as well as its capability to trace cross-sectional deformations. The research work proposed in this paper makes use of both Taylor Expansion (TE) and Lagrange Expansion (LE) [50]. The concept of TE is depicted in Fig. 2a, where the 3D displacement of a generic point “A” is a function of the 3D displacement of the structural node and its derivatives.

Fig. 2

Example of TE and LE for the cross-sectional deformation of thin-walled beam

The complete expression for an TE of order two (TE2) is reported in Eq. 6:

$$\begin{aligned} \begin{array}{ccccccccccccccccccccccc} u_{x_{A}}(x_{A},y_{A},z_{A}) &{} = &{} u_{x_{1}}(y_{A})&{}+&{}x_{A} \; u_{x_{2}}(y_{A})&{}+&{}z_{A} \; u_{x_{3}}(y_{A})&{}+&{}x_{A}^{2} \; u_{x_{4}}(y_{A})&{}+&{}x_{A}z_{A} \; u_{x_{5}}(y_{A})&{}+&{}z_{A}^{2} \; u_{x_{6}}(y), \\ u_{y_{A}}(x_{A},y_{A},z_{A}) &{} = &{} u_{y_{1}}(y_{A})&{}+&{}x_{A} \; u_{y_{2}}(y_{A})&{}+&{}z_{A} \; u_{y_{3}}(y_{A})&{}+&{}x_{A}^{2} \; u_{y_{4}}(y_{A})&{}+&{}x_{A}z_{A} \; u_{y_{5}}(y_{A})&{}+&{}z_{A}^{2} \; u_{y_{6}}(y), \\ u_{z_{A}}(x_{A},y_{A},z_{A}) &{} = &{} u_{z_{1}}(y_{A})&{}+&{}x_{A} \; u_{z_{2}}(y_{A})&{}+&{}z_{A} \; u_{z_{3}}(y_{A})&{}+&{}x_{A}^{2} \; u_{z_{4}}(y_{A})&{}+&{}x_{A}z_{A} \; u_{z_{5}}(y_{A})&{}+&{}z_{A}^{2} \; u_{z_{6}}(y) \\ \end{array}\nonumber \\ \end{aligned}$$

(6)

where $x_{A}$ and $z_{A}$ are the cross section coordinates of the point “A’, as reported in Fig. 2a. The number of the degrees of freedom (DOFs) is equal to the displacement and derivatives of the TE, and, for the case of a TE2, they are 18. Figure 2b shows an example of LE. Namely, the cross section is discretized with a set of Lagrange points (LP), opportunely subdivided into Lagrange polynomilals, locally defined. The 3D displacement of the point “A” is, then, an interpolation of the displacements evaluated at the LP of the belonging LE (in the Figure marked as red). The grade of the interpolation is defined by the number of the LP of the element, for instance, a 4 points LE (L4) ensures a linear interpolation, a 9 points LE (L9) a quadratic interpolation, and a 16 point LE (L16) a cubic interpolation. The LE used in this work is the L9, since it was demonstrated that it provides a high level of accuracy for a thin-walled structure (see [50]). The number of DOFs is equal to the number of the displacement for each LP, and, for the case of a 45LP cross section, it is 135. Equation (7) reports the expression of the 3D displacement of the generic point “A”, as follows:

$$\begin{aligned} \begin{array}{ccccccccccccccccccccccc} u_{x_{A}}(x_{A},y_{A},z_{A}) &{} = &{} F_{1}(x_{A},z_{A})u_{x_{1}}(y_{A})&{}+&{}F_{2}(x_{A},z_{A})u_{x_{2}}(y_{A})&{}+&{}\ldots &{}+&{}F_{9}(x_{A},z_{A})u_{x_{9}}(y_{A}), \\ u_{y_{A}}(x_{A},y_{A},z_{A}) &{} = &{} F_{1}(x_{A},z_{A})u_{y_{1}}(y_{A})&{}+&{}F_{2}(x_{A},z_{A})u_{y_{2}}(y_{A})&{}+&{}\ldots &{}+&{}F_{9}(x_{A},z_{A})u_{y_{9}}(y_{A}), \\ u_{z_{A}}(x_{A},y_{A},z_{A}) &{} = &{} F_{1}(x_{A},z_{A})u_{z_{1}}(y_{A})&{}+&{}F_{2}(x_{A},z_{A})u_{z_{2}}(y_{A})&{}+&{}\ldots &{}+&{}F_{9}(x_{A},z_{A})u_{z_{9}}(y_{A}) \\ \end{array} \end{aligned}$$

(7)

where the $F_{\tau }$ are not fully reported here for the sake of brevity, but can be found in [33]. The finite element method (FEM) is adopted to discretize the structure along the y axis. Thus, the generalized displacement vector $\varvec{\mathbf {u}}_{s}(y)$ and its variation are approximated as follows:

$$\begin{aligned} \begin{array}{ll} \varvec{\mathbf {u}}_{\tau }(y) =&{} N_{i}(y) \varvec{\mathbf {q}}_{\tau i} \quad i=1,2,\ldots ,N_n, \\ \delta \varvec{\mathbf {u}}_{s}(y) =&{} N_{j}(y) \delta \varvec{\mathbf {q}}_{s j} \quad j=1,2,\ldots ,N_n \end{array} \end{aligned}$$

(8)

where $N_{i}(y)$ and $N_j(y)$ stand for the i, jth 1D shape function, $N_n$ is the number of the structural nodes, $\varvec{\mathbf {q}}_{\tau i}$ is the vector of the FE nodal parameters, and i indicates summation. Among the shape functions $N_{i}$, which are described here (but a comprehensive description of the topic is reported in [48]), a cubic interpolation is assumed in this work. Basically, four nodes FEs (B4) are employed to approximate the 1D structures along the y axis Finally, the formula including Eqs. (5) and (8) reads as:

$$\begin{aligned} \begin{array}{ll} \varvec{\mathbf {u}}(x,y,z) = F_{\tau }(x,z)N_{i}(y) \varvec{\mathbf {q}}_{\tau i}, \\ \delta \varvec{\mathbf {u}}(x,y,z) = F_{\tau }(x,z)N_{i}(y) \delta \varvec{\mathbf {q}}_{\tau i}. \end{array} \end{aligned}$$

(9)

2.2 The node-dependent kinematics concept

As discussed in the Introduction, the possibility to couple local-to-global parts of the FE model can be done in many ways. NDK allows, by definition, to vary the kinematics node-by-node in the same element, by leading to the possibility to refine kinematics without any use of coupling mathematical methods. The refinement of the adopted kinematics can be taken a step further by assigning an own approximation to each node of an element. Thanks to the scalable nature of the CUF-based displacement models, it is possible to develop an FE with variable kinematic. The basic idea NDK is to describe the displacement field over each cross section of an element with different kinematics. In this way, one can refine the model only over the regions which require a higher-order theory to be accurately described, and associate lower-order theories in other domains, where particular phenomena do not occur, saving computational cost. That is, the $F_\tau , F_s$ functions are node-dependent, that is, the indices i, j of the nodes should distinguish the kinematics at each node. Then, Eq. (9) becomes:

$$\begin{aligned} \begin{array}{ll} \varvec{\mathbf {u}}(x,y,z) &{}= F^{i}_{\tau }(x,z)N_{i}(y) \varvec{\mathbf {q}}_{\tau i}, \quad \tau = 1,2,\ldots ,M^{i}, \\ \delta \varvec{\mathbf { u}}(x,y,z) &{}= F^{j}_{s}(x,z)N_{j}(y) \delta \varvec{\mathbf {q}}_{ s j}, \quad s = 1,2,\ldots ,M^{j} \end{array} \end{aligned}$$

(10)

where the indices i, j on $M^{i}, M^j$, and $F^{i}_{\tau }(x,z), F^{j}_{s}(x,z)$ underline that the expansion functions are associated with the i, jth node of the element, rather than with the entire element itself. An example of NDK approach is finally shown in Fig. 3, where a compact beam with four nodes is shown. The node 1 has a TE1 expansion, that is, Eq. (6) truncated at linear term and the number its DOFs is 9; node 2 is expanded with a TE2 [the Eq. (6)], with 18 DOFs.; node 3 and node 4 are approximated with 9 and 25 LP, respectively, so that their DOFs are 27 and 75, respectively. The number of the DOFs is the sum of the DOFs of every cross section, i.e.. $9 + 18 + 9\times 3 + 25\times 3 = 129$. Many examples will be discussed in the numerical results Section.

Fig. 3

NDK approach applied to a 1D element. Expansion functions (in red) associated individually with each node. The structure is modeled with 4 structural nodes, approximated with a TE1, TE2, 9LP and 25LP, respectively. Then, the number of the DOFs is $9 + 18 + 9\times 3 + 25\times 3 = 129$ (color figure online)

2.3 Nonlinear governing equations

The principle of virtual work is hereafter recalled for the derivation of the nonlinear FE governing equations. It states that the virtual work from the internal strain energy ($\delta L_\mathrm {int}$) is equal to the one made by the external loads ($\delta L_\mathrm {ext}$). Moreover, the geometrically nonlinear problem is obtained by introducing Eqs. (8) and (10) into Eq. (3), so that the strain vector can be written in algebraic form as follows:

$$\begin{aligned} \begin{array}{lll} \varvec{\mathbf {\epsilon }}= & {} (\varvec{\mathbf {B}}_l^{\tau i} + \varvec{\mathbf {B}}_{nl}^{\tau i}) \varvec{\mathbf {q}}_{\tau i} \end{array} \end{aligned}$$

(11)

where $\varvec{\mathbf {B}}_l^{\tau i}$ and $\varvec{\mathbf {B}}_{nl}^{\tau i}$ are the linear and nonlinear algebraic matrices with CUF and FEM formulations. For the sake of completeness, these operators are given below,

$$\begin{aligned} \varvec{\mathbf {B}}_l^{\tau i} = \left[ \begin{array}{ccc} F^{i}_{\tau _{,x}} N_{i} &{}\quad 0 &{}\quad 0 \\ 0 &{}\quad F^{i}_\tau N_{i_{,y}} &{}\quad 0 \\ 0 &{}\quad 0 &{}\quad F^{i}_{\tau _{,z}} N_{i} \\ F^{i}_{\tau _{,z}} N_{i} &{}\quad 0 &{}\quad F^{i}_{\tau _{,x}} N_{i} \\ 0 &{}\quad F^{i}_{\tau _{,z}} N_{i} &{}\quad F^{i}_{\tau } N_{i_{,y}} \\ F^{i}_{\tau } N_{i_{,y}} &{}\quad F^{i}_{\tau _{,x}} N_{i} &{}\quad 0 \end{array}\right] \end{aligned}$$

(12)

and

$$\begin{aligned} \varvec{\mathbf {B}}_{nl}^{\tau i} = \dfrac{1}{2} \left[ \begin{array}{ccc} u_{x_{,x}} F^{i}_{\tau _{,x}} N_i &{}\quad u_{y_{,x}} F^{i}_{\tau _{,x}} N_i &{}\quad u_{z_{,x}} F^{i}_{\tau _{,x}} N_i \\ u_{x_{,y}} F^{i}_\tau N_{i_{,y}} &{}\quad u_{y_{,y}} F^{i}_\tau N_{i_{,y}} &{}\quad u_{z_{,y}} F^{i}_\tau N_{i_{,y}} \\ u_{x_{,z}} F^{i}_{\tau _{,z}} N_i &{}\quad u_{y_{,z}} F^{i}_{\tau _{,z}} N_i &{}\quad u_{z_{,z}} F^{i}_{\tau _{,z}} N_i \\ u_{x_{,x}} F^{i}_{\tau _{,z}} N_i + u_{x_{,z}} F^{i}_{\tau _{,x}} N_i &{}\quad u_{y_{,x}} F^{i}_{\tau _{,z}} N_i + u_{y_{,z}} F^{i}_{\tau _{,x}} N_i &{}\quad u_{z_{,x}} F^{i}_{\tau _{,z}} N_i + u_{z_{,z}} F^{i}_{\tau _{,x}} N_i \\ u_{x_{,y}} F^{i}_{\tau _{,z}} N_i + u_{x_{,z}} F^{i}_{\tau } N_{i_{,y}} &{}\quad u_{y_{,y}} F^{i}_{\tau _{,z}} N_i + u_{y_{,z}} F^{i}_{\tau } N_{i_{,y}} &{}\quad u_{z_{,y}} F^{i}_{\tau _{,z}} N_i + u_{z_{,z}} F^{i}_{\tau } N_{i_{,y}} \\ u_{x_{,x}} F^{i}_{\tau } N_{i_{,y}} + u_{x_{,y}} F^{i}_{\tau _{,x}} N_i &{}\quad u_{y_{,x}} F^{i}_{\tau } N_{i_{,y}} + u_{y_{,y}} F^{i}_{\tau _{,x}} N_i &{}\quad u_{z_{,x}} F^{i}_{\tau } N_{i_{,y}} + u_{z_{,y}} F^{i}_{\tau _{,x}} N_i \end{array}\right] . \end{aligned}$$

(13)

The variation of the elastic internal work, considering constitutive (Eq. (4)) and geometrical relations (Eq. (8)), can be expressed as:

$$\begin{aligned} \begin{array}{lll} \delta L_\mathrm {int} &{} = &{} \displaystyle \int _{V} \delta \varvec{\mathbf {\epsilon }}^{T} \varvec{\mathbf {\sigma }} \;\mathrm {d}V \\ &{} = &{} \delta \varvec{\mathbf {q}}_{s j}^T \; \displaystyle \left( \displaystyle \int _{V} \left( \varvec{\mathbf {B}}_l^{s j} + 2 \, \varvec{\mathbf {B}}_{nl}^{s j} \right) ^T \varvec{\mathbf {C}} \left( \varvec{\mathbf {B}}_l^{\tau i} + \varvec{\mathbf {B}}_{nl}^{\tau i} \right) \;\mathrm {d}V\displaystyle \right) \; \varvec{\mathbf {q}}_{\tau i} \\ &{}=&{}\delta \varvec{\mathbf {q}}_{s j}^T \; \varvec{\mathbf {K}}_{S}^{i j \tau s} \; \varvec{\mathbf {q}}_{\tau i} \end{array} \end{aligned}$$

(14)

where $\varvec{\mathbf {B}}_l^{s j}$ and $\varvec{\mathbf {B}}_{nl}^{s j}$ come out from the variation of the strain components (Eq. (11)), and $\varvec{\mathbf {K}}_S^{i j \tau s}$ represents the so-called secant stiffness matrix. The complete form of the secant stiffness matrix $\varvec{\mathbf {K}}_S^{i j \tau s}$ can be found in [41, 51].

Omitting some mathematical steps, which interested readers can find in [52], the principle of virtual work, for the whole structure, becomes:

$$\begin{aligned} \varvec{\mathbf {K}}_{S} \; \varvec{\mathbf {q}} - \varvec{\mathbf {p}} = 0 \end{aligned}$$

(15)

where $\varvec{\mathbf {K}}_{S}$, $\varvec{\mathbf {q}}$, and $\varvec{\mathbf {p}}$ are the global FE arrays of the structure, with $\varvec{\mathbf {p}}$ corresponding to the loading vector.

For the resolution of the system of algebraic nonlinear equations (Eq. 15), in this work, the total Lagrangian formluation as the one presented in detail by Pagani and Carrera [42] is adopted. Nonlinear algebraic equations are solved via Newton–Raphson scheme, which requires the linearization of the equations. The related tangent stiffness matrix $\varvec{\mathbf {K}}_{T}$ can be obtained by the second variation of the strain energy at the equilibrium point, as follows:

$$\begin{aligned} \delta ^2 L_\mathrm {int} = \delta \varvec{\mathbf {q}}_{s j}^T K_T^{i j \tau s} \delta \varvec{\mathbf {q}}_{\tau i}. \end{aligned}$$

(16)

The explicit form of $\varvec{\mathbf {K}}_{T}$ is not given here, but it is derived in a unified form in [46]. Finally, a constraint is needed fo the resultant system of equations, and the arc-length method is employed (see the works by Carrera [53] and Crisfield [54, 55] for more details).

3 Numerical results

The process for the use of NDK methodology for nonlinear problems can be summarized in the following steps.

• At first, a mathematical model using B4 elements for the axis discretization and L9 for the cross section approximation (Fig. 4a) is built.

• A geometrical nonlinear analysis is performed, and the deformed configurations of the structure are investigated in the overall equilibrium path, from moderate to large displacements and rotations (Fig. 4b). Part of the structures could exhibit very high sectional deformations, other parts would shows quite rigid body motion for the sections.

• In this way, we can know which parts of the structure present a remarkable deformation (Fig. 4c in blue) and, therefore, require a refined kinematic for their description and which can be approximated by a lower-order theory (Fig. 4c in yellow).

• Finally, the node-dependent kinematic model is built, following the aforementioned assumptions (Fig. 4d), by introducing higher node kinematics in the zone related to the blue part and less refined ones in the yellow zone.

Fig. 4

Logical steps to perform the geometrically nonlinear analysis using the NDK approach

Various problems are addressed for demonstrating the application of NDK models in the large deflection field. For the analyzed cases, a reference solution is obtained by comparison with literature results, and proposed NDK models are compared, both in terms of displacement and stress distribution. B4 finite elements are employed in the analysis cases, and the notation used hereafter to identify the structural models involves the number of B4 elements approximated with LE or TE, followed by the number of FEs, so that, for example, “10LE” stands for “10B4 elements with Lagrange Expansion” and “5TE4” stands for “5B4 elements with Taylor Expansion of order 4”.

3.1 Thin-walled channel beam

The effects of NDK models in the large deflection field are shown and explained with the comparison with the results proposed by Carrera and Zappino [34], where a thin-walled channel beam was analyzed in the geometrically linear regime. The material is isotropic with Young’s modulus $E = 71.7$ GPa and Poisson’s ratio $\nu = 0.3$. The geometry is shown in Fig. 5a, with sides $h = b = 0.1$ m, length $L = 1$ m, and thickness $t = 0.005$ m. The structure was loaded at its tip with two different load configurations. Load case 1, described in Fig. 5b, involves two loads in the opposite direction and the same magnitude, whereas load case 2 (Fig. 5c) has one vertical load. The Figures also report the adopted theories to establish a reference solution. These theories were taken from the aforementioned paper, and they involve 10L9 for load case 1 (see Fig. 5b) and 8L9 for load case 2 (see Fig. 5c). For both load cases, eight B4 are adopted to approximate the displacement over the beam axis in the y direction.

Fig. 5

Geometric properties (a), loading condition and LE approximation of the cantilever thin-walled beam with symmetric C-section; 10L9 (b) and 8L9 (c)

Figure 6 reports the nonlinear static equilibrium curve for the load case 1, adopting a uniform refined kinematic model with 10LE (Fig. 5b). The Figure also shows two deformed configurations over the undeformed shape, for $P = 2$ kN and $P = 6.25$ kN. This solution was taken as the reference one for all subsequent NDK analyses for load case 1.

Fig. 6

Nonlinear static equilibrium curve of the cantilever thin-walled beam with symmetric C-section. Transverse displacement is evaluated at point “A” of the cross section (see Fig. 5b)

The NDK model of the beam structure is explained in Fig. 7. Basically, TE is employed in the structure gradually from the clamped area for every B4 element. Namely, in Fig. 7, the model described is 5TE–3LE. TE2 is adopted for load case 1, and TE4 is employed for load case 2.

Fig. 7

Description of the NDK model of the cantilever beam with symmetric C-section

The nonlinear equilibrium paths of analyzed NDK models are reported in Fig. 8. Clearly, until 4 FEs are approximated with TE2, the calculated displacement in the large deflection field is basically the same. A slight difference appears when employing 5TE, and it increases as the number of TE increases, going to the full TE model where the calculated displacement is completely wrong. The deformed tip cross sections for $P= 2$ kN are depicted for 6TE2–2LE, 7LE2–1LE and 8TE2 over the 8LE one, for completeness purpose.

Fig. 8

Nonlinear static equilibrium curve of the cantilever beam with symmetric C-section adopting the NDK approach with TE2 and LE (see Fig 5b). Transverse displacement is evaluated at point “A” of the cross section (see Fig. 5b)

In Table 1, the values of transverse displacement are reported. The linear solution for $P = 1$ kN is taken from the reference work [34], and the nonlinear investigation was performed for $P = 2$ kN, $P = 4$ kN, and $P = 6$ kN. The % error is evaluated with reference to the 8LE solution, and the trend is reported in Fig. 9. Evidently, the difference in the nonlinear regime between NDK models and the reference solution increases, but the effectiveness of the NDK model is independent of the load magnitude.

Table 1 Vertical displacement evaluated using different structural models

Fig. 9

$u_{zA}$ percentage error with respect to the reference solution of uniform LE refined model. Cantilever beam with symmetric C-section analysis case

Finally, an analysis of the stress condition is performed, and the obtained results are reported in Fig. 10. The components $\sigma _{xx}$ and $\sigma _{xz}$ are evaluated along the length of the beam for reference solution and three NDK models. The results demonstrate that even if the evaluated displacements do not show a relevant difference between reference and NDK models, the stress condition may not be accurately calculated. For instance, as reported in Fig. 8, the difference between 8LE and 5TE2–3LE displacements is almost null, whereas, in the range between $y/L = 0.5$ and $y/L = 0.8$, for $\sigma _{xx}$ and $\sigma _{xz}$ the difference is remarkable. From another point of view, it can be pointed out that the difference between the calculated stress distributions does not afflict the displacements of the structure in the nonlinear regime. This is the case for 6TE2–2LE, where the evaluated $\sigma _{xx}$ and $\sigma _{xz}$ in the range between $y/L = 0.5$ and $y/L = 0.9$ are relevant, but the difference between the displacement is smaller, reaching its maximum for 6 = kN with a 4.14% difference.

Fig. 10

Distribution of $\sigma _{xx}$ (a) and $\sigma _{xz}$ (b) components evaluated along the length of the beam at point “C” and “D”, respectively. Cantilever beam with symmetric C-section analysis case. $\sigma _{xx}^{*} = \dfrac{\sigma _{xx} \, h \, t}{P}$ and $\sigma _{xz}^{*} = \dfrac{\sigma _{xz} \, h \, t}{P}$

As far as the load case 2 is concerned, Fig. 11 reports the nonlinear static equilibrium curves for the transverse displacement of points “A” and “B” (see Fig. 5c). Deformed configurations are depicted, too, for $P = 6.25$ kN and 8.25 kN.

Fig. 11

Nonlinear static equilibrium curve of the cantilever thin-walled beam with symmetric C-section. Transverse displacements are evaluated at point “A” (red lines) and “B” (blue lines) of the cross section (see Fig. 5c) (color figure online)

Figure 12 reports the nonlinear curves with the NDK models. We underline that as the number of FEs with TE increases, the difference from the reference solution increases, except for the transverse displacement of the point “B” (see Fig. 5c), where the 8TE4 theory ensures a slightly more accurate result than the 7TE4–1LE model.

Fig. 12

Nonlinear static equilibrium curve of the cantilever beam with symmetric C-section adopting the NDK approach with TE4 and LE (see Fig 5c). Transverse displacements are evaluated at points “A” (a) and “B” (b) of the cross section (see Fig. 5c)

Finally, the values of the displacements are reported in Tables 2 and 3, and the % error distribution is reported in Fig. 13.

Table 2 Vertical displacement evaluated at point “A” (see Fig. 5c) using different structural models

Table 3 Vertical displacement evaluated at point “B” (see Fig. 5c) using different structural models

Fig. 13

$u_{zA}$ (a) and $u_{zB}$ (b) percentage error with respect to the reference solution of uniform LE (see Fig. 5c) refined model

3.2 Asymmetric C-section beam subjected to transverse loading

A cantilever beam undergoing large deflection due to a transverse loading was considered in the following analysis case. The considered material is aluminum alloy, as in the previous case. The cross-sectional shape and dimensions are shown in Fig. 14a. It is a asymmetric C-section with $b = 100$ mm, $h_{1} = 48$ mm, $h_{2} = 88$ mm, and $t = 10$ mm. Figure 14b reports the cross section discretization implementing 12L9. This polynomial pattern will be recalled as “LE” in the following analyses. 20 B4 elements are employed for the approximation of the beam axis.

Fig. 14

Geometric properties and loading condition (a) and Lagrange polynomial discretization (b) of the cantilever beam with asymmetric C-section undergoing large deflection due to a transverse loading

The static nonlinear analysis was performed using LE as the expansion function of every structural node. The result is shown in Fig. 15, along with some deformed configurations in the large deflection field. Clearly, due to the asymmetric shape of the cross section, a warping phenomenon occurs in the structure in the range between 25 and 50 kN. Then, it disappears as the load increases, as shown in the last depicted deformed configuration.

Fig. 15

Nonlinear static equilibrium curve of the cantilever beam with asymmetric C-section. Transverse displacement is evaluated at point “A” of the cross section (see Fig. 14)

Subsequently, the NDK approach was applied to the structure. The part of the beam near the clamp undergoes a remarkable cross section deformation (as clearly depicted in the deformed configurations reported in Fig. 15), and for this reason, a higher-order model is required to evaluate the displacements of this portion of the structure properly. Figure 16 shows the different kinematic assumptions used on the structural nodes along the beam axis.

Fig. 16

Description of the NDK model of the cantilever beam with asymmetric C-section undergoing large deflection due to a transverse loading

Four different patterns of refined and lower-order models are employed. The NDK models used in this analysis case are 7LE–13TE, 10LE–10TE, 10LE–5TE–5TE, and 15LE–5TE. A preliminary convergence analysis is carried out to check if the FE approximation along the beam axis adopted for the analysis of the uniform structure is appropriate also for NDK models. The first 7LE–13TE case is taken as example, and the results, involving 10 to 25B4, are reported in Fig. 17. Clearly, the best solution in terms of accuracy and computational cost is represented by 20B4. This approximation will be used for all the subsequent analyses, if not otherwise described. Then, nonlinear analyses are performed on the structure with the aforementioned theoretical assumptions. Figure 18 shows the results for the 7LE–13TE (Fig. 18a) and 10LE–10TE (Fig. 18b) NDK models. Every figure also shows the deformed configurations of the NDK structure at $P=25$ kN over the one from the reference configuration with 20LE. Moreover, sections at $y=400$ mm are highlighted to compare the results. Clearly, Fig. 18a shows how 7LE from the clamp are not enough to accurately evaluate the buckling-like phenomena which occur near $P = 25$ kN, as also highlighted by the sections at 400 mm. For this reason, Fig. 18b reports the nonlinear analysis with 10LE, and the other part of the beam approximated with TE1, TE3, and TE5. It demonstrates that the higher is the order of TE, the more accurate are the results, compared to the 20LE solution. Figure 19 reports more accurate results with variable kinematics over the beam length (10LE–5TE8–5TE1) and with 15LE–5TE. The higher accuracy of these models is proved by depicted sections. We have to point out that the DOF gain is lower than the previously shown results.

Fig. 17

Convergence analysis of the NDK 7LE–13TE model of the cantilever beam with asymmetric C-section undergoing large deflection due to a transverse loading. $e\% = \dfrac{u-u_{25B4}}{u_{25B4}}$

Fig. 18

Nonlinear static equilibrium curve of the cantilever beam with asymmetric C-section adopting the NDK approach with 7LE–13TE (a) and 10LE–10TE (b). Transverse displacement is evaluated at point “A” of the cross section (see Fig. 14). Cantilever beam with asymmetric C-section undergoing large deflection due to a transverse loading

Figure 19 reports the solution for the 10LE–5TE–5TE (Fig. 19a) and 15LE–5TE (Fig. 19b). We underline the increased accuracy of such models, in particular in catching the warping effects of the structure near $P=25$ kN.

Fig. 19

Nonlinear static equilibrium curve of the cantilever beam with asymmetric C-section adopting the NDK approach with 10LE–5TE–5TE (a) and 15LE–5TE (b). Transverse displacement is evaluated at point “A” of the cross section (see Fig. 14). Cantilever beam with asymmetric C-section undergoing large deflection due to a transverse loading

Finally, Tables 4 and 5 report the values of the transverse displacement for $P=25$, 30, and 40 kN, and $P=50$, 60, 120 kN, respectively.

3.3 Asymmetric C-section beam subjected to compression loading

The capability of the NDK approach to deal with the post-buckling of a thin-walled structure is tested in the following analysis case. The structure is the same as in the previous analysis case, with the same material and geometrical conditions. The cross-sectional discretization, used to assess a reference solution, is shown in Fig. 14, and this element pattern will be recalled as “LE” in the following reasonings. Figure 20 shows the post-buckling behavior of the structure, with the displacements over the x and z directions of points “A” and “B”. Two deformed configurations are reported, too, at $P^{*}= 5.9$ and $P^{*}= 7.9$, to appreciate the cross-sectional deformation in the middle part of the structure.

Table 4 Values of transverse displacement of the cantilever beam with asymmetric C-section undergoing large deflection due to a transverse loading for $P=25$, 30 and 40 kN

Table 5 Values of transverse displacement of the cantilever beam with asymmetric C-section undergoing large deflection due to a transverse loading for $P=50$, 60 and 120 kN

Fig. 20

Nonlinear static equilibrium curve of the cantilever beam with asymmetric C-section subjected to compression loading. A uniform LE is used to approximate the cross section displacement field (see Fig. 14). Transverse displacement is evaluated at points “A” and “B”. $P^{*} = \dfrac{P 4 L^{2}}{\pi ^{2} E I }$

Fig. 21

Nonlinear static equilibrium curve of the cantilever beam with asymmetric C-section subjected to compression loading adopting the NDK approach with LE and TE5. Transverse displacement is evaluated at points “A” (a) and “B” (b) of the tip cross section (see Fig. 20). Deformed sections refers to $y=\dfrac{L}{2}$

Then, the NDK was applied to this problem. Clearly, the middle part of the structure undergoes remarkable cross-sectional deformation. Thus, different and lower-order kinematics was applied to the elements near the clamp at the end of the structure. The results are shown in Fig. 21. Figure 21a shows the z component of the displacement of the point “A” (see Fig. 20). When NDK is applied only at the beginning of the structure, solutions similar to the reference one are achieved. The “3TE5–17LE” approximation provides very accurate results, although the computational cost-benefit is not so evident, with a gain of only 12% of DOFs. The “5TE5–15LE” configuration, with a DOF gain equals 19%, reports a higher error than the previous one, showing the importance of having a higher-order theory as we move closer to the middle portion of the beam. Interesting results are carried out by the “3TE5–14LE–3TE5” NDK. They provide slightly better results than the previous configuration, with a higher DOF gain, being 24%. The last two approximations, “5TE5–10LE–5TE5” and “7TE5–6LE–7TE5” provide not reliable results, as clearly shown also by the deformed sections, depicted in the Figure. Finally, Fig. 21b reports the displacement of the point “B” (see Fig. 20) in the x-direction. We can point out the same conclusions as in the previous Figure, although generally better results are achieved by the NDK assumptions, especially in the range between $P^{*}= 6.5$ and $P^{*}= 7.5$.

Fig. 22

LE approximation of the stiffener

Fig. 23

C-section beam with 3 (a) and 6 (b) stiffeners (in red) (color figure online)

Fig. 24

Nonlinear static equilibrium curve of the cantilever reinforced beam with asymmetric C-section. A uniform LE is used to approximate the cross section displacement field (see Fig. 14). Transverse displacement is evaluated at points “A” at $y = 1000$ mm (a) and “B” at $y = 250$ mm (b)

Fig. 25

Stiffeners’ deformation in the moderate (a, b) and large (c, d) displacement field

Fig. 26

Nonlinear static equilibrium curve of the cantilever reinforced beam with asymmetric C-section. NDK is employed to evaluate the displacement of the depicted points “A” (a) and “B” (b). Adopted notation for each structural bay between two stiffeners: $\circ$ = LE, $\blacksquare$ = TE5, $\square$ = TE1

Fig. 27

Nonlinear static equilibrium curve of the cantilever reinforced beam with asymmetric C-section. NDK is employed to evaluate the displacement of the depicted points “A” (a) and “B” (b). Adopted notation for each structural bay between two stiffeners: $\circ$ = LE, $\blacksquare$ = TE5, $\square$ = TE1

Fig. 28

Transverse displacement percentage error with respect to the reference solution of uniform LE refined model for bending case. $I= 1- \dfrac{{{\tilde{u}}}}{{D\tilde{O}F}}$, where ${{\tilde{u}}} = \dfrac{(u_{z}-u_{z\mathrm{Ref}})}{u_{z\mathrm{Ref}}}$ and ${D\tilde{O}F} = \dfrac{\hbox {DOF}}{\hbox {DOF}_{\mathrm{Ref}}}$

3.4 Asymmetric reinforced C-section beam subjected to transverse loading

As a final example, we want to investigate the influence of NDK models on a stiffened structure. For this reason, the same thin-walled structure present in the previous example was reinforced with 3 and 6 stiffeners. The Lagrange discretization of the stiffeners is reported in Fig. 22. The structures are reported in Fig. 23, where the stiffeners are highlighted in red. The geometric, material, boundary and loading conditions are the same as those presented in the analysis 3.2, but the discretization along the y-axis was made by 18B4, to have the same elements and length for each bay. Each stiffener is 5-mm thick and modeled with a linear interpolation FE. The reference solution was set by performing the geometrically nonlinear analysis with the cross-sectional discretization shown in Fig. 14. The results are shown in Fig. 24, along with some of the deformed shapes. As shown, the largest deformations occur near the clamp zone, so this part was kept approximated with the higher-order “LE” theory in the following NDK analyses. Particular attention must be paid to the stiffeners and their modeling. As depicted in Fig. 25a, b, in the moderate displacement field (around $P = 30$ kN), they do not show remarkable deformation, and a low-order kinematics could be enough to describe their behavior. However, in the large displacement field ($P = 70$ kN), their deformation becomes not negligible, and higher-order kinematics must be included (see Fig. 25c, d). However, to follow the main goal of this research activity, this aspect is not further analyzed, and the modeling of the stiffeners varies according to the NDK model of the closest bay.

The results of the NDK on the 3-stiffeners beam are shown in Fig. 26. Each bay of the structure has its own kinematics, described in the Figure by a $\circ$ for the LE, $\blacksquare$ for TE5, and $\square$ for TE1. By looking at Fig. 26a, the first three NDK sets lead to the same static results, with a decrease in terms of DOF equal to 33% and 25%. The fourth solution, i.e., LE–TE1–TE1, provides quite accurate results in the whole equilibrium curve but in the zone between P = 20 kN and $P = 40$ kN. This is more evident by looking at Fig. 26b, where the results in this load range are completely wrong. Finally, the last configuration, i.e., TE5-TE1-TE1, is reported to show what happens if a lower-order theory is employed for the first bay of the structure. It can be concluded that in the third bay of the structure a low-order theory can be employed, while in the first bay a higher-order theory is mandatory. In the middle bay, according to the desired output, one can use a low- or a higher-order. The same analyses were performed for the 6-stiffeners-beam, and the results are reported in Fig. 27. The second and third NDK approximations provide quite accurate results, with a DOF gain equal to 37% and 49%, respectively. From the fourth approximations, it was progressively increased the number of TE1, reaching the full TE1 in the sixth configuration. The DOF gain becomes larger, going from 65 to 80% and 96%, and the difference too, even if we can conclude that the fourth configuration leads to acceptable results.

4 Conclusions

The present research work was addressed to determine the effects and benefits of adopting node-dependent kinematics (NDK) in the geometrically nonlinear analysis of thin-walled structures. Both symmetric and asymmetric thin-walled structures were analyzed, including traverse and compressive loadings. The interesting case of thin-walled structures with transverse ribs was analyzed as well. A reference solution with a uniform refined kinematics was set for each case, and the capability of NDK models was evaluated. The results show the advantage of adopting NDK models in the proposed cases. The following main conclusions can be summarized:

• NDK is a powerful method to reduce computational cost in nonlinear problems, and this is confirmed clearly by the various sample problems reported in this paper;

• Higher/lower order kinematics are introduced easily in the regions of the structures which show higher/lower sectional deformations.

• NDK works well for geometrically nonlinear problems, and no drawbacks or numerical issues, with respect to linear analysis, were found.

• The mixing of various kinematics between two adjacent element is done simply by modifying the stiffness matrices, no mixing techniques or Lagrange multipliers are required;

For the case of the reinforced beam structure with 0, 3, and 6 stringers, Fig. 28 shows the trends of percentage error ${{\tilde{u}}}$ and I, which are defined in the following:

$$\begin{aligned} \begin{array}{ccccccccccccccccccccccc} {{\tilde{u}}} = \dfrac{u_{z}-u_{z\mathrm{Ref}}}{u_{z\mathrm{Ref}}}, \\ I= 1- \dfrac{{{\tilde{u}}}}{{D\tilde{O}F}}, \quad {D\tilde{O}F} = \dfrac{DOF}{DOF_{Ref}} \end{array} \end{aligned}$$

(17)

where the solution of a uniform, refined mathematical model is set as reference. Two different equilibrium points are considered. Clearly, the error distribution is not monotonous. It means that with a heavier model one can obtain a larger error than with a lighter one. In other words, the distribution and optimization of DOFs over the structure (namely, the NDK model), is a key step in the modeling process.