Skip to content

Two storey steel moment frame with W-sections for displacement-controlled sensitivity analysis¤

1
2
3
4
5
clc; clear;


opsMAT = OpenSeesMatlab();
ops = opsMAT.opensees;

 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
inch = 1.0;
kip = 1.0;
ft = 12 * inch;
ksi = kip / inch^2;


ops.wipe();
ops.model('basic', '-ndm', 2, '-ndf', 3);


ops.node(1, 0.0, 0.0);
ops.node(2, 30 * ft, 0.0);
ops.node(3, 0.0, 15 * ft);
ops.node(4, 30 * ft, 15 * ft);
ops.node(5, 0.0, 30 * ft);
ops.node(6, 30 * ft, 30 * ft);


ops.fix(1, 1, 1, 1);
ops.fix(2, 1, 1, 1);


matTag = 1;
Fy = 50.0 * ksi;
Es = 29000.0 * ksi;
b = 1 / 100;


ops.uniaxialMaterial('Steel01', matTag, Fy, Es, b);


colSecTag = 1;
beamSecTag = 2;


W18x76 = [18.2 * inch, 0.425 * inch, 11.04 * inch, 0.68 * inch];
W14x90 = [14.02 * inch, 0.44 * inch, 14.52 * inch, 0.71 * inch];


ops.section('WFSection2d', colSecTag, matTag, W14x90(1), W14x90(2), W14x90(3), W14x90(4), 20, 4);
ops.section('WFSection2d', beamSecTag, matTag, W18x76(1), W18x76(2), W18x76(3), W18x76(4), 20, 4);


colTransTag = 1;
beamTransTag = 2;
ops.geomTransf('Corotational', colTransTag);
ops.geomTransf('Linear', beamTransTag);


colIntTag = 1;
beamIntTag = 2;
nip = 5;
ops.beamIntegration('Lobatto', colIntTag, colSecTag, nip);
ops.beamIntegration('Lobatto', beamIntTag, beamSecTag, nip);


ops.element('forceBeamColumn', 10, 1, 3, colTransTag, colIntTag, '-mass', 0.0);
ops.element('forceBeamColumn', 11, 3, 5, colTransTag, colIntTag, '-mass', 0.0);
ops.element('forceBeamColumn', 12, 2, 4, colTransTag, colIntTag, '-mass', 0.0);
ops.element('forceBeamColumn', 13, 4, 6, colTransTag, colIntTag, '-mass', 0.0);
ops.element('forceBeamColumn', 14, 3, 4, beamTransTag, beamIntTag, '-mass', 0.0);
ops.element('forceBeamColumn', 15, 5, 6, beamTransTag, beamIntTag, '-mass', 0.0);


ops.timeSeries('Linear', 1);
ops.pattern('Plain', 1, 1);


ops.parameter(1);
ops.parameter(2);
ops.parameter(3);
for ele = [10, 11, 12, 13]
    ops.addToParameter(1, 'element', ele, 'E');
    ops.addToParameter(2, 'element', ele, 'fy');
    ops.addToParameter(3, 'element', ele, 'b');
end


figure;
opsMAT.vis.plotModel();
Output
[OpenSeesMatlab] Model summary Nodes: 6 Beam elements: 6
figure_0.png
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
run_gravity_analysis(10, ops);


ops.load(3, 1/3, 0.0, 0.0);
ops.load(5, 2/3, 0.0, 0.0);


max_disp = 20 * inch;
paramTags = ops.getParamTags();
pushover_output = run_sensitivity_pushover_analysis(5, [1, 2], 1, (1/25) * inch, max_disp, paramTags, false, ops);
Output
Analysis elapsed time is 0.377 seconds. 
 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
figure;
rows = numel(paramTags) * 2;
tiledlayout(rows, 1, 'TileSpacing', 'compact', 'Padding', 'compact');
ParamSym = {'E', 'F_y', 'b'};
ParamVars = [Es, Fy, b];


for s = 1:rows
    nexttile;
    plot_params(max_disp);


    if s <= numel(paramTags)
        sens = pushover_output.(sprintf('sensLambda_%d', s));
        plot(pushover_output.disp, pushover_output.force, '-k', 'LineWidth', 2.0); hold on;
        plot(pushover_output.disp, pushover_output.force + sens * ParamVars(s), '--k', 'LineWidth', 1.5);
        plot(pushover_output.disp, pushover_output.force - sens * ParamVars(s), '-.k', 'LineWidth', 1.5);
        fill([pushover_output.disp; flipud(pushover_output.disp)], ...
             [pushover_output.force + sens * ParamVars(s); flipud(pushover_output.force - sens * ParamVars(s))], ...
             [0.5 0.5 0.5], 'FaceAlpha', 0.15, 'EdgeColor', 'none');
        uistack(findobj(gca, 'Type', 'line'), 'top');
        if s == 1
            title('Displacement Control');
        end
        if s == 2
            ylabel('Base Shear, V_b [kip]');
        else
            ylabel('V_b [kip]');
        end
        legend({'U', ...
                sprintf('U + (\\partial V_b/\\partial %s)%s', ParamSym{s}, ParamSym{s}), ...
                sprintf('U - (\\partial V_b/\\partial %s)%s', ParamSym{s}, ParamSym{s})}, ...
                'Location', 'best');
    else
        k = s - numel(paramTags);
        sens = pushover_output.(sprintf('sensLambda_%d', k));
        plot(pushover_output.disp, sens * ParamVars(k), '-.k', 'LineWidth', 2.0); hold on;
        fill([pushover_output.disp; flipud(pushover_output.disp)], ...
             [sens * ParamVars(k); zeros(size(sens))], ...
             [0.5 0.5 0.5], 'FaceAlpha', 0.15, 'EdgeColor', 'none');
        uistack(findobj(gca, 'Type', 'line'), 'top');
        ylabel(sprintf('(\\partial V_b/\\partial %s)%s [kip]', ParamSym{k}, ParamSym{k}));
        if s == rows
            xlabel('Roof Displacement, U [in]');
        end
        legend({'DDM'}, 'Location', 'best');
    end
end


set(gcf, 'Position', [100, 100, 900, 1200]);
figure_1.png
  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
function run_gravity_analysis(steps, ops)
    ops.wipeAnalysis();
    ops.system('BandGeneral');
    ops.numberer('RCM');
    ops.constraints('Transformation');
    ops.test('NormDispIncr', 1.0E-12, 10, 3);
    ops.algorithm('Newton');
    ops.integrator('LoadControl', 1 / steps);
    ops.analysis('Static');
    ops.analyze(steps);
    ops.loadConst('-time', 0.0);
    ops.wipeAnalysis();
end


function outputs = run_sensitivity_pushover_analysis(ctrlNode, baseNodes, dof, Dincr, max_disp, SensParam, IOflag, ops)
    ops.wipeAnalysis();
    start_time = tic;
    ops.loadConst('-time', 0.0);


    testType = 'NormDispIncr';
    tolInit = 1.0E-8;
    iterInit = 10;
    algorithmType = 'Newton';


    ops.system('BandGeneral');
    ops.constraints('Transformation');
    ops.numberer('RCM');
    ops.test(testType, tolInit, iterInit);
    ops.algorithm(algorithmType);
    ops.integrator('DisplacementControl', ctrlNode, dof, Dincr);
    ops.analysis('Static');
    ops.sensitivityAlgorithm('-computeAtEachStep');


    if IOflag
        fprintf('Single Pushover: Push node %d to %.6g.\n', ctrlNode, max_disp);
    end


    outputs.time = [];
    outputs.disp = [];
    outputs.force = [];


    for i = 1:numel(SensParam)
        outputs.(sprintf('sensLambda_%d', SensParam(i))) = [];
    end


    currentDisp = ops.nodeDisp(ctrlNode, dof);
    ok = 0;


    while ok == 0 && currentDisp < max_disp
        ops.reactions();
        ok = ops.analyze(1);
        currentDisp = ops.nodeDisp(ctrlNode, dof);


        if IOflag
            fprintf('Current displacement ==> %.6g\n', currentDisp);
        end


        if ok ~= 0
            fprintf('\n==> Trying relaxed convergence...\n');
            ops.test(testType, tolInit / 0.01, iterInit * 50);
            ok = ops.analyze(1);
            ops.test(testType, tolInit, iterInit);
        end


        if ok ~= 0
            fprintf('\n==> Trying Newton with initial then current...\n');
            ops.test(testType, tolInit / 0.01, iterInit * 50);
            ops.algorithm('Newton', '-initialThenCurrent');
            ok = ops.analyze(1);
            ops.algorithm(algorithmType);
            ops.test(testType, tolInit, iterInit);
        end


        if ok ~= 0
            fprintf('\n==> Trying ModifiedNewton with initial...\n');
            ops.test(testType, tolInit / 0.01, iterInit * 50);
            ops.algorithm('ModifiedNewton', '-initial');
            ok = ops.analyze(1);
            ops.algorithm(algorithmType);
            ops.test(testType, tolInit, iterInit);
        end


        baseShear = 0.0;
        ops.reactions();
        for node = baseNodes
            baseShear = baseShear - ops.nodeReaction(node, dof);
        end


        outputs.time(end+1, 1) = ops.getTime();
        outputs.disp(end+1, 1) = ops.nodeDisp(ctrlNode, dof);
        outputs.force(end+1, 1) = baseShear;


        for i = 1:numel(SensParam)
            tag = SensParam(i);
            field = sprintf('sensLambda_%d', tag);
            outputs.(field)(end+1, 1) = ops.sensLambda(1, tag);
        end
    end


    elapsed = toc(start_time);
    fprintf('Analysis elapsed time is %.3f seconds.\n', elapsed);
end


function plot_params(max_disp)
    ax = gca;
    ax.Layer = 'top';
    xlim([0.0, max_disp]);
    grid on;
    ax.GridColor = [0.75 0.75 0.75];
    ax.GridAlpha = 0.75;
    ax.LineWidth = 0.75;
    ax.TickDir = 'in';
end