# Finite element method for stability and vibration analyses of thin-walled beams with arbitrary cross-section

Masoumeh Soltani, Behrouz Asgarian, Foudil Mohri

► To cite this version:

Masoumeh Soltani, Behrouz Asgarian, Foudil Mohri. Finite element method for stability and vibration analyses of thin-walled beams with arbitrary cross-section. The 2013 World Congress on Advances in Structural Engineering and Mechanics (ASEM2013), Sep 2013, Jeju (virtual), South Korea. pp.1403-1417. hal-04659523

HAL Id: hal-04659523

https://hal.science/hal-04659523v1

Submitted on 23 Jul 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Finite element method for stability and vibration analyses of thin-walled beams with arbitrary cross-section

*M. Soltani<sup>1)</sup>, B. Asgarian<sup>2)</sup>, and F. Mohri<sup>3)</sup>

1),2) Civil Engineering Faculty, K.N. Toosi university of Technology, Tehran, Iran

<sup>1)</sup> msoltani2@yahoo.com

3) Université de Lorraine, Laboratoire d'Etude des Microstructures et de Mécanique des Matériaux (LEM3), Metz. France

## ABSTRACT

In this paper, a generalized numerical and theoretical method is presented for the free vibration and stability analyses of tapered thin-walled beams with non-symmetric cross-section. The total potential energy is derived for an elastic behavior from the strain energy, the kinetic energy and work of the applied loads. Free vibration is considered in presence of harmonic excitations. The effects of the initial stresses and load eccentricities are also considered in the study. The governing equilibrium equations, motion equations and the associated boundary conditions are determined from the stationary condition. ......

## 1. INTRODUCTION

The use of thin-walled beams with I, Z and C cross sections has been increasing in many steel constructions due to their ability to utilize structural material more efficiently and optimize the distribution of weight. Flexural-torsional stability and vibration analyses are important topics in the design of these structures and research of accurate models has interested many researchers over the world from the last mid century until now. Since the early works of (Timoshenko 1961) and (Vlasov 1959) there have been extensive investigations on the elastic flexural-torsional buckling behavior of thin walled beam members. In these investigations, (Chen 1987) and (Bazant 1991), the buckling loads were evaluated either by closed-form solutions of the fourth-order differential equations governing the twisting and bending of the thin-walled beams or by using of complementary energy principle. Mentioned studies were focused on the stability and vibration analyses of prismatic thin-walled beams. However, in the last decades, thin-walled beams with variable cross-sections are extensively adopted in

<footnote n="1">¹ PhD Candidate</footnote>

<footnote n="2">² Associate Professor</footnote>

<footnote n="3">³ Assistant Professor</footnote>
different steel constructions. The investigations of elastic flexural-torsional buckling and vibration modes of tapered thin-walled with the use of efficient numerical techniques such as finite element method in the validation process has attracted many researchers. Among the first investigations on this topic, the most important ones are the studies of (Yang 1987) who formulated a finite element model for the beam that takes into account the effect of non-uniform torsion. More recently, (Kim 2000) investigated the linear stability and free vibration behavior of doubly symmetric I tapered thin-walled beams by a finite element approach. A theoretical model for large torsion context and equivalent beam finite element formulation in finite torsion amplitudes were developed by (Ronagh 2000). (Chen 2000) applied Hamilton's principle to the potential energy and obtained the motion equations governing of non-prismatic thin-walled beams. Based on the Rayleigh-Ritz method, a general variational formulation to analyze the lateral-torsional buckling behavior of tapered thin-walled beams with singly symmetric cross-sections was presented by (Andrade 2007). In (Asgarian 2013), the lateral torsional behavior of tapered beams with singly symmetric I cross-sections was studied. The previous work is extended to include eccentrically axial forces and dynamic loads for beams with non-symmetric cross-sections.

## 2. DERIVATION OF EQUILIBRIUM AND MOTION EQUATIONS FOR BUCKLING AND FREE VIBRATION

### 2.1. Kinematics

A tapered thin-walled beam with arbitrary cross-section is considered in the study (Fig.1a). A direct rectangular co-ordinate system is chosen, with x the initial longitudinal axis and y and z the first and second main bending axes. The origin of these axes is located at the centroid O. The shear centre with co-ordinates ($y_c, z_c$) in xyz is denoted C (Fig.1b). In the present study, the beam is subjected to arbitrary distributed force $p_x$ in structural domain in x direction. Under these conditions, the three displacement components of a point M on the section contour can be expressed as function of those of the shear center C:

$$U(x, y, z) = u_0(x) - y \frac{\partial(v(x) + z_c(x)\theta(x))}{\partial x} - z \frac{\partial(w(x) - y_c(x)\theta(x))}{\partial x} - \omega(y, z) \frac{\partial\theta(x)}{\partial x} \quad (1)$$

$$V(x, y, z) = v(x) - (z - z_c(x))\theta(x) \quad (2)$$

$$W(x, y, z) = w(x) + (y - y_c(x))\theta(x) \quad (3)$$

In which U represents the axial displacement. The displacement components V and W represent lateral and vertical displacements (in y and z directions). The torsion angle is denoted by $\theta$. The term $\omega(y, z)$ is the warping function, which can be defined based on Saint Venant's torsion theory. The Green's strain tensor components which incorporate the large displacements and including linear and non-linear strain part are given by:
<figure><img src="image_2.png" /><figcaption>Fig. 1 (a) Thin-walled beam with non-symmetric cross-section, (b) Coordinate system and notation for displacement parameters of the beam, (c) The nodal displacements of a thin-walled beam</figcaption></figure>

$$ \epsilon_{ij} = \frac{1}{2} \left( \frac{\partial U_i}{\partial x_j} + \frac{\partial U_j}{\partial x_i} \right) + \frac{1}{2} \left( \frac{\partial U_k}{\partial x_i} \frac{\partial U_k}{\partial x_j} \right) = \epsilon_{ij}^l + \epsilon_{ij}^* \quad i,j,k = x,y,z \quad (4) $$

The stress components for elastic homogeneous and isotropic materials are:

$$ \sigma_{xx} = E \epsilon_{xx}, \quad \tau_{xy} = G \gamma_{xy}, \quad \tau_{xz} = G \gamma_{xz}, \quad \sigma_{yy} = \sigma_{yy} = \tau_{zy} = 0 \quad (5) $$

The most general case of the normal and the shear stresses associated with the axial force **N**, the bending moments **M**<sub>y</sub> and **M**<sub>z</sub>, and shear forces **V**<sub>y</sub> and **V**<sub>z</sub> are considered as:

$$ \sigma_{xx}^0 = -\frac{M_y}{I_y}z - \frac{M_z}{I_z}y, \quad \tau_{xy}^0 = \frac{V_y}{A} = \frac{M_z'}{A}, \quad \tau_{xz}^0 = \frac{V_z}{A} = \frac{M_y'}{A} \quad (6) $$

$\sigma_{xx}^0$ denotes the initial stress and $\tau_{xy}^0$ and $\tau_{xz}^0$ are the mean value of the shear stresses.

## 2.2. Elastic equilibrium and motion equations

The equations of motion for non-prismatic thin-walled beam are derived by variation of total potential energy which is:

$$ \delta(U_l + U_0 + U_M - W) = 0 \quad (7) $$

In the last formulation, $\delta$ denotes a virtual variation. Where $U_l$ and $U_0$ are the elastic strain energy and the strain energy due to effects of the initial stresses, $U_M$ the kinetic energy, and $W$ the external load work. Their relationships for each term of the total potential energy are developed separately in the following:
$$
\begin{align}
\delta U_l &= \int_0^L \int_A (E(\varepsilon_{xx}^l \delta \varepsilon_{xx}^l) + G(\gamma_{xy}^l \delta \gamma_{xy}^l) + G(\gamma_{xz}^l \delta \gamma_{xy}^l)) dA dx \nonumber \\
\delta U_0 &= \int_0^L \int_A (\sigma_{xx}^0 \delta \varepsilon_{xx}^* + \tau_{xy}^0 \delta \gamma_{xy}^* + \tau_{xz}^0 \delta \gamma_{xz}^*) dA dx \tag{8} \\
\delta U_M &= \Omega^2 \int_0^L \int_A \rho (U \delta U + W \delta W + V \delta V) dA dx \nonumber \\
\delta W &= \int_L (p_x \delta u_p) dx \nonumber
\end{align}
$$

Where $\rho$ denotes the material density, $L$ and $A$ are the element length and area, respectively, and $\Omega$ is the natural frequency. For the external load work $W$ of the beam, an axial load $p_x$ is applied along a line ($PP'$) on the section contour (Fig.1a). In Eq. (8d), $u_p$ is the axial displacement of point $P$. The displacements of the point $P$ are:

$$ u_p = u_0 - y_p v_0' - y_p z_c' \theta - y_p z_c \theta' - z_p w_0' + z_p y_c' \theta + z_p y_c \theta' - \omega_p \theta' \quad (9) $$

Substituting the variation of the strain-displacement relations defined Eq. (4) and initial stresses Eq. (6) into Eq. (8a-d), and integration over the cross-section in the context of principal bending and warping axes conditions, the following equations are obtained:

$$
\begin{align}
\delta U_l &= \int_L (EAu_0' \delta u_0' + EI_z v'' \delta v'' + EI_y w'' \delta w'' + EI_\omega \theta'' \delta \theta'' + GJ \theta' \delta \theta') dx \nonumber \\
&\quad + \int_L EI_y (-y_c^2 \theta'' \delta \theta'' + y_c''^2 \theta \delta \theta + 4y_c'^2 \theta' \delta \theta' \theta - y_c'' \theta \delta w'' - y_c'' w'' \delta \theta \nonumber \\
&\quad \qquad - 2y_c' \theta' \delta w'' - 2y_c' w'' \delta \theta' + 2y_c' y_c'' \theta' \delta \theta + 2y_c' y_c'' \theta \delta \theta') dx \tag{10} \\
&\quad + \int_L EI_z (-z_c^2 \theta'' \delta \theta'' + z_c''^2 \theta \delta \theta + 4z_c'^2 \theta' \delta \theta' + z_c'' \theta \delta v'' + z_c'' v'' \delta \theta \nonumber \\
&\quad \qquad + 2z_c' \theta' \delta v'' + 2z_c' v'' \delta \theta' + 2z_c' z_c'' \theta' \delta \theta + 2z_c' z_c'' \theta \delta \theta') dx \nonumber
\end{align}
$$

$$
\begin{align}
\delta U_0 = \int_0^l & \left( N(v' \delta v' + w' \delta w' + r_c^2 \theta' \delta \theta' + z_c'^2 \theta \delta \theta + y_c'^2 \theta \delta \theta - (z_c z_c')' \theta \delta \theta \right. \nonumber \\
& \left. - (y_c y_c')' \theta \delta \theta + z_c \theta' \delta v' + z_c v' \delta \theta' + z_c' v' \delta \theta + z_c' \theta \delta v' - y_c \theta' \delta w' \right. \nonumber \\
& \left. - y_c w' \delta \theta' - y_c' w' \delta \theta - y_c' \theta \delta w' \right) - M_z w' \delta \theta' - M_z \theta' \delta w' + M_y v' \delta \theta' \nonumber \\
& + M_y \theta' \delta v' - M_z' w' \delta \theta - M_z' \theta \delta w' + M_y' v' \delta \theta + M_y' \theta \delta v' \nonumber \\
& - (M_z y_c'' + M_z'' y_c) \theta \delta \theta - (M_y z_c'' + M_y'' z_c) \theta \delta \theta - (\beta_y M_z + \beta_z M_y) \theta' \delta \theta' \right) dx \tag{11}
\end{align}
$$

$$
\begin{align}
\delta U_M = \Omega^2 \int_0^l & \left( \rho(Au_0 \delta u_0 + Av \delta v + Aw \delta w + I_S \theta \delta \theta + Az_c v \delta \theta + Az_c \theta \delta v \right. \nonumber \\
& \left. - Ay_c w \delta \theta - Ay_c \theta \delta w + I_\omega \theta' \delta \theta' - I_y y_c^2 \theta' \delta \theta' - I_z z_c^2 \theta' \delta \theta' \right. \nonumber \\
& \left. + I_z v' \delta v' + I_y w' \delta w' + I_z z_c'^2 \theta \delta \theta + I_y y_c'^2 \theta \delta \theta + z_c' I_z v' \delta \theta \right. \nonumber \\
& \left. + z_c' I_z \theta \delta v' - y_c' I_y w' \delta \theta - y_c' I_y \theta \delta w' \right) dx \tag{12}
\end{align}
$$
$$ W = \int_{0}^{l} (p_x \delta u_0 - M_{z0} \delta v' - M_{z0} z_c' \delta \theta - M_{z0} z_c \delta \theta' \\ - M_{y0} \delta w' + M_{y0} y_c' \delta \theta + M_{y0} y_c \delta \theta' + B_{\omega 0} \delta \theta') dx \quad (13) $$

In previous expressions, Eqs. (10)-(13), $I_\omega$ and $\theta$ denote warping and the Saint-Venant torsion constants. $I_y$ and $I_z$ are the second moments of area. $z_c$ and $y_c$ are the shear coordinates according to Vlasov's model. $\beta_y$ and $\beta_z$ are Wagner's coefficients which the exact formulations can be referred to Asgarian (2013). External bending moments and bimoment are:

$$ M_{z0} = p_x y_p, \quad M_{y0} = p_x z_p \quad \text{and} \quad B_{\omega 0} = -p_x \omega_p. $$

And the following geometrical relations are also used:

$$ r_c^2 = \frac{I_y}{A} + \frac{I_z}{A} + y_c^2 + z_c^2, \quad I_s = I_y + I_z + A(y_c^2 + z_c^2) \quad (14) $$

In the case of the constant axial load $N$ applied at the beam supports and according to Eq. (7) with respect to $u_0, v, w$ and $\theta$, the equations of motion for a non-prismatic thin-walled beam with non-symmetric cross-section are derived as:

$$ (EAu_0')' + \rho\omega^2 Au_0 = p_x \quad (15) $$

$$ (EI_z v'')'' + (EI_z z_c'' \theta)'' + 2(EI_z z_c' \theta')'' - (M_y \theta)'' - N(v + z_c \theta)'' \\ + \rho \Omega^2 A v + \rho \Omega^2 A z_c \theta - \Omega^2 (\rho I_z v')' - \Omega^2 (\rho z_c' I_z \theta)' = p_y + (M_{z0})' \quad (16) $$

$$ (EI_y w'')'' - (EI_y y_c'' \theta)'' - 2(EI_y y_c' \theta')'' + (M_z \theta)'' - N(w - y_c \theta)'' \\ + \rho \Omega^2 A w - \rho \Omega^2 A y_c \theta - \Omega^2 (\rho I_y w')' + \Omega^2 (\rho y_c' I_y \theta)' = p_z + (M_{y0})' \quad (17) $$

$$ \begin{aligned} & (EI_{\bar{\omega}} \theta'')'' - (GJ \theta')' - 4E(y_c'^2 I_y + z_c'^2 I_z)\theta' + E(y_c''^2 I_y + z_c''^2 I_z)\theta(x) + 2E(y_c' I_y w'')' \\ & - 2E(z_c' I_z v'')' - E(y_c'' I_y w'') + E(z_c'' I_z v'') - 2E(y_c' y_c'' I_y + z_c' z_c'' I_z)\theta(x) + M_z w'' - M_y v'' \\ & - (M_z y_c'' + M_z'' y_c)\theta(x) - (M_y z_c'' + M_y'' z_c)\theta(x) + ((\beta_z M_y + \beta_y M_z)\theta')' - N(r_c^2 \theta')' \\ & - N z_c'' + N y_c'' - N(z_c z_c'' + y_c y_c'')\theta(x) + \rho \Omega^2 I_s \theta(x) - \Omega^2 (\rho I_{\bar{\omega}} \theta')' + \rho \Omega^2 A z_c(x) v(x) \\ & - \rho \Omega^2 A y_c(x) w(x) + \rho \Omega^2 (y_c'^2 I_y + z_c'^2 I_z)\theta(x) + \rho \Omega^2 I_z z_c' v' - \rho \Omega^2 I_y y_c' w' \\ & = \hat{M}_t - M_t \theta(x) + (M_{z0} z_c - M_{y0} y_c) - (B_{\omega 0})' \end{aligned} \quad (18) $$

$$ \text{In Eq. (18), we put } (I_{\bar{\omega}} = I_{\omega} - y_c^2 I_y - z_c^2 I_z). $$

In what follows, a numerical method based on power series expansions is adopted for solution of the differential equations of motion (15-18). This procedure was used by different authors to evaluate the critical buckling loads and natural frequencies of the beam-column members (Asgarian 2013, Eisenberger 1994 and 1995 and Al-Sadder 2004).
### 3. NUMERICAL APPROACH

In the last expressions, due to presence of tapered elements, all the section properties including bending and torsion characteristics of the beam section, external and torsion moments are variable over the beam's length. For this reason, all these terms are taken as the following power series form:

$$ 
\begin{aligned}
I_z(x) &= \sum_{i=0}^{\infty} I_{zi} x^i, & I_y(x) &= \sum_{i=0}^{\infty} I_{yi} x^i, & I_{\bar{\omega}}(x) &= \sum_{i=0}^{\infty} I_{\bar{\omega}i} x^i, & I_s(x) &= \sum_{i=0}^{\infty} I_{si} x^i \\
A(x) &= \sum_{i=0}^{\infty} A_i x^i, & y_c(x) &= \sum_{i=0}^{\infty} y_{ci} x^i, & z_c(x) &= \sum_{i=0}^{\infty} z_{ci} x^i, & J(x) &= \sum_{i=0}^{\infty} J_i x^i \\
R_c &= r_c^2(x) = \sum_{i=0}^{\infty} R_{ci} x^i, & \rho(x) &= \sum_{i=0}^{\infty} \rho_i x^i, & \beta_z(x) &= \sum_{i=0}^{\infty} \beta_{zi} x^i, & \beta_y(x) &= \sum_{i=0}^{\infty} \beta_{yi} x^i
\end{aligned} 
\quad (19) $$

In order to solve the equations of motion, a non-dimensional variable ($\epsilon = \frac{x}{L}$) is introduced and three displacement parameters are also presented by the following infinite power series of the form:

$$ v(\epsilon) = \sum_{i=0}^{\infty} a_i \epsilon^i, \quad w(\epsilon) = \sum_{i=0}^{\infty} b_i \epsilon^i, \quad \theta(\epsilon) = \sum_{i=0}^{\infty} c_i \epsilon^i \quad (20) $$

By substituting Eq. (19) and Eq. (20) into Eqs. (16-18) and multiplying the two series in each terms of equations, the following recurrence formula for the coefficient $a_{k+4}$, $b_{k+4}$ and $c_{k+4}$ are obtained:

$$ 
\begin{aligned}
a_{k+4} = & \frac{-1}{EI_{z_0}^*(k+4)(k+3)(k+2)(k+1)} \\
& \times \left( E \sum_{i=1}^{k+2} I_{z_i}^* a_{k-i+4} (k-i+4)(k-i+3)(k+2)(k+1) \right. \\
& \left. + E \sum_{j=0}^{k+2} \sum_{i=0}^{j} I_{z_j}^* z_{c_{j-i+2}}^* c_{k-j+2} (j-i+2)(j-i+1)(k+2)(k+1) \right. \\
& \left. + 2E \sum_{j=0}^{k+2} \sum_{i=0}^{j} I_{z_j}^* z_{c_{j-i+1}}^* c_{k-j+3} (j-i+1)(k-j+3)(k+2)(k+1) \right. \\
& \left. - L^2 \sum_{i=0}^{k+2} M_{y_i}^* c_{k-i+2} (k+2)(k+1) - L^2 N a_{k+2} (k+1)(k+2) \right. \\
& \left. - L^2 N \sum_{i=0}^{k+2} z_{c_i}^* c_{k-i+2} (k+1)(k+2) + L^4 \Omega^2 \sum_{j=0}^{k} \sum_{i=0}^{j} \rho_i^* A_{j-i}^* c_{k-j} \right. \\
& \left. + L^4 \Omega^2 \sum_{m=0}^{k} \sum_{j=0}^{m} \sum_{i=0}^{j} \rho_i^* A_{j-i}^* z_{c_{m-j}}^* c_{k-m} - L^2 \Omega^2 \sum_{j=0}^{k+1} \sum_{i=0}^{j} \rho_i^* I_{z_{j-i}}^* a_{k-j+2} (k-j+2)(k+1) \right. \\
& \left. - L^2 \Omega^2 \sum_{m=0}^{k+1} \sum_{j=0}^{m} \sum_{i=0}^{j} \rho_i^* I_{z_{j-i}}^* z_{c_{m-j+1}}^* c_{k-m+1} (m-j+1)(k+1) - (L^4) p_{y_k}^* - (L^3)(k+1) M_{z_0 k+1}^* \right)
\end{aligned} 
\quad (21) $$
$$
\begin{equation}
\begin{aligned}
b_{k+4} = & \frac{-1}{EI_{y_0}^*(k+4)(k+3)(k+2)(k+1)} \\
& \times \left( E \sum_{i=1}^{k+2} I_{y_i}^* b_{k-i+4} (k-i+4)(k-i+3)(k+2)(k+1) \right. \\
& \qquad \left. - E \sum_{j=0}^{k+2} \sum_{i=0}^{j} I_{y_i}^* y_{c_{j-i+2}}^* c_{k-j+2} (j-i+2)(j-i+1)(k+2)(k+1) \right. \\
& \qquad \left. - 2E \sum_{j=0}^{k+2} \sum_{i=0}^{j} I_{y_i}^* y_{c_{j-i+1}}^* c_{k-j+3} (j-i+1)(k-j+3)(k+2)(k+1) \right. \\
& \qquad \left. + L^2 \sum_{i=0}^{k+2} M_{z_i}^* c_{k-i+2} (k+2)(k+1) - L^2 Nb_{k+2} (k+1)(k+2) \right. \\
& \qquad \left. + L^2 N \sum_{i=0}^{k+2} y_{c_i}^* c_{k-i+2} (k+1)(k+2) + L^4 \Omega^2 \sum_{j=0}^{k} \sum_{i=0}^{j} \rho_i^* A_{j-i}^* b_{k-j} \right. \\
& \qquad \left. - L^4 \Omega^2 \sum_{m=0}^{k} \sum_{j=0}^{m} \sum_{i=0}^{j} \rho_i^* A_{j-i}^* y_{c_{m-j}}^* c_{k-m} - L^2 \Omega^2 \sum_{j=0}^{k+1} \sum_{i=0}^{j} \rho_i^* I_{y_{j-i}}^* b_{k-j+2} (k-j+2)(k+1) \right. \\
& \qquad \left. + L^2 \Omega^2 \sum_{m=0}^{k+1} \sum_{j=0}^{m} \sum_{i=0}^{j} \rho_i^* I_{y_{j-i}}^* y_{c_{m-j+1}}^* c_{k-m+1} (m-j+1)(k+1) - (L^4) p_{z_k}^* - (L^3) (k+1) M_{y_0}^* \right)
\end{aligned}
\tag{22}
\end{equation}
$$

$$
\begin{equation}
\begin{aligned}
c_{k+4} = & \frac{-1}{EI_{\omega_0}^*(k+4)(k+3)(k+2)(k+1)} \\
& \times \left( E \sum_{i=1}^{k+2} I_{\omega_i}^* c_{k-i+4} (k-i+4)(k-i+3)(k+2)(k+1) - L^2 \sum_{i=0}^{k+1} GJ_i^* c_{k-i+2} (k+1)(k-i+2) \right. \\
& \qquad \left. - 4E \sum_{m=0}^{k+1} \sum_{j=0}^{m} \sum_{i=0}^{j} (I_{y_i}^* y_{c_{j-i+1}}^* y_{c_{m-j+1}}^* + I_{z_i}^* z_{c_{j-i+1}}^* z_{c_{m-j+1}}^*) c_{k-m+2} (j-i+1)(m-j+1)(k-m+2)(k+1) \right. \\
& \qquad \left. + E \sum_{m=0}^{k} \sum_{j=0}^{m} \sum_{i=0}^{j} (I_{y_i}^* y_{c_{j-i+2}}^* y_{c_{m-j+2}}^* + I_{z_i}^* z_{c_{j-i+2}}^* z_{c_{m-j+2}}^*) c_{k-m} (j-i+1)(j-i+2)(m-j+2)(m-j+1) \right. \\
& \qquad \left. + 2E \sum_{j=0}^{k+1} \sum_{i=0}^{j} (I_{y_i}^* y_{c_{j-i+1}}^* b_{k-j+3} - I_{z_i}^* z_{c_{j-i+1}}^* a_{k-j+3}) (j-i+1)(k-j+2)(k-j+3)(k+1) \right. \\
& \qquad \left. - E \sum_{j=0}^{k} \sum_{i=0}^{j} (I_{y_i}^* y_{c_{j-i+2}}^* b_{k-j+2} - I_{z_i}^* z_{c_{j-i+2}}^* a_{k-j+2}) (j-i+1)(j-i+2)(k-j+1)(k-j+2) \right. \\
& \qquad \left. - 2E \sum_{m=0}^{k} \sum_{j=0}^{m+1} \sum_{i=0}^{j} (I_{y_i}^* y_{c_{j-i+1}}^* y_{c_{m-j+3}}^* + I_{z_i}^* z_{c_{j-i+1}}^* z_{c_{m-j+3}}^*) c_{k-m} (j-i+1)(m-j+3)(m-j+2)(m+1) \right. \\
& \qquad \left. + L^2 \sum_{i=0}^{k} M_{z_i}^* b_{k-i+2} (k-i+2)(k-i+1) - L^2 \sum_{i=0}^{k} M_{y_i}^* a_{k-i+2} (k-i+2)(k-i+1) \right. \\
& \qquad \left. - L^2 \sum_{j=0}^{k} \sum_{i=0}^{j} M_{z_{j+2}}^* y_{c_{j-i}}^* (i+1)(i+2) c_{k-j} + L^2 \sum_{j=0}^{k+1} \sum_{i=0}^{j} (M_{y_i}^* \beta_{z_{j-i}}^* + M_{z_i}^* \beta_{y_{j-i}}^*) c_{k-j+2} (k+1)(k-j+2) \right)
\end{aligned}
\tag{23}
\end{equation}
$$
$$
\begin{align}
& -L^2 \sum_{j=0}^{k} \sum_{i=0}^{j} M_{y_i}^* z_{c_{j-i+2}}^* (j-i+1)(j-i+2)c_{k-j} - L^2 \sum_{j=0}^{k} \sum_{i=0}^{j} M_{z_i}^* y_{c_{j-i+2}}^* (j-i+1)(j-i+2)c_{k-j} \nonumber \\
& -L^2 \sum_{j=0}^{k} \sum_{i=0}^{j} M_{y_{i+2}}^* z_{c_{j-i}}^* (i+1)(i+2)c_{k-j} - L^2 N \sum_{i=0}^{k+1} R_{c_i}^* c_{k-i+2} (k-i+2)(k+1) \nonumber \\
& -L^2 N \sum_{i=0}^{k} z_{c_i}^* a_{k-i+2} (k-i+1)(k-i+2) + L^2 N \sum_{i=0}^{k} y_{c_i}^* b_{k-i+2} (k-i+1)(k-i+2) \nonumber \\
& -L^2 N \sum_{j=0}^{k} \sum_{i=0}^{j} (z_{c_i}^* z_{c_{j-i+2}}^* + y_{c_i}^* y_{c_{j-i+2}}^*) c_{k-j} (j-i+1)(j-i+2) + L^4 \Omega^2 \sum_{j=0}^{k} \sum_{i=0}^{j} \rho_i^* I_{s_{j-i}}^* c_{k-j} \nonumber \\
& -L^2 \Omega^2 \sum_{j=0}^{k+1} \sum_{i=0}^{j} \rho_i^* I_{\omega_{j-i}}^* c_{k-j+2} (k+1)(k-j+2) + L^4 \Omega^2 \sum_{m=0}^{k} \sum_{j=0}^{m} \sum_{i=0}^{j} \rho_i^* A_{j-i}^* z_{c_{m-j}}^* a_{k-m} \nonumber \\
& + L^2 \Omega^2 \sum_{n=0}^{k} \sum_{m=0}^{n} \sum_{j=0}^{m} \sum_{i=0}^{j} \rho_i^* (I_{y_{j-i}}^* y_{c_{m-j+1}}^* y_{c_{n-m+1}}^* + I_{z_{j-i}}^* z_{c_{m-j+1}}^* z_{c_{n-m+1}}^*) c_{k-n} (m-j+1)(n-m+1) \nonumber \\
& -L^4 \Omega^2 \sum_{m=0}^{k} \sum_{j=0}^{m} \sum_{i=0}^{j} \rho_i^* A_{j-i}^* y_{c_{m-j}}^* b_{k-m} + L^2 \Omega^2 \sum_{m=0}^{k} \sum_{j=0}^{m} \sum_{i=0}^{j} \rho_i^* I_{z_{j-i}}^* z_{c_{m-j+1}}^* a_{k-m+1} (m-j+1)(k-m+1) \nonumber \\
& -L^2 \Omega^2 \sum_{m=0}^{k} \sum_{j=0}^{m} \sum_{i=0}^{j} \rho_i^* I_{y_{j-i}}^* y_{c_{m-j+1}}^* b_{k-m+1} (m-j+1)(k-m+1) \nonumber \\
& -(L^3) \left( \sum_{i=0}^{k} z_{c_i}^* M_{z_{0-k+i+1}}^* (k-i+1) - \sum_{i=0}^{k} y_{c_i}^* M_{y_{0-k+i+1}}^* (k-i+1) \right) + (L^3) B_{\omega_{0k+1}}^* (k+1) \tag{23}
\end{align}
$$

Where

$$
\begin{align}
I_{z_i}^* &= I_{z_i} L^i, & I_{y_i}^* &= I_{y_i} L^i, & I_{\bar{\omega}_i}^* &= I_{\bar{\omega}_i} L^i, & I_{s_i}^* &= I_{s_i} L^i, & J_i^* &= J_i L^i \nonumber \\
& , & R_{c_i}^* &= R_{c_i} L^i, & A_i^* &= A_i L^i, & z_{c_i}^* &= z_{c_i} L^i, & y_{c_i}^* &= y_{c_i} L^i \nonumber \\
& , & \rho_i^* &= \rho_i L^i, & \beta_{z_i}^* &= \beta_{z_i} L^i, & \beta_{y_i}^* &= \beta_{y_i} L^i, & B_{\omega 0i}^* &= B_{\omega 0i} L^i \tag{24} \\
& , & M_{z_i}^* &= M_{z_i} L^i, & M_{y_i}^* &= M_{y_i} L^i, & M_{z_{0i}}^* &= M_{z_{0i}} L^i, & M_{y_{0i}}^* &= M_{y_{0i}} L^i \nonumber
\end{align}
$$

According to the recurrence formulas (21)-(23), the general solution of coupled system of differential equations (16)-(18) can be determined explicitly in terms of the twelve constants ($a_0, a_1, a_2, a_3$), ($b_0, b_1, b_2, b_3$) and ($c_0, c_1, c_2, c_3$). Then, the solutions of equilibrium equations can be expressed in the following matrix form:

$$
C = B \times A \tag{25}
$$

With

$$
A = \{a_0 \ a_1 \ a_2 \ a_3 \ b_0 \ b_1 \ b_2 \ b_3 \ c_0 \ c_1 \ c_2 \ c_3\}^T \quad (26)
$$

$$
B = \begin{bmatrix}
v_0(\varepsilon) & v_1(\varepsilon) & v_2(\varepsilon) & \cdots & v_{10}(\varepsilon) & v_{11}(\varepsilon) \\
w_0(\varepsilon) & w_1(\varepsilon) & w_2(\varepsilon) & \cdots & w_{10}(\varepsilon) & w_{11}(\varepsilon) \\
\theta_0(\varepsilon) & \theta_1(\varepsilon) & \theta_2(\varepsilon) & \cdots & \theta_{10}(\varepsilon) & \theta_{11}(\varepsilon)
\end{bmatrix} \quad (27)
$$
$$C = \begin{bmatrix} v(\varepsilon) \\ w(\varepsilon) \\ \theta(\varepsilon) \end{bmatrix} \qquad (28)$$

All terms of $v_i$, $w_i$ and $\theta_i$ are derived with the aid of the symbolic software MATLAB.

## 4. BOUDARY CONDITIONS AND SHAPE FUNCTION

Knowing that the all undefined coefficients $(a_0, a_1, a_2, a_3)$, $(b_0, b_1, b_2, b_3)$ and $(c_0, c_1, c_2, c_3)$ are functions of the displacements of degree of freedom (DOF) which can be derived by imposing four boundary conditions for each displacement parameters of a single-span beam (two at each end of the beam). Fig. 1c shows the nodal displacements of the beam element. Referring to this figure and using local coordinates at $x=0$ and $x=L$, the following boundary conditions must be satisfied:

$$D = V \times A \qquad (29)$$

Where

$$V = \begin{bmatrix} v_0(0) & v_1(0) & v_2(0) & \dots & v_{11}(0) \\ w_0(0) & w_1(0) & w_2(0) & \dots & w_{11}(0) \\ \theta_0(0) & \theta_1(0) & \theta_2(0) & \dots & \theta_{11}(0) \\ \frac{v_0'(0)}{L} & \frac{v_1'(0)}{L} & \frac{v_2'(0)}{L} & \dots & \frac{v_{11}'(0)}{L} \\ \frac{w_0'(0)}{L} & \frac{w_1'(0)}{L} & \frac{w_2'(0)}{L} & \dots & \frac{w_{11}'(0)}{L} \\ \frac{\theta_0'(0)}{L} & \frac{\theta_1'(0)}{L} & \frac{\theta_2'(0)}{L} & \dots & \frac{\theta_{11}'(0)}{L} \\ \frac{v_0(1)}{L} & \frac{v_1(1)}{L} & \frac{v_2(1)}{L} & \dots & \frac{v_{11}(1)}{L} \\ w_0(1) & w_1(1) & w_2(1) & \dots & w_{11}(1) \\ \theta_0(1) & \theta_1(1) & \theta_2(1) & \dots & \theta_{11}(1) \\ \frac{v_0'(1)}{L} & \frac{v_1'(1)}{L} & \frac{v_2'(1)}{L} & \dots & \frac{v_{11}'(1)}{L} \\ \frac{w_0'(1)}{L} & \frac{w_1'(1)}{L} & \frac{w_2'(1)}{L} & \dots & \frac{w_{11}'(1)}{L} \\ \frac{\theta_0'(1)}{L} & \frac{\theta_1'(1)}{L} & \frac{\theta_2'(1)}{L} & \dots & \frac{\theta_{11}'(1)}{L} \end{bmatrix} \qquad (30)$$

$$D = \begin{bmatrix} v^0 & w^0 & \theta^0 & v'^0 & w'^0 & \theta'^0 & v^1 & w^1 & \theta^1 & v'^1 & w'^1 & \theta'^1 \end{bmatrix} \qquad (31)$$

$$\Rightarrow C = B \times V \times D^{-1} \qquad (32)$$

In the linear finite element method of thin-walled beam, each element has six degrees of freedom at each node: vertical and lateral translations in the nodal z and y direction, respectively, and rotation about z-axis (the two nodes by which the finite element can be assembled into structure are located at its ends) (Fig. 1c). The
considered degrees of freedom at the left and right nods of each element respectively are: $v^0, v^1$ (the lateral displacement in the y direction), $w^0, w^1$ (the transverse displacement in the z direction), $\theta^0, \theta^1$ (the torsion angle), $v'^0, v'^1$ (the rotation about the z-axis, $\frac{\partial v}{\partial x}$), $w'^0, w'^1$ (the rotation about the y-axis, $\frac{\partial w}{\partial x}$), and $\theta'^0, \theta'^1$ (the rate of change twist, $\frac{\partial \theta}{\partial x}$).

The terms of the element mass and stiffness matrices can be found from the derivatives of the interpolation functions. The shape functions define the deformation shape of the element from applying unit translation, twist, warping or rotation at each of the 12 degrees of freedom, while constraining the other eleven nodal displacements. There are the following sets of boundary conditions:

$$ 
\begin{aligned}
& 1. v^0 = 1, w^0 = \theta^0 = v'^0 = w'^0 = \theta'^0 = v^1 = w^1 = \theta^1 = v'^1 = w'^1 = \theta'^1 = 0 \\
& 2. w^0 = 1, v^0 = \theta^0 = v'^0 = w'^0 = \theta'^0 = v^1 = w^1 = \theta^1 = v'^1 = w'^1 = \theta'^1 = 0 \\
& 3. \theta^0 = 1, w^0 = v^0 = v'^0 = w'^0 = \theta'^0 = v^1 = w^1 = \theta^1 = v'^1 = w'^1 = \theta'^1 = 0 \\
& 4. v'^0 = 1, v^0 = w^0 = \theta^0 = w'^0 = \theta'^0 = v^1 = w^1 = \theta^1 = v'^1 = w'^1 = \theta'^1 = 0 \\
& 5. w'^0 = 1, v^0 = w^0 = \theta^0 = v'^0 = \theta'^0 = v^1 = w^1 = \theta^1 = v'^1 = w'^1 = \theta'^1 = 0 \\
& 6. \theta'^0 = 1, v^0 = w^0 = \theta^0 = v'^0 = w'^0 = v^1 = w^1 = \theta^1 = v'^1 = w'^1 = \theta'^1 = 0 \\
& 7. v^1 = 1, v^0 = w^0 = \theta^0 = v'^0 = w'^0 = \theta'^0 = v^1 = w^1 = \theta^1 = v'^1 = w'^1 = \theta'^1 = 0 \\
& 8. w^1 = 1, v^0 = w^0 = \theta^0 = v'^0 = w'^0 = \theta'^0 = v^1 = w^1 = \theta^1 = v'^1 = w'^1 = \theta'^1 = 0 \\
& 9. \theta^1 = 1, v^0 = w^0 = v^0 = v'^0 = w'^0 = \theta'^0 = v^1 = w^1 = v'^1 = w'^1 = \theta'^1 = 0 \\
& 10. v'^1 = 1, v^0 = w^0 = \theta^0 = v'^0 = w'^0 = \theta'^0 = v^1 = w^1 = \theta^1 = v'^1 = w'^1 = \theta'^1 = 0 \\
& 11. w'^1 = 1, v^0 = w^0 = \theta^0 = v'^0 = w'^0 = \theta'^0 = v^1 = w^1 = \theta^1 = v'^1 = w'^1 = \theta'^1 = 0 \\
& 12. \theta'^1 = 1, v^0 = w^0 = \theta^0 = v'^0 = w'^0 = v^1 = w^1 = \theta^1 = v'^1 = w'^1 = \theta'^1 = 0
\end{aligned} 
\qquad (33) $$

According to these twelve boundary conditions corresponding to each shape function and using Eq. (32), 12 sets of interpolation functions can be derived.

## 5. LINEAR STIFFNESS MATRICES

The terms which lead to the element linear stiffness and the geometrical stiffness matrices can be extracted from Eq. (9) and Eq. (10) and the mass matrix can be derived by Eq. (11). The evaluated shape functions and the principle of internal virtual work can be used in formulating the element matrices for non-prismatic members. The generalized and geometrical stiffness matrices are derived as follows:
$$
\begin{equation}
\begin{aligned}
K_0 = E \int_L \{ &A[u'_0 \llbracket u'_0 \rrbracket]^T + I_z[v'' \llbracket v'' \rrbracket]^T + I_y[w'' \llbracket w'' \rrbracket]^T + I_{\bar{\omega}}[\theta'' \llbracket \theta'' \rrbracket]^T \\
& + I_z z_c''^2 [\theta \llbracket \theta \rrbracket]^T + 4I_z z_c'^2 [\theta' \llbracket \theta' \rrbracket]^T + I_z z_c'' ([v'' \llbracket \theta \rrbracket]^T + [\theta \llbracket v'' \rrbracket]^T) \\
& + 2I_z z_c' ([v'' \llbracket \theta' \rrbracket]^T + [\theta' \llbracket v'' \rrbracket]^T) + 2I_z z_c' z_c'' ([\theta \llbracket \theta' \rrbracket]^T + [\theta' \llbracket \theta \rrbracket]^T) \\
& + I_y y_c''^2 [\theta \llbracket \theta \rrbracket]^T + 4I_y y_c'^2 [\theta' \llbracket \theta' \rrbracket]^T - I_y y_c'' ([w'' \llbracket \theta \rrbracket]^T + [\theta \llbracket w'' \rrbracket]^T) \\
& - 2I_y y_c' ([\theta' \llbracket w'' \rrbracket]^T + [w'' \llbracket \theta' \rrbracket]^T) + 2I_y y_c' y_c'' ([\theta \llbracket \theta' \rrbracket]^T + [\theta' \llbracket \theta \rrbracket]^T) \} dx
\end{aligned}
\tag{34}
\end{equation}
$$

$$
\begin{equation}
\begin{aligned}
K_G = & \int_0^L \{ N([v'][v']^T + [w'][w']^T + R_c[\theta']\theta']^T - z_c z_c''[\theta]\theta]^T \\
& - y_c y_c''[\theta]\theta]^T + z_c ([\theta'][v']^T + [v']\theta']^T) \\
& + z_c' ([v']\theta]^T + [\theta][v']^T) - y_c ([w']\theta']^T + [\theta']\theta']^T) \\
& - y_c' ([w']\theta]^T + [\theta]\theta']^T) - M_z ([w']\theta']^T + [\theta']\theta']^T) \\
& + M_y ([\theta'][v']^T + [v']\theta']^T) - M_z' ([w']\theta]^T + [\theta]\theta']^T) \\
& + M_y' ([v']\theta]^T + [\theta]\theta']^T) - (M_z y_c'' + M_z'' y_c) [\theta]\theta]^T \\
& - (M_y z_c'' + M_y'' z_c) [\theta]\theta]^T - (\beta_y M_z + \beta_z M_y) [\theta']\theta']^T \} dx
\end{aligned}
\tag{35}
\end{equation}
$$

Where $K_{ij}^0$ is the first-order elastic stiffness matrix and $K_{Gij}$ is the geometric stiffness matrix, which accounts for the effect of the constant axial force ($N$) on the bending stiffness of the member. By applying the principle of virtual displacement along the column element with distributed mass, the mass matrix terms are given by:

$$
\begin{equation}
\begin{aligned}
M ={}& \rho\omega^2 \int_0^L \{ A[u_0 \llbracket u_0 \rrbracket]^T + A[v \llbracket v \rrbracket]^T + A[w \llbracket w \rrbracket]^T + I_s[\theta \llbracket \theta \rrbracket]^T \\
& + Az_c ([v \llbracket \theta \rrbracket]^T + [\theta \llbracket v \rrbracket]^T) - Ay_c ([w \llbracket \theta \rrbracket]^T + [\theta \llbracket w \rrbracket]^T) \\
& + I_{\bar{\omega}}[\theta' \llbracket \theta' \rrbracket]^T + I_z[v' \llbracket v' \rrbracket]^T + I_y[w' \llbracket w' \rrbracket]^T + I_z z_c'^2 [\theta \llbracket \theta \rrbracket]^T \\
& + I_y y_c'^2 [\theta \llbracket \theta \rrbracket]^T + z_c' I_z ([v' \llbracket \theta \rrbracket]^T + [\theta \llbracket v' \rrbracket]^T) \\
& - y_c' I_y ([w' \llbracket \theta \rrbracket]^T + [\theta \llbracket w' \rrbracket]^T) \} dx
\end{aligned}
\tag{36}
\end{equation}
$$

The structure stiffness and mass matrices can be obtained by assembling each element stiffness and mass matrices according to its nodal displacement. The process of assemblage is described in detail in most stability analysis textbook [1-4]. The critical buckling load and natural frequencies can be derived by solving the eigenvalue problems of following equations. In the buckling analysis, we have:

$$ (K^* + \lambda K_G)\phi = 0 \tag{37} $$

$\lambda$ are the eigenvalues and $\phi$ are the related eigenvectors. Under compressive loads, they lead to buckling loads and related eigenmodes. For the vibration, the
eigenvalue problem is put:

$$ (K^* - \lambda M)\phi = 0 \qquad (38) $$

Here $\lambda$ is related the eigenpulsation of the structure ($\lambda = \omega^2$). $\phi$ are the related vibration eigenmode.

## 6. APPLICATION

In this section, the accuracy and the efficiency of the power series method in stability and free vibration analyses of tapered thin-walled beams with arbitrary cross-sections and boundary conditions are investigated. In this regard, two comparative examples are studied. The buckling loads and natural frequencies provided by the numerical results of the aforementioned method are compared to finite element method by means of ANSYS code [15]. In this paper, each studied thin-walled beam is modeled with shell elements (SHELL63). Shell63 has both bending and membrane capabilities. Both in-plane and normal loads are permitted. The element has 6 degrees of freedom at each node, 3 translations in x, y and z directions and 3 rotations about the 3 axes. It is well known that for a system with n DOF, there exist n buckling modes and n vibration modes, but in practice only the lower ones are of interest.

### 6.1. Flexural-torsional buckling load of simply supported thin-walled beam-column

This section deals with the flexural-torsional and free vibration analyses of prismatic and web tapered simply supported beams. In the depicted beam, the top and bottom flanges are remained constant along the beam's length while, the web height is varying linearly from ($h_{max}$) at the left support to ($h_{min}$) at the other support. In this example, the tapering coefficient $\alpha$ has been varied from 0.4 to 1. In the case of stability analysis, the tip load is applied at the top flange of the bigger end section. In this example, two cross-sections are considered, a doubly and a singly symmetric I section with the bottom flange shorter than the top flange. The geometrical and material data are shown in Fig. 2.

<figure><img src="image_3.png" /><figcaption>Fig. 2 Simply supported tapered beam under axial compression load: geometry and material properties, (a) doubly symmetric section and (b) singly symmetric section.</figcaption></figure>

The values of the flexural-torsional critical buckling loads and the natural frequencies are derived according to the present study and finite element method by means of Ansys software using shell element are depicted in Table 1. The relative errors ($\Delta = (P_{cr} - P_{cr}^{FEM}) / P_{cr}^{FEM} \times 100$) between the present method and shell predictions are also provided.
<table><caption>Table 1: Critical buckling load $P_{cr}(kN)$ (stability analysis) and the natural frequency $\omega(\text{rad/s})$ (free-vibration analysis) for non-prismatic thin-walled beam.</caption><thead><tr><th rowspan="2">Analyses Type</th><th rowspan="2">α</th><th colspan="3">Equal Flanges</th><th colspan="3">Unequal Flanges</th></tr><tr><th>FEM (Ansys 5.4)</th><th>Present Method</th><th>Δ(%)</th><th>FEM (Ansys 5.4)</th><th>Present Method</th><th>Δ(%)</th></tr></thead><tbody><tr><td rowspan="4">Stability Analysis</td><td>0.4</td><td>164.84</td><td>163.72</td><td>0.68</td><td>103.65</td><td>103.416</td><td>0.23</td></tr><tr><td>0.6</td><td>165.58</td><td>165.137</td><td>0.27</td><td>103.54</td><td>103.486</td><td>0.05</td></tr><tr><td>0.8</td><td>166.57</td><td>166.188</td><td>0.23</td><td>103.07</td><td>102.78</td><td>0.28</td></tr><tr><td>1</td><td>167.54</td><td>168.813</td><td>0.76</td><td>102.29</td><td>101.77</td><td>0.51</td></tr><tr><td rowspan="4">Free Vibration Analysis</td><td>0.4</td><td>19.182</td><td>19.319</td><td>0.71</td><td>14.8</td><td>14.83</td><td>0.2</td></tr><tr><td>0.6</td><td>18.52</td><td>18.65</td><td>0.70</td><td>14.08</td><td>14.1</td><td>0.14</td></tr><tr><td>0.8</td><td>17.92</td><td>18.04</td><td>0.67</td><td>13.43</td><td>13.43</td><td>0.00</td></tr><tr><td>1</td><td>17.37</td><td>17.5</td><td>0.75</td><td>12.83</td><td>12.69</td><td>1.09</td></tr></tbody></table>

As shown in Table 1, the values of flexural-torsional buckling loads and natural frequencies derived from this method are in excellent agreement with those obtained by using Ansys code.

## 6.2. Natural frequencies and flexural-torsional buckling of cantilever thin-walled beam-columns

In this example, the stability and free vibration analyses of a non-prismatic cantilever thin-walled beam are investigated. The analyzed beam has two different cross-sections. The first one is a doubly I-section and the second is a singly symmetric I section, with the top flange shorter than the bottom flange. The web height (the distance between the flange mid-lines) of the cross-section is made to vary linearly from 600 mm at the clamped end to 300 mm at the free end. The beam also exhibits a linear flange tapering, as shown in Fig. 3. In the case of buckling analysis, a compressive load is applied at the top flange of the free end section (point B). Buckling analysis is performed for different beam slenderness $L=12m$. The modulus of elasticity is assumed $210GPa$.

<figure><img src="image_4.png" /><figcaption>Fig. 3 Tapered I-section cantilever thin-walled beam under axial compressive load: geometry and material properties.</figcaption></figure>
Table 2 presents the variation of the flexural-torsional buckling loads of the cantilever column with non-uniform cross-section computed by the proposed procedure and compared to those obtained by finite element method using Ansys software.

<table><caption>Table 2: Critical buckling load $P_{cr}(kN)$ for the non-prismatic beams. (Example 2)</caption><thead><tr><th rowspan="3">L(m)</th><th rowspan="3">Mode</th><th colspan="6">FTB of cantilever beam (KN), load on top flange</th></tr><tr><th colspan="3">Equal Flanges</th><th colspan="3">Unequal Flanges</th></tr><tr><th>FEM<br>(Ansys 5.4)</th><th>Present<br>Method</th><th>Δ(%)</th><th>FEM<br>(Ansys 5.4)</th><th>Present<br>Method</th><th>Δ(%)</th></tr></thead><tbody><tr><td rowspan="2">12</td><td>n=1</td><td>27.562</td><td>27.494</td><td>0.25</td><td>16.32</td><td>16.38</td><td>0.37</td></tr><tr><td>n=2</td><td>188.02</td><td>188.032</td><td>0.01</td><td>129.04</td><td>123.66</td><td>4.17</td></tr></tbody></table>

## 7. CONCLUSIONS

The total potential energy principle considering the effects of strain and initial stress, the kinematic energy and the work of the applied loads is used to derive the equilibrium equations of motion of thin-walled beams. Then, power series expansions method is used to derive the exact shape functions for the element with 12 degrees of freedom. Then, the stiffness and mass matrices are derived. The proposed method can be applied in various form of non-prismatic member under axially concentrated loads and variable mass per unit length. Furthermore; it can be used to evaluate both natural frequency and buckling load concurrently.

## REFERENCES

* Al-Sadder S.Z. (2004), "Exact expression for stability functions of a general non-prismatic beam-column member". Journal of Constructional Steel Research, Vol. 60, 1561-84.

* Andrade A, Camotim D, Borges Dinis P. (2007), "Lateral-torsional buckling of singly symmetric web-tapered thin-walled I-beams, 1D model vs. shell FEA". Computers and Structures, Vol. 85, 1343-1359.

* Asgarian B, Soltani M, Mohri F. (2013), "Lateral-torsional buckling of tapered thin-walled beams with arbitrary cross-sections". Thin-Walled Structures, Vol. 62, 96-108.

* Chen C-N. (2000), "Dynamic equilibrium of non-prismatic beams defined on an arbitrarily selected co-ordinate system". Journal of Sound and vibration, Vol. 230(2), 241-260.

* Eisenberger M. (1994), "Vibration frequencies for beams on variable one- and two-parameters elastic foundation". Journal of sound and vibration, Vol. 176(5), 577-584.

* Eisenberger M and Cohen R. (1995), "Flexural-torsional buckling of variable and open cross-section members". Journal of Engineering Mechanics, ASCE, Vol. 121(2), 244-254.

* Kim S-B, Kim M-Y. (2000), "Improved formulation for spatial stability and free vibration of thin-walled tapered beams and space frames". Engineering Structures, Vol. 22,
446–458.

Ronagh HR, Bradford MA and Attard MM. (2000), “Non-linear analysis of thin-walled members of variable cross-section. Part I: Theory”. *Computers and Structures*, Vol. 77, 285-299.

Ronagh HR, Bradford MA and Attard MM. (2000), “Non-linear analysis of thin-walled members of variable cross-section. Part II: Application”. *Computers and Structures*, Vol. 77, 301-313.

Yang YB, Yau JD. (1987), “Stability of beams with tapered I-sections”. *Journal of Engineering Mechanics, ASCE*, Vol. **113**(9), 1337–1357.

Timoshenko SP, Gere JM. (1961), *Theory of elastic stability*, 2nd ed. New York: McGraw-Hill.

Vlasov VZ. (1959), *Thin-walled elastic beams*, Moscow, French translation, Pièces longues en voiles minces, Eyrolles, Paris, 1962.

Chen WF, Lui EM. (1987), *Structural stability, theory and implementation*, New York, Elsevier.

Bazant ZP, Cedolin L. (1991), *Stability of structures Elastic, inelastic fracture and damage theories*, New York, Dover Publications.

MATLAB Version 7.6 .MathWorks Inc, USA, 2008.

ANSYS, Version 5.4, Swanson Analysis System, Inc, 2007.