A regularized higher-order beam elements for damage analysis of reinforced concrete beams

Abstract This work presents a numerical method for damage analysis of reinforced concrete beams using the higher-order beam theory based on Carrera unified formulation. The component-wise approach is employed to model the concrete and steel reinforcing bars as two independent one-dimensional finite elements. A modified Mazars damage model with tensile and compressive damage propagation laws is utilized for concrete, and an elastic-perfectly plastic law is used for steel rebars. To address the instability and mesh dependence caused by the strain-softening behavior of concrete, a fracture energy regularization technique based on the crack band model is developed, especially for the higher-order beam theory. The proposed method is validated by comparing its numerical results with three experimental benchmark results. The comparison indicates that the method accurately predicts the damage distribution of concrete and the flexural behavior of RC beams under quasi-static loading conditions while remaining computationally efficient.


Introduction
Reinforced concrete (RC) structures, which consist of concrete and steel, are widely employed in civil engineering.Many researchers utilized computational techniques such as the finite element method (FEM) for structural analysis and design instead of experimental campaigns.However, the elastic-plastic behavior of steel and the quasi-brittle nature of concrete poses challenges for non-linear numerical analysis of RC structures.Moreover, the accurate simulation of crack initiation and propagation is crucial in capturing the behavior of concrete.In RC structures, cracks may develop not only externally but also at the interfaces between steel rebars and the surrounding concrete, thereby adding complexity to the fracture analysis of RC structures.
The simulation of crack formation and growth in RC structures is crucial for evaluating their structural performance.Two commonly employed approaches to simulate the non-linear behavior of pure concrete structures or RC structures are the discrete model [1,2] and the continuum damage model [3,4].The discrete crack model allows displacement discontinuity by introducing interface elements to all element boundaries.However, the model depends on the mesh boundaries, and a re-meshing technique is required, leading to highly refined meshes [5].Consequently, the concept of strong discontinuity [6] has been developed to capture the arbitrary cracks in concrete, further contributing to the development of extended FEM (XFEM) [7,8], embedded FEM (EFEM) [9,10], and other related methods.
Alternatively, the continuum damage model is more popular due to its computational convenience.In this approach, concrete is treated as a continuum media, and the degradation of concrete stiffness represents the discontinuity caused by cracks.This method is easily combined with the analysis of rebars since physical discontinuous field deformation is explicitly required [11,12].Currently, the plastic deformation of steel can be simulated accurately through the FEM with Von Mises plasticity model [13].Several popular continuum damage models have been proposed and validated, including isotropic [14] and anisotropic damage models [15].For complex loadings, damage-plasticity models [16,17] have been proposed to consider both the inelastic and damage-dependent behavior of the material.
One issue associated with the continuum damage model is the occurrence of energy dissipation in localized regions rather than distributed zones when concrete strain softening occurs [18].This phenomenon can lead to the pathological sensitivity of the numerical results to the element characteristics.Various localization limiters have been proposed to mitigate the mesh-dependency issue, such as nonlocal models, including integral [19] and gradient models [20], where stress at a point depends on the strain at both the local and neighboring points.Another simpler limiter is the crack band model [21], which helps to restore the objectivity of numerical solutions by adjusting the concrete constitutive laws.This method introduces a characteristic element length to rescale the post-peak part of the concrete stress-strain law to avoid the tendency of fracture energy to approach zero.However, the correct way to determine the characteristic length for a given finite element remains an open question [22].
While traditional 1-dimensional (1D) models like the Euler-Bernoulli Beam Model (EBBM) and Timoshenko Beam Model (TBM) help capture the bending behavior of simple, slender structural members, their low computational costs are insufficient for capturing the material non-linearity due to cracking in RC heterogeneous structures.For detailed and accurate analysis of such structures, 3-dimensional (3D) models are required.Various 3D models have been developed in which 3D solid finite elements are adopted for concrete and beam or truss elements are employed for rebars [11,23,24].Interface models have also been added to account for bond-slip between concrete and reinforcing bars [25,26].However, the computational costs associated with these models are heavy, particularly for large-scale structures.Homogenization techniques [27,28] can be utilized to combine the behavior of concrete and steel to address these cost issues.
To balance the computational efficiency of 1D models and the elaboration of 3D models, a higher-order beam theory [29] based on Carrera Unified Formulation (CUF) [30] can be employed.According to CUF, a higher-order 1D beam model is utilized, and 3D displacement fields can be obtained by expanding the cross-section using various polynomials, such as Taylor [31] and Lagrange [32].In the framework of CUF, the orders of the approximation function for the beam element and the expanded function for the cross-section are considered as the input of analysis, eliminating the need for ad hoc assumptions.Lagrange expansion is preferred because it enables CUF to handle arbitrary geometries for both static and free-vibration analyses [33,34].Another developed method for analyzing composite structures with two or more parts, such as fiberreinforced structures, is the Component-Wise (CW) approach based on CUF [35].This approach divides the structures into different components based on the materials used, and each component can be modeled individually and simultaneously using Lagrange expansion cross-sectional elements.[36] has demonstrated the superior performance of the CW approach compared to other modeling approaches based on CUF, although only linear static analysis was conducted.Moreover, [37] has presented a damage analysis of notched RC beams based on 1D-CUF models under direct tension, which is not a typical loading condition in real engineering applications.
In this context, this study aims to develop a method for simulating the non-linear behavior of bend-dominated reinforced concrete (RC) structures using 1D-CUF models that accurately predict experimental results.To accomplish this, we combine a modified Mazars damage model, which is effective in pure concrete damage modeling [38,39], with the plastic behavior of steel rebars.To model cracks at the interface, we approximate them by locally initiating and evolving concrete damage around the steel rebars, avoiding the need for interface models.Additionally, we address the mesh dependence caused by concrete strain-softening by adopting a fracture energy regularization technique based on the crack band model and designing a consistent characteristic element length specifically for 1D-CUF models.Overall, the novelty of this work lies in developing an advanced numerical model based on a simple damage model for nonlinear analysis of RC structures, which can be easily implemented in actual engineering design.
This paper is organized as follows: First, some methodologies, including 1D higher-order beam theory, the CW approach, and a modified Mazars damage model, are presented.Next, three bending-dominated RC beams are simulated and compared with corresponding experimental results.Finally, some meaningful conclusions are drawn from the previous analysis.

Unified higher-order beam theory
The present work is based on a one-dimensional beam model derived from CUF and implemented using FEM [30].In the framework of CUF, 3D displacement of the 1D beam model can be enriched by expanding the corresponding cross-section, which is expressed as: uðx, y, zÞ ¼ F s ðx, zÞu s ðyÞ, s ¼ 1, 2, ::::, M where u s ðyÞ represents the generalized displacement vector of beam model; F s is the expansion function related to cross-section; s is an Einstein notation indicating summation and M is the number of terms in the expansion function.This work adopts Lagrange-type expansion (LE) for F s as it can handle any arbitrary cross-section.Quadrilateral elements such as four-node linear (L4), nine-node quadratic (L9), and sixteen-node cubic (L16) are typically used due to their higher accuracy.More detailed information on LE can be found in [30].
By discretizing the beam using the classical finite element method, the displacement field can be rewritten as: where u si is the nodal displacement vector, and N i is the beam shape function with N NE nodes per beam element.Typically, two-node linear (B2), three-node quadratic (B3), and four-node cubic (B4) beam elements are commonly used.The detailed shape functions for these elements can also be found in [30].However, it is essential to note that the choice of beam finite elements is independent of the selection of the class and order of the expansion function, which is one of the advantages of CUF models.
The governing equation for static problems can be obtained using the principle of virtual displacements.The detailed derivation process can be found in [39].For the sake of simplicity, the final equation is provided below: where F sj is the nodal external force vector, and K ssij is the fundamental nucleus of the stiffness matrix, which can be computed as: where l and X represent the length of beam element and area of cross-section, respectively; C is the material matrix; D is the differentiation operator; i, j, and s, s are indexes related to beam shape function and cross-sectional expansion function respectively.The integral in Equation 4 is obtained numerically using the Gauss quadrature technique, which is given explicitly in [30].

Component-wise approach
The reinforced concrete beam is a kind of composite structure that is made of steel and concrete (Figure 1a).An extension of 1D-CUF models, known as the componentwise approach, is utilized to model the concrete and steel individually and simultaneously.As shown in Figure 1b, the steel rebars and surrounding concrete are considered independent components using the same beam elements.The cross-section can then be discretized into multiple Lagrange elements with different materials, and the Lagrange points in the boundaries of other materials can ensure the continuity conditions among the various components, as shown in Figure 1c.Through this approach, the CUF model allows each component to maintain its independent material and geometrical properties.A refined model can also be determined for the interested component as needed.Some examples of the CW approach can be found in [36,37].

Modified Mazars damage model
This section discusses a modified Mazars damage model with fracture energy regularization technique to obtain the objectivity of RC damage analysis.Mazars damage model is an isotropic damage model for concrete described in [14].This model is based on the theory of elasticity coupled with damage mechanics.It employs a scalar damage variable to account for the uniform degradation of the stiffness properties in all directions.This variable depends only on the positive effective strains in the principle directions.
Compared to other damage models, such as concrete damage plasticity (CDP), the Mazars damage model does not consider permanent strains.However, it is adequate for damage analysis of RC structures under quasi-static load because the plastic behavior from rebars plays a more critical role in RC structures.Besides, when using fracture energy regularization for tension and compression, the structural response can get regularized and no longer be mesh dependent.Moreover, introducing this regularization can eliminate the need for fitting parameters.Therefore, the proposed model is ideal for practical structural design due to its simplicity and robustness.

Model formulation
Mazars [14] introduced a scalar variable for classical Hooke's law to consider the non-linear behavior of concrete material, which is written as: where r and e represent the stress and strain vector; d is the damage variable used to monitor the stiffness degradation.The positive principle strain governs the initiation of both tensile and compressive damage.Therefore, an equivalent strain e eq , as defined by Mazars, is calculated as: where hÁi þ is the Macauley bracket for picking out the positive value, and e i are the principle strains.Then the loading function is: f ðe, jÞ ¼ e eq ðeÞ À j where j is a threshold of damage growth that equals the ultimate tensile strain of material at first and then keeps updated as the e eq ðeÞ after first damage.

Damage propagation
Once f ðe, jÞ !0 in Equation 7, the damage is activated and determined by a combination of tensile and compressive damage responses.This can be expressed as: where a t and a c are weights that link the damage variable in tension d t and compression d c : The explicit calculation of a t and a c can be found in [38].
The damage variables are related to the internal variable j which equals to the equivalent strain under monotonic loading [22].While the original Mazars damage model [14] also provided damage propagation laws for d t and d c , these require fitting experimental results to obtain some parameters.In this work, modified damage evolution laws for tension and compression are adopted based on concrete constitutive laws (Figure 2) from fib MC2010 [40], which are more practical.Additionally, a fracture energy regularization technique based on crack band model is employed to regularize the softening behavior and prevent mesh dependency.
For simplicity and brevity, the damage evolution law of tension is shown as Eq 9, which is derived from a classical exponential softening constitutive law depicted in Figure 2a.
More detailed information on this derivation can be found [38].
where e d0 is the limit elastic strain, which is calculated by dividing the mean uniaxial tensile strength f ctm by Young's modulus E; p t is residual tensile stress ratio that is the ratio between the residual tensile stress and f ctm , which ensures the value of damage not equal 1.0 but infinitely close to 1.0; e tres is the corresponding residual tensile strain; e tu is the equivalent ultimate strain for bilinear softening, which is shown in Figure 2a and calculated as: where G ft is the fracture energy of mode I cracking.e tu can be adjusted to control the slope of the softening diagram through introducing the characteristic element length l c : Similarly, the damage evolution law of compression is derived from compressive stress-strain curve depicted in Figure 2b and the explicit formulation is expressed as: with where f cm is the mean compressive strength of the concrete; E cm is the secant Young's modulus; k, k 1 and k 2 are parameters from [41] to describe the softening part of constitutive laws; e c1 and e c2 are strain parameters adopted from [41], respectively; e c is a unidimensional strain ratio provided in [40].p c is residual compressive stress ratio that is the ratio between the residual compressive stress and f cm , which ensures stress never equals 0.0 to avoid convergence problem; e cres is the corresponding residual compressive strain; e cu is the extreme compressive strain, which is shown in Figure 2b and calculated by: where G fc represents the energy dissipation per unit area due to crushing.Similarly, the slope of softening part can be controlled by adjusting e cu through introducing the characteristic element length l c :

Characteristic element length
The characteristic element length l c depends on various aspects of mesh discretization such as element shapes, dimensions, interpolation functions, and so on [22].In the past, three main methods have been developed to estimate this length: (1) the method based on element area or volume [42], (2) the projection method [43], and (3) Oliver's method [44] which was extended to 3D linear elements by Govindjee [45].These methods are partially reviewed in [22].Recently, some modifications [46,47] have been made to Govindjee's method to improve the accuracy of l c estimation.However, these methods are limited to linear elements or only validated in 2D problems, making them unsuitable for higher-order beam elements based on CUF.
Given the aforementioned summary, this work employs a new method for estimating the characteristic element length, inspired by [47].The new estimation method is extensively discussed in [48].A simple example is employed here for illustration in Figure 3.One B3 element with three nodes is used to model the beam, while one L9 element with nine nodes is used for cross-section expansion.Then, A volume is assembled, as shown in the first two steps of Figure 3.According to the order of beam element and expansion element, the assembled volume is divided into eight small volumes in this case, each containing eight nodes similar to a linear solid element.Next, the middle points of each edge and the unit vector of the major principal strain at the target point are obtained and stored, assuming that the small gray shaded volume contains an interested point (Gaussian point).Finally, the projection method of [47] is conducted, as shown in the last step of Figure 3, and can be expressed as: with where x represents coordinates of the target point; x M means all coordinates of middle points from a small volume that containing the target point and nðxÞ is the unit vector of principal strain on the target point.

Numerical results
In this section, three experimental benchmarks are selected to assess the efficiency and accuracy of proposed method.One benchmark consists of an RC beam with only bottom rebars, while the other two consist of RC beams with stirrups for anti-shear failure.Quasi-static modeling is adopted for all numerical models through displacement-control method.

Leonhardt shear beam
The first experimental benchmark, named as Leohardt shear beam, is adopted from [49] to assess the proposed method.This beam is a reinforced concrete beam with a single bottom reinforcement, and its geometric information and boundary conditions are shown in Figure 4.The beam has two longitudinal rebars with a diameter of 26 mm, and the concrete cover is 37 mm.The loading process is displacement-controlled, with a maximum displacement value of 5 mm.The relevant material properties are listed in Table 1.
The hardening behavior of the steel rebar is assumed to be perfectly plastic.The half model, shown in Figure 5a is adopted in finite element analysis as the beam is symmetric across the y axes.The mesh discretizations on cross-section are shown in Figure 5b where there are three different sections, and the steel components are highlighted in red.Table 2 lists four models used to study the meshindependence of the proposed method.The discretizations of cross-sections are the same for all models, while Model C and D adopt linear elements and Models A and B use quadratic elements.Model A adopts higher-order beam elements with higher-order cross-section expansion, which should provide the highest accuracy due to highest number of degrees of freedom (DoFs).Compared to Model A, Model B and Model C investigate saving the computational costs by reducing the order of beam elements and cross-section expansion, respectively.Model D investigates convergence and mesh-independence by increasing the number of beam elements compared to Model C.
Figure 6 plots the displacement at the center of the beam versus the reaction force at one of the supports from four different models, including experimental one for comparison.Before the loading achieves around 36% of maximum displacement (1:8mm), all of the numerical curves match well with the experimental curve and are close to each other.
After this point, some small divergence occurs among the different numerical curves, but the stiffness of all numerical models remains similar.The experimental peak load is approximately 68:26kN with a corresponding displacement of approximately 3:1mm: All of the numerical models can achieve the peak value, but the corresponding displacements     are around 4:5mm: This difference is because the numerical stiffness of the latter part is clearly lower than the experimental stiffness.
To explain the above phenomenon, Figure 7 displays damage propagation during loading.For comparison, observed cracks in the experiment are plotted in Figure 8.Initially, damage occurs at the bottom and propagates vertically from Figure 7a-c due to beam bending.Following this, damage develops obliquely because of shear behavior, as demonstrated from Figure 7d-e.At this stage, the numerical stiffness exhibits a reduction, falling below the experimental value.Finally, horizontal damage is formed, connecting the support and the previously developed diagonal damage, as seen in Figure 7f.The final damage distribution is akin to the actual experimental cracks, except that a diagonal crack connecting the support and loading point was observed from the experimental campaign.However, it is essential to note that the proposed modified Mazars damage model only considers tensile and compressive damage based on fracture mode I and does not consider shear fracture.This drawback can account for the slight discrepancy between the numerical and experimental stiffness in the later stage.
Figure 9 illustrates the damage distributions at a displacement of 4 mm for different numerical models to investigate the objectivity of numerical models.It turns out that all models predict a similar damage distribution, with some vertical damage occurring in the mid-span and horizontal damage connected to one oblique damage in the side span.Combining this with the nonlinear performance shown in Figure 6, it can be concluded that the proposed method can provide mesh-independent results with any mesh configuration, even with Model C, which has less than 10,000 degrees of freedom.

Four-point bending RC beam
This experimental benchmark was reinforced with stirrups to resist the shear behavior.It was reported in [50,51] that the concrete collapse occurred due to the failure of concrete compression.Therefore, it is an important benchmark for validating regularized higher-order finite elements associated with compressive fracture energy.The geometric and boundary conditions are reported in Figure 10, with stirrups reinforced every 60 mm and a concrete cover thickness of 15 mm.The material properties are listed in Table 3.A bilinear model is used for steel initially, followed by an elastic-plastic model after ultimate strain.
Figure 11 illustrates the beam mesh and cross-sectional expansions used for the numerical analysis.A similar halfstructure is adopted as in the previous benchmark, but more beam elements are required due to the presence of stirrups, which result in a periodic change of the cross-section.To investigate the objectivity of numerical results obtained from different mesh configurations, four models listed in Table 4 are designed.The cross-sectional expansions are kept the same for four models to investigate the influence of the number and order of beam elements.
Figure 12 displays the load-displacement curves for various models, and the experimental curve is included for comparison.The numerical curves exhibit high similarity regardless of the mesh configurations adopted.The initial linear stiffness obtained from the numerical results agrees with the experimental value.However, the crack loads predicted by all the models are higher than the experimental value.This discrepancy is attributed to initial cracks or defects in the experimental specimens, also reported in [38].Subsequently, the stiffness at the inelastic stage with cracks predicted by the numerical results closely approximates the experimental campaign results.Finally, The numerical curves corresponding to phase III also show a good agreement with the experimental curve.Overall, these results demonstrate the robustness and reliability of the proposed numerical models, with mesh-independent results that accurately capture the material's behavior under both elastic and inelastic phases.Figure 13 depicts the damage propagation during the loading process, with all diagrams presented considering the deformed shape.Initially, multiple micro-cracks occur at the bottom of the beam, with the corresponding crack load being around 7:15kN, which matches well with the expected value from [38].As the loading continues, vertical damage propagates, and the steel rebars begin to yield, as shown in Figure 14, at 18% of the maximum displacement.The corresponding damage distribution is shown in Figure 13c.The concrete begins to crush at around 50% of the maximum displacement, as shown in Figure 13d, with damage also occurring at the top middle part.Last, Figure 13e presents the final vertical damage, which is similar to vertical cracks shown in Figure 15.

Three-point bending RC beam
The last benchmark is a three-point bending test of a concrete beam reinforced with longitudinal steel bars and stirrups conducted in [52].The geometric and reinforcement design of the tested specimen is shown in Figure 16.The beam is subjected to a single concentrated vertical load at   the top mid-span, with a maximum displacement control of 5 mm.The stirrups are spaced at 100 mm in the span and reduced to 50 mm near both ends (each along 500 mm).Only two longitudinal steel bars are placed at the bottom and top, and the concrete protective thickness is 20 mm.Therefore, the beam was characterized as under-reinforced and failed in ductile flexure, with the yielding of steel rebars followed by concrete crushing, as reported in [52].The beam is relatively large compared to the previous two benchmarks, making it ideal for validating the proposed method.Material properties of concrete and steel are listed in Table 5.
Similarly to the previous benchmarks, a half structure is adopted due to symmetry across the y axis.The beam element assignment and cross-section discretization are shown in Figure 17.Commonly, one beam element should be assigned along the thickness of the stirrups, and another one should be assigned along the space between two adjacent stirrups.Since more stirrups exist in this benchmark, more beam elements are required compared to previous benchmarks.
Due to the large number of stirrups in the structure, a minimum of 59 beam elements are required, as shown in Model I in Table 6.However, finer meshes with more beam elements are also presented, such as Model II and Model III in Table 6.Although only linear beam elements are used, the DoFs are still heavy.Therefore, higher-order elements are not investigated in this case.Figure 18 compares the load-displacement curves obtained from the three numerical models and the experimental result.The curves from the numerical models exhibit a high degree of similarity, providing evidence for the mesh independence of the proposed method.All numerical models capture the stiffness related to phase I and II well.However, the maximum loads obtained from all numerical models during phase II are slightly higher than the experimental value.The RC beam exhibited its weakest plane at the midspan in the experimental campaign.A stiff plate is added to the top of the FEM model to avoid stress concertation, which leads to the shift of the weakest plane away from the midspan, as shown in Figure 19.Additionally, a higher theoretical ultimate load of around 41:04kN, displayed in Figure 18, is expected if the weakest plane is assumed to occur 200 mm from the midspan.Therefore, the slight discrepancy between experimental and numerical results in the third phase is acceptable.
Figure 19 shows that the finer mesh can provide more detailed and accurate damage distributions compared to the coarser meshes.However, the overall damage patterns, such as the vertical damage caused by bending behavior, are similar for all numerical models.

Conclusions
This research work investigated the damage analysis of RC beams using 1D higher-order beam theory based on Carrera unified formulation.The proposed method utilized a component-wise approach to explicitly model the concrete and steel components.The modified Mazars damage model was adopted to evaluate concrete damage, while an elasticperfectly plastic material response was employed for steel rebars.A fracture energy regularization technique was adopted with a new estimation method for the characteristic element length to mitigate the mesh dependence of finite element models.The following conclusions were drawn based on the comparisons between numerical and experimental results from three benchmarks of bending-dominated RC beams: 1.In terms of load-displacement curves, the proposed 1D CUF model with the proposed material models can accurately predict both linear and non-linear behavior of RC beams under bending 2. Despite the lack of consideration for mode II shear fracture, the modified Mazars damage model can still capture the similar non-linear behavior of RC shear beams compared to experimental campaigns.3. The fracture energy regularization technique, with the new estimation method for the characteristic element length, effectively mitigates mesh dependence, allowing for computational costs to be saved by adopting coarser meshes.
Although accurate and robust results were obtained through the proposed model, a large number of DoFs was required if many stirrups were imposed, as seen in the third numerical example, which caused an increase in computational costs.Therefore, ongoing studies will explore adopting a node-dependent kinematic approach [53] to address this issue.

Figure 1 .
Figure 1.An illustration of CW approach for modeling RC structures: a) definition of RC structures, b) individual steel and concrete component, and c) assembled cross-section with Lagrange elements along with the beam.

Figure 2 .
Figure 2. Concrete constitutive laws for modified Mazars damage model: a) tension and b) compression.

Figure 3 .
Figure 3.An illustration of estimation method for l c :

Figure 6 .
Figure 6.Load-displacement curves of Leonhardt shear beams from different models.

Figure 8 .
Figure 8. Cracks observed in the experiment of Leonhardt shear beam.

Figure 14 .
Figure 14.Von Mises stress distribution of steel rebars at 18% of the maximum displacement.

Figure 12 .
Figure 12.Load-displacement curves of four-point bending RC beams from different models.

Figure 15 .
Figure 15.Crack distribution from the experiment of four-point bending RC beam.

Figure 18 .
Figure 18.Load-displacement curves of three-point bending RC beams from different models.

Figure 19 .
Figure 19.Final damage distribution of: a) Model I, b) Model II, and c) Model III.