Nonlinear analysis method of concrete structures under cyclic loading based on the generalized secant modulus

: A smeared crack model to represent cyclic concrete behavior is presented in this work. The model is based on analytical and experimental studies from the literature and proposes a numerical approach using a new concept, the generalized secant modulus. The monotonic formulation is described, followed by the changes to include the cyclic response, and the stress-strain laws to reproduce the hysteresis. Simulations adopting the proposed model were compared with experimental tests of cyclic tension and compression available in the literature, resulting in consistent load cycles. Three-point bending was simulated to display the structural response under non-elementary load. Finally, a reinforced concrete beam was studied to evaluate the model performance under usual loadings. The results show the model capacity to reproduce cyclic analyses and its potential to be extended to general loadings


INTRODUCTION
The complete description of a concrete structure subjected to mechanical loadings demands constitutive models capable of representing loading, unloading, and reloading regimes. Even under monotonic loading, the structural equilibrium can result in a stress redistribution related, e.g., to microcrack nucleation [1]. In other situations, the external load solicitation is cyclic [2]. Thus, cyclic constitutive models are required to represent the structural behavior.
where β r is the shear retention factor and G 0 is the shear modulus of the undamaged material.
Then, the secant moduli are calculated in function of the principal strains. At this point, stress-strain laws are required to obtain the elastic modulus degradation. Particular relations may be necessary to represent the concrete behavior under tension and compressive loads.
The total relation for the stress state calculation is symbolically written by Equation 3.
ℓ s . {ε ℓ }. (3) Or, in the matrix form, by Equation 4. (4) where , and are the stress components, while , and are the strain components in the local crack system. In an iterative-incremental solution based on the standard Newton-Raphson method, the tangential constitutive relation is also required, derived from the total relation, as expressed by Equation 5. The model was presented in the local system. The formulation in global systems is given by Equation 7 and Equation 8.
Where � ℊ � is the global stress vector, � ℊ � is the global strain vector, [ ] and [ ] are the stress and the strain transformation matrices, respectively. Thus, by replacing Equation 7 and Equation 8 in Equation 3, the Equation 9 is determined.
Since the global matrix is given by Equation 11 (11) the total relation in the global system is written by Equation 12.
The global tangent operator comes from the derivative of total relation in the global system to the strain vector, resulting in Equation 13. (13) where is the angle direction of the local crack system. If the local system is fixed during the crack nucleation, the transformation matrix does not present variation. Consequently, the global tangent matrix can be defined as Equation 14. (14) The secant modulus monitoring can establish the loading regime. Although, a more efficient control method based on loading functions is adopted. The loading functions are expressed in terms of strain components ( , ) by Equation 15 and Equation 16. (16) Where and are, respectively, the loading function in the normal and tangential directions to the crack plane; and are the strain components; and are historical parameters, i.e., the maximum values of and during the analysis. The initial values of the historical parameters ( 0 , 0 ) represents the strain in the elastic limit. When 0 or 0 are exceeded, the degradation starts, and the material behavior turns inelastic in this direction.
Regarding the loading regime, the loading function variation ( ̇) and the historical variable evolution () are calculated by Equation 17 and Equation 18, respectively, in a pseudo-time of the iterative-incremental process.
Since the regimes are known, the secant moduli ( , ) can be calculated by stress-strain laws. Such moduli reproduce the material degradation and are part of the constitutive matrix that provides the stress. In cyclic analyses, the secant moduli need to be reformulated for the unloading/reloading. Nevertheless, the historical variables remain the same during hysteresis.
A change in the constitutive matrix (Equation 1) has been proposed to perform the cycles. The generalized secant moduli , replaces the traditional secant moduli , in unloading/reloading stages. The stress-strain laws that provide the secant moduli must also be redefined to include cyclic paths and the calculation. Figure 1 summarizes the process.

Stress-strain laws for cyclic loading
The constitutive models are mechanical-mathematics representations that describe the behavior of a material media. The degradation is computed by the secant modulus variation for the proposed model, using stress-strain laws defined from approximations of experiments. Such laws can assume different formats: linear, polynomial, exponential, among others.
is the stress, is the strength limit, is the strain related with the elastic limit, is the current strain, and = represents tension, while = represents compression.
Where σ is the stress, f t is the tensile strength, ε is the current strain, ε t is the strain in the elastic tension limit, h is the characteristic length, and G f is the fracture energy. The stress-strain relations, although are restricted to monotonic loadings. In this work, four hypotheses have been incorporated in the original laws to reproduce the load cycles: the secant cycle, based on the traditional secant modulus; the linear cycle, using the undamaged elastic modulus; the linear cycle based on focal point concept; and an extended version presented by Bono [29] of the law initially proposed by Popopovics-Saenz. All cyclic strategies adopted the generalized secant modulus proposal detailed as follows.

Secant unloading/reloading
The secant hypothesis is the most common representation of concrete behavior because of its simplicity and precision in reproducing monotonic loadings. Another relevant feature is the maximum strain and the current material degradation monitoring. The degraded modulus can be calculated as the straight slope that starts in origin and goes until the unloading point. Figure 2 illustrates the secant hypothesis. The point with coordinates εr, σ r is a generic point where the reloading begins, that is the origin of the system σ-ε for a complete cycle.  Where the unloading strain (ε u ) is the maximum strain the material has ever experienced, the historical parameter. Despite the simplicity, the secant strategy loses representativity during cycles with an advanced strain stage since the permanent strains are not considered.

Linear elastic unloading/reloading
The linear elastic approach is based on the undamaged elastic modulus (E 0 ), as shown in Figure 3. The cycle is simplified by a straight line whose slope is E 0, which goes from the unloading point to the strain axis. The E GS is required during the unloading/reloading regimes. Such modulus is represented by a straight line that starts in origin and intersects the cyclic path. The E GS can be expressed as Equation 22. ) where ′ is the current strain in the cycle, and ′ is the stress associated with ′ ( Figure 3). The current stress ′ (Equation 23) is the unloading stress subtract from a stress variation Δσ and must be written as a function of the known parameters E 0 , , , ′ .
The unloading strain ( ) can be divided into elastic and permanent parcels, as indicated by Equation 24.
= + . (24) As the cyclic path follows the undamaged elastic modulus, the stress variation is calculated with Equation 25.
= 0 ( − ′ ). (25) Replacing in Equation 23, the stress at any point of the cycle is given by Equation 26.
The linear elastic approach represents permanent strains, neglecting stiffness degradation. The E SG is used in the cycles, but the degradation remains fixed as the historical variable.

Linear unloading/reloading based on the focal point
In experiments, Lee et al. [22] and Lee and Willam [27] observed that the material unloading is directed to a pole. As presented in Figure 4, the cyclic path is a straight line defined by the unloading point (ε u , σ u ) and the focal point (ε f , σ f ). In the proposed method, the first step is the linear modulus (E L ) calculus. This modulus represents the unloading/reloading straight slope given by Equation 27.
. (27) Then, Equation 27 is rewritten as Equation 28 to provide σ' At last, the generalized secant modulus is obtained as Equation 29.
The linear path based on a focal point improves cyclic material behavior representation since it can reproduce both stiffness degradation and permanent strain.

Nonlinear unloading/reloading
The previous focal point law has been extended to distinct unloading/reloading paths. Bono [29] adapted the original monotonic formulation of Popovics-Saez, introducing modifications in the origin and peak coordinates to make the cyclic analysis feasible. The Popovics-Saenz law proposed to describe the concrete monotonic compression behavior is calculated with Equation 30.
Where K = E 0 (ε i / f i ), R = K / (K -1); f i represents concrete strength, ε is the current strain, and ε i is the strain related to f i . Physically, R is the reason between the elastic stress related to ε i and the difference of this stress and f i . A, B, C, D are variables dependent on the branch curve and the adjustment parameters β ref and k ref .
Kwon and Spacone [41] defined the focal point for compressive loading as the tensile strength. Under tension loading, an adequate monotonic curve should be adopted. The unloading/reloading paths, although, can be used without restrictions. In this case, the focal coordinate is the compressive strength. The E GS is calculated as shown in Equation 29.

Cyclic monitoring
The iterative incremental method requires a technique to control the load cycles imposed on the numerical model. Therefore, two monitoring methods have been developed to inform the start and end of the loops: the step count method and the load limit method.
The step count method ( Figure 5) needs the step numbers where the unloading starts and ends. Despite its efficiency, the monotonic path must be previously known, and the load script depends on the increment size. On the contrary, the load limit method ( Figure 6) defines the cycles from the load values at the beginning of unloading and reloading. Tolerance is required since the equilibrium path points may not coincide with the specified limits. The load limit method is practical if the cycles loci are known from experimental data.  Both monitoring methods operate at the displacement increment. To perform the unloading, once the critical point is identified into the Newton-Raphson algorithm, such increment is multiplied by -1.0. The reloading occurs with a second sign inversion. At this moment, the loading is resumed, driving back to the monotonic curve. For material points under the elastic regime, the unloading/reloading path agrees with the monotonic curve. This strategy is applied to the direct displacement control method, and there is no cycle number limit. The algorithm presented in Figure 7 summarizes the cyclic monitoring procedure.

Direct tension
Golaparatnam and Sha [42] performed a series of experiments in plain concrete samples under cyclic tension with prismatic geometry (76 x 19 x 305 mm). The concrete properties are tensile strength f t = 3.53 MPa, tensile strain limit in the elastic regime ε t = 1.18 x 10 -4 , and fracture energy G f = 0.0564 N/mm.
Initially, a refinement test was performed in a monotonic analysis to define the mesh sensitivity. In the simulation, five meshes were analyzed, with 1, 4, 16, 60, and 320 elements. Since no evidence of localization effects was verified, the mesh with 60 four-node quadrilateral finite elements was adopted, considering an intermediate level of refinement.
The loading process assumed a reference load q = 19 N/mm, displacement increment of 5.0 x 10 -7 m, and tolerance for the convergence in displacement of 1.0 x 10 -4 . The cycles were introduced by the load limit monitoring method. Figure  8 illustrates the geometry configuration, the boundary conditions, the mesh, and the controlled degree of freedom. The material law by Boone et al. [25] and Boone and Ingraffea [26], and the Carreira and Chu [23], [24] were adopted to tension and compression, respectively. Besides the experimental parameters [42], other material properties were considered: Poisson ratio of 0.18 [43]; compressive strength of 40.4 MPa estimate from the relation [44], indicated in Equation 31.
⁄ ; (31) and the elastic limit of the compressive strain ε c = 0.002, from the NBR 6118 [45]. The characteristic length (h) is a parameter associated with the fracture process zone. Such value is usually related to the maximum aggregate size [46]. In the absence of information, (h) can be obtained from the region dimension used to measure the experimental strains. Here, it was considered h = 166 mm. The shear retention factor was admitted as βr = 0,05, a traditional value for plain concrete [47]. The maximum secant moduli degradation was limited to 92%. This number came from the monotonic path fitting to the experimental curve. For the linear unloading/reloading approach, the focal point is the concrete compressive strength. For the nonlinear cyclic strategy, the focal point of the stretch direction (εf1, σ f1 ) is also the compressive strength, while the compressed direction focus (ε f2, σ f2 ) is the tensile strength [41]. An extra simulation was performed for the nonlinear formulation considering the origin as the focus to both directions. Parametrization of the law by Popovic-Saez was done to fit Carreira and Chu [23], [24]  The numerical stress-strain relations are shown in Figure 9, contrasting with the experimental curve from Gopalaratnam e Shah [42].
Since the secant strategy is based on an elastic degradation model, all the cycles are directed to the origin, and its priority is stiffness deterioration reproduction. The agreement between the simulated cycle initiation and the unloading points from the experimental curve was used to achieve representative results. The permanent strain is represented in the elastic approach. Thus, the reloading points were prescribed to be as close as possible to the experimental residual strains. The stiffness degradation, however, is neglected.
The cyclic laws based on the focal point embrace both phenomena: stiffness degradation and permanent strain representation. The cyclic path remained linear for the simulation using the nonlinear law where the focus is concrete compressive strength, presenting a similar result to the elastic strategy. This fact is related to the distance between the unloading points and the focus. Although, this approach showed its potential to represent hysteresis with distinct unloading/reloading paths when the focal point was set in origin.
As observed in the comparison between experimental and numerical curves, the last cycles have presented a less accurate representation. This result is related to the high level of degradation [1], [12], [13]. Therefore, a possible solution includes a variable focal point, depending on the current material degradation.
The last analysis evaluated the potential of a variable focus. The focal point coordinates of each cycle coincide with the null stress and the plastic strain related to that loop in the experimental curve. Figure 9 shows that the focal point strongly affects the configuration of the cycles, permitting a better adjustment of the numerical response.
Despite the restrictions of the adopted stress-strain laws, described as a linear path in most cases, they have represented concrete cyclic behavior satisfactorily. The approaches based on the focal point are more efficient in reproducing experimental results since they couple stiffness degradation and residual strain. However, analyzing such phenomena individually and disregarding the configuration of the cycles, the elastic strategy has been the most adequate approach to describe permanent strain. The secant representation has shown numerical robustness using the actual secant modulus. Based on these issues, the cyclic law selection must be guided by the drawbacks and the advantages of each law.

Direct compression
Based on the experimental tests by Karsan and Jirsa [11], a cyclic compression was modeled. A refinement test was performed to define the mesh of 4 x 4 four-nodes quadrilateral ( Figure 10). A load of q = 19 N/mm, a displacement increment of -5.0 x 10 -7 m, and a tolerance convergence in load and displacement of 1.0 x 10 -4 drove the loading process. The cycles were controlled by the load limit method.
The material parameters from the experimental test are: Poisson ratio ν = 0.18; elastic modulus E = 31,700 MPa; fracture energy G f = 0.04 N/mm; concrete compressive strength f c = 27.6 MPa; and characteristic length h = 82.6 mm. The other parameters were estimated as: concrete tensile strength f t = 2.76 MPa (Equation 31); strain in the compressive elastic limit ε c = 0.0016 [45]; shear retention factor β r = 0,05 [47]. The concrete tensile strength parameter (f t ) is necessary to determine the focus of the linear law based on the focal point. Such focus was also adopted to the compressive direction for the law by Popovics-Saenz [28]- [31]. The concrete compressive strength, otherwise, was defined as the focal point for the tension direction. The parameters of the law proposed by Popovics-Saenz's were established to fit Carreira and Chu [23], [24] monotonic curve. The secant moduli maximum degradation was 90%. Table 2 summarizes all the parameters required by the stress-strain laws. The responses (Figure 11) highlight the conclusions obtained in the compression simulations. The secant law is restricted to degradation stiffness representation, while the elastic law can only reproduce residual strains. The results from the linear law based on the focal point show both phenomena, although the cycles remain simplified by a straight. Nevertheless, the cycles reproduced by the nonlinear law, considering the focus either in the concrete tensile strength or origin, showed the capacity to describe the experimental cyclic paths. Although, a discrepancy between the hysteresis loops from the test and the numerical simulations is noted. Such differences are related to the parametrization difficulty of the cyclic models. Since monotonic tests define the variables of the stressstrain laws, a relation among these parameters and the unloading/reloading regimes is unknown. The only variable that may change the cyclic configuration is the focal point.
In general, Popovics-Saenz [28]- [31] law with the focal point on the concrete tensile strength is a valid representation of the experimental curve. The proposed model could present a more precise reproduction of the cycles with an accurate definition of the focal point. Like the tension analysis, a numerical simulation was performed considering a variable focal point. Each cycle was directed to the intersection point between the empirical loop and the strain axis. The shapes of the load cycles fit the experimental curve better, although it is still necessary to determine a variable capable of controlling the hysteresis amplitude.

Three-point bending test
This simulation modeled a three-point bending to evaluate the model response in more general load conditions. The geometry of the beam is illustrated in Figure 12, as in the reference work of Hordijk [48]. The node and the direction controlled are also highlighted. Because of the geometric symmetry, only half of the beam was modeled. The finite element model consists of 180 four-nodes quadrilateral elements in a plane stress state.
The loading process has assumed a unit load P, a displacement increment of -2.0 x 10 -3 m, and a convergence in displacement of 1.0 x 10 -4 . The monitoring method by load limit was selected. The plain concrete is characterized by the parameters in Table 3, as given by [48]. The reference curve is the experimental test of Hordijk [48]. Variable focal points were adopted in the linear law based on the experimental reload points. The origin was the focus of the Popovics-Saenz law. The results are shown in Figure 13. The secant strategy showed robustness in cyclic reproduction, while the other approaches exhibited instabilities. The cyclic elastic approach reproduced only the first cycle, and not even the step size reduction allowed it to conclude the simulation. The laws based on the focal point reached a more advanced strain state in the analysis, but they failed after the third cycle when numerical instabilities interrupted the simulation by not meeting the convergence.
These instabilities must be related to the E GS in complex load states since the material points can present load sign inversion, even when an inversion is not observed in the structural behavior. A quadrant mapping of the stress-strain space is required to consider negative and positive quantities, improving the proposed model to embrace general load cases. The focal point is another subject that deserves attention in future works because of the dependency of the load loops on this variable. Nevertheless, the results show that all cycle hypotheses agree with the behavior observed in direct tension and direct compression from the literature.

Reinforced concrete beam subjected to uniform loading
A reinforced concrete (RC) beam subjected to a uniformly distributed load was modeled. Since RC is a material compound of concrete and steel, a combination of cyclic behaviors was adopted. The secant cycles were selected to represent the concrete elastic degradation, while the elastic cycles were adopted for steel to reproduce plastic strains during the yielding.
The structure dimensions were based on typical beams from projects ( Figure 14). The mesh has 420 four-node quadrilateral concrete elements, and two-node bar elements represent the reinforcing steel. The material parameters are summarized in Table 4. The steel behavior was considered elastic-perfectly plastic with a perfect bond between steel and concrete.
The vertical displacement in the midspan of the beam was controlled with an increment of -1.0 x 10 -4 m, convergence tolerance in displacement of 1.0 x 10 -5 , and a distributed loading of q = 0.33 MN/m. The monitoring method by step count was selected, introducing three arbitrary cycles in the equilibrium path after the steel yielding to reproduce an intense unloading/reloading regime. The secant cycles were performed considering two different compressive laws for concrete: Carreira and Chu [23], [24], and Popovics-Saenz [28]- [31].
The structural response is shown in Figure 15, and the critical degradation stages are highlighted along the equilibrium path. As a reference, a curve is presented based on the reinforced concrete stages representing the elastic behavior of the materials, the concrete cracking phase, and the yield of the reinforcement. The microcrack nucleation (1) is observed when the graphic becomes nonlinear, and concrete starts to damage. Then, the microcracks grow (2) until a critical point in which the nonlinearities are significant. Thus, the microcracks coalescence culminate in macrocracks (3). In the last stage, the macrocracks propagate in concrete, and the steel reaches the plastification (4).   The results from the simulation with the law by Carreira and Chu [23], [24] show an elastic cyclic behavior, while the simulation using the material description by Popovics-Saenz has a response between the elastic and the secant approaches. The selected concrete monotonic curve changed the cyclic configuration, but the steel conducted the unloading/reloading process in both analyses during the yielding stage.

CONCLUSIONS
A smeared crack model based on elastic degradation under cyclic loadings has been presented, aiming at a generalized secant modulus numerical approach. Four stress-strain laws with unloading/reloading descriptions have been adopted. The proposed method was validated by comparing numerical results of concrete under cyclic loading in compression and tension with experimental curves from the literature, resulting in good agreement. Three-point bending tests have been simulated for plain and reinforced concrete to evaluate the viability of applying this model to general cases. In this context, the main conclusions are: i. The numerical analyses have shown that the developed method based on the generalized secant modulus fulfill the proposed aim, simulating load cycles with distinct characteristics; ii. The examples validated by experimental curves prove the effectiveness of the strategies developed for concrete structures subjected to cyclic tension and cyclic compression; iii. The implementation of the cyclic stress-strain laws has been successfully realized, and the variety of material laws allows versatility in the hysteresis representation; iv. The approaches based on the focal point have the advantage of reproducing both stiffness degradation and residual strains; v. The bending analysis has evaluated the present formulation application to simulate structures under conventional loads. Since numerical instabilities have been observed, an investigation must be carried to improve the proposed model to embrace general loads; vi. The proposed methodology seems capable of simulating reinforced concrete structures, coupling the developed cyclic model for concrete with a plasticity model for steel.