Skip to content

Download MATLAB script

Nonlinear seismic response of a moment resisting frame¤

Original examples file: OpenSeesPy-Examples/NLSMRF.ipynb at master · AmirHosseinNamadchi/OpenSeesPy-Examples

The ground-motion file: RSN960_NORTHR_LOS270.AT2.

This is an OpenSeesMatlab simulation of a moment resisting frame subjected to seismic excitation. The model was introduced by C. Kolay & J. M. Ricles in their paper entitled Assessment of explicit and semi-explicit classes of model-based algorithms for direct integration in structural dynamics. The beams and columns of the MRF are modeled using dispBeamColumn fiber elements. The gravity load resisting system associated with the MRF is modeled using a lean-on column composed of linear elastic beam-column elements with 2nd Order \(P-\Delta\) effects [1].

image_0.png 

1
2
3
4
5
clear; clc; close all;

opsMAT = OpenSeesMatlab;
ops = opsMAT.opensees;
opsvis = opsMAT.vis;

Units¤

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
%% ------------------------------------------------------------------------
%  Units
% -------------------------------------------------------------------------
m    = 1.0;              % meter
KN   = 1.0;              % kilo-Newton
sec  = 1.0;              % second
inch = 0.0254*m;         % inch

kg = KN*(sec^2)/m;       % derived mass unit
g  = 9.81*(m/sec^2);     % gravitational acceleration

Earthquake record¤

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
%% ------------------------------------------------------------------------
%  Earthquake record
% -------------------------------------------------------------------------
dt = 0.01*sec;

recordFile = 'utils/RSN960_NORTHR_LOS270.AT2';
if ~isfile(recordFile)
    error('Ground-motion file not found: %s', recordFile);
end

rawMotion = readmatrix(recordFile, 'FileType', 'text', 'NumHeaderLines', 4);
rawMotion = rawMotion.';
rawMotion = rawMotion(:);
rawMotion = rawMotion(~isnan(rawMotion));

accel = rawMotion*3*g;
timeVec = (0:numel(accel)-1)'*dt;

figure;
plot(timeVec, accel, 'LineWidth', 1.5);
figure_0.png

Model definition¤

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
%% ------------------------------------------------------------------------
%  Model definition
% -------------------------------------------------------------------------
ops.wipe();
ops.model('basic', '-ndm', 2, '-ndf', 3);

%% ------------------------------------------------------------------------
%  Geometry
% -------------------------------------------------------------------------
node_coords = [ ...
    0,   0;
    6,   0;
    0,   1;
    6,   1;
    0,   2;
    6,   2;
    0,   2.5;
    6,   2.5;
    0,   3;
    0.5, 3;
    1,   3;
    3,   3;
    5,   3;
    5.5, 3;
    6,   3;
    0,   3.5;
    6,   3.5;
    0,   4;
    6,   4;
    0,   5;
    6,   5;
    0,   5.5;
    6,   5.5;
    0,   6;
    0.5, 6;
    1,   6;
    3,   6;
    5,   6;
    5.5, 6;
    6,   6] * m;

connectivity = [ ...
     1,  3;
     3,  5;
     5,  7;
     7,  9;
     9, 10;
    10, 11;
    11, 12;
    12, 13;
    13, 14;
    14, 15;
    15,  8;
     8,  6;
     6,  4;
     4,  2;
     9, 16;
    16, 18;
    18, 20;
    20, 22;
    22, 24;
    24, 25;
    25, 26;
    26, 27;
    27, 28;
    28, 29;
    29, 30;
    30, 23;
    23, 21;
    21, 19;
    19, 17;
    17, 15];

nel = size(connectivity, 1);

all_the_beams = [5:10, 20:25];
all_the_cols  = setdiff(1:nel, all_the_beams);

%% ------------------------------------------------------------------------
%  Sections and material
% -------------------------------------------------------------------------
sections.W24x55.d  = 23.57*inch;
sections.W24x55.tw = 0.395*inch;
sections.W24x55.bf = 7.005*inch;
sections.W24x55.tf = 0.505*inch;
sections.W24x55.A  = 16.2*(inch^2);
sections.W24x55.I1 = 1350*(inch^4);
sections.W24x55.I2 = 29.1*(inch^4);

sections.W14x120.d  = 14.48*inch;
sections.W14x120.tw = 0.590*inch;
sections.W14x120.bf = 14.670*inch;
sections.W14x120.tf = 0.940*inch;
sections.W14x120.A  = 35.3*(inch^2);
sections.W14x120.I1 = 1380*(inch^4);
sections.W14x120.I2 = 495*(inch^4);

leaning_col.A  = 9.76e-2*m^2;
leaning_col.I1 = 7.125e-4*m^4;

F_y = 345000*(KN/m^2);   % yield strength
E_0 = 2e8*(KN/m^2);      % initial elastic tangent
eta = 0.01;              % strain-hardening ratio
rho = 7850*(kg/m^3);     % mass density

lumped_mass = 50.97*kg;  % seismic floor mass
P_1 = 500*KN;            % nodal gravity load

%% ------------------------------------------------------------------------
%  Add nodes
% -------------------------------------------------------------------------
for n = 1:size(node_coords, 1)
    ops.node(n, node_coords(n, 1), node_coords(n, 2));
end

% Fictitious leaning-column nodes
ops.node(100, 7.0, 0.0);
ops.node(101, 7.0, 3.0, '-mass', lumped_mass, 1.0e-5, 1.0e-5);
ops.node(102, 7.0, 6.0, '-mass', lumped_mass, 1.0e-5, 1.0e-5);

%% ------------------------------------------------------------------------
%  Materials and sections
% -------------------------------------------------------------------------
ops.uniaxialMaterial('Steel01', 1, F_y, E_0, eta);

% Beam section
ops.section('WFSection2d', 1, 1, ...
    sections.W24x55.d, sections.W24x55.tw, ...
    sections.W24x55.bf, sections.W24x55.tf, 10, 3);

% Column section
ops.section('WFSection2d', 2, 1, ...
    sections.W14x120.d, sections.W14x120.tw, ...
    sections.W14x120.bf, sections.W14x120.tf, 10, 3);

%% ------------------------------------------------------------------------
%  Boundary conditions
% -------------------------------------------------------------------------
baseNodes = [1, 2, 100];
for i = 1:numel(baseNodes)
    ops.fix(baseNodes(i), 1, 1, 0);
end

% Rigid floor diaphragm in the horizontal direction
ops.equalDOF(12, 101, 1);
ops.equalDOF(27, 102, 1);

%% ------------------------------------------------------------------------
%  Transformations and beam integration
% -------------------------------------------------------------------------
ops.geomTransf('Linear', 1);   % beams and MRF members
ops.geomTransf('PDelta', 2);   % leaning columns

ops.beamIntegration('Lobatto', 1, 1, 5);   % beam section integration
ops.beamIntegration('Lobatto', 2, 2, 5);   % column section integration

%% ------------------------------------------------------------------------
%  Add elements
% -------------------------------------------------------------------------
% Beams
for i = 1:numel(all_the_beams)
    e = all_the_beams(i);
    iNode = connectivity(e, 1);
    jNode = connectivity(e, 2);
    ops.element('dispBeamColumn', e, iNode, jNode, 1, 1, ...
        '-cMass', rho*sections.W24x55.A);
end

% Columns
for i = 1:numel(all_the_cols)
    e = all_the_cols(i);
    iNode = connectivity(e, 1);
    jNode = connectivity(e, 2);
    ops.element('dispBeamColumn', e, iNode, jNode, 1, 2, ...
        '-cMass', rho*sections.W14x120.A);
end

% Leaning columns
ops.element('elasticBeamColumn', nel + 1, 100, 101, ...
    leaning_col.A, E_0, leaning_col.I1, 2);
ops.element('elasticBeamColumn', nel + 2, 101, 102, ...
    leaning_col.A, E_0, leaning_col.I1, 2);

1
opsvis.plotModel();
Output
[OpenSeesMatlab] Model summary Nodes: 33 Beam elements: 32 MP constraints: 2
figure_1.png

Rayleigh damping¤

1
2
3
4
5
6
7
8
%% ------------------------------------------------------------------------
%  Rayleigh damping: 2% damping in the first two modes
% -------------------------------------------------------------------------
lambda = ops.eigen('-fullGenLapack', 2);
omega = sqrt(lambda(:));
periods = 2*pi ./ omega;

disp('First two periods:');
Output
First two periods: 
1
disp(periods.');
Output
0.6169 0.1174 
1
2
3
4
5
6
7
8
9
zeta = [0.02; 0.02];
A_damp = [1/(2*omega(1)), omega(1)/2;
          1/(2*omega(2)), omega(2)/2];
coeff = A_damp \ zeta;
a_m = coeff(1);
b_k = coeff(2);

% Rayleigh damping based on initial stiffness
ops.rayleigh(a_m, 0.0, b_k, 0.0);

Gravity analysis¤

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
%% ------------------------------------------------------------------------
%  Gravity analysis
% -------------------------------------------------------------------------
ops.timeSeries('Linear', 1);
ops.pattern('Plain', 1, 1);
ops.load(101, 0.0, -P_1, 0.0);
ops.load(102, 0.0, -P_1, 0.0);

ops.constraints('Transformation');
ops.numberer('RCM');
ops.system('ProfileSPD');
ops.test('NormUnbalance', 1.0e-8, 100);
ops.algorithm('Newton');
ops.integrator('LoadControl', 0.1);
ops.analysis('Static');

okGravity = ops.analyze(10);
fprintf('Gravity analysis return code: %d\n', okGravity);
Output
Gravity analysis return code: 0

Time-history analysis¤

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
%% ------------------------------------------------------------------------
%  Time-history analysis
% -------------------------------------------------------------------------
ops.loadConst('-time', 0.0);
ops.wipeAnalysis();

% Uniform base excitation in the global X direction.
% The acceleration values are already scaled by 3g and written to file.
ops.timeSeries('Path', 2, '-dt', dt, '-values', accel);
ops.pattern('UniformExcitation', 2, 1, '-accel', 2);

ops.constraints('Plain');
ops.numberer('RCM');
ops.system('ProfileSPD');
ops.test('NormUnbalance', 1.0e-8, 100);
ops.algorithm('Newton');
ops.integrator('Newmark', 0.5, 0.25);
ops.analysis('Transient');

nSteps = numel(accel);
time_lst = zeros(nSteps, 1);
d_lst    = zeros(nSteps, 1);

for i = 1:nSteps
    ok = ops.analyze(1, dt);
    if ok ~= 0
        warning('Transient analysis failed at step %d, time %.4f sec, return code %d.', ...
            i, ops.getTime(), ok);
        time_lst = time_lst(1:i-1);
        d_lst = d_lst(1:i-1);
        break;
    end

    time_lst(i) = ops.getTime();
    d_lst(i) = ops.nodeDisp(27, 1);
end

Plot roof displacement response¤

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
%% ------------------------------------------------------------------------
%  Plot roof displacement response
% -------------------------------------------------------------------------
figure('Color', 'w', 'Position', [100, 100, 1000, 360]);
plot(time_lst, d_lst, 'LineWidth', 1.75);
grid on; box on;
xlabel('Time (sec)', 'FontName', 'Times New Roman', 'FontSize', 13);
ylabel('Horizontal displacement (m)', 'FontName', 'Times New Roman', 'FontSize', 13);
set(gca, 'FontName', 'Times New Roman', 'FontSize', 12);
title('Roof horizontal displacement response', 'FontName', 'Times New Roman', 'FontSize', 13);
figure_2.png
1
2
% Uncomment the following line if you want to clear the OpenSees domain at end.
% ops.wipe();