## Abstract

The state vector equation for lateral buckling in finite displacement theory is formulated using only the hypothesis of the Bernoulli–Euler beam. By using an appropriate orthogonalization of the warping functions, the normally complicated calculation has been processed systematically using only matrix notation. As a numerical analysis, the lateral buckling load on the cantilever receiving a concentrated end load on the upper flange was calculated using the coefficient matrix of the first-order increment; the post-buckling behaviour was investigated with increasing load. Since the state vector equation is a higher order nonlinear equation, the original coefficient matrix was fixed with an arbitrary initial value and the solution was provided by the Runge–Kutta transfer matrix method. Subsequent calculations were pursued in the same way with the solution obtained via Runge–Kutta methods as a new initial value and then shifted to the next load condition. This theory and analysis method does not employ an assumed displacement function, such as the Ritz's method; it is therefore useful for the finite displacement analysis of a beam with arbitrary boundary conditions and intermediate support conditions.

## 1. Introduction

The Bernoulli–Euler beam (hereafter, the B–E beam) is the basic and simplest model of beam bending. The micro- and nano-scaled beams, which have been the subject of much recent attention, must be developed from the B–E beam model (e.g. Park & Gao 2006) to obtain a higher order beam model. From this perspective, the lateral torsional buckling phenomena of the B–E beam will be taken up in this paper and a new finite displacement theory proposed.

As stated in the book of Timoshenko & Gere (1961), lateral torsional buckling phenomena are analysed with an assumed deformation function, such as trigonometric functions. As sinusoidal and/or cosine functions of rotation in the governing equation are developed in the Taylor series and approximated by only a few terms (e.g. Mohri *et al*. 2002), this treatment cannot trace the lateral torsional post-buckling behaviour accurately. Additionally, the effects of load height are practically important in lateral torsional buckling (e.g. Trahair 1993).

In this paper, based on the Bernoulli–Euler hypothesis, the state vector equation of the finite displacement theory was derived, beginning with the finite strain condition without any mathematical approximations. To be specific, the phenomena of lateral buckling of a beam, such as the elastica theory for columns, were described mathematically without carrying out the approximation that expands trigonometric functions by the Taylor series. To the author's knowledge, this is, to date, the first rigorous research to formulate state vector equations.

The following points were executed in the numerical calculation. The buckling load was determined by using the coefficient matrix that was derived from the first-order increment of the nonlinear state vector equation obtained. The load was further increased and then the post-buckling behaviour was observed. Owing to the higher order nonlinearity, the coefficient matrix of the original state vector equation includes unknown state quantity, which means that it cannot be used as is. To account for this, the coefficient matrix was fixed by substituting an appropriate initial value per load value, and the solution was obtained from the Runge–Kutta modified transfer matrix method (Friemann 1974). The calculation was similarly pursued with the obtained solution as the new initial value until the solution was obtained.

As an example of the numerical calculation, problems such as concentrated loads acting on the ends of thin-walled I-shaped (open section) and box-shaped (closed section) cantilevers were solved and it was shown how the state quantities change along with the increase of the load. The relationship between the acting height within the section and the buckling load with respect to the concentrated end load was specifically focused on.

Similar to the Rayleigh–Ritz method, this study is performed without the displacement approximation so that it is applicable to beams under all kinds of boundary conditions and intermediate support conditions.

## 2. The coordinate system and assumptions

In our coordinate system, *x* is the coordinate that is in the member axial direction passing through the centroid of the cross-section S, and *y* and *z* are the principal axes perpendicular to *x*. Moreover, in addition to the principal axes described above, in the cross-section perpendicular to the *x*-axis, *s* as the coordinate of the contour line that passes through the centre line of the thickness of each plate and *n* as the coordinate in the direction of thickness are assigned as shown in figure 1. The displacement components in the directions of *x*, *y*, *z*, *s* and *n* at a given point P of the structural element are *u*, *v*, *w*, *f*_{s} and *f*_{n}, respectively. However, the axial extension displacement of the cross-sectional centroid S is *u*_{0}.

The following assumptions are used for this study.

The cross-section that is made of plates is either open or closed, where the centre line of the thickness is linear.

The plate thickness

*t*inherent to the structural system is constant, and it is small in comparison with its width*b*.The member length is large in comparison with the size of the cross-section, and the derivative of the displacement with respect to the axial coordinate

*x*is small.No cross-sectional deformation.

Normal strain in the direction of the contour coordinate of the cross-section is zero:

*ϵ*_{s}=0Normal strain in the normal direction of the contour coordinate is zero:

*ϵ*_{n}=0Shear strain in the coordinate within the cross-section is zero:

*γ*_{sn}=0

Shear strain of the plate bending is zero:

*γ*_{nx}=0Secondary shear strain

*γ*_{xs,se}attributed to the secondary shear stress balanced with the normal stress*σ*_{x}is neglected:*γ*_{xs,se}=0Structural material is satisfied by Hooke's law.

The superscript asterisk in figure 1 implies the magnitude on the middle surface of the thickness.

## 3. Derivation of governing equation

### (a) Equations for the basic condition

The mathematical expression is simplified by combining the normal and tensor notations. Space coordinates, *x*, *s* and *n*, are expressed as *x*_{1}, *x*_{2} and *x*_{3}, respectively, i.e. *x*_{i}. Similarly, displacement components, *u*, *f*_{s} and *f*_{n} are *u*_{1}, *u*_{2} and *u*_{3}, respectively, i.e. *u*_{i}. Subscripts in Latin characters, *i*, *j* and *k* are expressed as 1, 2 and 3, respectively. According to the various rules of tensor mathematics, the general relationship between the strain tensor *ϵ*_{ij} and the displacement *u*_{i} of the finite displacement theory can be expressed as(3.1)From assumption (c), if the second-order derivation of the longitudinal displacement *u*(≡*u*_{1}) is ignored, by changing the dummy index of the Latin character of *k* to the Greek character of *γ*, the above equation can be written as(3.2)where the subscript *γ* can only be 2 or 3. If *n*=0 in this equation, then the strain displacement equation (3.3) on the centre line of the thickness is given as(3.3)From assumptions (d) and (e) *τ*_{nx} and *τ*_{sn} become zero. The normal stress *σ*_{n} in the thickness direction also becomes approximately zero, which means that the structural members are in the plane stress state. Given that the longitudinal elastic modulus is *E* and Poisson ratio is *ν*, the stress–strain relationship on the plane stress is(3.4)where, in this equation, the Greek characters *α*, *β* and *ρ* can only be 1 or 2.

Given that components in *x*- and *s*-directions of the body force are expressed as *p*_{x} and *p*_{s}, respectively, the stress equilibrium condition of the small elements of the plates composing the beam can be given as(3.5)where, in this equation, the Greek characters, *α* and *β* can only be 1 or 2. The shear stress *τ*_{xs} that forms the stress resultants of the warping torsion theory is composed of the primary component *τ*_{xs,pr} and the secondary component *τ*_{xs,se}. Then, the shear strain *γ*_{xs} that forms the strain resultant is composed of the primary and secondary components, *γ*_{xs,pr} and *γ*_{xs,se}, respectively.(3.6)where both sides of the equation (3.6)_{1} are divided by the shear modulus *G* and Hooke's law is applied. Since the secondary component is the only one equilibrated with the normal stress *σ*_{x}, the equation for the *x*-axis direction of the equilibrium condition (3.5) can be resolved as(3.7)If the secondary shearing strain *γ*_{xs,se} of the r.h.s. of the equation (3.6)_{2} is ignored as infinitesimal according to the assumption (f), *γ*_{xs}=*γ*_{xs,pr} is given. By adapting Hooke's law , the following equation is given:(3.8)where the primary shearing strain of the r.h.s. of the equation where the stress equation (3.7)_{1} was integrated by *s*. However, *T*(*x*) becomes Bredt shear flow that takes a constant value per plate composing the cross-section. The shear strain yields statically indeterminate shear flow per closed part of the cross-section. This is determined by the continuity condition of the warping displacement in the *x*-axis direction,(3.9)Here, *M* is the total number of cells enclosed by the midline of the section wall. Generally, this function is determined by solving a simultaneous equation regarding the unknown unit shear flow only with the total number of cells.

### (b) In-plane displacement functions of the cross-section

By expressing assumption (d) for the invariance of the cross-sectional shape as the strain expression, the following equation is given as:(3.10)By substituting the strain-displacement equation (3.2) in these equations, the in-plane displacement function is set as(3.11)The symbols (^{.}) and ( )′ imply differentiations with respect to the contour coordinate *s* and the axial coordinate *x*, respectively. In summary, the above equations become(3.12)Here, the in-plane function vectors *r*^{*}(*s*) and *b*^{*}(*s*) of the contour coordinate *s* of the cross-section are known functions and the displacement function vector ** V**(

*x*) of the axial coordinate

*x*is an unknown function. The content of those vector functions can be listed as follows:(3.13)where

*α*

^{*}is an angle that is formed by the principal axis

*y*of the cross-section and the contour coordinate

*s*, and and are, respectively, the vertical distance from the shear centre to the focus point on the middle surface of the thickness and the distance in the

*s*-direction. The beginning of the contour coordinate

*s*is the origin S that is the centre of the cross-section, and it is assumed that it is connected with an actual plane through a virtual plane of the thickness

*t*=0 (Nishino

*et al.*1974). Given that

*b*

_{T}is a value of the

*s*-coordinate, which is from the shear centre M to the foot of a perpendicular pulled down to the plate to which the focus point is related, the function can be expressed as (figure 2

*a*)(3.14)While and

*b*

_{T}are the constants indicating the point of the shear centre, is a linear function of contour coordinate

*s*that indicates the focus point. The elements

*ξ*,

*η*,

*ζ*and

*φ*of the vector function

**(**

*V**x*) that defines the in-plane displacement represent the integral quantity of the expansion displacement in the

*x*-axis direction, deflections in the

*y*- and

*z*-axes directions and the rigid rotation angle of the cross-section, respectively. The first element

*ξ*, which is the vector of the displacement function

**(**

*V**x*), is defined for descriptive purposes in order to equalize the number of vector elements, but it is eventually neglected due to its insignificant quantity.

### (c) Displacement function in the axial direction

Assumption (e), which is the shear strain of plate bending, is expressed as(3.15)By substituting the strain-displacement equation (3.2) and integrating it, the axial displacement function is given as(3.16)where the prime symbol at the upper right ( )′ implies differentiation with respect to *x*. Also, the symbol ** ϕ**(

*x*) is a matrix that gives the bidirectional in-plane displacement that is yielded in accordance with the distance

*n*from the centre of the thickness due to the cross-sectional rotation

*φ*, and it is defined as(3.17)where the symbol

**0**means that blank spaces are filled with all zero elements. The first term on the r.h.s. of equation (3.16) is the axial displacement on the centre line of the thickness. This is determined by the definition of the shear strain of the middle surface of the thickness, which is reflected in the strain-displacement equation (3.3),(3.18)If the secondary shear strain is ignored, since , the first term of the r.h.s. of the preceding equation becomes(3.19)where the primary shear strain of equation (3.8) is substituted for the shear strain at the l.h.s. of the preceding equation. However, Bredt shear flow, which has a constant value to each plate composing the cross-section, is(3.20)If equation (3.19) is integrated by

*s*and the axial expansion displacement that is at the centre of the cross-section is written as the constant of integration ‘−

*ξ*′(

*x*)’, the following equation is given as:(3.21)where a unit warping function and a unit primary shear flow function

*ψ*^{*}(

*s*) are defined as(3.22)The vector

*ψ*^{*}(

*s*) takes a constant value per plate, and it is determined by the continuity condition (3.9) of warping per closed cell

*m*. The equation (3.22)

_{2}is valid only for a single-cell cross-section.

In accordance with the technique of the virtual plate (Nishino *et al.* 1974), which connects the centre of the cross-section and the cross-section composing plates, the axial displacement is defined to be consistent here as well. If equation (3.21) at the centre line of the thickness is substituted into equation (3.16), the axial displacement of an arbitrary point can be expressed as(3.23)When the elements of *r*^{*}(*s*) are and , then the contour integration of the r.h.s. numerator of equation (3.22)_{2} becomes zero. Likewise, it can be easily expressed that the contour integration of the element becomes zero. From the above, the elements of the vectors *ψ*^{*}(*s*), and can be expressed as(3.24)where two types of the in-plane function of the torsion become(3.25)When the axial displacement equation (3.23)_{1} is expanded, the function of *x* that connects with the second in-plane function for torsion, , becomes zero. The following equation is the reason for this:(3.26)Therefore, this term does not appear in any traditional books or research papers. However, in order to orthogonalize effectively without changing the orders of these vectors, it shall be left here as a dummy placeholder. Later, the in-plane function will be changed to due to orthogonalization; but for reasons described above, the content of the axial displacement equation (3.23)_{1} becomes invariant and does not lose its consistency.

### (d) Strain

The displacement component is determined according to the condition that four of the six types of strain *ϵ*_{s}, *ϵ*_{n}, *γ*_{sn} and *γ*_{nx} become zero. The remaining strain components which do not become zero are *ε*_{x} and *γ*_{xs}.

Substituting equation (3.12)_{1} and (3.23)_{1} into equation (3.2), the normal strain *ϵ*_{x} is expressed as(3.27)The unit warping function ** w**(

*s*,

*n*) and the matrix

**′(**

*Ω**x*) are(3.28)(3.29)The coefficients

*c*

_{1}to

*c*

_{wT}are defined as the warping function to uncouple with others (see Usuki 1998).

Then, if only equation (3.28)_{1} is used for the unit warping function, the equation for the axial displacement *u*(*x*, *s*, *n*) of equation (3.23)_{1} can be rewritten as(3.30)

The displacement vector ** v**(

*x*) and the rotation vector

**(**

*ϑ**x*) of the beam, which do not couple the bending displacements and the torsional rotation, are denoted by the component expressions as(3.31)Next, by differentiating (3.13)

_{5}the derivative of the displacement vector

**(**

*V**x*) coupling with the large rotation

*φ*within the cross-section becomes(3.32)Therefore, the axial rotation vector

**(**

*θ**x*) of the cross-section that takes into account the large rotation within the cross-section

*φ*is defined as(3.33)Using the vector notation, these can be expressed as(3.34)where

*(*

**φ***x*) is defined as(3.35)

As it is generally known, the axial rotation vector ** θ**(

*x*) of the B–E beam is equal to the deflection angle

**′(**

*V**x*), as shown as(3.36)Therefore, equation (3.30) for the axial displacement

*u*(

*x*,

*s*,

*n*) can be also rewritten using the axial rotation angle

**(**

*θ**x*) as(3.37)

Along with the adoption of new vector symbol ** θ**(

*x*), equation (3.27) for the axial strain

*ϵ*

_{x}can be expressed as(3.38)As a result, the axial displacement

*u*(

*x*,

*s*,

*n*) and the axial strain

*ϵ*

_{x}(

*x*,

*s*,

*n*) are expressed as vectors based on the group of unit warping functions that are mutually perpendicular to each other.

By substituting equation (3.12)_{1} and (3.37) into equation (3.2),(3.39)is given, where the primary shear flow function at an arbitrary point of the cross-section is defined as(3.40)

### (e) Expression of the stress with the displacement quantities

The normal strain *ϵ*_{x} and the shear strain *γ*_{xs} are the two remaining strain components. By using them, the stress–strain equation (3.4) can be specified as(3.41)where is a modified modulus of elasticity, in which *ν* is the Poisson ratio, and *G* is the shear modulus. By substituting the strain equation (3.38) and (3.39) into the constitutive equations, they are shown as follows:(3.42)

(3.43)

Next, the secondary shear stress *τ*_{xs,se} is obtained, which is equilibrated with the normal stress *σ*_{x}. By substituting and transposing equation (3.42) into equation (3.7)_{2} and integrating it, the following equation is given:(3.44)However, the distributed load is omitted. Two types of vector functions that compose the unit secondary shear flow function (** S**(

*s*,

*n*)−

**) of the preceding equation are defined as(3.45)Similarly to the beam theory of the infinitesimal strain, the secondary statically indeterminate shear flow**

*Φ***is also determined, using the continuity condition (3.9) for the axial warping displacement. From equation (3.6)**

*Φ*_{1}, the total shear stress

*τ*

_{xs}can be given by the sum of equations (3.43) and (3.44)(3.46)The orthogonality between the primary shear flow function

**(**

*ψ**s*,

*n*) and the derived function (

*s*,

*n*) by the contour coordinate

*s*of the unit warping function, and the orthogonality between the function

**(**

*ψ**s*,

*n*) and the secondary shear flow function (

**−**

*S***) are formed as(3.47)which is the same as the infinitesimal displacement theory (Usuki & Sawada 1999).**

*Φ*### (f) The definition of the stress resultants

The stress resultant vectors of the warping moment *M*_{w}(*x*), the primary shear force *Q*_{pr}(*x*), the secondary shear force *Q*_{se}(*x*) and the total shear force ** Q**(

*x*), are defined as the cross-sectional integral that spans the entire cross-section

*A*of the beam members and written as follows:(3.48)The sum of the primary and the secondary shear forces is the total shear force. It can be expressed as the following equation:(3.49)To express these vectors in elements, they become(3.50)The first four elements in equation (3.50)

_{1}are the axial force

*N*

_{x}, the bending moment

*M*

_{z}about the

*z*-axis, the bending moment

*M*

_{y}about the

*y*-axis and the warping torsion moment

*M*

_{wT}, respectively, with respect to the shear centre. The secondary warping torsion moment

*M*

_{wt}of the fifth element is a new shear force that is yielded by using the Green–Lagrange strain. The fourth element of the primary shear force equation (3.50)

_{2}is the traditional primary torsion moment

*M*

_{T,pr}. Since the fifth element

*ψ*^{*}(

*s*) of equation (3.24)

_{1}is zero, the fifth element of the primary shear force equation (3.50)

_{2}also becomes zero. The shear forces

*Q*

_{y}and

*Q*

_{z}in the

*y*- and

*z*-directions of the bending equilibrated with normal stress, and the twisting moment

*M*

_{T,se}of the warping restraint, are attributed to the secondary shear force

*Q*_{se}of equation (3.50)

_{3}. The new warping restraint moment that is equilibrated with the secondary warping stress of torsion is denoted by

*M*

_{t}and placed at the fifth element. As above, the total shear force is given by equation (3.50)

_{4}as the sum of the primary and the secondary components. As generally known, the torque moment is .

### (g) Expression of the stress resultants with displacement components

By substituting the stress–displacement expressions (3.42) and (3.46) into the definition of stress resultants (3.48),(3.51)However, the matrices of the warping resistance ** F** and the torsion resistance

**of the cross-section are defined as(3.52)(3.53)Since the unit warping function**

*J***(**

*w**s*,

*n*) and the unit primary shear flow function

**(**

*ψ**s*,

*n*) are formed by a group of mutually linearly independent functions,

**and**

*F***become diagonal matrices and no coupling between deformations occurs.**

*J*By substituting the displacement vector of the r.h.s. of equation (3.51)_{1–3} into equation (3.42)–(3.44)*, the* stresses are given with the stress resultants as follows:(3.54)However, the torsion resistance ** J** is a singular matrix that has zero elements to the diagonal; therefore, the inverse of the matrix is used, for which only the element

*J*

_{T}is inverted.

### (h) Reissner variational principle

In order to derive the governing equation of the beam, the energy principle that uses the Reissner functional *ψ* is used (Vinson & Chou 1975). When the prescribed boundary conditions, equilibrium equations and the stress–displacement relationship are satisfied, the Reissner's functional(3.55)becomes a minimum. In other words, using the variational symbol *δ*, it can be expressed as(3.56)where the subscript *i* uses values 1–3 corresponding to *x*, *y* and *z*; *R* is the volume of the elastic body; *S*_{t} is the portion of the body surface over which tractions are prescribed; *F*_{i} is the *i*th component of a body force; *T*_{i} is the *i*th component of the surface traction; *u*_{i} is the *i*th component of the deformation; ; and *W*(*σ*_{ij}) is the strain energy function of the stress expression.

If the material is linearly elastic and isotropic, the strain energy function of the stress expression is represented as(3.57)

Neglecting the secondary shear deformation and the cross-sectional distortion, the Reissner functional is expressed as(3.58)where the beam with length *L* is receiving action from the distributed lateral load to the contour direction of the cross-section, *p*_{s}(*x*); the distributed lateral load to the normal line direction, *p*_{n}(*x*); the simultaneously concentrated lateral loads, *P*_{s} and *P*_{n}; the axially distributed load, *p*_{x}(*x*, *s*); and the axially concentrated load, *P*_{x}. By using the definitions of the stress resultants (3.48)_{1,3,4} and the warping and torsion resistance (3.52) and (3.53), the Reissner's functional *ψ* can be organized as(3.59)

It is necessary to make a few preparations before pursuing the variational calculation. When the first variation of the displacement function ** V**(

*x*),

**′(**

*V**x*) and

**(**

*ϕ**x*) are calculated, they can be expressed as(3.60)The new matrix is given as(3.61)The scalar variation quantity

*δφ*of the rotation angle can be written as(3.62)Also, in order to switch from the product of the variation of the matrix

**′(**

*Ω**x*) of equation (3.29) and the vector

**′(**

*V**x*) to the variation of

**′(**

*V**x*), the matrix is defined as(3.63)Then, matrix becomes equation (3.64)(3.64)

By substituting the Reissner's functional (3.59) into the variational principle (3.56), the relational equation of the stress resultants is given as(3.65)where *Γ* is the contour of the entire cross-section. The shear force *Q* of the preceding equation is used and defined as(3.66)If the shear force

**of the l.h.s. of the equation is transposed and equation (3.49) is used, the known relational expression of the warping moment equilibrium is written as(3.67)Now the distributed load , which is in the**

*Q**s*-direction of the contour and acting on the position coordinate of a given plate composing a cross-section, and the distributed load , which is in the

*n*-direction of the cross-sectional thickness, are united into the distributed load and used as a resultant force. Likewise, the concentrated lateral load is a resultant force component of the cross-section. Based on the above, the following is given:(3.68)The length of the perpendicular foot drawn from the shear centre M to the plate is

*r*

_{T}, and the position of the foot is

*b*

_{T}(

*s*)=0 of the plate. Therefore, if

*β*is the angle formed by the

*s*-coordinate of the plate's contour and the resultant force of the lateral load , it can be expressed as(3.69)Also, if the angle formed by the lateral load and the

*y*-axis of the cross-section is expressed as as shown in figure 3

*a*, then the second and third elements of vector become(3.70)They provide the components of the

*y*- and

*z*-axes of the lateral load

*p*

_{Q}(

*x*,

*s*). Likewise, the fourth and fifth elements of the vector become(3.71)where is the length of the perpendicular line drawn from the shear centre M to the lateral load , and is the distance from the shear centre M to the load acting point , which is parallel to the contour coordinate

*s*of the plate as shown in figure 3

*b*. The fifth element governs the influence on the buckling value of the lateral load where the point of action is in either the lower or upper side from the shear centre. Hence, the element components subjected to the product of the distributed lateral load and the vector are formulated as well as the product of the concentrated lateral load and the vector .

The boundary conditions are(3.72)where a line shown above the symbols implies the load term and are defined as(3.73)While the stress resultants and ** Q**(

*x*) are defined with respect to the cross-section that is perpendicular to the deformed member axis, the newly introduced warping moment and shear force are in the same direction with the vectors of the displacement components

**(**

*ϑ**x*) and

**(**

*v**x*), respectively. In other words, these are the bending moment and shear force of the finite displacement theory. Those vectors can be expressed as elements as follows:(3.74)However, in order to equalize the number of elements of a vector, a component that becomes zero is defined as the shear force vector for descriptive purposes but it does not have actual existence, and therefore has to be omitted from the actual calculation.

By substituting equation (3.51)_{1} into the r.h.s. of the new bending moment equation (3.73)_{1} and then substituting equation (3.34) into the resultant, those vectors are expressed as the displacement components such as vectors ** ϑ**(

*x*) of the cross-sectional rotation. Similarly, by substituting equation (3.51)

_{4}into the r.h.s. of the new shear force equation (3.73)

_{3}and then substituting equation (3.34) into the resultant, those vectors are expressed as the displacement components such as

**(**

*ϑ**x*) and

**(**

*v**x*), which means they become(3.75)The matrix used here is defined asThe coefficient matrices

**and**

*B***arewhereAlso, the components of the axially concentrated load**

*C**P*

_{x}and the laterally concentrated load

*P*

_{Q}, which are equilibrated with the stress resultants, are as follows:whereHere, of the external force vector equilibrated with the shear force is the angle formed by the laterally concentrated lateral load

*P*

_{Q}and the

*y*-axis. Therefore, the second element of the vector is the

*y*-axis component of the laterally concentrated lateral load

*P*

_{Q}, and the third element of the vector is the

*z*-axis component of the laterally concentrated lateral load

*P*

_{Q}. Also, the length of the perpendicular line drawn from the shear centre M of the beam cross-section to the concentrated lateral load is expressed as , and the distance from the shear centre M to the applied point measured parallel to the direction of the concentrated lateral load is expressed as . Using the in-plane functions and , the first half of the vector's fourth element represents the torsion effect of the lateral load

*P*

_{Q}that does not act directly on the shear centre.

Using the new stress resultants, the preceding shear force equilibrium condition (3.65) can be rewritten as(3.76)and the derived coefficient matrix is given by the following equation:whereAndwhere

### (i) The system of ordinary differential equations for state quantity

In order to facilitate easier prescription of the given boundary conditions, vectors such as displacement ** v**(

*x*), cross-sectional rotation

**(**

*ϑ**x*), warping moment and lateral force are employed as elements of the state quantity vector

**(**

*z**x*). In other words, they are defined as(3.77)The system of ordinary differential equations involving the preceding relational expression (3.36)

_{2}for the deflection and the cross-sectional rotation, the expression (3.75)

_{1}for the displacement quantity of warping moment, the relational expression (3.75)

_{2}for the lateral force and warping moment, and the equilibrium equation (3.76) for the lateral forces are given as(3.78)Or in matrix presentation, we have(3.79)Here,

**is an identity matrix. Among them, the equation, which is the first row of the first block, and the equation , which is the first row of the fourth block, are meaningless; therefore, they will be omitted from the actual calculation.**

*I*### (j) The state vector equation for the first-order increment

The state vector equation results in the following matrix expression if the first-order increment of the system of ordinary differential equations (3.78) is calculated,(3.80)The symbol (^{.}) refers to the first-order increment, which means that four types of state quantity elements , , , and the distributed loads, and , are the components of the first-order increment. The matrices used for the preceding equation can be specifically expressed aswhereBased on the reason mentioned at the end of (*i*), the omitted zero elements are indicated by the symbol .

According to the coefficient matrix of the r.h.s. of the obtained state vector equation (3.80) for the first-order increment, the blocks on each side of the cross-diagonal are symmetric. Also, the off-diagonal matrix with respect to the cross-diagonal appeared as a flagstone-like symmetry (superimposed property). In other words, each alternating symmetric block interchanges with each other (Pestel & Leckie 1963). The symmetric block subjected to alternation is a zero block in the infinitesimal displacement theory; therefore, this block is governing the finite displacement phenomena. Based on this coefficient matrix of the infinitesimal displacement theory, the cross-diagonally symmetric matrix is obtained when a field transfer matrix is created with Taylor expansion or Runge–Kutta formula. And, the cross-diagonally symmetric matrix reduces to a main diagonally symmetric stiffness matrix when the transfer matrix method is transformed to the displacement method. Then, the symmetric property of the principal diagonal of the stiffness matrix of the infinitesimal displacement theory is not guaranteed in the finite displacement theory owing to the alternatively symmetric blocks of the coefficient matrix (Simo & Vu-Quoc 1986).

### (k) The estimation method for buckling load

By taking the first-order increment of the state vector equation (3.79) for the finite displacement, the state vector equation (3.80) was obtained. In the equation, since there is no variation in terms of the load before and after buckling, the expression can be set as . The rotation angle *φ* is zero just before the cross-sectional rotation due to buckling, which means that if the rotation limit was *φ*→0 in equation (3.73)_{1,3}, the matrix becomes an identity matrix and the equation is given as

As for the obtained coefficient matrix of the r.h.s. of the state vector equation (3.80) for the first-order increment, under the limit of *φ*→0, each composing element is simplified.

The buckling condition matrix can be given by using the boundary conditions on the left and right ends, which is given by a transfer calculation with a field transfer matrix created in the coefficient matrix of the r.h.s. of the state vector equation for the first-order increment in accordance with the Runge–Kutta method (Schäfer 1970). This buckling condition becomes an eigenvalue equation, and the buckling load can be solved by using the smallest value as the determinant of the matrix becomes zero.

### (l) Analysis methods for post-buckling behaviour

The system of ordinary differential equation (3.79) is used. First, the applied load is set and it is substituted into the coefficient matrix. However, the coefficient matrix is an unknown matrix because the state quantity vector is involved as a function in the coefficient matrix. For that reason, an appropriate initial state quantity is assumed and substituted into the coefficient matrix. By solving the equation (3.79) of which the coefficient matrix becomes a known matrix with the Runge–Kutta method, the field transfer matrix can be derived. Second, the new state quantity can be given by computing the transfer calculation by the transfer matrix method and using the boundary conditions of the left and right ends. If the difference between the obtained and assumed state quantities is infinitesimal, those values satisfy the governing equation. If the gap is not infinitesimal, by substituting the obtained state quantity again into the coefficient matrix and repeating the same calculation, a new state quantity can be given. The state quantity of the post-buckling condition is calculated by the same calculation with a different load value.

## 4. Numerical examples

### (a) The buckling load of the cantilever receiving the concentrated end load at the acting height e

The buckling load on the cantilever beam of the cross-section that is biaxially symmetrical (*y _{M}*=0) is obtained, where the concentrated end load is received at the acting height

*e*(figure 4). The length of the cantilever beam is

*l*and the web height is

*h*. By using the web height

*h*as a reference length, the size (dimension) of the beam is determined. The acting height

*e*from the shear centre S is taken in the same direction as the

*z*-axis, and the lower direction below the

*x*-axis is positive. The cross-sections employed are I-shaped and narrow-box-shaped. In order to have the same width, height and area, the dimension of the cross-sections are set asFigure 5 shows the variation of the beam slenderness ratio

*λ*(≡

*l*/

*r*) and the buckling load

*P*

_{cr}with respect to two kinds of acting height,

*e*=0 (load acting at the shear centre S) and

*e*=−

*h*/2 (load acting at the upper flange). Here,

*r*indicates the radius of gyration. According to figure 5, the influence of the loading position becomes less as the slenderness ratio

*λ*increases. Also, the I-shaped cross-section provides a larger influence on the buckling load on the loading position than the box-shaped cross-section.

### (b) Post-buckling behaviour of the cantilever beam receiving the concentrated end load at the acting height e

The post-buckling behaviour of the cantilever beam (figure 4) whose cross-section (*y _{M}*=0) is biaxially symmetrical was obtained, where the concentrated end load is received at the acting height

*e*. The dimension of the cross-section is the same as the preceding one and the length of the member is also

*l*=5

*h*. Figures 6 and 7 show the relationship of the acting load

*P*and the displacement of the end of the cross-section. The vertical axis indicates the ratio of acting load

*P*to buckling load

*P*

_{cr}, and the horizontal axis is the displacement amount of the beam end. Figure 6

*a*shows

*φ*and 6

*b*, , which are the rotation and the rotation rate of the beam end of the I-shaped cross-section, respectively. Likewise, figure 7

*a*shows

*φ*and

*b*, , which are the rotation and the rotation rate of the beam end of the box-shaped cross-section, respectively.

Figure 8 shows the displacement of the cantilever beam end along with the load *P* increased after buckling. Figure 8*a* shows the dimensionless expression for the displacement of the I-shaped beam being divided by the web height *h*. To distinguish the difference of the two kinds of loading position, the cross-section images are shown with dashed and solid lines. The numbers in the figures indicate the magnification ratio with respect to the buckling load. Likewise, figure 8*b* shows those of the box-shaped cross-section. Next, figure 9 shows the three-dimensional view of the displacement when the acting load at the end of the I-shaped (figure 9*a*) and box-shaped (figure 9*b*) cross-sections is 1.2 and 1.1 times the buckling load, respectively. In order to clarify the variance of two loading positions, two different images are displayed together for both cross-sections.

## 5. Conclusion

The obtained governing state vector equation has been expressed with trigonometric functions whose argument is the cross-sectional rotation. Since this formulation was completed without any approximation such as Taylor expansion on the trigonometric function, it is therefore the elastica theory for the bending torsion phenomena. In order to confirm the validity of this theoretical equation, the first-order increment was derived from the obtained state vector equation, and the coefficient matrix of the increment state vector equation has indicated a perfect flagstone-like superimposed property to the cross-diagonal. In the linear infinitesimal displacement theory, there was a cross-diagonally symmetric property where a block with zero and non-zero values became a flagstone-like superimposed property. However, in this finite displacement theory, it has been confirmed that an alternative symmetric matrix appears at the location where the block is zero. This property was applied to the exponential function for the coefficient matrix of the state vector equation, or into the transfer matrix derived by the Runge–Kutta formula. Then, the symmetric property of the principal diagonal of the stiffness matrix of the beam, which can be obtained by converting the transfer matrix, is a mathematical reason for why the finite displacement theory breaks down.

The author validated previously that the buckling loads (moments) obtained are almost the same as the classical solution of Timoshenko & Gere (1961). Little difference arises from the fact that the present theory considers deformation before buckling (see Usuki 1998, pp. 262–263).

In order to numerically reconfirm the mathematical validity of this approach, a lateral torsional buckling problem was solved. First, a problem of buckling load determination and second, post-buckling behaviour were investigated. From the theoretically derived first-order increment of the state vector equation, the state vector equation for the increment has been obtained. But the buckling load was calculated from the transfer matrix of the coefficient matrix with respect to the incremental state quantity, and it is a highly accurate solution of the traditional theory. To be specific, the deformation that occurs immediately before buckling was accurately captured and therefore the influence that is given to the buckling load at the loading position can be calculated with high accuracy. Also, this solution method can be adapted to all kinds of boundary conditions and intermediate support conditions.

As for the analysis of the post-buckling behaviour, I-shaped (open section) and box-shaped (closed section) beams were used and the deformation behaviour of post-buckling was observed by comparing the influence of loading position. As a result, it was confirmed that the ratio of the buckling load becomes smaller when the loading position is higher, as well as the difference of the buckling load between open and closed cross-sections and the variance of the post-buckling behaviour.

As described, the Runge–Kutta transfer matrix method was applied to the post-buckling analysis, which is a higher order nonlinear problem. In the traditional transfer matrix method, the calculation cannot be properly completed due to the omittance of significant digits if the analysis length of the structural system is increased. Therefore, in controlling such aspects, the solutions were obtained based on the modified transfer matrix method that was previously developed by the author of this paper (Usuki & Nakamura 1986).

## Acknowledgments

The author gratefully acknowledges the Deutsche Akademische Austausch Dienst for providing the opportunity to study at Technische Hochschule Darmstadt. The author also thanks Yasushi Okada, Atsushi Kurosawa and Kiyoshi Yogo for their academic assistance.

## Footnotes

- Received October 9, 2007.
- Accepted February 7, 2008.

- © 2008 The Royal Society