Skip to content

Download MATLAB script

Moment-Curvature Analysis of a Reinforced Concrete Column Section¤

This live script is written as a guided walkthrough for a post-processing workflow. It focuses on retrieving, organizing, and visualizing model or response data after an OpenSees analysis. Read the text cells first, then run each code cell in order so that the variables, model state, and recorded results are available for the later sections.

This example demonstrates how to perform moment-curvature analysis on arbitrary OpenSees sections using anlys.MomentCurvature.

Section geometry : 400 mm x 600 mm rectangular RC column

Concrete cover : 40 mm (to stirrup outer face)

Longitudinal bars : 10 x D25 (rho_l ~ 1.23 %)

Axial load cases : 0, -500, -1000, -2000 kN (compression negative)

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

opsMAT = OpenSeesMatlab();
ops = opsMAT.opensees;
opsMC = opsMAT.anlys.MomentCurvature;  % Get MomentCurvature instance
disp(opsMC);
Output
MomentCurvature with properties: secTag: NaN P: NaN phi: [0x1 double] M: [0x1 double] FiberData: [] cyclePath: [0x1 double] NMMresult: [1x1 struct] 
1
2
ops.wipe();
ops.model("basic", "-ndm", 3, "-ndf", 6);

Define opensees section¤

This section defines the material or section properties. These choices control stiffness, strength, and the nonlinear behavior observed later.

 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
% ======================================================================
%  1.  Section geometry (SI: metres)
% ======================================================================
b     = 0.400;   % section width  [m]
h     = 0.600;   % section height [m]
cover = 0.040;   % clear cover to stirrup face [m]

% Clear distances from centroid to bar layers
d_top = h/2 - cover;   % distance from centroid to top/bottom bar row
d_mid = 0.0;            % distance from centroid to mid-height bar row

% Longitudinal reinforcement
%   Layout  : 3 bars top, 4 bars middle (2 per side), 3 bars bottom
%   Bar dia : D25  (diameter = 25 mm)
db  = 0.025;                 % bar diameter [m]
As1 = pi * db^2 / 4;        % single bar area [m^2]

nTop = 3;   % bars in top row
nBot = 3;   % bars in bottom row
nMid = 4;   % bars in middle rows (2 per side face)

% Total steel area and reinforcement ratio
As_total = (nTop + nBot + nMid) * As1;
rho_l    = As_total / (b * h);
% Fiber discretisation
nFibY_core  = 16;   % fibres along h in confined core
nFibZ_core  = 10;   % fibres along b in confined core
nFibY_cover = 3;    % fibres in top/bottom cover strips
nFibZ_cover = 10;   % fibres along b in cover strips
nFibY_side  = 10;   % fibres along h in side cover strips
nFibZ_side  = 2;    % fibres in side cover strips

  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
% ======================================================================
%  2.  Material models
% ======================================================================
% Material tags
matConc  = 1;   % confined concrete    (Concrete01 / Mander)
matCover = 2;   % unconfined concrete  (Concrete01)
matSteel = 3;   % reinforcing steel    (Steel02 / Menegotto-Pinto)

% --- Confined concrete (Mander model parameters) ---
%   Stirrup: D10 @ 100 mm spacing, 4-leg in both directions
%   (Parameters computed externally from Mander et al. 1988)
fpc_c  = -35.0e6;   % confined peak stress  [Pa]  (negative = compression)
epsc0  = -0.003;   % strain at peak stress
fpcu_c = -28.0e6;   % residual stress  (0.80 * fpc_c)
epscu  = -0.0160;   % ultimate confined strain (Mander formula)

ops.uniaxialMaterial('Concrete02', matConc, fpc_c, epsc0, fpcu_c, epscu);

% --- Unconfined (cover) concrete ---
fpc_u  = -30.0e6;   % cylinder strength  [Pa]
epsc0u = -0.0020;   % strain at peak (plain concrete)
fpcu_u =  -5.0;      % spalls at ultimate (brittle)
epscuu = -0.0040;   % ultimate strain of cover

ops.uniaxialMaterial('Concrete02', matCover, fpc_u, epsc0u, fpcu_u, epscuu);

% --- Reinforcing steel (Menegotto-Pinto / Steel02) ---
Fy  = 400.0e6;   % yield strength     [Pa]
Es  = 200.0e9;   % elastic modulus    [Pa]
b_s = 0.01;     % strain-hardening ratio (post-yield slope / Es)
%   Recommended Menegotto-Pinto constants (Filippou et al. 1983):
R0  = 18;  cR1 = 0.925;  cR2 = 0.15;

ops.uniaxialMaterial('Steel02', matSteel+10, Fy, Es, b_s, R0, cR1, cR2);
ops.uniaxialMaterial('MinMax', matSteel, matSteel+10, '-min', -0.15, '-max', 0.15);

eps_y = Fy / Es;   % yield strain (used later for yield-point search)

% ======================================================================
%  3.  Fiber section  (secTag = 1, bending about local y-axis)
% ======================================================================
%
%  Coordinate convention (OpenSees zeroLengthSection, local axes):
%    y : horizontal (width b)
%    z : vertical   (height h)
%  Bending about y-axis  => M acts about the z-direction on plan,
%  but OpenSees DOF 5 = My.
%
%  Patch corners: (z_min, y_min, z_max, y_max)  [for 'rect' patch]
%
%         z (height)
%         ^
%         |       b = 0.40 m
%    _____|_______________________
%   |     |       cover          |   <- z = +h/2
%   |     |  _______________     |
%   | cov |  |             | cov|
%   |     |  |  confined   |    |   <- core region
%   |     |  |   concrete  |    |
%   |     |  |_____________|    |
%   |     |       cover         |   <- z = -h/2
%   |_____|_____________________|-----> y (width)
%
secTag = 1;

opsMAT.pre.setSectionGeometryRecorder(true);
ops.section('Fiber', secTag, "-GJ", 1E12);

% Core concrete patch (confined)
ops.patch('rect', matConc, nFibY_core, nFibZ_core, ...
    -(h/2 - cover), -(b/2 - cover), ...   % z_min, y_min
     (h/2 - cover),  (b/2 - cover));      % z_max, y_max

% Cover patches (unconfined) - top strip
ops.patch('rect', matCover, nFibY_cover, nFibZ_cover, ...
     (h/2 - cover), -b/2, h/2, b/2);

% Cover patches - bottom strip
ops.patch('rect', matCover, nFibY_cover, nFibZ_cover, ...
    -h/2, -b/2, -(h/2 - cover), b/2);

% Cover patches - left side strip
ops.patch('rect', matCover, nFibY_side, nFibZ_side, ...
    -(h/2 - cover), -b/2, (h/2 - cover), -(b/2 - cover));

% Cover patches - right side strip
ops.patch('rect', matCover, nFibY_side, nFibZ_side, ...
    -(h/2 - cover), (b/2 - cover), (h/2 - cover), b/2);

% Reinforcement layers (layer 'straight': matTag, nBars, As, z1,y1, z2,y2)
%   Top row    (z = +d_top)
ops.layer('straight', matSteel, nTop, As1, ...
     d_top, -(b/2 - cover), d_top, (b/2 - cover));

%   Bottom row (z = -d_top)
ops.layer('straight', matSteel, nBot, As1, ...
    -d_top, -(b/2 - cover), -d_top, (b/2 - cover));

%   Middle rows on left face
ops.layer('straight', matSteel, 2, As1, ...
    -(h/4 - cover/2), -(b/2 - cover), (h/4 - cover/2), -(b/2 - cover));

%   Middle rows on right face
ops.layer('straight', matSteel, 2, As1, ...
    -(h/4 - cover/2), (b/2 - cover), (h/4 - cover/2), (b/2 - cover));
opsMAT.pre.plotSection(secTag);
figure_0.png
1
opsMAT.pre.setSectionGeometryRecorder(false);

Monotonic M-phi for four axial load levels¤

This section applies the actions on the model. The load pattern and scaling determine what response the analysis will try to reproduce.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
N_cases  = [0, -500e3, -1000e3, -2000e3];   % [N]
N_labels = {'N = 0', 'N = -500 kN', 'N = -1000 kN', 'N = -2000 kN'};
colors   = {'#1f77b4','#ff7f0e','#2ca02c','#d62728'};

maxPhi  = 0.15;   % maximum curvature to push to  [1/m]
incrPhi = 5e-4;   % base curvature increment       [1/m]

mc_list = cell(numel(N_cases), 1);

% fprintf('\n--- Monotonic analysis ---\n');
for k = 1:numel(N_cases)
    fprintf('  Case %d: %s\n', k, N_labels{k});
    mc_list{k} = opsMC.new(secTag, N_cases(k));
    mc_list{k}.analyze( ...
        'axis',           'y',   ...
        'maxPhi',         maxPhi, ...
        'incrPhi',        incrPhi, ...
        'limitPeakRatio', 0.80,  ...
        'smartAnalyze',   true);
end
Output
Case 1: N = 0 MomentCurvature: analysis complete. Case 2: N = -500 kN MomentCurvature: analysis complete. Case 3: N = -1000 kN MomentCurvature: analysis complete. Case 4: N = -2000 kN MomentCurvature: analysis complete.

Plot M-phi for all axial load levels

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
fig1 = figure('Position',[60 60 860 520], 'Color','w');
ax1  = axes(fig1);
hold(ax1, 'on');
for k = 1:numel(N_cases)
    [phi_k, M_k] = mc_list{k}.getMPhi();
    plot(ax1, phi_k, M_k/1e3, '-', ...
        'Color', colors{k}, 'LineWidth', 1.8, 'DisplayName', N_labels{k});
end
hold(ax1, 'off');
xlabel(ax1, 'Curvature  \phi  [1/m]',  'FontSize',14);
ylabel(ax1, 'Moment  M  [kN\cdotm]',   'FontSize',14);
title(ax1,  '400\times600 RC Column — Moment-Curvature', 'FontSize',15);
legend(ax1, 'Location','northeast', 'FontSize',12);
grid(ax1, 'on'); box(ax1, 'on');

Yield and ultimate curvature for each case:

1
fprintf('\n--- Limit-state identification ---\n');
Output
--- Limit-state identification --- 
1
2
fprintf('%-18s  %-10s %-12s  %-10s %-12s  %-8s\n', ...
    'Case', 'phi_y[1/m]', 'My[kN.m]', 'phi_u[1/m]', 'Mu[kN.m]', 'mu_phi');
Output
Case phi_y[1/m] My[kN.m] phi_u[1/m] Mu[kN.m] mu_phi 
 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
results = struct();
for k = 1:numel(N_cases)
    mc = mc_list{k};

    % --- Yield point: first yielding of tension steel (matTag=3, eps = +eps_y) ---
    [phi_y, M_y] = mc.getLimitState('matTag', matSteel, 'threshold', eps_y);

    % --- Ultimate point: cover concrete crushing (matTag=2, eps = epscuu) ---
    [phi_u1, M_u1] = mc.getLimitState('matTag', matCover, 'threshold', epscuu);

    % --- Ultimate point: confined concrete crushing (matTag=3, eps = epscu) ---
    [phi_u2, M_u2] = mc.getLimitState('matTag', matConc, 'threshold', epscu);

    % --- Ultimate point: 20 % post-peak drop ---
    [phi_u3, M_u3] = mc.getLimitState('peakDrop', 0.20);

    % Take the more critical (smaller) ultimate curvature
    if phi_u2 < phi_u3
        phi_u = phi_u2; M_u = M_u2;
        ult_ctrl = 'crushing';
    else
        phi_u = phi_u3; M_u = M_u3;
        ult_ctrl = 'peak-drop';
    end

    mu_phi = phi_u / phi_y;   % curvature ductility

    results(k).N     = N_cases(k);
    results(k).phi_y = phi_y;
    results(k).M_y   = M_y;
    results(k).phi_u = phi_u;
    results(k).M_u   = M_u;
    results(k).mu    = mu_phi;
    results(k).ctrl  = ult_ctrl;

    fprintf('%-18s  %10.4f %12.1f  %10.4f %12.1f  %8.2f  [%s]\n', ...
        N_labels{k}, phi_y, M_y/1e3, phi_u, M_u/1e3, mu_phi, ult_ctrl);
end
Output
Warning: matTag 1 strain -0.016 not reached; last step used. Warning: Moment has not dropped 20.0%; last step used. N = 0 0.0090 325.0 0.1500 362.6 16.67 [peak-drop] Warning: matTag 1 strain -0.016 not reached; last step used. Warning: Moment has not dropped 20.0%; last step used. N = -500 kN 0.0095 385.2 0.1500 412.7 15.79 [peak-drop] Warning: matTag 1 strain -0.016 not reached; last step used. Warning: Moment has not dropped 20.0%; last step used. N = -1000 kN 0.0105 446.3 0.1500 448.7 14.29 [peak-drop] Warning: Moment has not dropped 20.0%; last step used. N = -2000 kN 0.0125 543.1 0.1215 470.6 9.72 [crushing] 
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
% Mark yield and ultimate on the M-phi plot
hold(ax1, 'on');
for k = 1:numel(N_cases)
    r = results(k);
    plot(ax1, r.phi_y, r.M_y/1e3, 'o', ...
        'MarkerSize',8, 'MarkerEdgeColor', colors{k}, ...
        'MarkerFaceColor','w', 'LineWidth',1.5, 'HandleVisibility','off');
    plot(ax1, r.phi_u, r.M_u/1e3, 's', ...
        'MarkerSize',8, 'MarkerEdgeColor', colors{k}, ...
        'MarkerFaceColor', colors{k}, 'HandleVisibility','off');
end
% Legend entries for markers only (dummy plots)
plot(ax1, nan, nan, 'ko', 'MarkerSize',8, 'MarkerFaceColor','w', ...
    'LineWidth',1.5, 'DisplayName','Yield point');
plot(ax1, nan, nan, 'ks', 'MarkerSize',8, 'MarkerFaceColor','k', ...
    'DisplayName','Ultimate point');
legend(ax1, 'Location','northeast', 'FontSize',11);
hold(ax1, 'off');
figure_1.png

Bilinear approximation for each axial load case

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
fig2 = figure('Position',[50, 50, 860, 1000], 'Color','w');
for k = 1:numel(N_cases)
    mc = mc_list{k};
    r  = results(k);
    ax_k = subplot(numel(N_cases), 1, k, 'Parent', fig2);

    [phi_eq, M_eq] = mc.bilinearize( ...
        r.phi_y, r.M_y, r.phi_u, 'plot', true, 'ax', ax_k);
    results(k).phi_eq = phi_eq;
    results(k).M_eq   = M_eq;
    title(ax_k, N_labels{k}, 'FontSize',11);
    ax_k.XLabel.FontSize = 10;
    ax_k.YLabel.FontSize = 10;

    fprintf('%-18s  %10.4f %12.1f\n', N_labels{k}, phi_eq, M_eq/1e3);
end
Output
N = 0 0.0096 347.5 N = -500 kN 0.0100 403.8 N = -1000 kN 0.0105 447.1 N = -2000 kN 0.0113 489.9 
1
sgtitle(fig2, 'Equal-Area Bilinear Approximation', 'FontSize',14);
figure_2.png

Plot stress-strain history of fiber points:

1
mc_list{4}.plotFiberResponses();
figure_3.png

Ductility summary bar chart

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
fig3 = figure('Position',[940 60 560 420], 'Color','w');
ax3  = axes(fig3);
mu_vals = [results.mu];
N_kN    = [results.N] / 1e3;
bar(ax3, N_kN, mu_vals, 0.5, 'FaceColor','#4e79a7', 'EdgeColor','k');
xlabel(ax3, 'Axial Load  N  [kN]', 'FontSize',13);
ylabel(ax3, 'Curvature Ductility  \mu_\phi = \phi_u / \phi_y', 'FontSize',13);
title(ax3, 'Ductility vs Axial Load', 'FontSize',14);
text(N_kN, mu_vals + 0.15, arrayfun(@(x) sprintf('%.1f',x), mu_vals, ...
    'UniformOutput',false), 'HorizontalAlignment','center', 'FontSize',11);
ylim(ax3, [0, max(mu_vals)*1.2]);
grid(ax3, 'on'); box(ax3, 'on');
figure_4.png

N-My-Mz interaction surface (strong + weak axis)¤

The following commands carry out this step of the workflow. Run this cell after the previous sections so the required variables and model state already exist.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
% Use a fresh object (P handled internally by buildNMM)
mc_nmm = opsMC.new(secTag, 0);
mc_nmm.buildNMM( ...
    'NList',          linspace(-2400e3, 0, 21)', ...
    'axis1',          'y', ...
    'axis2',          'z', ...
    "capacityMode",   "strain", ...
    "matTag",         matSteel, ...
    "threshold",      eps_y, ...
    'phiMax',         maxPhi,  ...
    'incrPhi',        incrPhi, ...
    'nTheta',         36,      ...
    'limitPeakRatio', 0.80,    ...
    'smartAnalyze',   true,    ...
    'useParallel',    false);
Output
buildNMM: capacityMode=strain, 21 levels x 2 = 42 analyses. N = -2.400e+06 (1/21) N = -2.280e+06 (2/21) N = -2.160e+06 (3/21) N = -2.040e+06 (4/21) N = -1.920e+06 (5/21) N = -1.800e+06 (6/21) N = -1.680e+06 (7/21) N = -1.560e+06 (8/21) N = -1.440e+06 (9/21) N = -1.320e+06 (10/21) N = -1.200e+06 (11/21) N = -1.080e+06 (12/21) N = -9.600e+05 (13/21) N = -8.400e+05 (14/21) N = -7.200e+05 (15/21) N = -6.000e+05 (16/21) N = -4.800e+05 (17/21) N = -3.600e+05 (18/21) N = -2.400e+05 (19/21) N = -1.200e+05 (20/21) N = +0.000e+00 (21/21) buildNMM: done. (42 OpenSees analyses) 
1
mc_nmm.plotNMM('normalize', false, 'alpha', 0.45, 'colormap', 'turbo');
figure_5.png
figure_6.png
figure_7.png

Cyclic analysis at N = -1000 kN¤

This section configures and runs the analysis. The solver, constraints, convergence test, and step size should be read together because they control numerical robustness.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
N_cyc  = -1000e3;
mc_cyc = opsMC.new(secTag, N_cyc);

% Target drift: phi_u at N=-1000 kN, 12 amplitude levels, 2 repeats each
phi_u_ref = results(3).phi_u;   % index 3 = N=-1000 kN case
mc_cyc.setCyclePath(phi_u_ref, 'nCycle', 12, 'nHold', 2);

mc_cyc.analyze( ...
    'axis',           'y',          ...
    'maxPhi',         phi_u_ref,    ...
    'incrPhi',        incrPhi,      ...
    'limitPeakRatio', 0.80,         ...
    'cycleAnalyze',   true,         ...
    'smartAnalyze',   true);
Output
MomentCurvature: analysis complete. 
1
2
3
4
5
6
7
% Plot cyclic M-phi
fig4 = figure('Position',[940 540 620 440], 'Color','w');
ax4  = axes(fig4);
mc_cyc.plotMPhi(ax4);
title(ax4, sprintf('Cyclic M-\phi  (N = %.0f kN)', N_cyc/1e3), 'FontSize',14);
xlabel(ax4, 'Curvature  \phi  [1/m]', 'FontSize',13);
ylabel(ax4, 'Moment  M  [N\cdotm]',   'FontSize',13);
figure_8.png

Fiber stress-strain histories (N = -1000 kN, monotonic)

1
2
mc_cyc.plotFiberResponses();
sgtitle('Fibre Responses — N = -1000 kN', 'FontSize',13);
figure_9.png