Three-dimensional beam-columns

Pergamon

00457949(94)E0132-L

THREE-DIMENSIONAL

Compurcrs & Svucrures Vol. 52. No. 3. pp. 449460. 1994 Copyright \G 1994 Elsevier Science Ltd Printed in Great Britain. All rights reserved 0045-7949194 s7.00 + 0.00

BEAM-COLUMNS

W. R. SPILLERS and M. H. SHAMS Department of Civil and Environmental Engineering, New Jersey Institute of Technology, Newark, NJ 07012, U.S.A. (Received 7 April 1993)

Abstract-It is well known for two-dimensional beams that so-called beam-column effects predict a flexural stiffness which exhibits strain softening (hardening) in the presence of axial compression (tension). This paper discusses solutions of three-dimensional beamxolumns which exhibit coupled softening or hardening which possesses a rather general dependence upon the initial state of stress in the beam. Applications to the nonlinear analysis of three-dimensional frames are discussed as are some special cases of three-dimensional beam
INTRODUCTION Structural engineers have been well served by so-called

‘beam-columns’ over the years. Roughly speaking, the beam-column equations include the effect of member displacement upon bending moment to the end that the member stiffness becomes a function of the axial load [l]. The beam
The fact is that we are comfortable when considering the softening effect of axial compression upon member stiffness in two-dimensional beam+olumns, and at least willing to think about the softening effect of bending moment in lateral buckling problems in three-dimensional beams. It will be seen below that the entire initial state of stress in a member can interact with the member stiffness matrix in a rather complex manner. Since any incremental analysis has initial states of stress associated with it which are not simply compression, there is every reason to include the appropriate ‘beam-column effects’ described below. This paper describes a general approach which subsumes all the effects cited above. A procedure for computing the member stiffness for three-dimensional beams is presented in which four coupled differential equations are solved numerically using the method of finite differences. In order to find the 6 x 6 member stiffness matrix for a three-dimensional beam, this system of equations must be solved six times. Several special cases are also discussed. THE EQUATIONS OF THREE-DIMENSIONAL BEAM-COLUMNS

The three-dimensional beam column equations used here are now available in the literature [4] and are presented in Table 1. This section reviews their derivation. First of all, a perturbation method [S] is used to describe beam response. That is, a given straight beam, in equilibrium under given forces, is subjected to a small load perturbation. This method has the advantage describing what are typically thought of as nonlinear effects without recourse to theories of large strain. It also produces the tangent stiffness now common in nonlinear analysis and describes buckling

449

450

W. R. Table

SPILLERS and 1. First-order

M. H.

SHAMS

equations

Force equilibrium i camp:

P. +j?

j camp:

T:+p‘:+P:d,,,

k

- Pv6 /

,,, - P!S:,, - P$s ,,,\ - Py:,vv - ppcs,,, - py:,, = 0 - PY0, + P:6,,v, - py,,,

+ pyc5,,, - py,

= 0

camp: p; +p, + Pyd:,, + Py9, + P$? :,,, + Pp0,,, + pt.d,,, + p:H, = 0

Moment equilibrium i camp:

R’+Ci I

j camp:

mi;r;+ tf~\+ My6,,, - Mge, + M:6,:\., - A@?,,, + m:S,:, - rn;tI, - p: - P:B, = 0

k

\ -M’6 , ~,,-M~6,,,-M~S,.,,,-MP6,,,,-m~6,,,-m96,,,+PP6,,,-P~6,,,=0

camp: iiir_+ ii?:+ Mf’d,,, + M~f?v + M’$?,,,, + M;e ,,., + m’$S,, + m:t?, + P,. - P;e, = 0

as the response to load perturbation becomes singular. The starting point is the equilibrium of a threedimensional beam segment [6] described by force and moment equilibrium as P’+P=o,

M’fm+txP=O.

venient to introduce displacements in the following manner. Let w represent the rotation vector associated with any beam element. If in the initial configuration the beam is straight and lies along the x axis, for the case of small rotations w can be written as

(1) w = B,i - Sz,,j + 6,,,k,

Here P and M are the usual force and moment stress-resultants with p and m the applied forces and moments; the prime refers to the differentiation with respect to arc length. The vector t is of course the unit tangent vector which will also be referred to as i in the applications cited below. In component form

(6)

where 0, is the torsional rotation, 6: is the z component of the beam displacement vector, 6, is the y component of the beam displacement vector, and the comma is used to indicate differentiation. The rotation vector can, of course, be used to determine the changes of the base vectors as

P=P,i+P,.j+PTk M=Mxi+M~vj+Mzk

(2) j"-.+j'=jo+~j=jo+o

and the base vectors ijk selected so that P, and M, represent thrust and torque, respectively. Going from the initial to the perturbed configuration, changes are indicated as

xj"

k”-+k’=ko+tE=ko+w

xk”.

(7)

Clearly w x i” = 6,,,Xjo+ 61,,k0

M”+M’=Mo+&

w x i” = - 6,., j” + 0, k”

iO-+i’=iO+$

(3)

where L is the anticipated small parameter and the overbar is used to distinguish the perturbation term. The perturbed terms of eqn (2) are now inserted into the equilibrium equations of the perturbed configuration given in eqn (1). In component form these equations are Pii’ + Pli” + Pf,‘j’ + Pij” + Pfk’ + Pfk” +p$’ +pfj’ +pfk’ = 0

(4)

MLi’ + Myi + Mij’ + M$’ + Mf’k’ + Mfk” +

m,ki’ + m$’

+

mik’ + t’ x P’ = 0.

(5)

In order to complete the analysis, it is simply a matter of inserting the definition of the perturbations [eqn (3)] into eqns (4) and (5) and collecting terms in the coordinate directions. Before doing so it is con-

o x k” = -6,,.Yio - O,j”.

(8)

The first-order equations are collected in Table 1; the zero-order equations are assumed to be satisfied by the equilibrium of the initial configuration and will not be discussed here. The equations of Table 1 are obtained by writing the equilibrium equations in the perturbed configuration and keeping terms which are linear in the small parameter L. For convenience this parameter can then be set to one in which case the terms indicated by bars [eqn (3)] then represent the full perturbation. It should be noted that the assumptions of small displacement theory have been invoked above to allow arc length differentiation to be replaced by differentiation with respect to the space variable X. THE MEMBER

STIFFNESS

MATRIX

In the most simple case of elastic analysis, a straight uniform beam segment with no cross-sectional

Three-dimensional beamxolumns

warping can be described by four uncoupled events: axial deformation, torsional deformation, and bending deformation about two axes. If the two-dimensional beam--column equations are used to introduce nonlinear effects, there is a coupling between axial forces and bending. If the equations of Table I are used, they introduce a coupling between all the stress resultants present in a three-dimensional beam. This section describes how the method of finite differences can be used to compute the member stiffness matrix for a three-dimensional frame element which includes all these couplings. The resulting member stiffness matrix is a tangent stiffness matrix due to the use of the perturbation method. First of all, from considerations of equilibrium it can be argued that the member stiffness matrix is a 6 x 6 matrix. Second, using Newton’s method for nonlinear structural analysis, each step (iteration) of the nonlinear analysis becomes simply a linear analysis which uses the local tangent stiffness. Therefore, the terms in the member stiffness matrix are ‘forces’ due to unit ‘displacements’. Third, if the member forces are chosen properly the terms in the member stiffness matrix may be computed by introducing sequentially six discontinuities into the boundary conditions of the system of equations in Table 1: l

a unit axial discontinuity,

l

a unit torsional discontinuity,

l

four flexural discontinuities.

and

Note that the four flexural discontinuities are those which are used in moment distribution: a unit rotation about an axis of flexure is applied at one end of a beam while the other end is held fixed. Two beam ends and two axes of flexure then imply four flexural discontinuities. The next step will combine some of the equations in Table 1 to produce a system of four equations in the four displacements 6,, a,, 6:, and 0,. Quite generally, and following linear beam theory, the second equation can be combined with the sixth equation eliminating the shear term &; similarly the third equation can be combined with the fifth equation in the table to eliminate the shear term P:. Furthermore at this point all the member loads can be eliminated as not of interest

P”, = Pp = Pp = p, = tr = Fz = 0.

(9)

451

To complete the formulation, equations are appended p.y = k,&,,,

moments, M!, MP, can at most be linear functions. The six equations of Table 1 then reduce to the following four equations

constitutive

M, = kre,..r

MY= - k1&,,, , m: = k,$.n .

(11)

Here the KS are usual spring constants from considerations of strength of materials. Some general comments on this system of equations can now be made: l

The last three equations are coupled and must be solved simultaneously; then the first equation can be integrated to complete the solution.

l

The last two equations are fourth order in the beam displacements (like the linear, elastic beam equations); the other two equations are second order.

l

The equations themselves are linear in x since the initial moment diagrams M,!, Mf may be linear in X. (Timoshenko [2] remarks that equations of this type may be solved using Bessel functions.) FINITE DIFFERENCE SOLUTION

As a practical approach, each beam is divided into a number of segments and the solution of eqns (IO) represented by a finite number of points using the method of finite differences. In the case of the solutions described below, 20 spaces have been taken along each beam. Using central differences and introducing fictitious points to handle the boundary conditions, 75 simultaneous equations were solved on a PC for every case which involved the last three of eqns (10); once these equations have been solved together, the first of eqns (10) can be solved. The computed solutions were in agreement with the charts of Gere [I] when those solutions were appropriate. Other numerical checks are described below. A FORTRAN subroutine which produces a 6 x 6 member stiffness matrix is included in Appendix 1.

This implies that the initial axial thrust and torque, Pt, M$ must be constant and that the initial bending

four

SPECIAL CASES

The elastic beam

Where the initial forces are all zero, eqns (10) degenerate to the solution for an elastic beam in which axial compression, torsion, and y and z axis bending are all uncoupled

452

W. R.

SPILLERSand

M. H.

SHAMS

k,6’: = 0

k,6’: =0 k,6:

= 0

(12) Two -dimensional beam -columns

When only an initial axial load, P”,, is present, eqns (10) again uncouple as described above. In this case, the y and z axis bending is controlled by the familiar two-dimensional beam-column equations k,6:’ = 0 k#:) = 0 -k&‘+

(15) The solution Appendix 2.

of eqn (15) is described further

The effect of initial torsion

The discussion of special cases has thus far been motivated by practical applications and available solutions and therefore has been dominated by the effect of axial load. An alternative would have been to start with the case of initial torsion Mt. In this case the equations are

Po,6; = 0

k,6~{,’- Ppd.; = 0.

k,6’: = 0

(13)

This is the case which is most commonly used today in the nonlinear analysis of three-dimensional frames.

k,O:’ = 0 -k&“+ k,bf”+

Lateral buckling

If a constant initial moment Mt is added to the problem, the member twist 0, becomes coupled with 6, and 6, remain uncoupled

M0,6:’ = 0 M0,6; = 0.

(14) The third of these equations is particularly interesting and can be thought of as a beamcolumn equation in which the constant PO,/k,.is replaced by

charts can

A more complex case

If an initial constant moment M: is now added, (!I,, 6), and 6; become fully coupled and it is no longer possible to invoke charts for well-known solutions. In this case the equations take on the form

t Professor Maciej Bieniek of Columbia University suggests that this case is of interest to, and known to, those working on the mechanics of drill shafts.

(16)

Here only S, and 6, are coupled and it is possible to solve the system by hand. This case is discussed further in Appendix 3.t

CONCLUDING

With this correction Gere’s beamcolumn also be used in this case.

in

REMARKS

This paper has explored several special cases of the three-dimensional beam
453

Three-dimensional beam-columns REFERENCES 1.

2. 3. 4. 5.

6.

7. 8. 9.

J. J. Gere, Moment Distribution. Van Nostrand, New York (1963). S. Timoshenko, Theory of Elastic Stability. McGrawHill, New York (1936). M. A. Biot, Mechank of Incremental Deformation. John Wilev. New York (1965). W. R. Spillers, A note’ on some qualitative results for beams and rods. Quart. Appl. Math. 48, Y-580 ( 1990). A. E. Green. R. J. Knops and N. Laws, Large deformations, superimposed small deformations, and the stability of elastic rods. Inf. J. Solids Sfruct. 4. 555-577 (1968). E. Reissner, Variational considerations for elastic beams and shells. J. Engng Mech. Div., ASCE 88.23-57 (1962). T. See and R. E. McConnell, Large displacement elastic buckling of space structures. J. Struck Engng 112, 1052-1069 (1986). W. R. Spillers, Automated Structural Analysis: An Introduction. Pergamon Press, New York (1972). System/360 Scientific Subroutine Package, IBM Corp., White Plains, New York (1970).

APPENDIX I: PROGRAM

LISTING

This appendix lists the code of a FORTRAN subroutine which will compute the 6 x 6 member stiffness matrix for a three-dimensional frame given the member force matrix and the appropriate stiffness coefficients. It follows the notation of [8] and incorporates the linear equation solver SIMQ from the IBM scientific subroutine package[9] now in the public domain. The code is written to be easily read and consequently could be written more efficiently. The subroutine divides each beam into 20 segments and then uses the method of central differences to compute the displacements at the resulting nodes. As described above, the last three of eqns (10) are first solved; then the first equation is solved. Each of these steps requires the solution of simultaneous equations. In order to generate the 6 x 6 stiffness matrix it is necessary to solve eqns (10) six times. The input to the subroutine STIFF is the member force matrix FORCE and some stiffness parameters: FORCE (I)-axial

load

FORCE (2)-torsion FORCE (3)--y-axis bending moment at +ve end of the number FORCE (4)-z-axis bending moment at fve end of the member FORCE (5t-_y-axis bending moment at -ve end of the member FORCE (6)-z-axis bending moment at -ve end of the member ALEN-member AKL-axial

length

stiffness

AKT-torsional

stiffness

AKY-y-axis

bending stiffness

AKZ-z-axis

bending stiffness.

The subroutine returns the 6 x 6 member stiffness matrix AK(I, J). The first part of the subroutine solves for 8,. 6, and 6, at each node point in the former of a linear system Ax = b. The order of the variables in the matrix is that just given, i.e.

With 20 spaces (21 nodes) and two fictitious nodes at each end for boundary conditions, the system matrix A is 25 x 3 by 25 x 3. This implies an unused component 0, at each beam end since the torsional response is second order while the bending response is fourth order. When generating the system matrix A, each term in the differential equations eqns (10) ‘puts’ corresponding terms in A. In the program that is done by the PUT subroutines. For example, subroutine PUT 4 introduces the central difference coefficients of the fourth derivative operator into the matrix A. These equations are homogeneous (b = 0) except for boundary condition terms. For example, in order to introduce a unit torsional discontinuity it is necessary to specify that or, = B,,, = 0 fly, =0 0, _* = 1

(unused fictitious points)

(0, =0

at x =0)

(0, = 1

at x = /)

In terms of FORTRAN rise to O,, =O-+A(l,

code, these conditions

I)= 1.0,

B(1) = 0

Or,,= 0 -+ A(2, M - 2) = 1.O.

B(2) = 0

O,,=O+A(3,7)=

B(3) = 0

&,,_,=O,

1.0,

l+A(4,M-8)=1.0,

give

B(4)=0.1.

The other boundary conditions follow in a similar fashion. After the node displacements are computed by the subroutine SIMQ, the subroutines TORQ, BENDY, BENDZ are called to find the torque and the bending moments. The axial response [the first of eqns (IO)] is determined in a similar fashion. Finally, reactions are computed as elements of the member stiffness matrix. This, of course, requires corrections since the perturbation method gives results in the deformed coordinate system while the structural analysis program requires components in the member coordinate system. A listing of the FORTRAN code follows.

W. R.

454

C

SPILLERS and

M. H.

SHAMS

Dummy main program to access subroutine DIMENSION FORCE(6),AK(6,6) OPEN(UNIT=6,FILE="BMGEN.OUT") READ(5,l) AKL,AKT,AKY,AKZ,ALEN,FORCE WRITE(G,l)AKL,AKT,AKY,AKZ,ALEN,FORCE 1 FOFU+AT(5F5.0/6F5.2) CALL STIFF(ALEN,AKL,AKT,AKY,AKZ,FORCE,AK) WRITE(6,75)((AK(I,J),J=1,6),1=1,6) 75 FORMAT(6E18.8) STOP END

C C

SUBROUTINE STIFF(ALEN,AKL,AKT,AKY,AKZ,FORCE,AK) GENERALIZED MEMBER STIFFNESS COMMON A,B,C,Al COMMON /PROP/ H DIMENSION A(75,75),B(75),C(25),.A1(25,25),FORCE(6) 1 ,AK(6,6),DISC(6) *X-FORCE(5) AMYO(X)=PZO AMZO(X)=-PYO *X-FORCE (6) PY~=-(F~RCE(~)+F~RCE(G))/ALEN PZO=(FORCE(3)tFORCE(5))/ALEN P=FORCE(l) TOR=FORCE(2) NSPACES=24 N=NSPACES+l M=3*N N2=N-2 H=ALEN/FLOAT(NSPACES-4) DO 22 1=1,6 22 DISC(I)=O. Solve for bending and torsion DO 999 ICOL=1,6 DISC(ICOL)=l. DO 1 I=l,m B(I)=O. DO 1 J=l,M 1 A(I,J)=O. DO 2 1=3,N2 X=(1-3) X=X*H CALL PUT4(1,3,3,AKZ) CALL PuT1(1,3,1,2.*PZO) CALL PUT3(1,3,2,TOR) CALL PUT2(1,3,1,AMYO(X)) CALL PuT2(1,3,3,-P) CALL PUT4(1,2,2,-AKY) CALL PUT1(1,2,1,2.*PY0) CALL PUT3(1,2,3,TOR) CALL PUT2(1,2,1,-AMZO(X)) CALL PUT2(1,2,2,P) CALL PUT2(I,l,l,AKT) CALL PUTl(I,1,3,-PZ0) CALL PUT1(1,1,2,PYO) CALL PUT2(1,1,3,-AMYO(X)) CALL PUT2(1,1,2,-AMZO(X)) CALL PUT1(1,1,3,PZO) CALL PUT1(1,1,2,-PYO) 2 CONTINUE

Three-dimensional beam-columns

BC FOR THX A(l,l)=l. A(2,M-2)=1. A(3,7)=1. A(4,M-8)=1. B(4)=DISC(Z) BC FOR DZ A(5,8)=1. A(6,M-7)=1. A(M,11)=.5/H A(M,5)=-.5/H B(M)=-DISC(5) A(M-l,M-4)=.5/H A(M-l,M-10)=-.5/H B(M-l)=-DISC(3)

66

67

11

4

69

BC FOR DY A(M-2,M-6)=1. A(M-3,9 )=l. A(M-4,12 )=.5/H A(M-4,6 )=-.5/H B(M-4)=DISC(6) A(M-5,M-3)=.5/H a(m-5,M-9)=-.5/H B(M-5)=DISC(4) CALL SIMQ(A,B,M,KS) CALL TORQ(AKT,3,l,ATORM) CALL TORQ(AKT,N_2,l,ATORP) CALL BENDY(AKY,3,2,AMYMl) AMYMl=-AMYMl CALL BENDZ(AKZ,3,3,AMZMl) AMZMl=-AMZMl CALL BENDY(AKY,N_2,2,AMYPl) CALL BENDZ(AKZ,N_2,3,AMZP1) WRITE(6,66) (I,(B(3*1-3tJ),J=l,3),I=l,N) FORMAT(I5r3E20.8) . Solve for axial compression WRITE(6,67)ATORM,ATORP,AMYMl,AMZMl,AMYPl,AMZPl FORMAT(6E20.8) DO 11 I=l,N C(I)=O. DO 11 J=l,N Al(I,J)=O. DO 4 1=3,N2 CALL PUT22(I,l,l,AKL) CALL GET2(1,3,VAL) C(I)=C(I)tVAL*PYO CALL GET2 (1,2,VAL) C(I)=C(I)tVAL*PZO Al(l,l)=l. Al(N,N)=l. A1(2,3)=1. Al(N-l,N-2)=1. c(n-l)=DISC(l) CALL SIMQ(Al,C,N,KS) WRITE(6,69) (I,C(I),I=l,N) FORMAT(I5,E20.8) CALL THR(AKL,3,PM)

455

W. R. SPILLERS and

456

M. H. SHAMS

CALL THR(AKL,N-2,PP) WRITE(6,68)PM,PP 68 FORMAT(2E20.8) L

C

C

555

556

557

558

559

999

Components for stiffness matrix Corrections for deformed geometry IF(ICOL.NE.2) GO TO 555 AMYPl= AMYPl-FORCE(4) AMZPl= AMZPltForce(3) CONTINUE IF(ICOL.NE.3) GO TO 556 AMZPl=-FORCE(2)tAMZPl PP=PM ATORP=ATORM CONTINUE IF(ICOL.NE.4) GO TO 557 AMYPl=AMYPltFORCE(2) PP=PM ATORP=ATORM CONTINUE IF(ICOL.NE.5) GO TO 558 AMZMl=AMZMltFORCE(2) CONTINUE IF(ICOL.NE.6) GO TO 559 AMYMl=AMYMl_FORCE(2) CONTINUE AK(l,ICOL)=PP AK(2,ICOL)=ATORP AK(3,ICOL)= AMYPl AK(4,ICOL)=AMZPl AK(5,ICOL)= AMYMl AK(G,ICOL)= AMZMl DISC(ICOL)=O. CONTINUE RETURN END SUBROUTINE PUTl(IROW,IEQ,IVAR,COEFF) COMMON A COMMON /PRO /H DIMENSION A(75,75) 1=3*IROW-3tIEQ J=3*IROW-3tIVAR A(I,J+3)=A(I,Jt3)tCOEFF*.5/H A(I,J-3)=A(I,J-3)-COEFF*.5/H RETURN END

C

SUBROUTINE PUT2(IROW,IEQ,IVAR,COEFF) COMMON A COMMON /PROP/ H DIMENSION A(75,75) 1=3*IROW-3tIEQ J=3*IROW-3+IVAR A(IIJt3)=A(I,Jt3)tCOEFF/H**2 A(I,J-3)=A(I,J-3)+COEFF/H**2 A(I,J )=A(I,J )-COEFF*2./H**2 RETURN END

Three-dimensional beamxolumns c

SUBROUTINE PUT22(IROW,IEQ,IVAR,COEFF) COMMON A,B,C,Al COMMON /PROP/ H DIMENSION A(75,75),B(75),C(25),Al(25,25) I= IROW-1tIEQ J= IROW-1tIVAR A1(I,Jt1)=A1(I,Jt1)tCOEFF/H**2 A1(I,J-1)=Al(I,J-1)tCOEFF/H**2 Al(I,J )=Al(I,J )-COEFF*2./H**2 RETURN END SUBROUTINE PUT3(IROW,IEQ,IVAR,COEFF) COMMON A COMMON /PROP/ H DIMENSION A(75,75) 1=3*IROW-3tIEQ J=3*IROW-3tIVAR A(I,Jt3)=A(I,Jt3)-COEFF/H**3 A(I,J-3)=A(I,J_3)tCOEFF/H**3 A(I,Jt6)=A(I,Jt6)tCOEFF*.5/H**3 A(I,J-6)=A(I,J-6)-COEFF*.5/H**3 RETURN END L

SUBROUTINE PUT4(IROW,IEQ,IVAR,COEFF) COMMON A COMMON /PROP/ H DIMENSION A(75,75) 1=3*IROW-3tIEQ J=3*IROW-3tIVAR A(I,Jt3)=A(I,Jt3)-COEFF*4./H**4 A(I,J-3)=A(I,J-3)-COEFF*4./H**4 A(I,Jt6)=A(I,Jt6)tCOEFF /H**4 A(I,J-6)=A(I,J_6)tCOEFF /H**4 A(I,J )=A(I,J )tCOEFF*6./H**4 RETURN END SUBROUTINE TORQ(AKT,NODE,IVAR,TOR)

COMMON /PROP/H

COMMON A,B DIMENSION A(75,75),B(75) 1=3*NODE-3tIVAR TOR=AKT*(B(It3)-B(I-3))*.5/H RETURN END SUBROUTINE BENDY(AKY,NODE,IVAR,BEND)

COMMON /PROP/H

COMMON A,B DIMENSION A(75,75),B(75) 1=3*NODE-3tIVAR BEND=-AKY*(B(It3)tB(I-3)-2.*B(I))/H**2 RETURN END C SUBROUTINE BENDZ(AKZ,NODE,IVAR,BEND)

como~ /PROP/H

W.

458

R.SPILLERSandM. H. SHAMS

COMMON A,B DIMENSION A(75,75),B(75) 1=3*NODE-3tIVAR BEND=AKZ*(B(It3)+B(I-3)-2_*B(I))/H**2 RETURN END

C SUBROUTINE GET2 (NODE,IVAR,VAL ) COMMON /PROP/H COMMON A,B DIMENSION A(75,75),B(75) 1=3*NODE-3tIVAR VAL = (B(It3)tB(I-3)-2.*B(I))/H**2 RETURN END C SUBROUTINE THR(AKL,NODE,P) COMMON /PROP/H COMMON A,B,C DIMENSION A(75,75),B(75),C(25) I=NODE P =AKL*(C(Itl)-C(I-1))*.5/H RETURN END C SUBROUTINE SIMQ(A,B,N,KS) DIMENSION A(l),B(l) C C C

FORWARD SOLUTION TOL=O.O KS=0 JJ=-N DO 65 J=l,N JY=J+l JJ=JJtNtl BIGA=O IT=JJ-J DO 30 I=J,N

C C C

SEARCH FOR MAXIMUM COEFFICIENT IN COLUMN

IJ=ITtI IF(ABS(BIGA)-ABS(A(IJ))) 20,30,30 20 BIGA=A(IJ) IMAX=I 30 CONTINUE TEST FOR PIVOT LESS THAN TOLERANCE (SINGULAR MATRIX) C C IF(ABS(BIGA)-TOL) 35,35,40 35 KS=1 RETURN C INTERCHANGE ROWS IF NECESSARY C C 40 Il=JtN*(J-2) IT=IMAX-J DO 50 K=J,N Il=IltN 12=11tIT

Three-dimensional SAVE=A

beam+columns

459

(11)

A(Il)=A(I2) A (12)

=SAVE

C

C

DIVIDE

EQUATION

BY

LEADING

COEFFICIENT

C 50

A(Il)=A(Il)

/BIGA

SAVE=B

(IMAX)

B(IMAX)=B(J) B

C C C

(J)=SAVE/BIGA

ELIMINATE IF

55

NEXT

(J-N)

VARIABLE

55,70,55

IQS=N*(J-1)

DO 65

IX=JY,N

IXJ=IQS+IX IT=J-IX DO

60

JX=JY,N

IXJX=N*

(JX-1)

+1X

JJX=IXJX+IT

60 65

A(IXJX)

=A(IXJX)

- (A(IXJ)

*A(JJX)

)

B(IX)=B(IX)-(B(J)*A(IXJ))

C C C

BACK

70

SOLUTION

NY=N-1 IT=N*N 80 J=l,NY IA=IT-J IB=N-J IC=N DO 80 K=l, J B(IB)=B(IB)-A(IA)*B(IC) IA=IA-N

DO

80

IC=IC-1 RETURN END APPENDIX 2

This Appendix discusses the solution of a beam subject to an initial state constant biaxial bending My and My. is trivial. The remaining equations can form as

A4: 0 k,

of eqns (1 S), the case of axial load Pt and The first of eqns (15) be written in operator

+ Bs sin k,x + B,, cos kzx

-Mf

MY

k,D’-P’:

The determinant

equation

+B,~~+B,sink,x+E,cosk,x

6, =

c, + czx + C,x’ + C,d + C,x”

0 +C,x’+C,sink,.x+C,cosk,.r

for this system

is + C, sin kg

+ M”‘(k,D’

- P:‘, + M;‘(k, D’-

which is basically quadratic in D’ and regarded to be accessible. The solutions generally have the form 0, = A, + A2x + A+‘+

A,\-‘+

+A,r’+A,sink,s+A,cosk,s + A,, sin k2_u + A,,, cos k2s

P:)) = 0

thus should be of this system

A,.u4

+ C,, cos kg,

where k, and kz are the roots of the quadratic described above. Coupling of these coefficients requires that the polynomial terms of order quadratic and higher must vanish and also provides eight relationships between the coefficients of the trigonometric terms, The IO boundary conditions then complete the problem statement, For one particular case the solution of this system has been completed. In this case there was at most a 2% disagreement with the numerical solution generated using the finite difference method.

W. R. SPILLERS and M. H. SHAMS

460 APPENDIX

3: A PROBLEM

Euler buckling

OF TORSION

Y

This appendix discusses the case of a beam subjected to an initial state of torsion, Mt. First the solution of the special case of eqns (16) is discussed. Then, starting from first principles, the case of flexural bucking under torsion is discussed. For eqns (I 6) the first two equations are trivial but the last two give a coupling of 6, and 6, as

Here D represents the derivative d/dx. If the determinant this system is set to zero it follows that

and that 6, and 6, have the form 6, = A, + A,.~ + A$

+ A,x’

Mz = P:Sy

of

(deformed

quilibrium)

kiz = kzSj (constitutive equation)

Torsional

buckling

Y

+ A$

+ A,.~’ + A, sin kx + A, cos kx

-I- -Mf

6, = B, + B2x + B,x’ + B4xj + 8,x4

,/------.\

M,"

)I

+ B,x5 + E, sin kx + BR cos kx with k’ = My'/(k, , k:). Coupling of these coefficients that A, = A, = A, = E4 = E, = Bh = 0 and that -K,.A,k

+ B,M:

= 0,

K,.A,k

+ B,M:

implies

=O.

Given the eight boundary conditions (flexure about two axes) of this system, the As and Bs can be computed explicitly by hand. For the particular case investigated here there was at most a 3% disagreement between the formal solution and the results of the finite difference solution. There is a problem, a subset of eqns (16), in which torsion can produce flexural buckling which is accessible to hand calculation (Fig. Al). For the torsional buckling problem of the figure 6, and S, couple as

tii,=

liiy = -M;,;

(deformed

equilibrium,

xy plane)

tiz = -M;S;

(deformed

equilibrium,

xz plane)

, MY=-ky6;

K.8;

Fig. Al. Simple buckling

with k2 = My/k, implies that A, = B, = 0,

-k,.kA,

+ M:B4 = 0.

S,.=6: equation

implies solutions

of the coefficients

k,.kA, + MPB, = 0.

is complemented

at x =0

and

with the four

x = I.

in this case is

D2(k,.klD1 + My) Which

=0

cqs.)

cases.

kl. In this case coupling

The problem statement boundary conditions

The determinant

(constittnive

Setting the determinant of the matrix of the coefficients equal to zero in this case implies the characteristic equation

= 0.

cos KL(cos

in the form

S, = A, + Az.x + A, sink,

+ A,cos

k,

6,= B,+ &.Y + B,sin k, + B,cos k,

The smallest nontrivial buckling torque is

root

KL - I) = 0. is KL = ?r/2 or the critical