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

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.


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 [1], 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 [1]. 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 [2]; Papanikolaou [1]; Vaz Rodrigues [3]; Simão et al. [4]. Tsao and Hsu [5] presented a general approach for building moment-curvature relations of reinforced concrete sections that was latter named the fiber approach (Spacone et al. [6]; Sfakianakis [7]). Fafitis [8] presented an approach based on Green's Theorem that only requires boundary integrals, evaluated with Gauss quadrature. The work by Kwak and Kim [9] 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. [10] presented a general approach for building the moment curvature relation or reinforced concrete beams, employing Gaussian quadrature and integration cells. Papanikolaou [1] 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 [1]. Vaz Rodrigues [3] 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 [2]; Simão et al. [4] 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. [11] studied the momentcurvature 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 [5]- [7], 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 [1], 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 [8]; Bonet et al. [10]; Papanikolaou [1]; Vaz Rodrigues [3]. In this case, integration of the stress field is made using standard quadrature rules from numerical analysis (Atkinson [12]; Burden and Faires [13]). 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 [12], [13]. 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 [14].
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 momentrotation 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.

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.
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 . 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. The curvature κ at any point of a beam is (Timoshenko [15]; Hibbeler [16]) 2 2 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 For small slopes we have d / dz 0 → and, from Equation 4, κ ≈ θ . This indicates that the curvature is approximately equal to the rotation under small slopes.

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 [5]- [7]. 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 that contains the entire cross section, i.e. that satisfy The integration domain Ω and the cross section Ω are represented in Figure 4. Note that the integration domain is a box that contains the entire cross section. ∈ Ω  =  ∉ Ω  (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 In this case, integration over the cross section Ω is substituted by integration over the entire integration domain Ω . 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 [17], [18]. 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. If n x points are employed in each direction, then the total number of grid points is n = n x 2 . These grid points are collected into the sample set 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 Applying the quadrature rule from Equation 11 gives 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 Substitution of the quadrature weights into Equation 11 gives the quadrature rule  (15) From Equation 15 we observe that the integral is given by multiplication of the area of the integration domain A by the sample mean of the integrand, and thus the approach is similar to Monte Carlo Simulation [17], [18] where σ c (x i , y i ) is the stress on the concrete at position (x i , y i ).

MOMENT-ROTATION RELATIONSHIP
To evaluate the moment-rotation diagram, the resulting moment is evaluated for a set or prescribed rotations θ 0 , θ 1 , θ 2 , ... in a given interval [0, θ f ]. Assuming m rotations increments are applied we have 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 [12], [13]. 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.

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 [19]; CEB-FIP [20]; Bentz [14] 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 where ε c is the strain in the concrete, f c is the resistance in compression, ε c2 = 2.0/1000 and ε cu = 3.5/1000.
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.

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 2 . 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 4 . 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.

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 2 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 4 . 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.  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.

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 2 . 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 4 . 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.

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.
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 2 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. 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 Mu is the maximum bending moment resisted by the cross section and is evaluated in line 4. 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