A design tool for the analysis of reinforced and prestressed concrete columns

This article aims to describe an automatic calculation algorithm developed for the analysis of reinforced concrete and / or prestressed concrete columns. The article describes two analysis types: sectional analysis and analysis on the element. The sectional analysis consists in the generation of three characteristic curves: bending moment - curvature (M-1/r), axial force - bending moment (N-M) and bending moment – bending stiffness diagrams. The analysis on the element provides the height – displacement (H-w) diagram, obtained by using the finite element method (FEM). For the two analysis, the calculation approach and the implementation of the algorithm is briefly described by the mathematical equations used. To ease the usage of the algorithm an application was developed in MATLAB appDesigner. For validating the algorithm, the results were compared with those obtained from experimental tests in the laboratory but also with numerical simulations using commercial computational software based on FEM. MATLAB programming language was used to implement the algorithm, while results of the analysis proved to be conservative when compared to the experimental ones. Moreover, the bending moment resistance for prestressed concrete columns using bonded tendons is better estimated with at least 10% when compared with a commercial software, based on fiber-based cross-section nonlinear analysis.


Introduction
Facing a fast urban development in Romania, the need for buildings in industrial regions increased significantly.The most widely used buildings in these regions are single-story warehouses with long slender columns fully fixed in foundation and hinged at top.The main problem of such structures is the significant lateral displacement at the roof level caused by horizontal actions such as wind or earthquake [1] [2].This large drift can be diminished by increasing the lateral stiffness of the superstructure [3].Such improvement can be reached by the using columns with larger cross-section or by using materials with higher modulus of elasticity.However, it has been observed that by switching from reinforced concrete to prestressed concrete we get the same effect but with significant material savings [4] [5].
The aim of this publication is to present a new design tool developed to calculate the lateral displacement of reinforced and/or prestressed concrete long columns.The algorithm is composed of two parts, first the generation of the characteristic curves on the cross-section: bending moment -curvature (M -ϕ) and axial force -bending moment (N -M); and second, the variation of bending stiffness in accordance with the magnitude of the bending moment (M -E•I), the variation of the lateral displacement along the column length (H -w), and generation of the force diagrams: axial loading, shear force, and bending moment.Computation on the cross-section is iterative, while the calculations on the entire element are done using the Finite Element Method (FEM).Investigations for determining the accuracy of the algorithm were done by comparing the computational results with the ones from experimental tests done by other researchers [4] [6], as well as by running in parallel numerical analysis in a FEM based commercial software, SAP2000, for the same input values.
The programming language used for writing the algorithm was MATLAB, which afterwards was transformed into an app using MATLAB appDesigner, resulting into a design tool with a friendly user interface.

Algorithm development
From a functional point of view, the algorithm can be divided in two main parts: the starting point preceding the determination of the curve and the computation of the actual curve.The input data for computation include: cross-section geometry, properties of materials, loss of prestress and material degradation parameters, followed by processing the input data: cross-section transformation and the determination of the actual prestressing force.Afterwards, the following diagrams are determined: plastic bending moment -axial force (Mpl -N), bending moment -curvature (M -ϕ), bending moment -flexural stiffness (M -E•I), height -displacement (H -w) and force diagrams.

Input
In this section the description of the considered column is made.As a departure point, the column crosssection definition is required.The cross-section is given in two steps, first step: the height (h) and the width (b) of the column cross-section, and second step: position of the longitudinal reinforcement [2].The reinforcement is defined by coordinates and geometrical characteristics (diameter, area).A visual representation of the input data is shown in figure 1 and in figure 2.  For the rest of column`s geometry, the position of the nodes is given [7].In this purpose a three by m matrix is used, where: m is the number of nodes defined, the first column consists of the node index, the second and the third column comprise the node coordinates.A second matrix is needed to specify which degree of freedom should be blocked for each node, therefor a four by k matrix is used where k represents the number of nodes, the first column is the index of the node which is affected by restraints and the last three columns represent the following: translation on the horizontal direction, translation on the vertical direction and rotation by the z axis.Is assigned 1 if the degree of freedom is blocked and 0 if it's free.Configuration of notations for commonly used supports are presented in table 1.The nodal loads are defined also by using matrixes.
The elements that connect the nodes are defined in a matrix, indexing the element in the first column, mentioning the starting (i) and the end (j) nodes in the second and third column, and the fixation type in the last column.The fixation type is needed to correctly determine the stiffness matrix.
The definition of the used materials and the parameters from losses of prestress are presented in [2].The constitutive curves are determined according to [8] and [9].
Lastly, the parameters of material degradation are needed.By determining the effect of the shrinkage and creeping of the concrete, the material degradation can be numerically quantified [10].Each effect results in a value with which the ultimate compressive strain of the concrete is reduced, equation (1).The data required is the following: the age of concrete when loaded, t0, coefficient dependent of used cement, s, the coefficients, kh, , relative humidity, RH, the age of the concrete at the beginning of drying shrinkage, ts.
where εt is the value with which εcu is reduced due to concrete shrinkage and εcc_inf_t0 is the reduction value from the concrete creeping.

Determining the transformed (idealized) cross-section
After the input section is filled, the transformation of the cross-section is made [11].The idea behind the transformation (idealization) of the cross-section is that the axial stiffness of passive and active reinforcement is taken into consideration.The reinforcement is transformed into rectangular concrete evasions which are attached to the cross-section at the corresponding coordinates in the y direction, as shown in figure 3. The process in which this transformation is reached is made in the following order: the calculation of the transformation coefficient presented in equations (2), the computation of the idealized area with equation (3), the determination of the geometrical center of the idealized cross-section with equation (4) and the calculation of the moment of inertia for the transformed cross-section, using the parallel axis theorem, presented in equation (5).
where ΣAi is the sum of the equivalent concrete areas, yi is the distance from a chosen end of the cross-section to each equivalent concrete area.
where Ai are the equivalent concrete areas, Ii are the moment of inertia of each area.

Plotting the axial force -plastic bending moment interaction curve (N -Mpl)
Succeeding the preliminary calculations, the first interaction curve determined is the axial forcebending moment (N -Mpl).To determine the curve, the behavior of the column is divided in different stages as presented in [2].The theoretical basis for this section is presented in [12].The step-by-step calculation is presented in [10].An iterative calculation is realized for three out of the six stages: 2a, 2b, 2c.The reason for is that usually the behavior of the concrete columns is in these areas.The extreme points: Ia and IIIb are determined with equations ( 6) and (7).The variation between these end points and the iterated segment are considered linear.
At pure compression (point IIIb), only concrete is considered to carry load.At this stage the steel reinforcement is prone to buckling and its contribution is neglected in the calculation.At pure tension (point Ia), only the reinforcement is considered to carry tensile stress, consequently, the concrete tensile strength is neglected throughout the computation.For prestressed concrete, the compressed active reinforcement is disregarded in calculations [13].
The three iterated segments are determined using a linear strain diagram.The strain diagram is determined as follows: for each segment the strain value in the extreme compressed concrete fiber is considered equal to εcu, the strain value of the tensioned reinforcement is chosen with regards to each segment's interspace, a line that connects the two values is drawn and the other fibers strain values are determined geometrically, with equation (8).
where ϕ is the curvature of the section, x is the compressed concrete height and yfib is the height of each fiber.
If the cross-section includes active reinforcement, then the iterated strain value is that of the tendon, otherwise the strain of the passive reinforcement is considered.The difference between the three iterated segments is the interspace in which the value of reinforcement strain is chosen for each calculation cycle.These interspaces are presented in figure 4. The next step is the determination of the stress diagram based on the strain diagram, using Hook`s law.For each fiber the stress value is determined by equation (9).After obtaining the stress diagram, the force value is determined for each fiber with equation (10).
F fib =σ fib •A fib (10) where Efib is the modulus of elasticity of each analyzed fiber, Afib is the area of each analyzed fiber.
Last, the axial force and the plastic bending moment are determined.The axial force is calculated by summing up the forces that are present in the fibers, while the plastic bending moment is determined by summing the moments that each force generates around the cross-section`s geometrical center.The analytical formulation is presented in equation (11).

Plotting the bending moment -curvature interaction curve (M -Φ)
The theoretical background for the determination of the M -ϕ interaction curve is like that of the N -Mpl curve.The strain-stress-force diagram determination method is identical to the one presented in the previous subsection.The determination of the M -ϕ curve is done in three different stages as presented in [2].The difference between each stage is the chosen strain that's going to be iterated.These stages are presented in figure 5.The theoretical basis is presented in [12] and [14].Figures 6 to 9 present an example of graphical outputs obtained with the algorithm by before described computational approach.The three stages are color coded: blue represents stage one, green stage two, and red stage three.

Plotting the bending stiffness -bending moment interaction curve (M -E•I)
The bending stiffness -bending moment curve is determined from the previously obtained M -ϕ curve.
The theoretical basis of the determination of this curve is presented in figure 10.For each point of the M -ϕ curve a line is drawn from the origin through that point.The angle between the horizontal axis and this theoretical line is the value of the stiffness in that point.For each determined plotted point of the moment -curvature diagram, the algorithm determines the stiffness value by using equation (12).Then the value obtained is capped by values obtained in equations 13 and 14.If the cross-section is uncracked then equation ( 13) is used as the capping value, otherwise equation ( 14) takes its place.The way the algorithm determines if the section is cracked or uncracked is by comparing the stress value at the extreme tensioned concrete fiber with fctk.If the obtained value exceeds fctk, then the section is considered cracked [15].
where (E•I)id is the stiffness value determined for the transformed cross-section, Atot is the summed area of the reinforcements present in the cross-section, dtot represents the effective height (the distance between the center of mass of the tensioned reinforcement and the top of the section), β is a reduction coefficient.

Plotting the lateral displacement -height interaction curve (w -H)
To determine the height -displacement curve the finite element method (FEM) is used, for linear elements only.The benefit of such method is that the geometry is easily customizable.The drawback of the method is that the accuracy of the results is highly dependent on the discretization used.The theoretical basis is presented in [14] and [7], while the implementation method is detailed in [10].
The implementation of the method consists in matrix calculation, having the first step the column input.This aspect was discussed in subsection 2.1.The data needed is the description of the nodes and the elements that connect the nodes.Second, the degrees of freedom for each inputted element are indexed and gathered in a matrix called 'The matrix of degrees of freedom, M_GL'.The content of the M_GL matrix is shown in expression (15).

M_GL= �
FE index u x i u y i w z i u x j u y j w z j ⋮ � where FE index is the index of the current element and ux, uy and wz are the degrees of freedom for each node of the element.Afterwards the 'Localization' matrix, L, is constructed, which helps collect each specific matrix element into a single matrix.The localization matrix is presented in expression (16).
Next a rigidity matrix is selected for each element.These matrixes are attributed to each element using the fixation type defined in the input.In addition, matrixes are defined in local coordinates for each element so a transformation matrix, R, is needed for the conversion to global coordinates, expression (17).After transformation, each rigidity matrix is gathered in a matrix using L, as in equation K=L T •k g •L (18) The next step consists in the generation of the nodal load vector, P. The vector cumulates the values of the nodal loads for each node.Afterwards, the rigidity matrix and the nodal load vector are reduced with regards to the nodes that suffered losses of degrees of freedom.Then the displacement, U, and the force, S, matrix is determined.The expression that helps determine these vectors are presented in equations (19) and (20).After obtaining the displacement values in the direction of each degree of freedom, a displacement matrix is filled with the blocked displacements that were previously reduced.
U free =K reduced -1 •P reduced (19) where Kreduced is the reduced rigidity matrix and Preduced is the reduced load vector.S=K•U (20) After calculating the value for the moment for each finite element, the stiffness is altered in the inputted data and the calculation is repeated.The stiffness value is picked using the M -E•I interaction curve.In the input, the stiffness is defined using the value for moment of inertia, I, so the expression that considers the rigidity reduction is shown in expression (21).

App development
To make the algorithm more intuitive to use, an application was developed in MATLAB appDesigner [10].The application interface can be seen in figure 11.The interface is divided in two main parts: input on the left-side and output on the right-side, both are further divided in subsections.For input these subsections are named as: Material, Section, Element, Prestress and time effect data, and Results, while for output three subsections are created: Geometry, Section curves, Element curves.The data required for the input section is identical to that listed in subsection 2.1 `Input`.The 'Material' subsection has three main parts: the definition of concrete, reinforcement and tendon.The layout of each part is alike.They include the following elements: data input area (top left), informative graphic (top right), input data description (bottom left) and graph area for the added data (bottom right) activated by the `Plot` button.The layout is presented in figure 11.
The subsection 'Section' similarly to the previous one has three main parts: the definition of the concrete cross-section and the positioning of the passive and active reinforcement.The layout of the concrete cross-section definition is different compared to the reinforcement positioning.The crosssection is defined by its height (h) and with (b), values which can be added in the two input boxes (top left).By pressing the 'Save' button the application will calculate the area and the moment of inertia of the given cross-section.The layout is presented in figure 12.The reinforcement is added by mentioning the position in the cross-section and the diameter.For the active reinforcement, additionally the area of the tendon is also required.The reinforcement is added one by one, and it can be verified in the table shown at the bottom side.
The 'Element' subsection is also divided in three parts: the definition of nodes and elements, and the attribution of nodal restraints and loads.Nodes can be added by giving them an index and coordinates.The interface is presented in figure 13.Restraints and loads can be attributed to the nodes by mentioning the index of the node and the desired values, in regards of the prescriptions from the subsection 2.1 'Input'.Elements are added by defining an index, the beginning and end nodes ('i' respectively 'j') and the type of fixation.The data will be inserted to a table at the bottom of the subsection, along with the transformed area and moment of inertia, the length of the element and the angle between the horizontal axis and longitudinal axis of the element.The 'Prestressing and time effect data' subsection is divided in two parts.The loss of prestress calculation is made on the left-side of the subsection, and the material degradation on the right-side.In addition, each coefficient is explained in the text boxes, as in figure 14.
Through the last subsection, 'Results', the user can control de increments in which the diagrams are determined and can also obtain some numerical results as well, figure 15.This subsection is divided in four parts.The first part transforms the cross-section in an idealized one, utilizing the above-mentioned steps.If elements where defined, then this step is automatically made when the first element is added.The second part determines the number of iterations in which the N -Mpl interaction curve is split.In this section, the coordinates of the balance point and the corresponding plastic moment to an arbitrarily chosen axial force is given.The third part determines for different values of axial force and draws the curves M -ϕ and the M -E•I, through an iterative process.The last section determines the number of finite elements in which the elements are discretized to determine the H -w interaction curve.

Results
To test the accuracy of the results provided by the algorithm, eight types of cross-sections were modelled.Three cross-sections were taken from [4], figure 19, and the remaining five from [6], figure 20.The first three cross-sections were tested experimentally and modelled in a commercial software for structural analysis, SAP2000.The generated N -Mpl and M -ϕ curves were compared with the ones obtained from the algorithm.The detailed results are presented in [2].The maximum displacement was determined for the maximum lateral load and compared with the values presented in [4], table 2. The M -E•I curves are compared in figure 21.
Table 2: Displacement value comparison between experiment (from [4]) and the computed ones.

Conclusions
From the bending moment -bending stiffness (M -E•I) diagrams we can deduce the following conclusions: the axial loading has a positive effect on the stiffness of the element.For higher axial loading (for example columns S02 as well as for S03) due to prestressing the bending stiffness remains high on a larger interval of bending moments.The prestressed columns started losing stiffness after the bending moment exceeded 25kNm, while the columns that were not prestressed had a lower threshold of 10kNm.The results of the analysis when compared to the ones from experiments [4] and other published scientific papers [6], show that the approach developed initially by Faur [12] [16] [17] for calculating reinforced concrete members in normal cross-sections can be successfully further developed and applied to prestressed concrete members as well and through programming transformed into an app.The comparison between the M -E•I diagram generated by the algorithm and the ones taken from [4] show that the value of the initial stiffness it has a lower calculated value that the one taken from references.This difference in initial stiffness can be observed in the comparison made between the algorithm and the cross-sections taken from [6].
In the case of cross-sections SJ2, SJ3, SJ4, and SJ5, we can see that the allure of the curves is similar, but the maximal bending stiffness value differs.The algorithm generates conservative bending stiffness values (between 60 and 70 MNm 2 ) when compared to the ones from [6] (between 75 and 85 MNm 2 ).

Figure 1 :
Figure 1: Size and position of passive reinforcement.

Figure 2 :
Figure 2: Size and position of pretensioned reinforcement.

Figure 3 :
Figure 3: Graphical representation of the transformed cross-section.

Figure 4 :
Figure 4: Stress state and force diagram for each segment.The reinforcement strain value interspace for each segment is the following: [εyp, εup] for 2a, [0, εyp] for 2b, and for 2c the neutral axis positioning interval is [h, dp] from [10].

Figure 9 :
Figure 9: σ -ε curve for reinforcement in tension generated by the algorithm.

Figure 10 :
Figure 10: Bending stiffness value for the bending moment value at iteration "i".

Figure 14 :
Figure 14: Loss of prestress and material degradation data input.

Figure 18 :
Figure 18: Output of displacement and force diagrams on the element.

Figure 21 :
Figure 21: Comparison between algorithm generated M -E•I curve (red) and the one obtained by [4] (blue).For the last five cross-sections the M -E•I curve were compared.The results are presented in figure 22.

Figure 22 :
Figure 22: Comparison between algorithm generated M -E•I curve (red) and the one obtained by [6] (blue).

Table 1 :
Notations for common supports.