Shell element for constrained finite element analysis of thin-walled structural members

Thin-Walled Structures 105 (2016) 135–146

Contents lists available at ScienceDirect

Thin-Walled Structures journal homepage: www.elsevier.com/locate/tws

Full length article

Shell element for constrained finite element analysis of thin-walled structural members Sándor Ádány Budapest University of Technology and Economics, Department of Structural Mechanics, 1111 Budapest, Műegyetem rkp. 3, Hungary

art ic l e i nf o

a b s t r a c t

Article history: Received 7 January 2016 Received in revised form 6 April 2016 Accepted 10 April 2016

In this paper a novel shell finite element is introduced, specifically proposed for constrained shell finite element analysis. The proposed element is derived from the finite strips used in the semi-analytical finite strip method. The new finite element shares the most fundamental feature of the finite strips, namely: transverse and longitudinal directions are distinguished. Moreover, the new element keeps the transverse interpolation functions of finite strips, however, the longitudinal interpolation functions are changed from trigonometric functions (or function series) to classic polynomials. It is found that the proper selection of the polynomial longitudinal interpolation functions makes it possible to perform modal decomposition similarly as in the constrained finite strip method (cFSM). This requires an unusual combination of otherwise well-known shape functions. If the so-constructed shell finite elements are used to model a thin-walled member, (hence, with using discretization in both the transverse and the longitudinal directions,) modal decomposition can be done essentially identically as in cFSM, whilst the practical applicability of the method is significantly extended (e.g., various restraints, holes, certain crosssection changes can easily be handled). In this paper the focus is on the derivation of the novel shell finite element. Constraining capability is illustrated by some basic examples. Practical application of the novel element will be presented in subsequent papers. & 2016 Elsevier Ltd. All rights reserved.

Keywords: Constrained finite element method Shell finite element

1. Introduction The finite strip method (FSM) can be regarded as a special version of finite element method (FEM) in which special “finite element”-s are used. The most essential feature of FSM is that there are two pre-defined directions, and the base functions (or: interpolation functions) are different in the two directions. (Typically, though not necessarily and not always, the two characteristic directions are perpendicular to each other). In classical (sometimes also referred as to semi-analytical) FSM, as in [1–4] the structural member to be analyzed is discretized only in one direction (say: transverse direction), while in the other direction (say: longitudinal direction) there is no discretization, i.e., in this direction there is only one element along the member. This is why the dimensions of this “finite element” are typically distinctly different in the two directions, and that is why such an element is called “finite strip” rather than “finite element”. This special finite strip discretization has various consequences. An advantageous consequence is that the total number of elements, therefore the total number of degrees of freedom is much smaller than in case of a classic FEM, which means faster E-mail address: [email protected] http://dx.doi.org/10.1016/j.tws.2016.04.012 0263-8231/& 2016 Elsevier Ltd. All rights reserved.

calculation, as well as simpler pre- and post-processing. The price of the calculation efficiency is restricted generality: the analyzed member must be highly regular (e.g., typically it must be straight, prismatic, etc.). A major restriction of the classical semi-analytical FSM is that a certain longitudinal shape function can properly be applied for only a certain problem with given boundary conditions, since in the lack of longitudinal discretization accurate solution can be expected only if the longitudinal interpolation function well represent the real behavior (i.e., if the applied shape function satisfies the differential equation and boundary conditions of the problem).This restriction can (partially) be released if either trigonometric series or splines are used for the longitudinal interpolation. In either case the problem size is significantly increased (compared to FSM), while practical applicability is still limited (compared to FEM). The original idea of constraining a shell-model is proposed in [5,6] then in [7–10]. The idea is to define special constraints, based on some pre-defined mechanical criteria, the introduction of which enables modal decomposition. Modal decomposition transforms the original displacement field into a set of modal displacements that can solve two basic problems: calculation in a reduced but practically meaningful space (e.g., calculating global buckling directly, by using only a few degrees of freedom), and modal identification (e.g., assigning participation percentages from

136

S. Ádány / Thin-Walled Structures 105 (2016) 135–146

the modal deformations to a general deformation field). The constrained finite strip method (cFSM) was first proposed and developed for the semi-analytical FSM with sine-cosine longitudinal shape functions that correspond to locally and globally pinned-pinned end restraints of the thin-walled beam or column. Later other end conditions have also been considered [11,12], but still within the semi-analytical FSM. The method was then generalized to be able to handle general cross-sections [13,14], which also required a more systematic definition of the deformation modes [15]. An attempt to constraining a spline FSM is presented in Ref. [16]. In Refs. [17,18] the constraining technique is applied to shell FEM in the context of a commercial finite element code. Though all these above-mentioned constrained methods share the same basic mechanical background, there are distinct differences. The advantageous features of cFSM are as follows: (a) the method provides a full modal decomposition (i.e., the whole displacement filed is transformed into a modal system), (b) the mechanical criteria of the modes are exactly satisfied, and (c) adding constraints reduces the DOF number of the problem. However, since cFSM is based on FSM, it has all the restrictions of FSM: for example the member has to be regular, or, FSM is efficient only if the longitudinal shape function is defined specific to the end restraints, which means that arbitrary boundary conditions cannot be handled in an efficient way. The so-far proposed constrained FEM has potentially more general applicability than cFSM, but (a) it provides only partial modal decomposition (that is why modal identification is not readily be handled), (b) the mechanical criteria are satisfied only approximately, and (c) adding constraints increases the DOF number of the problem. It can be concluded that all of the existing constrained methods possess limitations and/or disadvantages. A better constrained method should provide full decomposition, should satisfy the mechanical criteria exactly, should be generally applicable as much as possible (as far as loading, boundary conditions, etc. are concerned). In this paper the first step toward such a method is presented. A novel shell finite element is introduced. The application of the novel shell element leads to a method which shares (most of) the advantageous features of the original cFSM, while providing significantly extended general practical applicability. The proposed element has a major geometrical limitation of being rectangular and that its local coordinate system must match the longitudinal and transverse direction of the member. Otherwise, no further restrictions are included, that is, the proposed element can be used similarly to any other shell finite elements. By using the proposed novel shell element various problems can readily be handled: various end and intermediate restraints, nearly arbitrary loading, arbitrary buckling modes including shear buckling or web crippling, certain cross-section changes along the

member length, holes. Some of these problems can be solved by other methods, some not. For example, various loading and restraints can be handled by the generalized beam theory (GBT) together with modal decomposition, see e.g. [19–22]. Though attempts to extend GBT for sections with holes are recently reported [23,24], it is still believed that the cFEM based on the here proposed finite element is powerful, since it integrates the advantageous features of all the available modal decomposition methods, and it is based on the shell finite element method which is widely used in research and even in design practice. In this paper the derivation of the novel shell finite element is presented in detail. Then the constraining capability of the element is illustrated by some basic semi-analytical examples. The examples justify the applicability. The details of the constrained finite element method (cFEM) will be presented in subsequent papers together with various examples to illustrate the advantageous features of the new method.

2. Derivation of the proposed finite element 2.1. Short overview on existing semi-analytical FSM and cFSM In finite strip method a member is discretized into longitudinal strips, instead of finite element method, which applies discretization in both the longitudinal and transverse directions. In Fig. 1a single strip is shown, along with the typically used local coordinate system and the degrees of freedom (DOF) for the strip, the dimensions of the strip, and the applied end tractions. By using the simplest longitudinal trigonometric functions, the displacements are approximated as follows.

⎡⎛ x ⎞ ⎛ x ⎞⎤ ⎡ u1⎤ mπy u (x, y) = ⎢ ⎜ 1 − ⎟ ⎜ ⎟⎥ ⎢ ⎥ sin ⎣⎝ b ⎠ ⎝ b ⎠⎦ ⎣ u2 ⎦ a

(1)

⎡⎛ x ⎞ ⎛ x ⎞⎤ ⎡ v1⎤ mπy v (x, y) = ⎢ ⎜ 1 − ⎟ ⎜ ⎟⎥ ⎢ ⎥ cos ⎣⎝ b ⎠ ⎝ b ⎠⎦ ⎣ v2 ⎦ a

(2)

⎡⎛ ⎤ 3x 2 2x 3 ⎞ ⎛ 2x 2 x 3 ⎞ ⎛ 3x 2 2x 3 ⎞ ⎛ x 2 x 3 ⎞ ⎥ − 2⎟ ⎜ 2 − 3 ⎟ ⎜ − 2⎟ w (x, y ) = ⎢ ⎜ 1 − 2 + 3 ⎟ ⎜ −x+ ⎢⎣ ⎝ b b b ⎠ ⎝ b ⎠ ⎝ b b ⎠ ⎝ b b ⎠ ⎥⎦ ⎡ w1 ⎤ ⎢Θ ⎥ ⎢ 1 ⎥ sin mπy ⎢ w2⎥ a ⎢⎣ Θ 2 ⎥⎦

(3)

The above formulae represent pinned-pinned boundary conditions. Two important features of the longitudinal shape functions are that the same longitudinal function is used for u and w, and the longitudinal function for v is the derivative of that used for

Fig. 1. Coordinates and DOF in finite strip method.

S. Ádány / Thin-Walled Structures 105 (2016) 135–146

137

u and w. Therefore, the above formulae can be generalized as follows:

base vectors for the given M space. Applying RM for the intended space (M ¼G, D, L, S, and/or T) Eq. (7) becomes:

⎡⎛ x ⎞ ⎛ x ⎞⎤ ⎡ u1⎤ u (x, y) = ⎢ ⎜ 1 − ⎟ ⎜ ⎟⎥ ⎢ ⎥ f (y) ⎣⎝ b ⎠ ⎝ b ⎠⎦ ⎣ u2 ⎦

R M TK e R M Φ M − Λ M R M TK g R M Φ M = 0

(10)

(4)

⎡⎛ x ⎞ ⎛ x ⎞⎤ ⎡ v1⎤ f ′(y) v (x, y) = ⎢ ⎜ 1 − ⎟ ⎜ ⎟⎥ ⎢ ⎥ ⎣⎝ b ⎠ ⎝ b ⎠⎦ ⎣ v2 ⎦ g (y)

K eM Φ M − ΛM K gM Φ M = 0

(11)

(5)

⎡⎛ 3x2 2x3 ⎞ ⎛ 2x2 x3 ⎞ ⎛ 3x2 2x3 ⎞ ⎛ x2 x3 ⎞⎤ w (x, y) = ⎢ ⎜ 1 − 2 + 3 ⎟ ⎜ −x+ − ⎟ ⎜ − 3 ⎟ ⎜ − 2 ⎟⎥ ⎢⎣ ⎝ b b ⎠ ⎝ b b2 ⎠ ⎝ b2 b ⎠ ⎝ b b ⎠⎥⎦ ⎡ w1⎤ ⎢Θ ⎥ ⎢ 1 ⎥ f (y) ⎢ w2 ⎥ ⎢⎣ Θ2 ⎥⎦

(6)

By appropriately selecting the functions, various end conditions can be described, as in [2,11,12]. The proposed functions are all trigonometric functions. Depending on the boundary condition, the functions might be given in series form. Note, in Eq. (5) f(y) is the longitudinal base function used for the transverse in-plane and out-of-plane displacements, while g(y) is usually defined as the first derivative of the internal function of f(y), see e.g., Eq. (2). The local elastic and geometric stiffness matrices can be constructed by following conventional FEM steps, by considering the 2D generalized Hooke's law (for the elastic stiffness matrix) and by considering the second-order strain terms (for the geometric stiffness matrix). The stiffness matrices can be determined analytically. From the local stiffness matrices the member's (global) stiffness matrices (elastic and geometric, Ke and Kg) can be compiled as in FEM, by transformation to global coordinates and assembly. For a given distribution of edge tractions on a member the geometric stiffness matrix scales linearly, resulting in the classic eigen-buckling problem, namely:

K eΦ − ΛK g Φ = 0

(7)

with

Λ = diag ⎡⎣ λ1 λ2 ... λnDOF ⎤⎦ and Φ = ⎡⎣ φ1 φ2 ... φnDOF ⎤⎦

(8)

where λi is the critical load multiplier and φi is the mode shape vector. The constrained FSM (cFSM) is an extension to FSM that uses mechanical assumptions to enforce or classify deformations to be consistent with a desired set of criteria. The method is originally presented in Refs. [5–10], and implemented in Refs. [25,26]. The cFSM constraints are mechanically defined, and are utilized to formally categorize deformations into global (G), distortional (D), local (L), and other (i.e., shear and transverse extension, S þT) deformations. The mechanical criteria are mostly given by setting certain displacement and displacement derivatives to zero. Once the mechanical criteria are transformed into constraint matrices, any FSM displacement field d (e.g. an eigen-buckling mode φ is an important special case) may be constrained to any modal dM deformation space via:

d = R M dM

where KeM and KgM are reduced-size elastic and geometric stiffness matrices for the eigen-buckling solution constrained to space M. Modal identification, i.e. categorization of a general deformation into the M spaces, is also possible, due to the fact that G þD þL þSþT spans the entire FSM space. As such, the RGDLST constraint matrix represents an alternative basis for the FSM space, in which deformations are categorized. This basis transformation of displacement vector d may be expressed as:

d = ⎡⎣ R G R D R L R ST ⎤⎦ c

where c now provides the deformations within each class: cG, cD, cL, cST. The values of c are dependent on the normalization of the base vectors within R. A full discussion of the normalization selection for R is provided in Ref. [27]. Once c is determined, pi participation of an individual mode or pM participation of a deformation space M can also be determined as follows:

pi = ‖ci ‖/‖c‖ or pM = ‖cM ‖/‖c‖

(13)

The cFSM method has recently been generalized [13,14]. An important feature of the generalization is the more refined definition of the modes, (e.g., introduction of primary and secondary mode, as well a practically meaningful handling of the shear modes), and the extension of the method to arbitrary crosssections. 2.2. New longitudinal shape functions Our goal here is to transform the “finite strip” into a shell “finite element”. Since the above-summarized semi-analytical FSM uses classic polynomials in the transverse direction, the new shell element can inherit the transverse interpolation functions from FSM. The longitudinal interpolation function should be changed, however, by keeping some important characteristics of the functions of FSM. Namely, the new longitudinal shape functions must have the following features:

 they must be able to exactly satisfy the constraining criteria for 

mode decomposition (no-shear criterion, no-transverse-extension criterion, etc.), the transverse in-plane displacements must be interpolated by using the same shape functions as used for the out-of-plane displacements,

(9)

where RM is a constraint matrix, the derivation of which can be in found Refs. [5–10] for open cross-sections, and M might be G, D, L, S and/or T. Though modal decomposition is not restricted to eigen-buckling solution, this is the problem where modal decomposition is mostly used. It can be completed by introducing the desired constraint matrix RM, the columns of which can be interpreted as

(12)

Fig. 2. Coordinates and DOF in finite strip method.

138

S. Ádány / Thin-Walled Structures 105 (2016) 135–146

 they must provide C(1) continuous interpolation for the out-ofplane displacements (which is practically useful for defining various end restraints). Moreover, we change the local coordinate system, see Fig. 2, in order to be more conform with the common practice:

(14)

(20)

(21)

(22)

(3) Nx,2 =x−

3x2 2x3 + 3 a2 a

2x2 x3 + a a2

(23)

(24)

3x2 2x3 − 3 a2 a

(25)

x2 x3 + a a2

(26)

First-order shape functions in y: (1) N y,1 =1

(1) N y,2 =

−

y b

y b

(27)

(28)

Third-order shape functions in y: (3) N y,1 =1 −

(15) ⎛ 2y 2 y 3 ⎞ ⎟⎟ ⎜⎜ y − + b ⎝ b2 ⎠

⎛ 3y 2 2y 3 ⎞ ⎟⎟ ⎜⎜ − ⎝ b2 b3 ⎠

⎛ 2x 2 x 3 ⎞ ⎜ −x+ − ⎟ a ⎝ a2 ⎠

(3) N y,2 =y −

⎛ y2 y3 ⎞ ⎤ ⎟⎟ ⎥ ⎜⎜ − + ⎝ b b 2 ⎠ ⎥⎦

⎛ 3x 2 2x 3 ⎞ ⎟ ⎜ − ⎝ a2 a3 ⎠

⎛ x2 x 3 ⎞ ⎜ − ⎟ ⎝ a a2 ⎠

(3) N y,3 =

⎤ ⎥ ⎥⎦

(16)

⎡ cu1⎤ ⎡ u1⎤ u (x, y) = ⎡⎣ N y(1,1) N y(1,2) ⎤⎦ ⎢ ⎥× ⎡⎣ Nx(2,1) Nx(2,2) Nx(2,3) ⎤⎦ ⎢ cu2 ⎥ ⎢ ⎥ ⎣ u2 ⎦ ⎣ cu3 ⎦

(17)

⎡ cv1⎤ ⎢ ⎥ ⎡ (1) (1) ⎤ ⎡ v1⎤ ⎡ (3) (3) (3) (3) ⎤ ⎢ cv2 ⎥ v (x, y) = ⎣ N y,1 N y,2 ⎦ ⎢ ⎥× ⎣ Nx,1 Nx,2 Nx,3 Nx,4 ⎦ ⎣ v2 ⎦ ⎢ cv3 ⎥ ⎣ cv4 ⎦

(18)

with the shape functions as follows. Second-order shape functions in x:

2y 2 y 3 + 2 b b

(29)

(30)

(31)

y2 y 3 + b b2

(32)

The proposed longitudinal shape functions and its coefficients are illustrated in Fig. 3. 2.3. Finite element interpretation The above formulae include separate sets of coefficients for the transverse and longitudinal directions. However, these coefficients can easily be exchanged by classic finite element nodal displacements. As an example, the in-plane longitudinal displacement is expressed as follows, see Eq. (17):

u (x, y)=u1cu1N y(1,1) Nx(2,1)+u2 cu1N y(1,2) Nx(2,1)+u1cu2 N y(1,1) Nx(2,2)+u2 cu2 N y(1,2) Nx(2,2)+u1cu3 N y(1,1) Nx(2,3)+u2 cu3 N y(1,2) Nx(2,3)

⎡ w1 ⎤ ⎥ ⎢ (3 ) ⎤ ⎢ Θ 1 ⎥ N y,4 ⎦ ⎢ w2⎥ ⎢⎣ Θ 2 ⎥⎦

⎡ cw1 ⎤ ⎥ ⎢ ⎡ (3 ) (3 ) (3 ) (3) ⎤ ⎢ cw 2 ⎥ × ⎣ N x,1 − N x,2 N x,3 − N x,4 ⎦ ⎢ cw 3 ⎥ ⎢⎣ cw 4 ⎥⎦

3y 2 2y 3 + 3 b2 b

3y 2 2y 3 − 3 b2 b

(3) N y,4 =−

In a shorter form:

w (x, y ) = ⎡⎣ N (y3,1) N (y3,2) N (y3,3)

x 2x2 + a a2

(3) Nx,4 =−

⎡ cv1⎤ ⎢c ⎥ ⎢ v2 ⎥ ⎢ cv3 ⎥ ⎣ cv4 ⎦

⎡ ⎛ 3x 2 2x 3 ⎞ ⎟ ×⎢ ⎜ 1 − + ⎢⎣ ⎝ a2 a3 ⎠ ⎡ cw1 ⎤ ⎢c ⎥ ⎢ w2 ⎥ ⎢ cw 3 ⎥ ⎢⎣ cw 4 ⎥⎦

4x 4x2 − 2 a a

(3) Nx,1 =1 −

(3) Nx,3 =

⎛ x 2x2 ⎞ ⎤ ⎜− + 2 ⎟ ⎥ ⎝ a a ⎠ ⎦

⎡⎛ y ⎞ ⎛ y ⎞⎤ ⎡ v1⎤ v (x, y) = ⎢ ⎜ 1 − ⎟ ⎜ ⎟⎥ ⎢ ⎥ ⎝ ⎣ b ⎠ ⎝ b ⎠⎦ ⎣ v2 ⎦ ⎡⎛ 3x2 2x3 ⎞ ⎛ 2x2 x3 ⎞ ⎛ 3x2 2x3 ⎞ ⎛ x2 x3 ⎞⎤ × ⎢ ⎜ 1 − 2 + 3 ⎟ ⎜ x− + ⎟ ⎜ − 3 ⎟ ⎜ − + 2 ⎟⎥ ⎢⎣ ⎝ a a2 ⎠ ⎝ a2 a ⎠ ⎝ a a ⎠⎥⎦ a a ⎠ ⎝

⎡⎛ 3y 2 2y 3 ⎞ ⎟⎟ w (x , y ) = ⎢ ⎜⎜ 1 − + ⎢⎣ ⎝ b2 b3 ⎠

3x 2x2 + 2 a a

Third-order shape functions in x:

The notations for the displacement functions and nodal displacements will be changed accordingly, as detailed below. The proposed interpolation functions are summarized as follows:

⎡ cu1⎤ ⎢c ⎥ ⎢ u2 ⎥ ⎣ cu3 ⎦

(2) Nx,2 =

(2) Nx,3 =−

 The longitudinal axis will be the global X and local x axis,  the coordinate system will be right-handed.

⎡⎛ y ⎞ ⎛ y ⎞⎤ ⎡ u1⎤ u (x, y) = ⎢ ⎜ 1 − ⎟ ⎜ ⎟⎥ ⎢ ⎥ ⎣⎝ b ⎠ ⎝ b ⎠⎦ ⎣ u2 ⎦ ⎡ ⎛ ⎛ 4x 4x2 ⎞ 3x 2x2 ⎞ × ⎢ ⎜1 − + 2 ⎟ − 2 ⎟ ⎜ ⎝ a a a ⎠ a ⎠ ⎣ ⎝

(2) Nx,1 =1 −

(33)

The in-plane longitudinal DOF are the ui cuj constants, where i¼1...2 and j¼1...3. Thus, finally, there are 6 such DOF, all of them are translational, and will be denoted here as in Fig. 4. The interpolation with the finite element DOF:

(19)

u (x, y)=u11N y(1,1) Nx(2,1)+u13 N y(1,2) Nx(2,1)+u21N y(1,1) Nx(2,2)+u23 N y(1,2) Nx(2,2)+u31N y(1,1) Nx(2,3)+u33 N y(1,2) Nx(2,3)

(34)

Similarly, the in-plane transverse DOF are the vi cvj constants,

S. Ádány / Thin-Walled Structures 105 (2016) 135–146

see Eq. (18), where i ¼1...2 and j¼1...4. Thus, finally, there are 8 such DOF, which will be denoted here as in Fig. 4. Therefore, the v displacement function is interpolated as follows:

v (x, y)=v11N y(1,1) Nx(3,1)+v13 N y(1,2) Nx(3,1)+v31N y(1,1) Nx(3,3)+v33 N y(1,2) Nx(3,3)+ϑz11N y(1,1) Nx(3,2)+ϑz13N y(1,2) Nx(3,2)+ϑz 31N y(1,1) Nx(3,4) +ϑz 33N y(1,2) Nx(3,4)

(35)

The out-of-plane displacement function can be expressed similarly from Eq. (19), by using finite element nodal displacement DOF, as follows:

w (x,

dNx(2,1) dNx(2,1) dNx(2,2) dNx(2,2) ∂u =u11N y(1,1) +u13 N y(1,2) +u21N y(1,1) +u23 N y(1,2) +u31 ∂x dx dx dx dx (2) (2) dNx,3 dNx,3 +u33 N y(1,2) N y(1,1) dx dx dNx(2,1) dx

+

dNx(2,3) dx

dNx(2,2) dx

(u21N y(1,1)+u23 N y(1,2) )

(u31N y(1,1)+u33 N y(1,2) ) = 0

(39)

Considering the shape functions and its derivatives, it is easy to conclude that the actual strain function is linear both in x and y. Therefore, the function can be expressed in the form:

Nx(3,4) −ϑxy11N y(3,2) Nx(3,2)−ϑxy13N y(3,4) Nx(3,2)−ϑxy31N y(3,2)

C11xy +C10 x+C01y +C00=0 (36)

Therefore, the proposed element has 30 DOF: 6 for u, 8 for v. and 16 for w. Each corner node has 7 DOF (1 for u, 2 for v, and 4 for w), while there are two additional nodes at (x,y)¼(a/2,0) and (x, y) ¼(a/2,b) with one DOF per node for the u displacement. The DOF are illustrated in Fig. 4.

3. Constraints The constraints that are embedded in cFSM are given in [13,14]. It can be observed that the constraints are formulated by setting various displacement derivatives to zero. It can also be observed that the criteria are practically independent of the longitudinal shape functions, therefore, the same criteria can be used for the here-proposed finite element that are used for finite strips. It is also important that the introduction of the mechanical criteria must lead to simple relationships in between the nodal displacement DOF, since this is a necessary condition if the mechanical criteria are intended to be exactly satisfied. In case of the no-longitudinal-extension criterion, the criterion is:

∂u =0 ∂x

(u11N y(1,1)+u13 N y(1,2) ) +

(38)

Nx(3,3)−ϑy11N y(3,1) Nx(3,2)−ϑy13N y(3,3) Nx(3,2)−ϑy31N y(3,1) Nx(3,4) −ϑy33N y(3,3)

Nx(3,4) −ϑxy33N y(3,4) Nx(3,4)

εx=

notations as in Fig. 4):

y)=w11N y(3,1) Nx(3,1)+w13 N y(3,3) Nx(3,1)+w31N y(3,1) Nx(3,3)+w33 N y(3,3) Nx(3,3)+ϑx11N y(3,2) Nx(3,1)+ϑx13N y(3,4) Nx(3,1)+ϑx31N y(3,2) Nx(3,3)+ϑx33N y(3,4)

139

with the C coefficients as follows:

C11= − 4 (u11−u13−2u21+2u23+u31−u33 )/b/a2

(41)

C10=4 (u11−2u21+u31)/a2

(42)

C01=(3u11−3u13−4u21+4u23+u31−u33 )/b/a

(43)

C00= − (3u11−4u21+u31)/a

(44)

The longitudinal strain is zero for any x-y if (and only if) all the C coefficients are zero. This is satisfied if:

u11 = u21 = u31

(45)

u13=u23=u33

(46)

This means that rigid-body longitudinal displacement is allowed, otherwise the longitudinal displacement must be zero. The implementation of the various criteria requires the repetition of essentially the same steps, therefore, the other criteria will not be discussed in detail. Nevertheless, the resulting relationships are summarized as follows. In case of the no-transverse-extension criterion, the criterion is:

(37) εy=

Using the assumed shape functions (with using the general

(40)

∂v =0 ∂y

Fig. 3. Longitudinal shape functions.

(47)

140

S. Ádány / Thin-Walled Structures 105 (2016) 135–146

Fig. 4. Nodal displacement DOF of the finite element.

Using the assumed shape functions, the criterion is satisfied for any x-y if:

v11=v13

ϑx13 − ϑx33 =ϑxy13 (or=ϑxy33) a

(48)

(64)

In case of the no- transverse-curvature criterion, the criterion is:

v31=v33

(49) κy=

ϑz11=ϑz13

(50)

ϑz 31=ϑz 33

(51)

In case of the no-shear criterion, the criterion is:

γxy=

∂u ∂v + =0 ∂y ∂x

(52)

Since in the constraining procedure the no-shear criterion is applied always together with the no-transverse extension criterion, it is more straightforward (and leads to simpler conclusions) to immediately consider the results of no-transverse-extension criterion, which finally leads to the following equations:

u11 − u13 =ϑz11(or = ϑz13) b u21 − u23 3 v31 − v11 ϑz11 + ϑz 31 = − b 2 a 4 u31 − u33 =ϑz 31(or = ϑz 33) b

(53)

(54)

is:

∂ 2w =0 ∂x2

Using the assumed shape functions, the criterion is satisfied for any x-y if:

ϑx11=ϑx13

(66)

w13 − w11 =ϑx11(or=ϑx13) b

(67)

ϑx31=ϑx33

(68)

w33 − w31 =ϑx31(or=ϑx33) b

(69)

ϑxy11=ϑxy13

(70)

ϑy13 − ϑy11 =ϑxy11(or=ϑxy13) b

(71)

ϑxy31=ϑxy33

(72)

ϑy33 − ϑy31 =ϑxy31(or=ϑxy33) b

(73)

In case of the no- mixed-curvature criterion, the criterion is:

(56)

Using the assumed shape functions, the criterion is satisfied for any x-y if:

ϑy11=ϑy31

(65)

(55)

In case of the no- longitudinal-curvature criterion, the criterion

κ x=

∂ 2w =0 ∂y2

(57)

κxy=

∂ 2w =0 ∂x∂y

(74)

Using the assumed shape functions, the criterion is satisfied for any x-y if:

ϑxy11=ϑxy13=ϑxy31=ϑxy33=0

(75)

ϑx11=ϑx31

(76)

ϑx13=ϑx33

(77)

(60)

ϑy11=ϑy13

(78)

ϑxy11=ϑxy31

(61)

ϑy31=ϑy33

(79)

ϑx11 − ϑx31 =ϑxy11(or=ϑxy31) a

w31−w11=w33−w13

(80)

(62)

ϑxy13=ϑxy33

(63)

w11 − w31 =ϑy11(or=ϑy31) a

(58)

ϑy13=ϑy33

(59)

w13 − w33 =ϑy13 (or=ϑy33) a

In case of no-warping-shear criterion the shear strain is allowed to be non-zero, but only due to transverse displacements. The criterion therefore is:

S. Ádány / Thin-Walled Structures 105 (2016) 135–146

∂u =0 ∂y

(81)

Using the assumed shape functions, the criterion is satisfied for any x-y if:

u11=u13

(82)

u21=u23

(83)

u31=u33

(84)

In case of no-transverse-shear criterion the shear strain is allowed to be non-zero, but only due to longitudinal displacements. The criterion therefore is:

∂v =0 ∂x

(85)

Using the assumed shape functions, the criterion is satisfied for any x-y if:

ϑz11=ϑz13=ϑz 31=ϑz 33

(86)

v11=v31

(87)

v13=v33

(88)

141

above proposed finite elements and constraining equations, but the problems are solved analytically in order to provide a better comparison to exact analytical solutions. The problem is a simple straight column or beam member with a (narrow) rectangular cross-section with width and thickness b and t, respectively. The length is L. It is discretized by two (equal) finite elements along the length, as shown in Fig. 5. The applied FE model thus has 10 nodes, 6 nodes having 7 DOF and 4 nodes having 1 DOF. Altogether the model has 46 DOF, which is reduced due to constraining and end restraints, as discussed as follows. Since our goal is to calculate critical forces to global modes, the constraints characterizing global modes must be applied. Namely: (i) transverse (membrane) extension is zero, (ii) (membrane) shear deformations are zero, and (iii) there is no transverse curvature. No-transverse-extension criteria means that the v and ϑz , DOF are halved, thus, instead of having 6 different v and 6 different ϑz , there are only 3 independent v and 3 independent ϑz . Furthermore, no-shear criterion means that the u (warping) displacements are dependent on v-s and ϑz -s. Therefore, instead of having 10 u-s, we have only 5 independent u DOF. Finally, no-transverse-curvature criterion means:

   

instead instead instead instead

of of of of

having having having having

6, 6, 6, 6,

there there there there

are are are are

only only only only

3 3 3 3

independent independent independent independent

ϑx ϑxy , w ϑy

4. Demonstrative examples 4.1. Model, restraints, constraints A few examples are presented here. The examples are simple classical buckling problems to which analytical solution (for the critical load) is known. The examples are solved by using the

There are 23 independent DOF, shown in Fig. 6, The constraining equations are given in Table 1: Two classic end restraint conditions will be considered here: S-S, i.e. simple-simple (or: pinned-pinned) and C-C, i.e., clampedclamped. In case of S-S, “fork” supports are assumed at both ends, i.e.,

Fig. 5. Beam or column member and its discretization.

142

S. Ádány / Thin-Walled Structures 105 (2016) 135–146

u (x, y)=u1N y(1,1) Nx(2,1)+u2 N y(1,2) Nx(2,1)+u3 N y(1,1) Nx(2,2)+u4 N y(1,2) Nx(2,2)+u5 N y(1,1) Nx(2,3)+u6 N y(1,2) Nx(2,3)

(89)

Substituting the constraint equations (for the notations, see Fig. 6):

u (x, y)=(ua +ϑz, a b/2) N y(1,1) Nx(2,1)+(ua −ϑz, a b/2) ⎛ ⎛ 1. 5va ϑz, a 1. 5vc ϑz, c ⎞ b ⎞ ⎟ ⎟ + − + N y(1,2) Nx(2,1)+⎜ ub −⎜ ⎝ a ⎝ a 4 4 ⎠ 2⎠ ⎛ ⎛ 1. 5va ϑz, a 1. 5vc ϑz, c ⎞ b ⎞ ⎟ ⎟ + − + N y(1,1) Nx(2,2)+⎜ ub +⎜ ⎝ a ⎝ 4 a 4 ⎠ 2⎠ N y(1,2) Nx(2,2)+(uc +ϑz, c b/2) N y(1,1) Nx(2,3)+(uc −ϑz, c b/2) N y(1,2) Nx(2,3)

(90)

In case of C-C restraints, uc =0, va=0 and ϑz, a=0, therefore:

⎛ ⎛ 1. 5vc ϑz, c ⎞ b ⎞ ⎟ ⎟ u (x, y)=(ua ) N y(1,1) Nx(2,1)+(ua ) N y(1,2) Nx(2,1)+⎜ ub −⎜ − + ⎝ ⎝ a 4 ⎠ 2⎠ ⎛ ⎛ 1. 5vc ϑz, c ⎞ b ⎞ (1) (2) ⎟ ⎟ N N +(ϑz, c b/2) N y(1,1) Nx(2,2)+⎜ ub +⎜ − + ⎝ ⎝ 4 ⎠ 2 ⎠ y,2 x,2 a Fig. 6. Nodal DOF of the constrained column or beam member.

va=0, ve=0, wa=0, we=0, ϑx, a=0 and ϑx, e=0. In addition, one crosssection in the longitudinal direction should be restrained, e.g.,: uc =0. Thus, in case of S-S end restraints there are 16 independent DOF. In case of C-C, the following DOF are known to be zero: va=0, ve=0, wa=0, we=0, ϑx, a=0, ϑx, e=0, ϑy, a=0, ϑy, e=0, ϑz, a=0, ϑz, e=0, ϑxy, a=0 and ϑxy, e=0. In addition, one cross-section in the longitudinal direction should be restrained, e.g.: uc =0. Thus, in case of C-C end restraints there are 10 independent DOF. Since we are aiming analytical solution, the number of DOF must further be reduced. It is known that the buckling shapes (for the considered end restraints) are either symmetrical or pointsymmetrical. If they are symmetrical, this means further reduction of the independent DOF as follows: ϑy, c=0, ϑz, c=0, ϑxy, c=0, ue=−ua , ud=−ub , ϑy, e=−ϑy, a , ϑxy, e=−ϑxy, a and ϑz, e=−ϑz, a . If they are point-symmetrical, the following equations must be satisfied: vc =0, wc =0 and ϑx, c=0, ue=ua , ud=ub , ϑy, e=ϑy, a , ϑxy, e=ϑxy, a and ϑz, e=ϑz, a . (Note, in case of C-C, the last 3 conditions are automatically satisfied, since the involved displacement components are zero.). Thus, utilizing symmetry or point-symmetry condition, the independent (non-zero) DOF number reduces to 8 in case of S-S, while to 5 in case of C-C. With these 8 or 5 parameters the whole displacement field is determined. These parameters are listed for the considered cases, as follows:

   

S-S, symmetric: ua , ϑy, a , ϑz, a , ϑxy, a , ub , vc , wc , ϑx, c . S-S, point-symmetric: ua , ϑy, a , ϑz, a , ϑxy, a , ub , ϑy, c , ϑz, c , ϑxy, c . C-C, symmetric: ua , ub , vc , wc , ϑx, c . C-C, point-symmetric: ua , ub , ϑy, c , ϑz, c , ϑxy, c .

4.2. Solution procedure To have the buckling solution, the energy method is followed. The main steps are summarized as follows. The u(x,y), v(x,y) and w(x,y) displacement functions are given above. As an illustration, let us express v(x,y) for the first element, in case of C-C restraints and symmetrical displacements. The displacement function for one finite element is given in general by Eq. (34). By using the notations of Fig. 5, it can be written as follows:

N y(1,1) Nx(2,3)+(−ϑz, c b/2) N y(1,2) Nx(2,3)

(91)

⎛ ⎛ 1. 5vc ϑz, c ⎞ b ⎞ ⎟ ⎟ u (x, y)=ua N y(1,1) Nx(2,1)+ua N y(1,2) Nx(2,1)+⎜ ub +⎜ − ⎝ a ⎝ 4 ⎠ 2⎠ ⎛ ⎛ 1. 5vc ϑz, c ⎞ b ⎞ (1) (2) b ⎟ ⎟ N N +ϑz, c N y(1,1) N y(1,1) Nx(2,2)+⎜ ub −⎜ − ⎝ a ⎝ 4 ⎠ 2 ⎠ y,2 x,2 2 b Nx(2,3)−ϑz, c N y(1,2) Nx(2,3) 2

(92)

Finally, in case of symmetrical displacements ϑz, c=0, therefore:

⎛ 3b ⎞ (1) (2) ⎛ 3b ⎞ u (x, y)=ua N y(1,1) Nx(2,1)+ua N y(1,2) Nx(2,1)+⎜ ub +vc ⎟ ⎟ N N +⎜ ub −vc ⎝ 4a ⎠ y,1 x,2 ⎝ 4a ⎠ N y(1,2) Nx(2,2)

(93)

As it can be seen, the u(x,y) displacement function of the first finite element is expressed by three displacement parameters, namely: ua, ub and vc. The other displacement functions can be expressed similarly. Dependency on z can be considered as follows:

u (x, y , z )=u (x, y)−z

∂w (x, y) ∂x

(94)

v (x, y , z )=v (x, y)−z

∂w (x, y) ∂y

(95)

w (x, y , z )=w (x, y)

(96)

The first-order strains can be determined as follows:

⎤ ⎡ ∂u (x, y , z ) ⎥ ⎢ x ∂ ⎡ εx (x, y , z ) ⎤ ⎢ ⎥ ⎥ ⎢ ⎢ ∂v (x, y , z ) ⎥ ε=⎢ εy (x, y , z ) ⎥=⎢ ⎥ ∂y ⎢⎣ γ (x, y , z )⎥⎦ ⎢ ⎥ xy u x , y , z v x , y , z ∂ ( ) ∂ ( ) ⎥ ⎢ + ⎥⎦ ⎢⎣ ∂y ∂x

(97)

The material follows Hooke's law:

⎤ ⎡ E νE 0⎥ ⎢ 2 2 ⎥ ⎢ 1−ν 1−ν D=⎢ νE ⎥ E ⎢ 1−ν 2 1−ν 2 0 ⎥ ⎥ ⎢ ⎣ 0 G⎦ 0

(98)

S. Ádány / Thin-Walled Structures 105 (2016) 135–146

143

Table 1 Summary of constraint equations for the example. Actual constraint equations

Reference to eqs.

u2=ua −ϑz, a b/2 v2=va w2=wa+ϑx, a b/2 ϑx,2=ϑx, a

ϑ y,1=ϑ y, a−ϑxy, a b/2

ϑ y,2=ϑ y, a+ϑxy, a b/2

Eq. (66) Eq. (71)

ϑz,1=ϑz, a ϑxy,1=ϑxy, a

ϑz,2=ϑz, a ϑxy,2=ϑxy, a

Eq. (50) Eq. (70)

u3=ub −

(

1.5va ϑz , a 1.5vc ϑz , c + − + a 4 a 4

)

b 2

u4=ub +

2

t /2 L /2 b /2

i = 1 −t /2 0

0

ϑz,6=ϑz, c ϑxy,6=ϑxy, c

Eq. (50) þ(51) Eq. (70)þ(72)

(

1.5vc ϑz , c 1.5ve ϑz , e + − + a 4 a 4

)

b 2

u8=ud +

(

Eq. (48) þ (49) Eq. (67)þ(69)

1.5vc ϑz , c 1.5ve ϑz , e + − + a 4 a 4

)

b 2

Eq. (54) Eq. (55)

u9=ue +ϑz, e b/2 v9=ve w9=we−ϑx, e b/2 ϑx,9=ϑx, e

u10=ue −ϑz, e b/2 v10=ve w10=we+ϑx, e b/2 ϑx,10=ϑx, e

ϑ y,9=ϑ y, e−ϑxy, e b/2

ϑ y,10=ϑ y, e+ϑxy, e b/2

Eq. (68) Eq. (73)

ϑz,9=ϑz, e ϑxy,9=ϑxy, e

ϑz,9=ϑz, e ϑxy,10=ϑxy, e

Eq. (51) Eq. (72)

Eq. (49) Eq. (69)

which the values of critical load can be determined. 4.3. Flexural and torsional buckling of columns

1 T ε i Dε i dxdydz 2

(99)

(100)

(101)

The first formula, i.e. option a) is used in classical (beammodel-based) stability solutions, while the second formula is typically used in shell models (such as FSM or shell FEM). These options are discussed in detail in [28,29]. Here, both options will be considered. The external potential thus: b L /2

∑∫ ∫ ∫

II

εx, i px dxdydz

0

(102)

In equilibrium the total potential is stationary, thus:

∂(Πint + Π ext ) =0 ∂

Eq. (53) þ(55)

ϑz,5=ϑz, c ϑxy,5=ϑxy, c

1 ⎡ ⎛ ∂u (x, y , z ) ⎟⎞2 ⎜⎛ ∂v (x, y , z ) ⎟⎞2 ⎜⎛ ∂w (x, y , z ) ⎟⎞2⎤ εx II = ⎢ ⎜ + + ⎥ ⎠ ⎝ ⎠ ⎝ ⎠⎦ ∂x ∂x ∂x 2⎣⎝

i = 1 −t /2 0

Eq. (54)

Eq. (66) þ (68) Eq. (71) þ (73)

or

Π ext = −

b 2

ϑ y,6=ϑ y, c +ϑxy, c b/2

1 ⎡ ⎛ ∂v (x, y , z ) ⎟⎞2 ⎜⎛ ∂w (x, y , z ) ⎟⎞2⎤ εx II = ⎢ ⎜ + ⎥ ⎠ ⎝ ⎠⎦ ∂x ∂x 2⎣⎝

t /2

)

ϑ y,5=ϑ y, c −ϑxy, c b/2

The external energy is the negative of the work done by the external loading on the second-order displacements. The external loading is either uniformly distributed over the cross-section or linearly changing with y. The second-order strain can be expressed as:

2

1.5va ϑz , a 1.5vc ϑz , c + − + a 4 a 4

u6=uc −ϑz, c b/2 v6=vc w6=wc +ϑx, c b/2 ϑx,6=ϑx, c

The internal potential, that is the strain energy is expressed as follows:

∑∫ ∫ ∫

(

Eq. (48) Eq. (67)

u5=uc +ϑz, c b/2 v5=vc w5=wc −ϑx, c b/2 ϑx,5=ϑx, c

u7=ud −

Πint =

Eq. (53)

u1=ua +ϑz, a b/2 v1=va w1=wa−ϑx, a b/2 ϑx,1=ϑx, a

(103)

which leads to a homogeneous system of linear equations. Nontrivial solution exists if the coefficient matrix is stationary, from

The loading is a concentric axial force, which is assumed to distribute uniformly over the end sections.

px =

F bt

(104)

In case of S-S restraints and symmetric displacements, there are 8 independent displacements parameters. This means that Eq. (103) leads to a (homogeneous) system of linear equations consisting of 8 equations, to which maximum 8 eigen-values, therefore, maximum 8 buckling modes might belong. However, if the coefficient matrix is rank-deficient, the number of existing modes might be less than 8, too. Similarly, there are maximum 8 pointsymmetric buckling modes for the S-S problem. Since the symmetric and point-symmetric displacements are certainly linearly independent of one another, this finally means that maximum 16 buckling modes might exist. In case of C-C restraints there are 5 independent displacement parameters, which means that Eq. (103) leads to 5 equations for both the symmetric and point-symmetric cases, therefore, the maximum number of buckling modes is 10. 4.3.1. Option a) In case of S-S with option a) the equation system determined by Eq. (103) is rank-deficient, and there are only 6 solutions for the symmetrical modes and 6 solutions for the point-symmetrical modes as follows:

 2 symmetrical and 2 point-symmetrical solutions in minor-axis buckling mode

 2 symmetrical and 2 point-symmetrical solutions in major-axis buckling mode

144

S. Ádány / Thin-Walled Structures 105 (2016) 135–146

 2 symmetrical and 2 point-symmetrical solutions in pure tor-

  

sional buckling mode. Therefore, 2 symmetric and 2 point-symmetric modes are nonexistent in S-S with option a). It can be observed from the solutions, that the existing modes are independent of the longitudinal displacements. In other words, in option a) the longitudinal displacements do not lead to buckling. In case of C-C with option a) there are 6 solutions as follows: 1 symmetrical and 1 point-symmetrical solution in minor-axis buckling mode 1 symmetrical and 1 point-symmetrical solution in major-axis buckling mode 1 symmetrical and 1 point-symmetrical solution in pure torsional buckling mode.

Therefore, 2 symmetric and 2 point-symmetric modes are nonexistent in C-C with option a). Again, the non-existing modes are associated with the longitudinal displacements. The modes are illustrated in Fig. 7, where the distribution of the characteristic displacement is shown along the member length. The characteristic displacement is v in case of major-axis buckling, w in case of minor-axis buckling, and ϑx in case of torsional buckling. In case of flexural buckling the solution (i.e., the expression for the critical force) can be written as:

Fcr, F =c

EI L2

relevant c constants are summarized in Table 2, where the exact c values are also given from the well-known classic analytical solutions. Here are some comments.

 The applied discretization provides only the 4 lowest critical









(105)

where c is a constant depending on the end restraint conditions and the buckling shape, while EI is the flexural rigidity about the relevant axis. In case of pure torsional buckling the solution (i.e., the expression for the critical force) can be written as:

Fcr, T =

1 ⎛ cEI ⎞ ⎜ GIT + 2w ⎟ r02 ⎝ L ⎠

(106)

where c is a constant depending on the end restraint conditions and the buckling shape, while EIw is the warping rigidity, GIt is the (St Venant) torsion rigidity and r0 is the polar radius of gyration. Considering both the symmetric and point-symmetric cases, there are altogether 4 solutions for minor-axis and 4 solutions for major-axis flexural buckling, and 4 solutions for pure torsional buckling in case of S-S, and 2 minor-axis solutions, 2 major-axis solutions, and 2 pure torsional solutions in case of C-C. The



forces (and corresponding buckled shapes) for a certain mode (e.g. minor-axis mode) in case of S-S, while only the lowest 2 critical forces in case of C-C end restraint conditions. Obviously, higher-order modes exist and could be found by applying mode elements (i.e., a finer discretization). The lowest critical forces can reasonably well approximated even by the applied two finite elements, both in S-S and C-C. The error of approximation of the lowest critical values is less than 1%. For critical forces associated with higher modes the tendency is obvious: the higher the mode, the larger the error. The error is caused solely by the longitudinal shape functions. It is known that the sine-cosine functions belong to the exact solutions, which are approximated by the Hermite cubic polynomials in case of FEM. It is to note that solution for the lowest critical force for flexural buckling in case of S-S with using exactly the same approximate longitudinal function can be found in the literature, with exactly the same results. The EI flexural rigidity in classical solutions is equal to Ebt3/12 and Eb3t/12 for minor-axis and major-axis buckling, respectively. In case of the presented FEM solution the rigidities are Ebt3/12/(1-ν2) and Eb3t/12(1-ν2), respectively, where ν is the Poisson's ratio. The 1/(1-ν2) term is due to the fact that the presented FEM solution is based on shell model, and the global modes are defined with no transverse (membrane) extension, which is resulted in a (virtual) increase of axial and bending rigidities. The same phenomenon has already been reported and discussed e.g. in [28,30]. The EIw warping rigidity in classical solutions is equal to Eb3t3/144 (for the given rectangular cross-section, taking the second-order warping into account). In case of the presented FEM solution the rigidity is Eb3t3/144/(1-ν2). The 1/(1-ν2) term is due to the restrained (membrane) transverse strain, as discussed above.

4.3.2. Option b) The consideration of the (∂u (x, y, z ) /∂x )2 term in the secondorder longitudinal strain will change the equations systems from

Fig. 7. Distribution of the characteristic displacements for flexural and torsional buckling.

S. Ádány / Thin-Walled Structures 105 (2016) 135–146

145

Table 2 Critical forces for flexural and pure torsional buckling. End restr. Symmetry nr of half waves Exact c S-S S-S S-S S-S C-C C-C

Sym Point-sym Sym Point-sym Sym Point-sym

1 2 3 4 1 2

π2 ¼9.8696 4π2 ¼39.478 9π2 ¼88.826 16π2 ¼157.91 4π2 ¼39.478 80.763

FEM approximation c 9.9438 48.00 128.72 240.00 40.00 120.00

Eq. (103). It turns out that the system is not rank-deficient any more, therefore, the number of buckling modes is equal to the number of equations: 16 in case of S-S, and 10 in case of C-C. The 4 additional buckling modes (e.g., two symmetrical, two pointsymmetrical) involve longitudinal displacements only, that is why they can readily be referred as to axial modes. These axial modes are independent of the end restraints, due to the fact “C” support is defined so that the longitudinal displacement is not restrained, therefore, S-S and C-C cases are practically identical from the longitudinal displacements point-of-view. The axial modes are illustrated in Fig. 8, where the distribution of the characteristic displacement, i.e., u displacement is shown along the member length. The critical forces that belong to the additional (axial) modes are identical for all the 4 modes, and can be given as follows:

EA Fcr, A= 1−ν 2

(107)

where A is the cross-sectional area (i.e., in this case ¼bt), and ν is the Poisson's ratio. The consideration of the (∂u (x, y, z ) /∂x )2 term will have some effect on the critical force formulae to flexural and pure torsional modes, too, as follows. The flexural modes are changed in accordance with the following formula:

Fcr, F (b)=

1 1 Fcr , F (a)

+

Fcr , T (a)

+

Fcr , A

(108)

Fw Fcr , A F w + Ft

Ft =GIt

(110)

(111)

where Iw is the warping constant (with considering secondary warping) and It is the torsion inertia and c is the same constant listed in Table 2. Some comments:

 The additional modes are characterized by no flexural dis-





The loading is a concentric end moment (about the strong axis of the beam), which is assumed to distribute over the end sections, linearly changing with y.

(109)

EIw

( 1−ν2) L2



px =

1

where

Fw=c



[28], the longitudinal shape function is assumed to be cosine (with one or multiple waves). In case of FEM, there are two identical critical forces, which means that the buckled shape can be any arbitrary linear combination of the corresponding shape functions (with the ua and ub and ud and ue displacement parameters). The option b) critical force formulae for F and T buckling are in full accordance with the earlier findings of [28]. Note, the here presented option b) formulae correspond to the yyy case of [28]. However, it is to highlight that the earlier formulae have been derived by using the exact longitudinal shape functions, while here the same formulae are found for approximate longitudinal shape functions. Again, in [28] only S-S end restraint has been considered, however, here the formulae are found to be valid for other (namely C-C) cases, too.

4.4. Lateral-torsional buckling of beams

1 1

 In case of the analytical solution for the axial mode presented in

1

The pure torsional modes are changed in accordance with the following formula:

Fcr, T (b)=

Fig. 8. Distribution of the characteristic u displacement for axial buckling.

placement, that is the axis of the member remains straight during the buckling, while the longitudinal displacements are non-zero. This type of mode is discussed in [28, 29], called axial mode. In case of axial mode the value of the critical force is independent of whether the end restraint is S-S or C-C.

M ⎛ b⎞ ⎜ y− ⎟ b3t /12 ⎝ 2⎠

(112)

For the sake of simplicity, only option a) is discussed here. In case of S-S there are only 8 existing solutions, 4 plus/minus pairs. In case of C-C there are only 4 existing solutions, 2 plus/ minus pairs. The critical moments can be written as follows:

Mcr = ±

Fcr, F ( Ft + Fw ) = ±

⎛ ⎞ cEI ⎜ GIT + cEIw ⎟ ⎜ 2 2 2 2 ( 1−ν ) L ⎝ ( 1−ν ) L ⎟⎠

(113)

where Fcr , F is the flexural critical force to minor-axis buckling, and the c constants are exactly identical to those given in Table 2. Here are some comments.

 The applied discretization provides only the 4 lowest critical



moments (and corresponding buckled shapes) in case of S-S, while only the lowest 2 critical moments in case of C-C end restraint conditions. Obviously, higher-order modes exist and could be found by applying more elements (i.e., a finer discretization). The lowest critical moments can reasonably well be approximated even by the applied two finite elements, both in S-S and

146





S. Ádány / Thin-Walled Structures 105 (2016) 135–146

C-C. The error of approximation of the lowest critical values is less than 1%. For critical moments associated with higher modes the tendency is obvious: the higher the mode, the larger the error. The error is caused solely by the longitudinal shape functions. It is known that the sine-cosine functions belong to the exact solutions, which are approximated by the Hermite cubic polynomials in case of FEM. The appearance of the 1/(1-ν2) term is due to the restrained (membrane) transverse strain, as discussed above.

5. Concluding remarks In this paper a novel shell finite element is introduced. The introduced element is specifically intended to be used in constrained finite element method. The novel finite element is derived from the semi-analytical finite strip method, by keeping the transverse interpolation functions of FSM, but changing the longitudinal interpolation functions into classic polynomials. As it is proved, the proper selection of the longitudinal polynomial interpolation functions makes it possible to apply essentially the same constraining technique that is used in constrained finite strip method. This paper focuses on the novel finite element. The derivation of the new element is shown in detail. It is proved that the important mechanical criteria of the constraining procedure can directly and easily be implemented into the newly proposed element. Some basic semi-analytical examples are provided to illustrate the applicability of the proposed shell finite element. The new finite element can be used to model thin-walled structural members. If a highly regular mesh is used for the discretization, the constraining can be done, that is modal decomposition of the deformation modes can be completed, essentially identically to how cFSM solves the modal decomposition. Hence, the method that uses the here-introduced finite element can properly be referred as to constrained finite element method (cFEM). Since the so-defined cFEM method is essentially a shell FEM, it can handle a wide range of practical problems together with modal decomposition: members with holes, various loading and restraint conditions, shear buckling, web crippling, certain crosssection changes along the length. To all these problems first- or second-order static analysis, linear buckling analysis, dynamic analysis, etc. can be completed, with or without modal decomposition. Some preliminary results of such cFEM calculations have already been reported in [31]. More details and various applications are going to be published in subsequent papers in the near future.

Acknowledgments The work was conducted with the financial support of the OTKA K108912 project of the Hungarian Scientific Research Fund.

References [1] Y.K. Cheung, Finite strip method in the analysis of elastic paltes with two opposite ends simply supported, Proc. Inst. Civ. Eng. 40 (1968) 1–7. [2] Y.K. Cheung, Finite Strip Method in Structural Analysis, Pergamon Press, United Kingdom, 1976. [3] G.J. Local Hancock, Distortional, and lateral buckling of I-beams, ASCE J. Struct. Eng. 104 (11) (1978) 1787–1798. [4] B.W. Schafer, Cold-Formed Steel Behavior and Design: Analytical and Numerical Modeling of Elements and Members with Longitudinal Stiffeners (PhD

Dissertation), Cornell University, Ithaca, NY, 1997. [5] S. Ádány, Buckling mode classification of members with open thin-walled crosssections by using the Finite Strip Method, Research Report, Johns Hopkins University, 2004. (available at) 〈http://www.ce.jhu.edu/bschafer〉. [6] S. Ádány, B.W. Schafer, Buckling mode classification of members with open thin-walled cross-sections, Proceedings of the Fourth International Conference on Coupled Instabilities in Metal Structures, (Rome,Italy); Sept 27–29, 2004. pp. 467–476. [7] S. Ádány, B.W. Schafer, Buckling mode decomposition of single-branched open cross-section members via Finite Strip Method: derivation, Thin-Walled Struct. 44 (5) (2006) 563–584. [8] S. Ádány, B.W. Schafer, Buckling mode decomposition of single-branched open cross-section members via Finite Strip Method: application and examples, Thin-Walled Struct. 44 (5) (2006) 585–600. [9] B.W. Schafer, S. Ádány, Buckling analysis of cold-formed steel members using CUFSM: conventional and constrained finite strip methods, Proceedings of 18th International Specialty Conference on Cold-Formed Steel Structures, October 26–28, 2006, Orlando, USA, pp. 39–54. [10] S. Ádány, B.W. Schafer, A full modal decomposition of thin-walled, singlebranched open cross-section members via the constrained finite strip method, J. Constr. Steel Res. 64 (1) (2008) 12–29. [11] Z. Li, B.W. Schafer, Finite Strip Stability Solutions for General Boundary Conditions and the Extension of the Constrained Finite Strip Method, in Trends in Civil and Structural Engineering Computing.. Stirlingshire, UK: Saxe-Coburg Publications, 2009. [12] Z. Li, B.W. Schafer, Buckling analysis of cold-formed steel members with general boundary conditions using CUFSM: Conventional and constrained finite strip methods, Proceedings of the 20th International Specialty Conference on Cold-Formed Steel Structures - Recent Research and Developments in ColdFormed Steel Design and Construction, p. 17, 2010. [13] S. Ádány, B.W. Schafer, Generalized constrained finite strip method for thinwalled members with arbitrary cross-section: primary modes, Thin-Walled Struct. 84 (2014) 150–169. [14] S. Ádány, B.W. Schafer, Generalized constrained finite strip method for thinwalled members with arbitrary cross-section: secondary modes, orthogonality, examples, Thin-Walled Struct. 84 (2014) 123–133. [15] S. Ádány, Decomposition of in-plane shear in thin-walled members, ThinWalled Struct. 73 (2013) 27–38. [16] M. Djafour, H. Dib, M. Djelli, N. Djafour, M. Matallah, D. Zendagui, Buckling Mode Decomposition of Thin-Walled Members using a Constrained Spline Finite Strip Method, Proceeding of the 6th International Conference on Coupled Instabilities in Metal Structures, (Glasgow, Scotland), Dec 3–5, 2012. (Eds: J Loughlan, D H Nash, J Rhodes), pp 467–474. [17] M. Casafont, F. Marimon, M.M. Pastor, Calculation of pure distortional elastic buckling loads of members subjected to compression via the finite element method (June-July), Thin-Walled Struct. 47 (6–7) (2009) 701–729. [18] M. Casafont, F. Marimon, M.M. Pastor, M. Ferrer, Linear buckling analysis of thin-walled members combining the Generalised Beam Theory and the Finite Element Method, Comput. Struct. 89 (21–22) (2011) 1982–2000. [19] R. Schardt, Verallgemeinerte Technische Biegetheorie, Springer Verlag, Berlin, Heidelberg, 1989. [20] P.B. Dinis, D. Camotim, N. Silvestre, GBT formulation to analyse the buckling behaviour of thin-walled members with arbitrarily ‘branched’ open crosssections (January), Thin-Walled Struct. 44 (1) (2006) 20–38, ISSN 0263–8231. [21] R. Gonçalves, M. Ritto-Corrêa, D. Camotim, A new approach to the calculation of cross-section deformation modes in the framework of generalized beam theory, Comput. Mech. 46 (5) (2010) 759–781. [22] R. Bebiano, R. Gonçalves, D. Camotim, A cross-section analysis procedure to rationalise and automate the performance of GBT-based structural analyses (July), Thin-Walled Struct. 92 (2015) 29–47, ISSN 0263–8231. [23] M. Casafont, J. Bonada, M.M. Pastor, F. Roure, GBT calculation of distortional and global buckling of cold-formed steel channel columns with multiple perforations, Eighth International Conference on Advances in Steel Structures, Lisbon, Portugal, July 22–24, 2015, paper no: 87, p. 18. [24] J. Cai, C.D. Moen, Generalized beam theory buckling analysis for members with holes, Eighth International Conference on Advances in Steel Structures, Lisbon, Portugal, July 22–24, 2015, paper no: 86, p. 18. [25] CUFSM: Elastic Buckling Analysis of Thin-Walled Members by Finite Strip) Analysis. CUFSM v3.12, 2006. 〈http://www.ce.jhu.edu/bschafer/cufsm〉. [26] CUFSM: Elastic Buckling Analysis of Thin-Walled Members by Finite Strip) Analysis. CUFSM v4.05, 2012. 〈http://www.ce.jhu.edu/bschafer/cufsm〉. [27] Z. Li, M.T. Hanna, S. Ádány, B.W. Schafer, Impact of basis, orthogonalization, and normalization on the constrained finite strip method for stability solutions of open thin-walled members, Thin-Walled Struct. 49 (9) (2011) 1108–1122. [28] S. Ádány, Global buckling of thin-walled columns: analytical solutions based on shell model, Thin-Walled Struct. 55 (2012) 64–75. [29] S. Ádány, D. Visy, Global buckling of thin-walled columns: numerical studies, Thin-Walled Struct. 54 (2012) 82–93. [30] S. Ádány, N. Silvestre, B.W. Schafer, D. Camotim, GBT and cFSM: two modal approaches to the buckling analysis of unbranched thin-walled members, Int. J. Adv. Steel Constr. 5 (2) (2009) 195–223. [31] S. Ádány, Constrained finite element method: demonstrative examples on the global modes of thin-walled members, Proceedings of the 22nd International Speciality Conference on Cold-Formed Steel Structures, Nov 5–6, 2014, (St. Louis, USA (eds: R.A. LaBoube, W.W. Yu), 2014, pp. 67–82.