1.Introduction

As is well known, cracks can occur in tension zones of reinforced concrete structures even under relatively light loads. These cracks, when further loaded, significantly affect the stress distribution and strength of the structure. The crack width affects the rate of reinforcement corrosion and the durability of reinforced concrete structures. Numerous scientific studies have been devoted to modeling the process of crack formation and opening in the analysis of reinforced concrete structures. In [1], a continuous-discontinuous method was used to model cracking deformation and rock contact, taking into account the internal cracking of quadrangular elements. This method divides the computational domain into quadrilateral elements with high accuracy. It combines the advantages of the continuous and discontinuous methods. Before crack formation, this method does not require the introduction of normal and tangential contact stiffness. After formation, cracks can propagate along the element boundaries and diagonal lines of the quadrilateral elements. In [2], a crack model was proposed that focuses on calculating the crack spacing based on the reinforcement strain profile. For small beams, a linear crack width profile was assumed. For primary cracks in large beams, a bilinear width profile was used. The crack spacing and width predictions obtained using the proposed model are in good agreement with test results.

One of the most widely adopted approaches for modeling crack formation is Extended Finite Element Method (XFEM) [3–6]. Within this framework, additional degrees of freedom along with corresponding enrichment functions are introduced at the nodes of the cracked element, utilizing the Heaviside step function to represent the displacement jump across the crack faces. The strain fields within the finite element are subsequently determined through the derivatives of these displacement approximation functions. In [5], this methodology was successfully employed to analyze massive structures while accounting for crack propagation. Furthermore, [7] integrated XFEM with the substructure technique for crack modeling. By partitioning the structure into multiple substructures, XFEM could be efficiently applied to capture localized cracking within each subdomain.

Papers [8, 9] focus on experimental investigations of crack propagation in reinforced concrete beams, evaluating crack widths and spacing in comparison with standardized design values. In [10], a novel family of finite elements based on deformation analysis was introduced to capture material softening in planar frame structures. Within this framework, the softening zone is represented by a discrete crack treated as an ‘excluded’ point of the finite element. Consequently, crack strains are computed separately from the bulk element strains, being linked to the element solely via kinematic variables that define the crack opening. To further advance structural modeling, a rectangular finite element formulated through a complementary energy functional was proposed in, where nodal equilibrium equations were integrated using the Lagrange multiplier method and derived via the principle of virtual displacements. Additionally, a new approach for modeling frictional contact within XFEM was developed in. This work introduced an innovative nodal scheme incorporating a coupled Coulomb friction model. The primary novelty of this technique lies in its unprecedented ability to account for stress redistribution within localized regions of the crack surface, offering a distinct advantage in mitigating the blocking problem during numerical simulations of frictional contact.

[14] introduces a finite element formulation for the analysis of reinforced concrete slabs, accounting for both concrete cracking and reinforcement yielding. Developed within the framework of localized damage mechanics, this formulation assumes that all inelastic effects are concentrated along discrete lines, directly linking all numerical parameters to the mechanical properties of the slab. Furthermore, this finite element approach has been numerically validated to yield mesh-independent results. Alternatively, an approach for modeling discrete cracks in three-dimensional space using a particle method was described in. Under this framework, cracks can be arbitrarily oriented, and their propagation is discretely represented by activating crack surfaces on individual particles, thereby eliminating the need for explicit crack topology tracking.

In [16], an adaptive crack propagation algorithm utilizing the multifunctional super singular element  method was proposed. This approach yields highly accurate numerical solutions for singular stress fields and stress intensity factors, thereby enabling precise trajectory predictions. The implementation of adaptive mesh refinement technology allows for the automatic generation of background meshes for complex geometries based on three distinct intersection strategies applied within the update zone, significantly reducing the complexity of mesh subdivision. Similarly, the primary objective of study [17] was to develop a local multiscale mesh refinement strategy. To capture singularities accurately during crack growth, a localized refinement scheme around the crack tip region was implemented. The formulation integrates elements with a variable number of nodes and their corresponding shape functions to effectively resolve the mismatch problem arising from adjacent elements of different sizes.

Various plastic flow theories are frequently employed to model the elastic-plastic behavior of concrete under loading [18–20]. In these theories, the relationship between stress increments and plastic strains is governed by specifying yield and loading surfaces. This approach enables the development of effective incremental-iterative algorithms for analyzing reinforced concrete components under diverse loading scenarios. Under active loading, however, a simpler orthotropic concrete model can be utilized, such as the one outlined in the CEB-FIP Model Code 1990.

In [21], an innovative approach to describe crack evolution in quasi-brittle materials, modeled in the framework of isotropic damage mechanics, has been presented. Based on innovative Finite Element Method – Virtual Element Method coupled discretization of the body, this technique gives satisfactory results. They are comparable to the results obtained by well established approaches in the literature, where crack propagation is correctly caught by usually employing remeshing procedures or adding degrees of freedom. The aim of the work [22] was to reduce the computing efforts in the framework of the Global Cracking Elements Method, which involves finite elements with nonlinear interpolation of the displacement field. This aim was achieved with the help of hybrid linear and nonlinear interpolation. The numerical examples have shown that the proposed approach saves a lot of computing time as compared to the original global cracking elements. Besides, with this approach, curved boundaries can be more accurately modeled by adding edge nodes. In [23], to solve the mesh dependency problems of the cracking elements method with a 6-node triangular element (CEM-T6), the numerical formulation and procedures of CEM-T6 were investigated. The characteristic length used in the original CEM-T6 was found to be inaccurate and inconsistent. Hence, a new form of the characteristic length that depends only on the size and shape of the T6 element was presented. With the new form, the mesh dependency of CEM-T6 was greatly mitigated.

The study aimed to construct a finite element model of bending reinforced concrete beams taking into account crack formation. An incremental orthotropic model was used for modeling the elastic-plastic deformation of concrete. A stepwise iterative algorithm was used to solve the nonlinear problem.

2.Methods

To determine the location and opening width of cracks in a bending reinforced concrete beam, the plane elasticity problem was formulated under physically nonlinear conditions. An orthotropic model was applied for the stress-strain relationship, following the CEB Model Code 1990 guidelines. Under this model, the concrete stresses under uniaxial compression (Fig. 1a) were calculated:

                                            (1)

a)

b)

Figure 1. The stress-strain diagrams for concrete: a) under compression; b) under tension.

 is the initial modulus of elasticity of concrete;  = –0.0022;    is the compressive strength of concrete,  is the tensile strength of concrete,  is the ultimate deformation under compression.

                                     (2)

Figure 2. Biaxial strength of concrete (fcm = 37 MPa).

Under biaxial stress condition, the concrete strength was evaluated using Equations (3)–(5), which represent failure criterion curves (Fig. 2). In case of biaxial compression or tension-compression at

                                                       (3)

under biaxial tension

                                                                   (4)

under biaxial tension-compression at

                                                              (5)

     are the greatest, mean and least principal stresses.

In the coordinate system aligned with the principal stress axes, the relationship between strain and stress increments was determined as follows:

                                           (6)

The tangent deformation moduli    expressed as functions of the principal stresses, were determined using the concept of equivalent uniaxial strains

                                                                (7)

                                                                     (8)

 is the number of stress increments.

Tangent modulus of elasticity:

                            (9)

Equation (7) can be expressed in matrix form as follows:

                                                     (10)

To determine the tangent moduli  based on equivalent uniaxial deformations  it is necessary to replace  in Equation (1) with the corresponding value of the ultimate concrete strength

By substituting  for  in Equation (3), the concrete strength value was obtained.

                                               (11)

Similarly, by substituting  for  in Equation (3), the following relation was obtained:

                                              (12)

This approach implies that at each iteration step, the concrete strength in the principal direction under consideration is determined based on the current principal stress in the orthogonal direction.

For the biaxial tension-compression state  the following relations were derived in accordance with Equation (5):

                                                           (13)

                                                            (14)

The strain value  must be replaced by the equivalent strain  accompanied by the corresponding Poisson’s ratio  In the case of biaxial compression or tension-compression at  the following relation was obtained:

                        (15)

Under biaxial tension

                                                   (16)

Under biaxial tension-compression at  the following relation was obtained:

     (17)

                                   (18)

The shear modulus was determined using Equation (19).

                                                          (19)

In Equation (6), the relationship between the shear stress and shear strain increments was incorporated as follows:

                           (20)

The following notations were introduced for the vectors and the concrete stiffness matrix:

       (21)

Equation (20) can be expressed in matrix form as follows:

                                                                (22)

The prime symbol indicates that the equation is formulated in the principal axes

The principal stresses and the inclination angle of the crack normal are evaluated as follows:

                                                (23)

                                                            (24)

A crack forms in a finite element when the principal stress satisfies the following condition:

The vector of stress increments  in the global coordinate system  was expressed via the matrix

                        (25)

Using Equations (22) and (25), the strain vector  in the global coordinate system can be expressed as follows:

                                               (26)

Within the finite elements where a crack has formed, the concrete deforms elastically, and the matrix  is given by:

                             (27)

The stress-strain relationship for the reinforcement is represented by a bilinear diagram as follows:

                                          (28)

where  is the yield strength of the reinforcement,  is the modulus of elasticity of the reinforcement steel, and  is the hardening modulus.

3.Results and Discussion

Figure 3 Finite element model.

The stiffness matrix of the rectangular finite element was derived by approximating the stress fields within the element domain using linear functions [12, 24, 25]:

                                            (29)

The local coordinate system was centered at the finite element origin. The stress functions (29) satisfy homogeneous differential equilibrium equations in the finite element domain for plane elasticity theory. In [24], such approximations were used to obtain the stiffness matrix of a rectangular finite element with a linear stress-strain diagram of the material.

Using Equation (29), the tangent stiffness matrix of the finite element is derived, accounting for physical nonlinearity:

                           (30)

After a crack had formed in the finite element, additional unknown displacements were introduced into its nodes (Fig. 4b), and the stresses on each side of the crack were approximated independently, yielding:

                    (31)

                                              (32)

The superscript in Equation (31) denotes that the corresponding parameters are associated with the finite element subdomains  or  (Fig. 4a).

a)

b)

Figure 4. a) The variants of the crack cd location; b) the unknowns’ number in a node, depending on the crack location [24].

Increments of finite element complementary strain energy:

                                     (33)

where  is the element thickness.

The matrix  (26) for a finite element without a crack is calculated at its geometric center point; its elements depend solely on the parameters      and are assumed to be constant. The tangent matrix  is calculated numerically using the 3-point Gauss quadrature rule:

                              (34)

   are the weighting factors. The expressions for the matrix components of a cracked finite element were provided in [24].

To this end, the equilibrium equations derived via the principle of virtual displacements were incorporated into the complementary energy increment (Fig. 4) following the approach in [24]:

                                          (35)

where  is the element nodal displacement increment vector,  is the element nodal reaction vector. Expressions for the matrix  elements, which depends only on the geometric parameters of the finite element, are given in [24]. Each node has two primary and two additional displacements. Based on these nodal displacements, the crack opening width within the finite element was evaluated (Fig. 5). The expressions for evaluating the crack opening width were provided in [24].

a)

b)

c)

d)

Figure 5. The basic δu and additional δucr possible displacements of nodes after the crack cd forming: а) δu4,x(y); b) δucr4,x(y); c) δu3,x(y); d) δucr3,x(y). The figure is given in [24].

For the iterative scheme, the iteration index  was incorporated into Equation (35):

                                   (36)

Expression (36) was transformed as follows:

                         (37)

where  is the vector of unknowns from the previous iteration. Equating the derivatives of  to zero, yields the following equations:

                                (38)

Solving system of Equation (36) yields:

                                             (39)

                                                         (40)

                                                             (41)

where  is the finite element tangent stiffness matrix.

Figure 6. Reinforcement modeling within the finite element.

If the reinforcement passes through the finite element (Fig. 6), its stiffness matrix must be incorporated into the element stiffness matrix defined by Equation (41). The bond slip of the reinforcement relative to the concrete was neglected, therefore, the end displacements of the reinforcement were expressed in terms of the nodal displacements of the rectangular finite element. Thus, the following relation was obtained:

                                                                (42)

                                   (43)

The reinforcement stresses are constant along the element length and  is the only unknown parameter for the reinforcement.  is the tangent modulus of elasticity of the reinforcement. Accounting for the reinforcement, equation (40) is expressed as follows:

                                               (44)

The global system of linear equations for the entire domain is formulated as follows:

                                   (45)

where  is the global tangent stiffness matrix for the entire region,  is the global nodal displacement increment vector,  is global external load vector,  is the total number of finite elements,  is the number of reinforcement-yielding elements.

At each step of the iterative process, the finite element  with the maximum principal stress  was determined. If  a crack initiates within the element. Subsequently, the vector  was updated using the following formula:

                                                                  (46)

Reducing the step size is permissible because the tangent stiffness matrix is constant. Thus, at each iteration, crack initiation is restricted to a single finite element exhibiting the maximum principal stress. This algorithm allows one to track the sequence of crack formation. Upon crack initiation, additional unknown nodal displacements are introduced into the finite element structure [24], and a new iteration is executed.

If no cracks initiation occurs within any finite element, the iteration termination condition is checked. The Euclidean norm was used to calculate the solution residual, and the iteration termination criterion was the defined by the following inequality:

                                                                 (47)

where  is an acceptable computational tolerance. If the vector  approaches zero during the iteration, the nodal equilibrium conditions of the finite element mesh nodes are satisfied, and the iterative process terminates. Fig. 6 shows the step-by-step algorithm of iterations.

1. , , , , .

2. For  = 1 to number of load step do

3.    .

4.    . Crack=false.

5.     For  = 1 to number of iterations do

6.       Compute for all finite elements:

            (9),

           .

           If (j > 1) .

           If (j > 1) and (Crack=false) and ( )

             then

                , .

                 Terminate the current load step.

7.       Compute: .

8.       Solve the system: .

9.      ,

          Compute for all finite elements:

          , ,

         , ,

          principal stresses (23),

10.     Find out , .

11.     If  

            then

                 Crack=true.

                  Add the new unknowns to the nodes of finite element .

                  Compute for all finite elements:

                  , ,

                  , ,

                   ,

                  , ,

                   principal stresses  (23).

              else

                   Crack=false.

           Compute the finite element’s equivalent deformations

          .

 12.   End of iterations loop.

        Form the load step results.

 13. End of load step loop.

Figure 7. Step-by-step algorithm of iterations.

To validate the proposed model, a numerical analysis of a reinforced concrete beam subjected to a uniformly distributed load was performed. The beam parameters are listed in Table 1. The reinforcement consists of four 24-mm diameter bars Fig. 3.

Table 1. The beam parameters.

Length, m

h, m

b, m

fcm, MPa

fctm, MPa

Eci, MPa

Es, MPa

Rs, MPa

As, m2

ys, m from the bottom

6

0.6

0.3

22

1.805

32500

200000

235

0.00181

0.03

Due to symmetry, the numerical analysis was performed for only half of the beam. To model the support, vertical displacements were restricted at the nodes of the leftmost cross-section. The beam was analyzed under a uniformly distributed vertical load  = 10, 20, 30, 34 kN/m.

Fig. 8 shows the evolution of the solution residual (47) during the iterations. The loading step size was set to  = 2 kN/m. The spikes in the residual norm correspond to the initiation of new cracks during the loading process. Once no new crack initiation occurs during a load step, the residual norm approaches zero within just two iterations (as indicated by the nearly vertical lines in Fig. 8). Thus, the number of required iterations is determined chiefly by the total number of initiated cracks.

Figure 8. Convergence history of the solution residual. qz = 24 kN/m, grid 20 by 80.

Figs. 9–12 show the numerical results for the beam analysis under various loads. The scale of the crack width visualization is 50 times larger than that of the finite element mesh. The maximum crack width is 0.128 mm under a load of  = 34 kN/m. Because the rightmost crack (numbered 1 in Fig. 12) lies directly on the symmetry plane, its calculated width was doubled to reflect the actual physical crack opening.

Figure 9. qz = 10 kN/m. Grid 20 by 80.

Figure 10. qz = 20 kN/m. Grid 20 by 80.

Figure 11. qz = 30 kN/m. Grid 20 by 80.

Figure 12. qz = 34 kN/m. Grid 20 by 80.

The concrete stresses in the X-direction are presented in Fig. 13. Fig. 14 presents the mid-span load-deflection curves obtained using the proposed method in comparison with the LIRA-SAPR simulation results. In LIRA-SAPR, an exponential stress-strain diagram for class B40 concrete was adopted. According to LIRA-SAPR calculations, the maximum tensile stress in the reinforcement is 107182 kN/m2, and the peak compressive stress in the concrete along the X-direction is –10529 kN/m2. These findings are in good agreement with the numerical predictions obtained via the proposed method. Maximum vertical displacements obtained using LIRA-SAPR calculations are fifteen percent smaller than those obtained using the proposed method (Fig. 14).

Figure 13. The stresses along the X-axis. LIRA-SAPR.

Figure 14. The blue line represents the results of the proposed method; the red line shows the LIRA-SAPR results.

It should be noted that the proposed method does not take into account possible slippage of the reinforcement relative to the concrete, which can lead to increased crack width and displacements. To account for reinforcement bond slip, independent degrees of freedom must be defined for the rebar nodes. The reinforcement-to-concrete interface can be captured via shear spring elements, governed by a prescribed bond stress-slip diagram.

4.Conclusion

A finite element model for the non-linear analysis of reinforced concrete beams incorporating crack initiation and propagation has been presented. Plane stress rectangular finite elements were employed to discretize the concrete beam. The finite element stiffness matrix was formulated based on a complementary energy functional. By employing Lagrange multipliers, the functional was augmented with nodal equilibrium equations derived via the principle of virtual displacements. An incremental orthotropic model was applied for the stress-strain relationship of concrete.

Crack formation was incorporated by introducing additional unknown nodal displacements and using independent stress approximation functions within the element subdomains separated by the crack. The stiffness matrix of the rectangular finite element was derived by approximating the stress fields within the element domain using linear functions. To solve the non-linear problem, an incremental-iterative algorithm was proposed to track the progressive sequence of crack initiation within the finite elements. The Euclidean force norm was used to evaluate the solution residual. The number of required iterations is determined chiefly by the total number of initiated cracks. Once no new crack initiation occurs during a load step, the residual norm approaches zero within just two iterations.

To validate the proposed model, a numerical analysis of a reinforced concrete beam subjected to a uniformly distributed load was performed. The numerical results derived via the proposed method were validated against the LIRA-SAPR simulation data.