A simple numerical approach for the pushover analysis of slender cantilever bridge piers taking into account geometric nonlinearity

The response of slender bridge piers to horizontal actions may be significantly influenced by geometric nonlinearities. In such conditions, the use of sophisticated models implemented in complex structural analysis software can be economically disadvantageous, especially in the preliminary design phases. This paper proposes a simple numerical procedure to compute the nonlinear pushover response of cantilever bridge piers subject to horizontal loads. The procedure is based on an iterative approach to enforce the element equilibrium under large displacements, efficiently accounting for P-Delta effects induced by vertical loads. Evaluation of the bending moment–curvature response of the pier base cross section is required and used as basic input data. For fast preliminary analyses, sectional response can be manually computed in simplified linearized form, thus completely eliminating the need to use structural analysis software. Indeed, the entire procedure can be implemented in standard programming codes, such as PythonTM or Matlab®, and used to evaluate the pushover response of piers with arbitrary cross section. Comparison with experimental test results and solutions based on Finite Element models shows that proposed procedure can be used to get a fast, yet accurate, estimate of the entire force–displacement curve and, in particular, of the pier ultimate displacement.


Introduction
Cantilever Reinforced Concrete (RC) piers are a common solution to realize substructures for simply supported as well as continuous bridge decks. Thanks to advanced casting techniques such as slipforms, this type of piers may reach considerable heights and aspect ratios such that geometric nonlinearities become relevant.
Nonlinear structural modeling of bridges under earthquake actions is often carried out by modeling this type of piers as link elements. This approach simplifies the study of bridge behavior and reduces the computational efforts. The constitutive relations for such zero-dimensional elements are expressed in terms of base shear-top displacement curves that can be computed either by simplified methods based on sectional response and plastic hinge length assumptions or by explicit nonlinear pushover analysis of isolated piers (Priestley et al., 1996).
For slender geometries, nonlinear geometric effects must be considered, as horizontal loads usually produce large transverse deformations and induce significant P-Delta Paolo Di Re, Davide Bernardini, Daniela Ruta and Achille Paolone contributed equally to this work. effects that can significantly change the shape of the capacity curves with respect to the case of linear geometry (Priestley et al., 1996;Fenwick et al., 1992;Dhakal & Maekawa, 2000).
Several methods have been proposed to determine the response of cantilever beams subject to concentrated and distributed loads under nonlinear geometric effects. However, most of them assumes linear elastic material response, neglecting the effects of constitutive nonlinearities. A comprehensive review of the relevant scientific literature was discussed by Al-Saddler and co-authors (Dado & Al-Sadder, 2005;AL-Sadder & Al-Rawi, 2006), including both analytical and numerical formulations. Earlier proposals adopted analytical methods based on the solution of elliptic integrals to analyze linear elastic structures (Bisshopp & Drucker, 1945;Timoshenko & Gere, 1963), usually performed by imposing some limitations on the element geometry, loading conditions and maximum displacements. Only in few cases, e.g., Ohtsuki & Ellyin (2001), the approach was extended to the analysis of more complex structures, such as square frames. Recently, Zhang & Chen (2013), Tari et al. (2015) and Cammarata et al. (2019) extended the method to a wider range of geometric and loading configurations. Finally, Wang et al., (2018) introduced a new procedure based on the analysis of rocking effects to reproduce the behavior of post-tensioned precast piers.
Although analytical approaches are attractive, they are limited to simple material constitutive behavior. Hence, numerical methods have been preferred and widely developed in the last decades, due to their higher flexibility and applicability to complex geometries and loading conditions. In linear elastic material range, many authors, e.g., Wang & Kitipornchai (1992), Lee (2001) and Banerjee et al. (2008), proposed numerical procedures that iteratively solve the governing equations by assuming starting values of specific parameters. These are known as numerical integration methods based on iterative shooting techniques. Similar approach was also used by Dado & Al-Sadder (2005). Other numerical approaches were based on Finite Difference (Wang et al., 1961;AL-Sadder & Al-Rawi, 2006) and Boundary Element Method (Katsikadelis & Tsiatas, 2003) and applied to the analysis of curved and non-prismatic beams.
For slender RC piers, accounting for both geometric and material nonlinearities, is essential to perform correct simulation of structural response. Babazadeh et al. (2017) focused on nonlinear response of RC bridge piers and developed a new analytical model based on the evaluation of elastic and inelastic cross-sectional rotations along the pier height. The model provides closed form expressions for transverse displacement, bending moment and curvature distributions, as well as length of the plastic region, for given values of the base bending moment induced by top concentrated loads. However, this is limited to bilinear monotonically increasing moment-curvature diagram of the element cross section and neglects pier self-weight.
For general cross-sectional response due to more realistic, but complex, nonlinear material behavior, Finite Element (FE) method is usually employed, largely adopted in both professional practice and research works. In this case, piers are modeled with nonlinear beam-column FEs, using either lumped (Scott & Fenves, 2006) or distributed (Lee and Filippou, 2009; plasticity approaches. Within this framework, common approach for inclusion of material nonlinearities is the adoption of cross-sectional fiber discretization (Spacone et al., 1996;Kostic & Filippou, 2012;Poston, 1986;Sessa et al., 2019), where beam cross section is divided into multiple sub-areas (fibers) characterized by nonlinear constitutive laws. When Euler-Bernoulli beams are considered, onedimensional (1D) material constitutive relationships suffice to describe the fiber response and can be also exploited to account for corrosion of the material in existing structures .
Among beam-column FE models, force-based and mixed formulations (Ciampi & Carlesimo, 1986;Taylor et al., 2003) have proven superior performances because of the advantages given by the exact interpolation of the internal forces (Spacone et al., 1996;Alemdar & White, 2005). Hence, modeling of slender bridge piers with forcebased beam-column FEs is nowadays a consolidate procedure (Mackie & Scott, 2019) that has been widely adopted in many recent studies, e.g., Su et al., (2020), Yang et al., (2021), Pozo et al., (2022) and Bernardini et al., (2021). However, in addition to pier strength assessment, bridge seismic design and safety check usually request the evaluation of pier ductility capacity, which necessarily entails the correct evaluation of the element ultimate displacement. As force-based FEs rely on the exact interpolation of beam internal forces, they enforce the kinematic compatibility conditions in weak form and, therefore, may yield to approximate results on strain and displacement values. Hence, accuracy on ultimate drift displacement can be obtained at the computational price of increasing the number of integration points or that of the FEs composing the mesh.
Moreover, nonlinear geometric effects are usually included in beam-column FEs through corotational approaches (Crisfield & Moita, 1996;Felippa & Haugen, 2005; assuming the element able to undergo large nodal displacements, yet considering small strains at the sectional level. Similar simplified computational schemes are also adopted when the analysis only requires to include the main nonlinear contribution associated to P-Delta effects (Gaiotti & Smith, 1989;Barros et al., 2010;Burgueño et al., 2016;. In these cases, mesh refinement is required, i.e., single structural member is discretized into multiple beam FEs to approximately account for secondorder effects due to large element deflection (P-effects), in addition to those caused by large global displacements of the element ends (P-Δ effects) . This paper proposes a simple numerical procedure to evaluate nonlinear pushover behavior of slender cantilevered RC bridges piers under horizontal top loads. The procedure is based on the assumption that moment-curvature response ( M-K ) of the pier base cross section is known for a constant value of the vertical compressive stress and can be used to uniquely determine the force-displacement curve response of the pier. Indeed, typical bridge pier M-K responses are characterized by a first monotonically increasing branch, up to the bending peak strength, M p , followed by rapidly decreasing branch (Fig. 1), due to the element collapse (Babazadeh et al., , 2017Mackie & Scott, 2019).
Focusing on the first monotonically increasing part of the M-K curve, the procedure assigns bending curvature value at the pier base and iteratively determines the corresponding value of the top horizontal force, together with transverse cross-sectional displacements along the entire element height.
The proposed procedure is also applicable by assuming linearized M-K curves ( Fig. 1) as starting data for the base cross-sectional response. This can be convenient for preliminary fast assessment of the pier seismic capacity and is one of the main advantages of the approach.
After describing the main assumptions and fundamental equations of the model (Sect. 2), the proposed procedure is described in Sect. 3. Finally, Sect. 4 shows three numerical applications. In the first one, the behavior of a circular pier is investigated to study the ability of the proposed model in capturing the nonlinear geometric effects and compare them with those of a fiber beam model. In the second one, the effects of variable axial forces are discussed. Finally, third test is used to validate the proposed procedure against experimental test results.

Basic equations for a geometrically nonlinear cantilever
A straight RC cantilever pier with height H is considered (Fig. 2). The pier is subject to its self-weight and, at the top, to a constant vertical concentrated force, P, and a variable horizontal concentrated force, F. Pier cross section is uniform along H and has slenderness ratio sufficiently high to neglect the effects associated to transverse shear deformations. Indeed, Euler-Bernoulli beam theory is considered to model the pier in the reference system (x, y, z), being x the vertical axis and y and z the axes lying on the cross-sectional plane, with y parallel to the direction of force F. Self-weight is defined as q = A ,  where A is the cross-sectional area and the weight per unit volume.
Equilibrium is defined in the element deformed configuration, described by the transverse displacement profile, v(x), of the general cross section at x. As a result, full P-Delta analysis is performed, i.e., both P-Δ and P-effects  are considered. Pier is assumed as infinitely rigid in the axial direction, i.e., vertical displacements and deformations are neglected in the equilibrium definition.
Proposed procedure is based on the hypothesis that moment-curvature response ( M-K ) of the cross section does not vary along H and is always equal to that of the pier base cross section. In other words, the procedure supposes that linear variation of the axial force due to the element self-weight produces negligible variation of the bending behavior of the cross sections along the height. Moreover, as usually occur in typical RC elements, M-K curve ( Fig. 1) is assumed as characterized by a monotonically increasing branch, going from zero to peak point (K p , M p ) . The descending collapsing branch is neglected.
The equations governing the proposed numerical procedure are derived referring to the general deformed configuration of the pier, Fig. 2b. The equilibrium conditions of the pier portion above the general cross section at x permits to compute the variation of bending moment, M(x), along the height of the pier. This results in the following local equilibrium conditions where is used as substitute variable of x to integrate the contribution of q along the portion of the pier above x. For the cross section at the base of the pier, bending moment results as which represents the global element equilibrium equation.
Constitutive condition relates M(x) to the corresponding curvature K(x). This is preliminary computed for the base cross section of the pier and is equivalently used for all cross section along H. Thus, the biunivocal expression giving M(x) as function of K(x), and vice versa, is written in closed form, that is Finally, by imposing the element compatibility conditions, bending rotation, (x) , and transverse displacement, v(x), of (1) and the general cross section at x are derived through integration of curvatures K(x), that is and Numerical procedure for the computation of force-displacement curve For cantilever piers under linear geometry, numerical determination of F-v curve is straightforward when M-K response curve is monotonic. In this case, a step-by-step procedure can be followed, assuming force F as independent controlling variable. Increasing values of transverse force, F, are assigned and moment distribution along the pier, and, thus, by solving the integrals in Eqs. 4 and 5, transverse displacements, v(x), and, thus, top displacement, v , are easily computed. Under these assumptions, peak point of the F-v response is reached when base cross section reaches its peak bending moment M p , (Fig. 3) which, in this case, also coincides with the maximum pier capacity. However, when P-Delta effects are included, same approach cannot be followed. In fact, for a given deformation state of the pier, nonlinear geometric effects produce an increase of bending moments along the pier. As a result, the peak of the global force, F, is reached when bending moment at the base cross section is lower than M p (Fig. 3). In the post-peak part of the global response, force F is expected to decrease, while bending moments along the pier keep increasing. Hence, if force F is used as independent controlling variable of the computational algorithm, objectivity of the solution is lost.
To avoid this issue, a new approach is proposed in the following, where M-K response curve of the pier base cross section is followed and used as starting point of the analysis. In other words, the proposed numerical procedure permits the determination of pier behavior in terms of pushover force-displacement response curve, F-v , based on the curve parametrization with respect to the curvature, K B , of the base cross section, which is assumed as given independent controlling variable of the analysis. This approach gives a significant advantage for the computational algorithm, when P-Delta effects are included.
In the proposed numerical procedure, bending curvature, K B , at the base of the pier is used as independent controlling quantity. This is increased in a step-by-step approach, up to the peak value, K p . At each step, corresponding values of force, F, and displacement, v , are determined. Hence, the following algorithm is applied, where superscripts n − 1 and n denote the previous and current step, respectively. However, when P-Delta effects are included in the analysis, moment distribution along the pier is influenced by the pier deformation, that is by the transverse displacement, v(x), to be determined at each step (Eq. 1). Hence, an iterative approach must be followed to find, at step n , the values of v(x) that satisfy the equilibrium conditions in the deformed configuration. In the following, subscript i − 1 and i denote the previous and current iteration, respectively, performed at current step n.
The algorithm works as follows and is summarized in 3. Iteration counter is set to one, i = 1 , and the initial values of the transverse displacements are recovered from From the element equilibrium condition, horizontal force, F n i , is obtained by inverting Eq. 2, and assuming the deformation state of the pier as frozen at the previous iteration, that is

11)
If r > tolerance → i = i + 1 and GO TO point 4) else → State at n set equal to state of iteration i and CONTINUE with step n + 1 5. Bending moment distribution along the pier is determined from Eq. 1): 6. Curvature distribution along the pier is determined for current iteration from the constitutive relationship, Eq. 3, as: 7. Bending rotation variation along the pier is determined for current iteration from Eqs. 4, that is 8. Updated transverse displacement variation along the pier is determined for current iteration from Eqs. 5, that is 9. Top displacement is computed for current iteration as: 10. The following residual quantity is computed: represents the quantity increment at iteration i with respect to iteration i − 1; 11. If r is greater than a tolerance value, iteration counter increases, i = i + 1 , and procedure returns to point 4); otherwise, response quantities of current step are stored and the procedure continues with next step, n + 1 , of the analysis.
To be noted is that, when nonlinear geometric effects are neglected, displacement-dependent terms in the Eqs. 1 and 2, used at point 4) and 5), vanish, that means transverse displacement variation along H does not affect the equilibrium conditions. As a result, one iteration is required in the proposed procedure at each step of the analysis.

Applications of the proposed procedure
To validate the proposed procedure and investigate its performance in determining the response of RC bridge piers, three numerical tests are presented in the following. First test is a benchmark conducted on the cylindrical pier numerically studied by Mackie & Scott (2019) under linear geometry assumption and neglecting the effects due to pier self-weight, i.e., under uniform axial force. Here, contribution of P-Delta effects is included and the performance of the proposed model in capturing the actual pier response is investigated, comparing the results with those obtained with a beam FE model with fiber cross-sectional discretization. The second test considers a rectangular hollow pier.
As proposed procedure assumes that M-K response does not vary along H, regardless of the axial force variation due to pier self-weight, this test discusses the validity of this assumption by comparing the responses obtained by including and neglecting the pier self-weight. Finally, third test reproduces the behavior of a circular pier experimentally tested during the campaign described in , Burgueño et al. (2016) and  and aims at validating the proposed procedure against real specimen results. Analyses are performed using a home-made Python TM -based implementation of the procedure, where Simpson quadrature rule is used to solve all integrals defined along the pier axis.

Response under uniform axial force and comparison with fiber beam model
The RC circular pier depicted in Fig. 4 and originally studied in Mackie & Scott (2019) is considered first. This application Fig. 4 Circular pier test: pier geometry, reference system and loading scheme is used to compare the results of the proposed procedure with those obtained with force-based FE fiber models, aiming at reaching the accurate prediction of pier ultimate displacement capacity. Pier cross section has diameter and clear cover equal to D = 50.8 cm and c = 2.54 cm , respectively, and is reinforced with eight equally spaced #8 bars (Ø = 25.4 mm ). Pier height is not given in Mackie & Scott (2019), but H = 6.1 m is assumed, in accordance with the other reported data, corresponding to a aspect ratio H/D equal to 12.
A constant vertical load, P = 600.51 kN , acts at the top of the pier, while variable transverse force, F, is parallel to reference axis y, defining with z the cross-sectional plane. Pier self-weight, q, is neglected.
Same constitutive assumptions as in Mackie & Scott (2019) are adopted, i.e., concrete is modeled as an elastic-no-tension material with compressive stiffness equal to E c = 24855.61 MPa and reinforcing steel is modeled as an elastic-perfectly-plastic material with elastic stiffness and yielding stress equal to E s = 199948.04 MPa and f y = 344.74 MPa , respectively.
Analyses are conducted by considering both actual and linearized M-K cross-sectional responses shown in Fig. 5. Actual response (solid blue curve) is obtained through sectional analysis performed in OpenSees software (2021), for a value of the axial compressive stress equal to P. The response coincides with that reported in Mackie & Scott (2019). In particular, a midpoint fiber cross-sectional discretization is used with a polar grid of 50 and 30 fibers placed along the circumferential and radial directions, respectively, for the core and 50 and 5 fibers placed along the circumferential and radial directions, respectively, for the cover. Reinforcing bars are modeled as additional fibers.
In Fig. 5, points ( K y , M y ) and ( K p , M p ) indicate the curvature-moment pairs at first steel yielding and at peak strength, respectively. These are used to construct the linearized M-K response (dashed black curve), by assuming equal area under the post-peak part of the curve. Peak strength curvature is assumed equal to K p = 0.0880 m −1 .
Following the study in Mackie & Scott (2019), F-v response curve is first computed for linear geometry hypothesis and, then, nonlinear geometric effects are included.
Linear geometry analysis Figure 6 shows the solution obtained with the proposed procedure (dotted black curve) by considering the actual M -K cross-sectional behavior. The results are compared with those obtained by adopting the force-based Euler-Bernoulli beam-column FE with distributed plasticity available in OpenSees and same fiber cross-sectional discretization used for sectional analysis. FE model details coincide with those used in Mackie & Scott (2019), except for the number of Gauss-Lobatto quadrature points placed along the axis of the FE modeling the pier. A 6-points rule is considered in Mackie & Scott (2019), giving the results represented by the dashed red curve, while Fig. 6 also shows the those obtained with 10 (solid blue curve) and 7 (solid orange curve) points. For OpenSees FE model, analyses are carried out by monitoring the bending curvature attained at the pier base cross section and execution is stopped when this reaches the value K B = K p = 0.0880 m −1 . In this way, ultimate point of F-v curve directly corresponds to ( K p , M p ) at the pier base, as happens in the proposed procedure. Although all model solutions agree very well with each other, OpenSees model underestimates pier ultimate deformation when only 6 quadrature points are used and provides a lower ultimate displacement value ( 6% error) with respect to the solution obtained with the proposed procedure, which gives a more correct result. A higher number of quadrature points must be used in the beam FE model to capture the correct displacement capacity of the pier.

Nonlinear geometry analysis
When nonlinear geometric effects are included, force-displacement curves are significantly different. Figure 7 compares the F-v response curve obtained with the proposed procedure (dotted black curve) with those evaluated through similar OpenSees beam FE model used for linear geometry analysis, yet accounting for P-Delta effects. Given the results in Fig. 6, this assumes a 7-points Gauss-Lobatto quadrature rule. However, mesh discretization considers 20 (solid blue curve), 10 (solid orange curve) or 5 (dashed red curve) FEs, to better capture the nonlinear phenomena due to P-effects.
For this level of slenderness ( H∕D = 12.0 ), it turns out that geometric nonlinearities have strong influence on the force-displacement curve with respect to the case of linear geometry.
All solutions agree with each other also in this case, but FE models with coarser mesh discretization give stiffer nonlinear behavior with respect to the correct result obtained with the proposed procedure, up to 9% error for the ultimate displacement. However, all models provide identical peak strength. This is indicated in the figure by the green point ( v * , F * ).
To be noted is that, under P-Delta effects, peak force, F * , is reached when bending moment at the base cross section is lower than its maximum value, M p . This is confirmed by the plots in Fig. 8 Finally, Fig. 9 shows the F-v response curves obtained under both linear and nonlinear geometry, by assuming the linearized M-K cross-sectional behavior (dashed red curves). These are compared with previous results evaluated by As expected, the model underestimates the initial stiffness of the pier and provides a sharper transition from the elastic to the inelastic material behavior, in accordance with the piecewise linear variation of the linearized M-K response. However, global behavior of the pier is well captured, providing similar peak strength, yet higher ultimate displacement values.

Response under variable axial force: influence of pier self-weight
The second application is the analysis of a rectangular hollow RC bridge pier. To investigate the influence of the simplifying assumption of uniform M-K response made in the proposed procedure, the analysis is carried out by including and neglecting the pier self-weight. Specimen geometry is depicted in Fig. 10, where cross-sectional reinforcement arrangement and element loading scheme are also indicated. Cross-sectional dimensions are equal to L y = 450 cm and L z = 250 cm , along the strong, y, and weak, z, direction, respectively. Thicknesses are equal to s y = 50 cm and s z = 40 cm . Concrete clear cover is uniform for both external and internal edges and equal to c = 3 cm , measured from the external edge of the stirrups. Longitudinal and transverse reinforcements are symmetrically arranged. Four rows of eighteen bars with diameter equal to 14 mm are located along the small dimension, L z , and four rows of 24 bars with diameter equal to 28 mm are located along the large dimension, L y . Stirrups have diameter and spacing equal to 14 mm and 20 mm , respectively. Pier height is equal to H = 30 m , corresponding to a aspect ratio of 6.67 and 12.0 for the strong, y, and weak, z, bending direction, respectively.
A constant vertical load, P = 7703.1 kN , acts at the top of the pier, while variable transverse force is alternatively applied in the strong, F y , and weak, F z , direction, respectively, to study the response of the pier for two different slenderness ratios. Pier self-weight is assumed as equal to q = 132.5 kN/m. M-K cross-sectional response is obtained through sectional analysis performed in OpenSees software by adopting a detailed fiber discretization model. Cross section is divided into quadrilateral patches, as depicted in Fig. 11 and each patch is discretized in a regular grid of fibers, according to midpoint integration rule. Yellow patches indicate the unconfined concrete regions, having thickness equal to 5 cm , i.e., approximately the distance between the crosssectional edge to the center of the longitudinal bars, while blue patches indicate confined concrete. For the latter, confinement criteria adopted to compute material mechanical properties are based on Mander's approach (Mander et al.   The number of fibers considered in the different directions is indicated by the numbers in Fig. 11. One fiber is always placed across the unconfined concrete thickness, while six and three fibers are used in the confined regions across the cross-sectional thicknesses in the directions parallel to s y and s z , respectively. Seventy-two and ten fibers are placed along the central parts of the cross section in the directions parallel to L y and L z , respectively. Reinforcing bars are modeled as additional fibers.
Kent-Scott-Park model ('Concrete02') is used to simulate the behavior of concrete fibers. Adopted mechanical properties are reported in Table 2 for unconfined and confined regions, being f pc and cp the compressive peak strength and strain, respectively, f pu and cu the compressive ultimate strength and strain, respectively, c the ratio between unloading slope at compression peak and initial slope, f ct the tensile strength and E cts the tension softening stiffness. Table 3 reports the mechanical parameters assumed to reproduce the behavior of longitudinal bars, for which 'Rein-forcingSteel' model is used, being f sy and f su the yield and ultimate stress, respectively, E s and E sh the initial elastic stiffness and the tangent stiffness at initial strain hardening, respectively, sh the strain at initial strain hardening and su the strain at peak stress. Figure 12 shows the M-K response curves corresponding to the axial forces acting on the top and base cross sections when self-weight is considered, i.e., N 0 = 7703.1 kN (black curve-top cross section) and N 30 = 11678.1 kN (blue curve-base cross section). If self-weight is neglected, axial force is uniform and equal to N 0 . Corresponding axial compression ratios are equal to N 0 ∕(f pc A) = 0.0354 and N 30 ∕(f pc A) = 0.0536 for the two cases, respectively. M-K response curves are plotted for both strong y and weak z bending directions. Figures 13 and 14 show the F-v response curves for the strong and weak bending behavior, respectively. These are computed both neglecting ( q = 0 ) and accounting for ( q ≠ 0 ) the pier self-weight, to highlight the effects due to axial force variation. Proposed model solution (dotted OpenSees -q = 0 OpenSees -q ≠ 0 Proposed -q = 0 Proposed -q ≠ 0 OpenSees -q = 0 Proposed -q = 0 OpenSees -q ≠ 0 Proposed -q ≠ 0 Fig. 14 Global response F-v curves obtained for the rectangular hollow pier for bending in the weak z direction curves) are compared with those obtained using the forcebased Euler-Bernoulli beam-column FE with distributed plasticity available in OpenSees and same fiber crosssectional discretization used for sectional analysis. For this model, a fine mesh is used considering 15 FEs with 7 Gauss-Lobatto quadrature cross sections, to obtain accurate comparison results also in terms of pier ultimate displacement. As done in the previous test of the circular pier, F -v curves are computed by monitoring the bending curvature attained at the pier base cross section and execution is stopped when this reaches the peak value, as plotted in Fig. 12, i.e., K p = 0.0267 m −1 for N 0 and K p = 0.0214 m −1 for N 30 . As observed for the circular pier, when self-weight is neglected ( q = 0 ), proposed model results perfectly match those obtained with the FE model. However, a small overestimate of the nonlinear pier stiffness is observed for q ≠ 0 , i.e., lower values are obtained for the ultimate displacement ( 5% error for the strong direction and 7% error for the weak direction). In fact, as axial force reduces moving from the pier base to the top, cross section above the base should follow a more flexible M-K curve than that resulting at x = 0 . However, in the proposed procedure, same M-K response is used for all cross section along H. To better explain this issue, Fig. 15 plots the variation along pier height of bending curvatures, K(x), and moments, M(x), at ultimate point of the F-v response, for bending in the weak direction, where the error is larger. Solid blue and black dotted curves represent the solution obtained with the FE and proposed model, respectively.
As observed, although the two models consider exactly the same moment distribution along the pier, in the plastic zone above the pier base, curvatures obtained with the proposed procedure are slightly smaller than those resulting from the FE software.
To eliminate the error, different M-K responses can be used to determine the bending curvature values of the various cross sections along the pier, at point 6) of the proposed procedure (see Table 1). However, this is more costly, as M -K curve must be preliminary computed for many values of the axial force. Alternatively, an amplification factor could be used to correct curvature variation along the pier at each iteration of the procedure (point 6)). These solutions are not explored in this work, since the error resulting by the assumption of uniform response is in any case small.

Comparison with experimental results
As last validation analysis, one of the four circular RC piers tested during the experimental campaign conducted in , Burgueño et al. (2016) and  and labeled as pier M121505 is considered. The pier has height H = 7.315 m and circular cross section with diameter D = 0.61 m (Fig. 16). Steel reinforcements consist of sixteen equally spaced #6 longitudinal bars (Ø = 19 mm ) and a #4 spiral (Ø = 12.7 mm ) with spacing equal to 76 mm . Variable horizontal load is applied through a hydraulic actuator that acts at the center of the top loading block, while constant vertical load is applied using two external post-tensioning bars. Nominal Babazadeh's circular pier: geometry and reference system vertical load, i.e., axial force in the rods, is controlled by two hydraulic jacks located on the top of the pier and is equal to P * = 712 kN , corresponding to axial compression ratio of P * ∕(f pc A) = 0.05 . Top loading block has height equal to 0.457 m and pier self-weight is assumed as q = 7.09 kN/m. Experimental test was conducted by applying the horizontal force under quasi-static control of the pier top displacement and following a fully reversed cyclic loading pattern consisting in 21 cycles with increasing level of deformation. Monotonic pier response was preliminary computed through fiber cross-sectional beam modeling, as reported in Burgueño et al. (2016). Hence, analyses reported in the following focus on the numerical evaluation of the monotonic response of the pier obtained through the proposed procedure. This is compared with both monotonic and cyclic response behaviors reported in the reference works.
Due to the particular configuration of the test setup, loading scheme of the pier slightly differs from that considered in Fig. 2 to develop the proposed procedure. In fact, the direction of the axial load is dictated by that of the posttensioning rods, that changes when pier deforms according to angle , as depicted in Fig. 17. This result as : Hence, vertical force P varies during the loading process, resulting equal to being L b = 0.483 m the height of the footing. Moreover, the effective horizontal force F eff given by the actuator differs from the net horizontal force acting in the pier, i.e., the horizontal base reaction F. Indeed, it results: To account for the correct orientation of the loading forces, for this test, the proposed procedure scheme is slightly modified. Referring to Table 1, at each iteration, before evaluating force F n i at point 4), angle is computed according to Eq. 6 and using the last value, v n i , of the top displacement obtained at the previous iteration. Thus, vertical force P to be used at points 4) and 5) is computed according to Eq. 7. Finally, when convergence of the iterative process is reached for the current loading step (point 11)), effective force F eff is obtained from Eq. 8. Moment-curvature cross-sectional response is also reported in Babazadeh's work, evaluated by means of fiber cross-sectional analysis. This is plotted in Fig. 18 and used to apply the proposed procedure. Figure 19 shows the F-v response curve obtained with the proposed procedure, distinguishing between the evolution of the effective force, F eff , applied at the top of the pier (dotted gray curve) and the base resisting force, F, (black dotte curve). These are compared with the corresponding monotonic responses (orange and blue solid curve, respectively) reported in Burgueño et al. (2016). Moreover, the experimental cyclic behavior of the pier is depicted as red solid curve.
Response curves obtained with the proposed procedure are in very good agreement with reference results, confirming that monotonic evolution of the horizontal resisting force well approximates the envelope of the cyclic behavior, although a very slight increase in pier strength can be observed, if compared with the results by Babazadeh.
Good match with the experimental results is also observed by comparing the variation along pier height of bending curvatures, K(x), proving the accuracy of the proposed procedure in reproducing the correct behavior of real specimen. Variation of K(x) is plotted in Fig. 20 for the ultimate loading step, corresponding to a value of v = 0.48 m.

Conclusions
A simple numerical procedure has been proposed to determine the structural response of RC bridge piers subject to constant vertical loads and variable transverse forces applied at the top, taking into account geometric nonlinearities arising under large displacements. The procedure considers the moment-curvature response curve of the cross section at pier base as input data and provides the force-displacement response curve of the pier as primary output. This is parametrized with respect to the base bending curvature, K B . At each step of the analysis, a simple and fast iterative approach is used to enforce the element equilibrium in the deformed configuration, so that transverse displacement variation along the pier height is univocally determined.
If axial force is uniform over the height, the procedure converges to the exact solution. This allows the user to get a quick and reliable computation of both pier strength and ultimate displacement capacity. Comparison with forcebased FE fiber models shows that, to obtain the same level of accuracy in the determination of ultimate displacement, a fine pier discretization is required, either in terms of integration points or FE number.
If axial force varies along the height, the assumption of uniform sectional response leads to an approximate solution. The comparison with force-based FE fiber model results shows that, at least for axial force variations usually associated with the weight of typical piers, the error induced by this approximation is practically negligible. This is confirmed by the comparison with the experimental tests by Babazedeh, that showed very good agreement between measurements and numerical results obtained by the proposed procedure.
The proposed procedure can be easily implemented without resorting to full FE codes and is suitable for accurate prediction of the entire force-displacement curve of slender RC piers, with, in particular, a reliable estimate of the ultimate displacement.
In this work, the analysis was limited to monotonically increasing M-K response curves, i.e., computation is performed until base pier cross section reaches its maximum bending capacity. Extension to more complex behaviors of the cross section is currently ongoing and will be reported in future publications. Funding Open access funding provided by Università degli Studi di Roma La Sapienza within the CRUI-CARE Agreement. No funding was received to assist with the preparation of this manuscript.

Conflicts of interest
On behalf of all authors, the corresponding author states that there is no conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.