High flexibility problem

[eroz]

Getting higher flexibility (larger period and larger lateral deflection under unit lateral load) from a model in OpenSees in comparison to same one developed in SAP and result of hand calculations.

The model is relatively simple: 2-dimensional steel frame, 1 story, 1 bay, with chevron bracing. All members elasticBeamColumn. Fixed base. Brace elements moment released at both ends.

Tcl script is below. What seems to be causing the discrepancy; was I not able to reflect what I described above inside the script correctly (error in defining moment releases perhaps or etc.)? Thanks in advance.

p.s. certain by now that section properties and dimension have been entered in correctly

_____________________________________________________
model BasicBuilder -ndm 2 -ndf 3

#DEFINE NODES
node 1 0 0
node 111 0 0
equalDOF 1 111 1 2

node 2 0 156

node 6 360 180

node 5 360 0
node 551 360 0
equalDOF 5 551 1 2

node 9 180 156
node 991 180 156
node 992 180 156
equalDOF 9 991 1 2
equalDOF 9 992 1 2

#ESTABLISH FIXED SUPPORTS
fix 1 1 1 1
fix 5 1 1 1

#ASSIGN JOINT MASS AT ROOF
mass 2 0.593636833 0.0 0.0
mass 6 0.593636833 0.0 0.0

geomTransf Linear 1

#COLUMNS
element elasticBeamColumn 1 1 2 28.2 29000 833 1
element elasticBeamColumn 4 5 6 28.2 29000 833 1

#BEAM
element elasticBeamColumn 7 2 9 14.1 29000 484 1
element elasticBeamColumn 8 9 6 14.1 29000 484 1

#BRACES
element elasticBeamColumn 13 111 991 9 45595 6.75 1
element elasticBeamColumn 14 992 551 9 45595 6.75 1

#CHECK PERIOD
set xDamp 0.05; # damping ratio (0.02-0.05-typical)
set lambda [eigen 1]
set omega [expr pow($lambda,0.5)]
set Tperiod [expr 2*3.14159265358979323846/$omega]; # period (sec.)
puts "T = $Tperiod"

#CHECK LATERAL DEFLECTION AT ROOF
recorder Node d6.xls disp -node 6 -dof 1

#APPLY UNIT LATERAL LOAD AT ROOF
pattern Plain 1 "Linear" {
load 6 1.0 0.0 0.0
}

#STATIC ANALYSIS
constraints Plain
numberer RCM
system UmfPack
test NormDispIncr 1.0e-8 10
algorithm Newton
integrator LoadControl 1
analysis Static
analyze 1

[silvia]

add a mass equal to 1e-9 to a secondary dof, not sure if that will change anything.

you have to figure out what is different between the models, and how the programs handle the differences.

[eroz]

The mass addition did not change anything.

The models are the same. Such a simple and elastic model probably is not handled differently between programs. I was able to verify my result (lateral deflection) from SAP with hand calculations. But getting something different from the OpenSees model. Obviously I am doing something wrong in the script for OpenSees - stg. missing or unintended from what I described in my previous post.

I have been looking at this for days now and it is beyond me what is causing the discrepancy.

[silvia]

how are the braces handled?
why aren't you just using trusses?

[fmk]

murat,

there is nothing wrong with the results you are getting from OpenSees. Also you should use trusses as silvia says. also based on the latest script you sent node 6 is too high.

for those interested here is the corrected scriptfollowed by a matlab script that does the same thing and gives the same results:

OpenSees Script:
[code]
model BasicBuilder -ndm 2 -ndf 3

# Constants
set PI 3.14159265358979323846;
set g 386.6; # gravitational acceleration
set nodalmass 0.593636833937; # mass per node

# Nodes
node 1 0 0
node 2 0 156
node 5 360 0
node 6 360 156
node 9 180 156

# Boundary Conditions
fix 1 1 1 1
fix 5 1 1 1

# Mass
mass 2 $nodalmass 0.0 0.0
mass 6 $nodalmass 0.0 0.0

geomTransf Linear 1

# Columns
element elasticBeamColumn 3 1 2 28.2 29000 833 1
element elasticBeamColumn 4 5 6 28.2 29000 833 1

# Beams
element elasticBeamColumn 5 2 9 14.1 29000 484 1
element elasticBeamColumn 6 9 6 14.1 29000 484 1

# BRACES
uniaxialMaterial Elastic 2 45595
element truss 1 1 9 9.0 2
element truss 2 9 5 9.0 2

#======================================================================================
# FIND STRUCTURAL PERIOD
#======================================================================================
set lambda [eigen 1]
set omega [expr pow($lambda,0.5)]

set Tperiod [expr 2*$PI/$omega]; # period (sec.)
puts "T = $Tperiod"
[/code]

Matlab Script:

[code]

% IDs for mapping nodal quantaties to structure dof
% node 1 - fixed -> -1 -1 -1
% node 2 - free -> 1 2 3
% node 5 - fixed -> -1 -1 -1
% node 6 - free -> 4 5 6
% node 9 - free -> 7 8 9
%

%
% Structure matrices (9 free dof in model)

Kstruct = zeros(9,9);
Mstruct = zeros(9,9);

%
% Add Column Stifnesses to K matrix
%

A=28.2;
E=29000;
I=833;
L=156;

ncol = 2;
Kcol=[[4*E*I/L 2*E*I/L 0];[2*E*I/L 4*E*I/L 0]; [0 0 E*A/L]];
Tcol = [[-1/L 0 1 1/L 0 0];[-1/L 0 0 1/L 0 1];[0 -1 0 0 1 0]];
Kcol = transpose(Tcol)*Kcol*Tcol

IDcols=[[-1 -1 -1 1 2 3];[-1 -1 -1 4 5 6]];

for i = 1:ncol,
for j = 1:6,
for k = 1:6,
row = IDcols(i,j);
col = IDcols(i,k);
if (row ~= -1) && (col ~= -1)
Kstruct(row,col) = Kstruct(row,col) + Kcol(j,k);
end
end
end
end

%
% Add Beam Stifnesses to K matrix
%

A=14.1;
E=29000;
I=484;
L=180;

nbeam = 2;
Kbeam=[[4*E*I/L 2*E*I/L 0];[2*E*I/L 4*E*I/L 0]; [0 0 E*A/L]];
Tbeam = [[0 1/L 1 0 -1/L 0];[0 1/L 0 0 -1/L 1];[-1 0 0 1 0 0]];
Kbeam = transpose(Tbeam) * Kbeam * Tbeam

IDbeams=[[1 2 3 7 8 9];[7 8 9 4 5 6]];


for i = 1:nbeam,
for j = 1:6,
for k = 1:6,
row = IDbeams(i,j);
col = IDbeams(i,k);
if (row ~= -1) && (col ~= -1)
Kstruct(row,col) = Kstruct(row,col) + Kbeam(j,k);
end
end
end
end

%
% Add brace stifnesses
%

nbrace = 2;
Lbrace = sqrt(180*180+156*156);
theta = atan(156.0/180.0);
Abrace = [9 9];
Ebrace = [45595 45595];
thetaBrace = [theta -theta theta -theta theta -theta];
IDbrace = [[-1 -1 7 8]; [7 8 -1 -1]];

for i = 1:nbrace,
A = Abrace(i);
E = Ebrace(i);
cs = cos(thetaBrace(i));
sn = sin(thetaBrace(i));
Tbrace = [-cs -sn cs sn];
kbrace = E * A / Lbrace;
kbrace = transpose(Tbrace) * kbrace * Tbrace;

for j = 1:4,
for k = 1:4,
row = IDbrace(i,j);
col = IDbrace(i,k);
if (row ~= -1) && (col ~= -1)
Kstruct(row,col) = Kstruct(row,col) + kbrace(j,k);
end

end
end
end

%
% Add Masses to mass matrix
%

mass = 0.593636833937;
nmass = 2;
IDm = [1 4];
for i = 1:nmass,
loc = IDm(i);
Mstruct(loc,loc) = mass;
end

%
% compute eigenvalues & print non-zero periods
%

sprintf('Eigenvalues Generalized Problem');
a=eig(Kstruct, Mstruct);
for i = 1:numel(a),
if isinf(a(i)) == 0
period = 2*pi()/sqrt(a(i))
end
end
[/code]