MATLAB computational routines for moment-curvature relation of reinforced concrete cross sections

Rotinas computacionais em MATLAB para relação momento-curvatura de seções transversais de concreto armado

ABSTRACT:

In this paper we present a set of MATLAB routines for evaluation of the moment-curvature relationship of reinforced concrete cross sections. This is a topic of major importance for both academic and practical design purposes in the context of structural engineering. The computational routines were developed to be simple, general and flexible. This allows wide practical application and future improvements and modifications. The well-known fibers approach is employed, but an alternative development of the method is also presented. This is interesting from the conceptual point of view. Finally, numerical comparisons are presented to validate the routines.

Keywords:
reinforced concrete; moment-curvature; computational routines; non-linear analysis

RESUMO:

Neste trabalho apresenta-se um conjunto de rotinas do MATLAB para calcular a relação momento-curvatura de seções transversais de concreto armado. Este é um assunto de grande importância para fins acadêmicos e práticos de projeto no contexto da engenharia estrutural. As rotinas computacionais foram desenvolvidas de forma a serem simples, gerais e flexíveis. Isso permite ampla aplicação prática e futuras melhorias e modificações. A conhecida abordagem das lamelas foi utilizada, no entanto também é apresentado um desenvolvimento alternativo do método, que é interessante do ponto de vista conceitual. Por fim, são apresentadas comparações numéricas para validação das rotinas.

Palavras-chave:
concreto armado; momento-curvatura; rotinas computacionais; análise não linear

1 INTRODUCTION

Evaluation of the moment-curvature relation of reinforced concrete cross sections is a problem of significant importance in structural engineering, with wide application to both scientific and design purposes. For this reason, it has been subject of intense investigation over the last decades. The main goal of this work is to present a set of simple, general and flexible computational routines for the problem at hand. The routines were developed in MATLAB programing language and are intended mainly for engineering education and academic purposes.

Several approaches to obtain the moment-curvature relation have been proposed over the years. In general, the most significant difference between them concerns how integration of the stress field is made. Conceptually, evaluation of the moment-curvature relation requires integration of the stress field over the cross section to evaluate force and moment resultants. According to Papanikolaou [¹], most approaches can be broadly classified into methods that employ: a) analytical integration; b) fiber integration and c) numerical integration. An interesting conceptual comparison between these approaches was presented by Papanikolaou [¹]. Here we present only a brief account of works from the past three decades that are relevant to the present study. For an extensive review on the subject the reader may refer to Kwak and Kim [²]; Papanikolaou [¹]; Vaz Rodrigues [³]; Simão et al. [⁴].

Tsao and Hsu [⁵] presented a general approach for building moment-curvature relations of reinforced concrete sections that was latter named the fiber approach (Spacone et al. [⁶]; Sfakianakis [⁷]). Fafitis [⁸] presented an approach based on Green’s Theorem that only requires boundary integrals, evaluated with Gauss quadrature. The work by Kwak and Kim [⁹] employed moment-curvature relations for non-linear structural analysis of reinforced concrete beams, but no emphasis was given to how moment-curvature relations for arbitrary cross sections can be obtained. Bonet et al. [¹⁰] presented a general approach for building the moment curvature relation or reinforced concrete beams, employing Gaussian quadrature and integration cells. Papanikolaou [¹] presented a general approach for arbitrary cross sections that employs an adaptive Gaussian quadrature scheme. A particularly good review on the subject was also presented in the work by Papanikolaou [¹]. Vaz Rodrigues [³] adopted a novel technique based on computer graphics for subdivision of the cross-section domain together with Gaussian quadrature. The works by Kwak and Kim [²]; Simão et al. [⁴] also addressed the application of moment-curvature relations for structural analysis, but no emphasis was given to how moment-curvature relations can be obtained. Finally, even though Liew et al. [¹¹] studied the moment-curvature relation for metal cross section, the work illustrates how these concepts are important in a wider sense in the context of structural engineering.

Unfortunately, analytical integration of the stress field is possible only in the case of very simple cross sections and simple constitutive laws. For this reason, approaches based on analytical integration are generally limited to simple cross sections and polynomial constitutive laws. Otherwise, some computational scheme is necessary.

In the case of the fiber approach [⁵]–[⁷], the cross section is divided into rectangles (fibers). Force and moment resultants are then evaluated considering the contribution of each fiber. The main advantage of this approach is its simplicity from both the conceptual and computational point of view, since it employs the intuitive notion of summation of the fiber’s contribution instead of abstract numerical quadrature schemes. Consequently, it is also very flexible and general. However, some authors argue that the resulting computational routines may be slow [¹], since a large number of fibers may be required to obtain sufficient accuracy in some cases.

To improve computational efficiency, later works focused on improved quadrature schemes, such as the works by Fafitis [⁸]; Bonet et al. [¹⁰]; Papanikolaou [¹]; Vaz Rodrigues [³]. In this case, integration of the stress field is made using standard quadrature rules from numerical analysis (Atkinson [¹²]; Burden and Faires [¹³]). The resulting methods are generally more efficient than the fiber approach, but the use of standard quadrature rules leads to some issues. First, in several cases it is necessary to employ some subdivision of the cross section into integration cells, a procedure that increases the overall complexity of the computational routines. Besides, most available quadrature rules are inefficient for integration of discontinuous functions [¹²], [¹³]. Consequently, numerical quadrature schemes require a fine integration mesh when discontinuous constitutive laws are employed. In some cases, the resulting computational effort may be increased so much that no significant gain in efficiency is observed in comparison to fiber based approaches.

In this work we present a set of general, flexible, and simple computational routines for building the moment-curvature of reinforced concrete sections. The routines are intended for engineering education and academic purposes. Integration of the stress field is made using the fiber approach. An alternative development of the fibers approach is presented, that is interesting from the conceptual point of view. The resulting algorithm does not require meshing of the cross section into integration cells and has no issues with discontinuous constitutive laws. The results are compared to those obtained with the software Response 2000, that is described in detail by Bentz [¹⁴].

This paper is organized as follows. In Section 2 we present a brief review on the behavior of RC cross sections. The quadrature scheme used for evaluation of the axial and moment resultants is described in Section 3. A procedure for obtaining the moment-rotation diagram is described in Section 4. In Section 5 the constitutive laws employed are described. Numerical examples are presented in Section 6 to validate the routines developed. The conclusions are presented in Section 7. The computational routines developed are detailed in the ^{Appendix A}.

2 STRUCTURAL MODEL OF THE CROSS SECTION

We assume the z axis is aligned with the beam axis, the x axis is the direction of the applied bending moment and the y axis is that orthogonal to both x and y. In this case the cross section is defined in the xy plane, with bending moment applied on the direction of x. The axes are represented in Figure 1.

Figure 1
Axes.

If cross sections remain plane after deformation, the strain field on the cross section can be written as

$$\epsilon (y,\theta ,y_c)=(y-y_c)\tan \theta$$ (1)

where θ is the rotation of the cross section and y_c is the position of the neutral axis measured in the direction of y. The resulting strain field is represented in Figure 2. Compression strains and stresses are assumed to be positive. Finally, note that that the strain field varies linearly over the cross-section height y and depends on θ and y_c.

Figure 2
Strain Field.

The axial force and the bending moment resultants are given by

$$N(\theta ,y_c)={\displaystyle {\int }_{\Omega }{\sigma }_{c}dxdy+{\displaystyle \sum _{j=1}^{n_r}{A}_{sj}{\sigma }_{sj}}}$$ (2)

$$M(\theta ,y_c)={\displaystyle {\int }_{\Omega }(y_c-y){\sigma }_{c}dxdy+{\displaystyle \sum _{j=1}^{n_r}(y_c-y_{sj})A_{sj}{\sigma }_{sj}}}$$ (3)

where Ω represents the cross section domain (see Figure 3), σ_c is the stress field in the concrete, A_sj are the cross sectional areas of the rebars (reinforcement bars), σ_sj are the stresses in the rebars, y_sj are the position of the rebars and n_r is the number of rebars. Note that both N and M depend on θ and y_c, since these two variables define the strain field and, consequently, the stress field.

Figure 3
Cross section.

The curvature κ at any point of a beam is (Timoshenko [¹⁵]; Hibbeler [¹⁶])

$$\kappa =\frac{1}{\rho }=\frac{d^2\upsilon /d{z}^2}{{[1+{(d\upsilon /dz)}^2]}^{3/2}}$$ (4)

where ρ is the curvature radius, υ is the transversal displacement and z is the longitudinal axis. If sections orthogonal to the neutral axis remains so after displacements take place, we can write

$$\theta =\frac{d^2\upsilon }{d{z}^2}$$ (5)

For small slopes we have $\text{d}/\text{dz}\to 0$ and, from Equation 4, $\kappa \approx \theta$. This indicates that the curvature is approximately equal to the rotation under small slopes.

3 QUADRATURE RULE: FIBERS APPROACH

For given values of θ and y_c, the contributions of the rebars for N and M (Equations 2 and 3) can be easily evaluated by summation. Evaluation of the contribution from concrete, on the other hand, requires integration of the stress field over the cross section. As discussed in the introduction, this is the main difficulty in evaluation of the moment-curvature relation for arbitrary cross sections. If the geometry of the cross section and the constitutive law are simple, the required integrals can be evaluated analytically. However, when the cross-section geometry is complex, the above integrals require numerical quadrature.

The fiber approach is employed for integration of the stress field [⁵]–[⁷]. However, we present an alternative development of the scheme that puts in evidence some interesting conceptual properties. For this purpose, we first define a rectangular integration domain

$$\overline{\Omega }=\left\{\left(x,y\right)\in {ℝ}^2|x_l\le x\le x_u,y_l\le y\le y_u\right\}$$ (6)

that contains the entire cross section, i.e. that satisfy

$$\Omega \subseteq \overline{\Omega }$$.(7)

The integration domain $\overline{\Omega }$ and the cross section Ω are represented in Figure 4. Note that the integration domain is a box that contains the entire cross section.

Figure 4
Cross section Ω and integration domain $\overline{\Omega }$.

The geometry of the cross section is then represented using an indicator function defined as

$$I(x,y)=\left\{\begin{array}{c}\mathrm{1,}(x,y)\in \Omega \\ \mathrm{0,}(x,y)\notin \Omega \end{array}\right.$$ (8)

The indicator function simply returns 1 if the point (x, y) is inside the cross-section domain Ω and 0 otherwise.

Integration of some function f (x, y) over the cross section Ω can then be written as

$${\displaystyle {\int }_{\Omega }f\left(x,y\right)dxdy}={\displaystyle {\int }_{\overline{\Omega }}I\left(x,y\right)f\left(x,y\right)dxdy}$$ (9)

In this case, integration over the cross section Ω is substituted by integration over the entire integration domain $\overline{\Omega }$. Multiplication by the indicator function is employed to take into account that only points inside Ω should be considered.

We then generate a sample of points with uniform distribution on the integration domain Ω. Several schemes can be employed for this purpose, such as random sample generation [¹⁷], [¹⁸]. A uniform grid of equally spaced points in each direction is employed, as illustrated in Figure 5. This uniform grid is built as described in the computational routines in ^{Appendix A}.

Figure 5
Uniform quadrature grid.

If n_x points are employed in each direction, then the total number of grid points is n = n_x². These grid points are collected into the sample set

$$S=\left\{\left(x_1,y_1\right),\left(x_2,y_2\right)\mathrm{,...,}\left(x_n,y_n\right)\right\}$$ (10)

where n is the sample size. The quadrature rule for the integral from Equation 9 can then be written as

$${\displaystyle {\int }_{\Omega }f(x,y)dxdy={\displaystyle {\int }_{\overline{\Omega }}I(x,y)f(x,y)dxdy\approx {\displaystyle \sum _{i=1}^n{w}_{i}I(x_i,y_i)f(x_i,y_i)}}}$$ (11)

where w_i are quadrature weights and the grid points (x_i, y_i) are taken as the quadrature nodes.

To obtain the quadrature weights, note that the area of the integration domain is given by

$$\overline{A}=(x_u-x_l)(y_u-y_l)$$ (12)

Applying the quadrature rule from Equation 11 gives

$$\overline{A}={\displaystyle {\int }_{\overline{\Omega }}1dxdy\approx {\displaystyle \sum _{i=1}^n{w}_i}}$$ (13)

Since the sample points follow a uniform distribution, we can assume that all quadrature points have the same weight w_i and from Equation 13

$$w_i=\frac{\overline{A}}{n}$$ (14)

Substitution of the quadrature weights into Equation 11 gives the quadrature rule

$${\displaystyle {\int }_{\Omega }f(x,y)dxdy\approx \overline{A}\left[\frac{1}{n}{\displaystyle \sum _{i=1}^{n}I(x_i,y_i)f(x_i,y_i)}\right]}$$ (15)

From Equation 15 we observe that the integral is given by multiplication of the area of the integration domain $\overline{A}$ by the sample mean of the integrand, and thus the approach is similar to Monte Carlo Simulation [¹⁷], [¹⁸] from the computational point of view. Finally, multiplication by I (x_i, y_i) in Equation 15 implies that only quadrature points inside the cross section are actually considered. Consequently, the summation from Equation 15 can be evaluated by discarding quadrature points outside the cross section.

Employment of the quadrature rule from Equation 15 to Equations 2 and 3 gives

$$N(\theta ,y_c)\approx \overline{A}\left[\frac{1}{n}{\displaystyle \sum _{i=1}^{n}I(x_i,y_i){\sigma }_c(x_i,y_i)}\right]+{\displaystyle \sum _{j=1}^{n_r}{A}_{sj}{\sigma }_{sj}}$$ (16)

$$M(\theta ,y_c)\approx \overline{A}\left[\frac{1}{n}{\displaystyle \sum _{i=1}^{n}I(x_i,y_i)(y_c-y_i){\sigma }_c(x_i,y_i)}\right]+{\displaystyle \sum _{j=1}^{n_r}(y_c-y_{sj})A_{sj}{\sigma }_{sj}}$$ (17)

where σ_c (x_i, y_i) is the stress on the concrete at position (x_i, y_i).

4 MOMENT-ROTATION RELATIONSHIP

To evaluate the moment-rotation diagram, the resulting moment is evaluated for a set or prescribed rotations θ₀, θ₁, θ₂, ... in a given interval [0, θ_f]. Assuming m rotations increments are applied we have

$$\left\{\begin{array}{l}{\theta }_0=0\hfill \\ {\theta }_1={\theta }_0+\Delta \theta \hfill \\ {\theta }_2={\theta }_1+\Delta \theta \hfill \\ ⋮\hfill \\ {\theta }_m={\theta }_{m-1}+\Delta \theta \hfill \end{array}\right.$$ (18)

where

$$\Delta \theta =\frac{{\theta }_f}{m}$$ (19)

is the rotation increment.

For a given rotation θ_i, the position of the neutral axis y_c is the root of the equation

$$N({\theta }_i,y_c)=0$$ (20)

where we consider, without loss of generality, that no axial forces are applied. The root of the above equation is obtained with the Bisection Method [¹²], [¹³]. The axial resultant is evaluated using Equation 16. After the position of the neutral axis y_c is found, the resulting moment can be evaluated using Equation 17.

5 CONSTITUTIVE LAWS

The proposed approach is general and can be employed together with any constitutive law for concrete and steel. In this work very simple constitutive laws are employed in order to focus on the approach rather than on the mechanical properties of the materials. More refined constitutive laws can be found in Chen [¹⁹]; CEB-FIP [²⁰]; Bentz [¹⁴] and other more recent references on the subject.

In the case of concrete, the resistance in tension is neglected and the parabolic model is employed in compression. The resulting stress-strain relationship is

$${\sigma }_c({\epsilon }_c)=\left\{\begin{array}{ll}0\hfill & ,{\epsilon }_c<0\hfill \\ {\mathrm{f}}_c\left[1-{\left(1-\frac{{\epsilon }_c}{{\epsilon }_{c2}}\right)}^2\right]\hfill & \mathrm{,0}\le {\epsilon }_c<{\epsilon }_{c2}\hfill \\ {\mathrm{f}}_c\hfill & ,{\epsilon }_{c2}\le {\epsilon }_c<{\epsilon }_{cu}\hfill \\ 0\hfill & ,{\epsilon }_c>{\epsilon }_{cu}\hfill \end{array}\right.$$ (21)

where ε_c is the strain in the concrete, f_c is the resistance in compression, ε_c2 = 2.0/1000 and ε_cu = 3.5/1000.

For steel, an elastic-perfectly plastic model is employed. We also assume perfect bond between steel and concrete, resulting in the stress-strain relationship

$${\sigma }_s({\epsilon }_s)=\left\{\begin{array}{ll}{\epsilon }_s{E}_s\hfill & {\mathrm{,0}}_{}{}_{}\le \left|{\epsilon }_s\right|\le {\epsilon }_{sy}\hfill \\ {\mathrm{f}}_y{\epsilon }_s/\left|{\epsilon }_s\right|\hfill & ,{\epsilon }_{sy}\le \left|{\epsilon }_s\right|\le {\epsilon }_{su}\hfill \\ 0\hfill & ,\left|{\epsilon }_s\right|>{\epsilon }_{su}\hfill \end{array}\right.$$ (22)

where ε_s is the strain in the steel, f_y is the yielding stress, ε_su=10/1000, ε_sy=f_y /E_s and E_s=210GPa. In the above equation the quantity ε_s /|ε_s| is only employed to consider the sign of the strain.

Note that positive signs are employed for compression. Besides, the resistances must be weighted by partial factors when the model is employed for design purposes. The stress-strain relation for concrete with f_c = 30MPa and for steel with f_y = 400MPa are presented in Figures 6 and 7, for illustration purposes.

Figure 6
Stress-strain relation for concrete with f_c = 30MPa.

Figure 7
Stress-strain relation for steel with f_y = 400MPa.

6 EXAMPLES

6.1 Rectangular cross section

In the first example we study the cross section from Figure 8. The cross section is a rectangle with width b=200mm and height h=500mm. The material properties were taken as f_c=40MPa and f_y=500MPa. The centroids of the three bottom and two top rebars are positioned at 40mm and 460mm from the bottom of the cross section, respectively. Each rebar has a cross sectional area equal to A_sj=123mm². The limits of the integration domain were set as x_l=y_l=0mm, x_u=200mm and y_u=500mm. The sample size employed for numerical quadrature was n=10⁴.

Figure 8
Rectangular cross section.

The results obtained with the approach proposed in the present work are compared to the results obtained with Response 2000 in Figure 9. As can be observed, the results obtained are very similar, which validates the proposed approach. The minor differences observed are mainly due to divergences in the constitutive laws, as there are only a few predefined constitutive laws that can be chosen in Response 2000 and some parameters cannot be modified by the user.

Figure 9
Moment-curvature of rectangular cross section.

6.2 Ellipsoidal cross section

In the second example we study the cross section from Figure 10. The input data of this example is detailed in the ^{Appendix A}. The cross section is an ellipse with radius a=100mm and b=200mm about axis x and y, respectively. The material properties were taken as f_c=30MPa and f_y=400MPa. The rebars have cross sectional areas A_s=123mm² with centroids positioned at (−30, −150) mm, (+30, −150) mm, (−70, 0) mm, (+70, 0) mm and (0, +150) mm. The limits of the integration domain were set as x_l=y_l=−200mm and x_u=y_u=+200mm. The sample size employed for numerical quadrature was n=10⁴. The sample points located inside the cross section and the rebars are shown in Figure 11. We observe that the sample reproduces the geometry of the cross section, as expected.

Figure 10
Ellipsoidal cross section.

Figure 11
Ellipsoidal cross section sample points.

The moment-rotation diagram is presented in Figure 12. The results are compared to those obtained by Response 2000. Again, the results agree, indicating that the proposed approach can obtain accurate results.

Figure 12
Moment-rotation diagram for the ellipsoidal cross section.

6.3 Bridge cross section

In the last example we study the cross section from Figure 13, where the dimensions are presented in mm. The coordinates of the vertices of the cross section are detailed in the ^{Appendix A}. The material properties were taken as f_c=50MPa and f_y=500MPa. The rebars have cross sectional areas A_s=491mm². The reinforcement layers are distant 50 mm, 100 mm, 400 mm, 575 mm, 750 mm, 925 mm, 1100 mm and 1450 mm from the bottom of the cross section, respectively. The limits of the integration domain were set as x_l=−1300mm, x_u=1300mm, y_l=0mm and y_u=1500mm. The sample size employed for numerical quadrature was n=10⁴. The sample points located inside the cross section and the rebars are shown in Figure 14. The moment-rotation diagram is presented in Figure 15. We observe that the results are very similar to those obtained with Response 2000.

Figure 13
Bridge cross section.

Figure 14
Bridge cross section sample points.

Figure 15
Moment-rotation diagram for the bridge cross section.

7 CONCLUDING REMARKS

In this work we presented a set of MATLAB computational routines for obtaining the moment-curvature relation of reinforced concrete cross sections. The routines are simple, general, and flexible, allowing wide practical application, future improvements and modifications. The routines were validated by numerical examples. An alternative development of the well-known fibers approach is also demonstrated, that is interesting from the conceptual point of view.

We emphasize that several improvements and modifications can be easily implemented in the presented routines. Axial forces can be included by modification of Equation 20 and more advanced constitutive laws can be modeled in Equations 21 and 22. It is also possible to replace Bisection Method for finding the position of the neutral axis by more advanced root finding methods. Finally, biaxial bending can be addressed by rotation of the cross-section axes.

These routines are valuable for engineering education and academic purposes. From the academic point of view, the routines are valuable since they are open source. There are many computer programs that are able to obtain the moment-curvature relation of reinforced concrete cross sections. However, no open source codes on the subject are known by the authors. Researchers wishing to test some novel constitutive law in the context of beams, for example, will likely need to implement a new computational code or employ some commercial software, what is not always desirable. The provided routines try to fill this gap on the field.

It is worth mentioning that the proposed method differs from other methods (ex: Gauss Numerical Integration) because in addition to being less computationally efficient, it adequately solves continuous functions with discontinuities. Hence, the method proposed in this work aims to solve general problems in which the constitutive law is continuous but may contain some discontinuity.

The computational routines developed in this paper can be improved for instance by adopting a different number of integration points in x and y directions, using more points on one direction than on the other. However, since our aim is to implement a didactic and elementary strategy that allows for improvements, such improvements should be carried out in future works.

From the education point of view, the routines can be employed in a wide range of situations. First, they can be used to evaluate the ultimate resistance and the moment-curvature relation in several interesting cases. The routines can also be employed to show how different constitutive laws will impact the results obtained, for example. Finally, more advanced students may modify the routines to take into improvements and modifications.

APPENDIX A COMPUTATIONAL ROUTINES

In this section we present the computational routines in MATLAB programming language. Interested readers can contact the corresponding author to obtain these routines. The flowchart that explains the operation of the Moment-Curvature routine is shown in Figure 16. In addition, the authors are willing to place the computer routines on the FileExchange platform of Matlab or on another open source platform, besides keeping them in the paper.

Figure 16
Moment-Curvature routines flowchart.

The input data is defined in function data. The input data for the example of the ellipsoidal cross section is presented in Algorithm 1. In lines 2-3 the resistances of concrete and steel are defined.

Algorithm 1 data.m (ellipsoidal)

1	function [ fc, fy, As, xl, xu, yl, yu, nx, tol, thetaf, m] = data
2	fc = 30;
3	fy = 400;
4	As = [123 -30 -150;
5	123 +30 -150;
6	123 -70 0;
7	123 +70 0;
8	123 0 +150];
9	xl = -200;
10	xu = 200;
11	yl = -200;
12	yu = 200;
13	nx = 100;
14	tol = 1e-3;
15	thetaf = 0.08 / 1000;
16	m = 100;
17	end

In lines 4-8 the rebars are described, using a matrix where each line represents a single rebar. The columns of matrix As contain the cross sectional areas, the x coordinate and the y coordinate of the rebars, respectively. The data at line 5 then represents a rebar with cross sectional area 123mm² positioned at x = +30mm y = -150mm, for example.

The limits of the rectangular integration domain are defined in lines 9-12. In line 13 we define the number of grid points in each direction. In line 14 the stopping criterion of the Bisection Method for finding the position of the neutral axis is defined. The total rotation to be applied to the section and the number of increments are defined in lines 15-16.

The indicator function of the cross section is defined in function geometry, that receives as input the position of some point and returns 1 if the point is inside the cross section and 0 otherwise. The function geometry for the examples studied is presented in Algorithms 2-4.

Algorithm 2 geometry.m (rectangular)

1	function I = geometry (x)
2	xv = [ 0 200 200 0 ] ;
3	yv = [ 0 0 600 600 ] ;
4	I = inpolygon (x (1), x (2), xv, yv) ;
5	end

Algorithm 3 geometry.m (ellipsoidal)

1	function I = geometry (x)
2	a = 100 ;
3	b = 200 ;
4	f = (x (1)/ a)^2 + (x (2)/ b)^2 - 1 ;
5	I = f <0;
6	end

Algorithm 4 geometry.m (bridge)

1	function I = geometry (x)
2	xv = [-300 -300 -150 -150 -300 -1300 -1300 1300 1300 300 150 150 300 300 ] ;
3	yv = [ 0 200 400 1100 1200 1200 1500 1500 1200 1200 1100 400 200 0 ] ;
4	I = inpolygon (x (1), x (2), xv, yv) ;
5	end

In these routines, the x and y coordinates of some point are stored in the components x(1) and x(2) of vector x, respectively. In Algorithms 2 and 4, the indicator function is evaluated using the native MATLAB function inpolygon. This function returns 1 when some point with coordinates x(1) and x(2) is inside the polygon defined by vertices xv and yv, that indicate the x and y coordinates of the vertices, respectively. In Algorithm 3, the indicator function is taken as 1 when the ellipse equation gives f < 0.

The moment-rotation diagram can be plotted using the commands from Algorithm 5. In line 1 the input data is read using data. In line 2 the function diagram is called with the data from the input file. The diagram is plotted in line 3. The ultimate moment M_u is the maximum bending moment resisted by the cross section and is evaluated in line 4.

Algorithm 5 Commands for plotting the moment-curvature diagram

1	[ fc, fy, As, xl, xu, yl, yu, nx, tol, thetaf, m ]=data ;
2	[ theta, M ] = diagram(fc, fy, As, xl, xu, yl, yu, nx, tol, thetaf, m) ;
3	plot (theta, M)
4	max(M)

The moment-rotation diagram is built with the function diagram from Algorithm 6. In lines 2-4 the quadrature weights are obtained. The number of rebars is evaluated at line 5. The uniform quadrature grid is built in lines 6-16. In lines 6-14 a uniform grid is built in the square [0, 0] × [1, 1]. Coordinate transformation to the rectangular integration domain is then made in lines 15-16. The resulting grid is stored in matrix S, where each line stores a given point and each column stores a given coordinate. The indicator function of the sample points is evaluated in lines 17-19, using the function geometry defined previously. Sample points outside the cross section are then deleted from the sample and the new sample size ns is evaluated in lines 20-21. In lines 22-26 the moment-rotation diagram is built, by finding the resulting bending moment for prescribed rotations with function momrot.

Algorithm 6 diagram.m

1	function [ theta, M ] = diagram (fc, fy, As, xl, xu, yl, yu, nx, tol, thetaf, m)
2	A0 = (xu - xl) * (yu - yl) ;
3	n = nx ^ 2 ;
4	w = A0 / n ;
5	nr = size (As, 1) ;
6	du = 1 / nx ;
7	u = (du / 2): du: (1-du / 2) ;
8	k = 0 ;
9	for i =1: nx
10	for j = 1: nx
11	k = k + 1 ;
12	S(k,:) = [ u(i) u(j) ] ;
13	end
14	end
15	S (:, 1) = xl + S (:, 1) * (xu - xl) ;
16	S (:, 2) = yl + S (:, 2) * (yu - yl) ;
17	for i =1:n
18	I (i) = geometry (S(i,:)) ;
19	end
20	S(I ==0,:) = [ ] ;
21	ns = size (S, 1) ;
22	dtheta = thetaf / m;
23	for i =2: m
24	theta (i) = i * dtheta ;
25	M(i) = momrot (theta(i), ns, S, w, nr, As, fc, fy, yl, yu, tol) ;
26	end
27	end
28	function M = momrot (theta, ns, S, w, nr, As, fc, fy, yl, yu, tol)
29	a = yl ;
30	b = yu ;
31	Na = resultants (theta, ns, S, w, nr, As, a, fc, fy) ;
32	h0 = b - a;
33	it = ceil (log2 (h0 / tol)) ;
34	for i = 1: it
35	c = 0.5 * (a + b) ;
36	Nc = resultants (theta, ns, S, w, nr, As, c, fc, fy) ;
37	if Na * Nc < 0
38	b = c ;
39	else
40	a = c ;
41	Na = Nc ;
42	end
43	end
44	yc = 0.5 * (a + b) ;
45	[ N, M ] = resultants (theta, ns, S, w, nr, As, yc, fc, fy) ;
46	end
47	function [ N, M ] = resultants (theta, ns, S, w, nr, As, yc, fc, fy)
48	N = 0 ;
49	M = 0 ;
50	for i =1: ns
51	y = S(i, 2) ;
52	e = tan(theta) * (y - yc) ;
53	sigma = sigmac (e, fc) ;
54	N = N + w * sigma ;
55	M = M + w * sigma * (y - yc) ;
56	end
57	for i =1: nr
58	y = As (i, 3) ;
59	e = tan(theta) * (y - yc) ;
60	sigma = sigmas (e, fy) ;
61	N = N + sigma * As (i, 1) ;
62	M = M + sigma * As (i, 1) * (y - yc) ;
63	end
64	end
65	function sigma= sigmac (e, fc)
66	ec2 = 2.0 / 1000 ;
67	ecu = 3.5 / 1000 ;
68	n = 2 ;
69	if e < 0
70	sigma = 0 ;
71	else
72	if e > ecu
73	sigma = 0 ;
74	elseif e > ec2
75	sigma = fc ;
76	else
77	sigma = fc * (1 - (1 - e / ec2)^n) ;
78	end
79	end
80	end
81	function sigma = sigmas (e, fy)
82	Es = 210000;
83	esy = fy / Es ;
84	esu = 10 / 1000 ;
85	if abs (e) > esu
86	sigma = 0 ;
87	elseif abs (e) > esy
88	sigma = fy * e / abs (e) ;
89	else
90	sigma = e * Es ;
91	end
92	end

Evaluation of the resulting moment for a give rotation is made using the function momrot, defined at line 28. In lines 29-44 the position of the neutral axis is obtained using the Bisection Method in the interval [y_l, y_u], where function resultants evaluate the axial force resultant. The number of iterations to be performed is evaluated in line 33, considering the size of the original interval h0 and the stopping criterion tol, since each iteration of the method halves the search interval (Atkinson [¹²]; Burden and Faires [¹³]). After the position of the neutral axis is found, the resulting moment is evaluated in line 45.

The axial and moment resultants are evaluated with the function resultants, defined at line 47. Lines 50-56 are related to quadrature of the stress in the concrete, while lines 57-63 are related to the contribution of the axial forces on the rebars. The concrete and steel constitutive laws are modeled in the functions sigmac and sigmas, respectively, in lines 65-92.

• Financial support: None.
• How to cite: G. F. Melo, A. J. Torii, E. M. Medeiros, and A. K. L. Kzam, “MATLAB computational routines for moment-curvature relation of reinforced concrete cross sections,” Rev. IBRACON Estrut. Mater., vol. 14, no. 2, e14202, 2021, https://doi.org/10.1590/S1983-41952021000200002

REFERENCES

• ¹
  V. K. Papanikolaou, "Analysis of arbitrary composite sections in biaxial bending and axial load," Comput. Struc., vol. 98-99, pp. 33–54, 2012, http://dx.doi.org/10.1016/j.compstruc.2012.02.004
  » http://dx.doi.org/10.1016/j.compstruc.2012.02.004
• ²
  H.-G. Kwak and S.-P. Kim, "Simplified monotonic moment-curvature relation considering fixed-end rotation and axial force effect," Eng. Struct., vol. 32, no. 1, pp. 69–79, 2010, http://dx.doi.org/10.1016/j.engstruct.2009.08.017
  » http://dx.doi.org/10.1016/j.engstruct.2009.08.017
• ³
  R. Vaz Rodrigues, "A new technique for ultimate limit state design of arbitrary shape RC sections under biaxial bending," Eng. Struct., vol. 104, pp. 1–17, 2015, http://dx.doi.org/10.1016/j.engstruct.2015.09.016
  » http://dx.doi.org/10.1016/j.engstruct.2015.09.016
• ⁴
  P. D. Simão, H. Barros, C. C. Ferreira, and T. Marques, "Closed-form moment-curvature relations for reinforced concrete cross sections under bending moment and axial force," Eng. Struct., vol. 129, pp. 67–80, 2016, http://dx.doi.org/10.1016/j.engstruct.2016.09.033
  » http://dx.doi.org/10.1016/j.engstruct.2016.09.033
• ⁵
  W. H. Tsao and C. T. T. Hsu, "A nonlinear computer analysis of biaxially loaded L-shaped slender reinforced concrete columns," Comput. Struc., vol. 49, no. 4, pp. 579–588, 1993, http://dx.doi.org/10.1016/0045-7949(93)90062-I
  » http://dx.doi.org/10.1016/0045-7949(93)90062-I
• ⁶
  E. Spacone, F. C. Filippou, and F. F. Taucer, "Fibre beam-column model for non-linear analysis of R/C frames: part I. Formulation," Earthquake Eng. Struct. Dynam., vol. 25, no. 7, pp. 711–725, 1996, http://dx.doi.org/10.1002/(SICI)1096-9845(199607)25:7<711:AID-EQE576>3.0.CO;2-9
  » http://dx.doi.org/10.1002/(SICI)1096-9845(199607)25:7<711:AID-EQE576>3.0.CO;2-9
• ⁷
  M. G. Sfakianakis, "Biaxial bending with axial force of reinforced, composite and repaired concrete sections of arbitrary shape by fiber model and computer graphics," Adv. Eng. Softw., vol. 33, no. 4, pp. 227–242, 2002, http://dx.doi.org/10.1016/S0965-9978(02)00002-9
  » http://dx.doi.org/10.1016/S0965-9978(02)00002-9
• ⁸
  A. Fafitis, "Interaction surfaces of reinforced-concrete sections in biaxial bending," J. Struct. Eng., vol. 127, no. 7, pp. 840–846, 2001., http://dx.doi.org/10.1061/(ASCE)0733-9445(2001)127:7(840)
  » http://dx.doi.org/10.1061/(ASCE)0733-9445(2001)127:7(840)
• ⁹
  H.-G. Kwak and S.-P. Kim, "Nonlinear analysis of RC beams based on moment–curvature relation," Comput. Struc., vol. 80, no. 7-8, pp. 615–628, 2002, http://dx.doi.org/10.1016/S0045-7949(02)00030-5
  » http://dx.doi.org/10.1016/S0045-7949(02)00030-5
• ¹⁰
  J. L. Bonet, M. L. Romero, P. F. Miguel, and M. A. Fernandez, "A fast stress integration algorithm for reinforced concrete sections with axial loads and biaxial bending," Comput. Struc., vol. 82, no. 2-3, pp. 213–225, 2004, http://dx.doi.org/10.1016/j.compstruc.2003.10.009
  » http://dx.doi.org/10.1016/j.compstruc.2003.10.009
• ¹¹
  A. Liew, L. Gardner, and P. Block, "Moment-curvature-thrust relationships for beam-columns," Structures, vol. 11, pp. 146–154, 2017, http://dx.doi.org/10.1016/j.istruc.2017.05.005
  » http://dx.doi.org/10.1016/j.istruc.2017.05.005
• ¹²
  K. E. Atkinson, An Introduction to Numerical Analysis, 2nd ed. New York, NY: Wiley, 1989.
• ¹³
  R. L. Burden and J. D. Faires, Numerical Analysis, 9th ed. Belmont, CA: Thomson Brooks/Cole, 2011.
• ¹⁴
  E. C. Bentz, “Sectional analysis of reinforced concrete members,” Ph.D. thesis, Graduate Dept. Civil Eng., Univ. Toronto, Toronto, Canada, 2000.
• ¹⁵
  S. Timoshenko, Strengh of Materials, 2nd ed. New York, NY: D. Van Nostrand Company, Inc., 1940.
• ¹⁶
  R. C. Hibbeler, Mechanics of Materials, 8th ed. Upper Saddle River, NJ: Prentice Hall, 2011.
• ¹⁷
  S. M. Ross, Simulation, 4th ed. New York, NY: Elsevier, 2006.
• ¹⁸
  R. Y. Rubinstein and D. P. Kroese, Simulation and the Monte Carlo Method, 2nd ed. Hoboken, NJ: Wiley-Interscience, 2008.
• ¹⁹
  W. F. Chen, Plasticity in Reinforced Concrete. New York, NY: McGraw-Hill, 1982.
• ²⁰
  Comité Euro-International du Béton, CEB-FIP model code 1990. Lausanne, Switzerland: Thomas Telford, 1991.

Edited by

Editors: Osvaldo Luís Manzoli, José Luiz Antunes de Oliveira e Sousa, Guilherme Aris Parsekian.

Publication Dates

• Publication in this collection
  05 Feb 2021
• Date of issue
  2021

History

• Received
  23 May 2019
• Accepted
  18 June 2020