Inelastic Lateral Torsional Buckling of a I-Shaped Steel Beam

I-shaped beams which are widely used in steel design; under transverse load (or any flexural effect) loaded by gradient normal stress and partially acts like a compression member which tends to buckle. Compressive stress leads a portion of the member to shift from the centerline while a portion of the member under tensile stress tends to keep its initial position. Eventually, resulting motion is a combination of sectional torsion and lateral flexural buckling. 

When the problem is elastic -phenomenon can be elastic or inelastic dependent on the resulting stresses on the member- it is possible to approximate to the solution by a combination of linear static analysis and an buckling eigenvalue analysis. 

But in the cases that plastic strains occur, updating geometry to the deformed shape to create an initial imperfection and gradually increase the static load until the system loses its stability is necessary.

Following model does both;

Pre-processing commands

1. Open pro-processing phase
/PREP7

2. Define a set of parameters
H = 120E-3
B = 58E-3     
TF =7.7E-3       
TW =5.1E-3
L = 2000E-3       

It's always good to check console output if you're using command line. 

3. Define a new element type
ET,1,185

4. Material Data Input
MPTEMP,,,,,,,, 
MPTEMP,1,0 
MPDATA,EX,1,,200E9     !Pa - N/m2
MPDATA,PRXY,1,,0.3     
MPDATA,DENS,1,,7850    !kg/m3

TB,BISO,1,1,2, 
TBTEMP,0
TBDATA,,235E6,200E7,,,,

5. Derivation of keypoint locations and geometric properties of the section from given parameters

X1 = TW / 2
X2 = B / 2
X3 = B / 2
X4 = -B / 2
X5 = -B / 2
X6 = -TW / 2
X7 = X6
X8 = X5
X9 = X4
X10 = X3
X11 = X2
X12 = X1

Y1 = H / 2 - TF
Y2 = H / 2 - TF
Y3 = H / 2
Y4 = H / 2
Y5 = H / 2 - TF
Y6 = H / 2 - TF
Y7 = -Y6
Y8 = -Y5
Y9 = -Y4
Y10 = -Y3
Y11 = -Y2
Y12 = -Y1

IX1 = TW*(H-2*TF)*(H-2*TF)*(H-2*TF)/12
IX2 = TF*TF*TF*B/12+B*TF*(H/2-TF/2)*(H/2-TF/2)
IX = IX1 + 2*IX2

WX = 2*IX/H
MX = 1
SIGX = MX/WX
SLOPE = 2*SIGX/H

6. Adding keypoints forming the sections contours

K,1,X1,Y1
K,2,X2,Y2
K,3,X3,Y3
K,4,X4,Y4
K,5,X5,Y5
K,6,X6,Y6
K,7,X7,Y7
K,8,X8,Y8
K,9,X9,Y9
K,10,X10,Y10
K,11,X11,Y11
K,12,X12,Y12

7. Creating a cross-sectional area (A1) using keypoints 


A,1,1,2,3,4,5,6,7,8,9,10,11,12

8. Extrude the area (A1) along its normal extursion of L and create a volume (V1)

VEXT,1,,,,,L,,,,

9. Specify a global element size 

ESIZE,10E-3,0,

10. Mesh the volume (V1) by sweeping method 

VSWEEP,1

11. Apply translational restraint to the specified lines (L21 and L9)

DL,21,,UX,0
DL,21,,UY,0
DL,21,,UZ,0

DL,9,,UX,0
DL,9,,UY,0

12. Load the surface with gradient pressure implying the unit moment on cross-section

SFGRAD,PRES,0,Y,H/2,SLOPE
SFA,2,,PRES,SIGX
SFA,1,,PRES,SIGX

13. Linear static analysis

/SOL
SOLVE,1
FINISH

14. Buckling analysis which extracts buckling mode shapes and frequecies

/SOL
ANTYPE,1
BUCOPT,LANB,1,1E-6,1E6,RANGE 
SOLVE
FINISH

15. Post-processing commands

/POST1
/EDGE,1,0,45
/GLINE,1,-1

SET,LAST11

*GET,FACTOR,ACTIVE,0,SET,FREQ


16. Go back to pre-processing

/PREP7

17. Update initial geometry with the deformed shape (buckling mode shape)

UPGEOM,1,,,FILE, RST
CDWRITE,DB,FILE,CDB

18. Calculate the ultimate bending stress (from elastic case) and apply them to sections

SIGX_2 = FACTOR/WX
SLOPE_2 = 2*SIGX_2/H

SFGRAD,PRES,0,Y,H/2,SLOPE_2
SFADELE,1,1,PRES
SFADELE,2,1,PRES

SFA,2,,PRES,SIGX_2
SFA,1,,PRES,SIGX_2

19. Nonlinear analysis settings;
-Analysis Type: Static,
-Geometric Nonlinearities: On,
-Substeps (Current, Max, Min) : 50,100,25,
-Save Results: Each Step
-Auto Time-step: On,
-Line Search: On,
-Maximum iterations at each timestep: 50
-End Time. 1 s.

/SOLU
ANTYPE,0
NLGEOM,1
NSUBST,50,100,25
OUTRES,ERASE
OUTRES,ALL,ALL 
AUTOTS,1
LNSRCH,1
NEQIT,50
TIME,1 

20. Solve the case
SOLVE  
FINISH

21. Since a nonlinear buckling case is simulated here using static analysis type (buckling eigenvalue is only valid for linear cases) solution will be ceased after some point that the system lost its stability. At this point, you can either integrate the axial normal stresses on cross-sections to obtain ultimate moment capacity to LTB case, or just check the time (0.72 out of 1 s) to see the ratio between elastic and inelastic capacities.

/POST1
SET,LAST
SET,PREV
*GET,FACTOR_2,ACTIVE,0,SET,TIME

22. Export the result to a text file

ADIV = ' | '
*CFOPEN,OUTPUT.OUT
*VWRITE,'ELAS:',FACTOR,ADIV,'INELASTIC:',FACTOR*FACTOR_2,
(A6,g16.8,A3,A6,g16.8)
*CFCLOSE

Resulting text file content is;

ELAS:    17280.442     | INELAS   12269.114   

to validate the simulation using analytical approach, I've written the following MATLAB script that uses AISC360-10 equations to calculate ultimate moment capacities under constant moment value;

E = 20000;
Fy = 23.5;

b = 5.8;
h = 12.0;
tf = .77;
tw = .51;

h0 = h - tf;

[Ix,Iy,A,ix,iy,Wx,Wy,Zx,J,its]=calculateIProperties(h,b,tf,tw);


AA = J/Wx/h0;
BB = E/Fy/0.7;

Lr = 1.95*its*BB*sqrt(AA+sqrt(AA^2+6.76*(1/BB)^2));
Lp = 1.76*iy*sqrt(E/Fy);

% Yielding Strength

Mp = Zx*Fy;


% Inelastic LTB Strength

Mn = Mp - (Mp - 0.7*Fy*Wx)*((Lb - Lp)/(Lr - Lp));

% Elastic LTB Strength

Fcr = pi^2*E/(Lb/its)^2*sqrt(1+0.078*(J/(Wx*h0))*(Lb/its)^2);

Mne = Fcr*Wx;

where the function [..] = calculateProperties(..) is defined as;

function [Ix,Iy,A,ix,iy,Wx,Wy,Zx,J,its] = calculateIProperties(h,b,tf,tw)
    A = 2*b*tf + (h-2*tf)*tw;
    Ix = tw*(h-2*tf)^3/12+2*(tf^3*b/12+b*tf*(h/2-tf/2)*(h/2-tf/2));
    Iy = tw^3*(h-2*tf)/12+2*b^3*tf/12;
    ix = sqrt(Ix/A);
    iy = sqrt(Iy/A);
    Wx = Ix/(h*0.5);
    Wy = Iy/(h*0.5);
    Zx = (b*tf*(h-tf)*0.5+(h/2-tf)^2*tw*0.5)*2;
    Cw = (h-tf)^2*b^3*tf/24;
    its = sqrt(sqrt(Iy*Cw)/Wx);
    J = (2*b*tf^3 + (h-tf)*tw^3)/3;
end


Matlab script yields strength values for inelastic and elastic buckling modes;
- Mn = 1211.9829
- Mne = 1624.8328