Incorporation of pre-existing cracks in finite element analyses of reinforced concrete beams without transverse reinforcement

Cracking in reinforced concrete (RC) bridges and other structures is common and not necessarily detrimental. However, some cracks may grow past specified limits and, aside from aesthetic and durability aspects, may influence the ductility and structural capacity of an RC member. This is not generally reflected in current assessment methods and, therefore, improved methods are needed. The aim of the current work was to develop a modelling methodology to incorporate pre-existing cracks into finite (FE) analysis for improved structural as- sessments. Two different approaches were investigated: (1) weakening the continuum elements at the position of a crack and (2) introducing discrete crack elements with weakened properties. In both approaches, a total-strain based model was used in the continuum elements. These modelling approaches were applied to analyses of experiments, in which concrete beams had been pre-cracked and tested in four-point bending. The pre-existing cracks led to differing failures limiting the deformation capacity, plus varying ultimate capacity and ductility. In the current study the weakened-elements approach captured the failure characteristics, ultimate capacity and ductility more accurately than a standard FE analysis without cracks included; the discrete-crack approach, on the other hand, did not. Furthermore, the bending stiffness differed between the experimental tests and the FE analyses. Damaged bond properties and closure of cracks in the compressive zone were identified as probable causes. Moreover, the choice of shear retention used for the weakened elements was shown to noticeably affect the predicted capacity and ductility. In conclusion, the weakened-elements approach was the most straightforward to implement. It was less time-consuming and led to better agreement with experimental results, compared to the discrete-crack approach. Based on this study, the weakened-elements approach is regarded as a promising approach for the structural assessments of tomorrow.


Introduction
Countries in Europe and worldwide rely on their road networks for economic and social development. Investment in road networks is vast, with bridges as the most expensive element. Although bridges only contribute some 2% of the length of the road network, bridge infrastructure represents 30% of total investment [1]. Many bridges are composed entirely of structural concrete, or designed to function in combination with other materials; steel girders supporting a concrete deck, for example [2].
Cracking in reinforced concrete (RC) bridges and other structures is common but not necessarily detrimental. Many structures are designed to crack under service loads. However, cracking needs to be controlled if it is to meet durability and aesthetic requirements. Current codes specify crack width limitations and minimum reinforcement levels, see [3] for example. The causes of cracking in concrete members are numerous; from plastic settlement or plastic shrinkage cracks in young concrete, to cracks from live loads acting on the structure [4]. Cracks may also occur in concrete structures due to internal or external restraints. For example, temperature variations between different parts of a structure can lead to varying deformation needs. If those deformations are prevented, cracking may occur [5,6]. In practice though, some cracks may grow past specified limits even if design codes are followed. More importantly, not only does cracking negatively affect durability and aesthetics, it may also influence the ductility and structural capacity of an RC member [7]. It has been experimentally shown that the influence on the ductility and ultimate capacity can be both positive and negative, depending on the position and characteristics of the crack as well as the failure mode [8,9]. Cracks should be incorporated in structural assessment to consider their impact on ductility, capacity and failure mode and there is, therefore, a need for improved assessment of reinforced concrete infrastructure.
One way to improve assessment is by including updated information on the structure; information on pre-existing cracks, for example. This may be facilitated by Digital Twin (DT) models [10] which, in the future, are anticipated to play a critical role in optimised management of critical infrastructures [11]. A DT is a virtual replica of the structure. It stores information collected during the structure's lifetime (through different types of sensors, for example). It may provide insights into the structural capacity via finite element (FE) analysis. It may also function as a decision-making support tool [11]. Given the latest advancements in inspection techniques (such as inspection drones and fibreoptic sensing, see [12,13;14,15] respectively), data on cracks will be collected more accurately and more frequently than ever before. Consequently, the FE simulations used in structural assessments need to incorporate information on pre-existing cracks.
There are several ways to model cracking in nonlinear FE analyses of reinforced concrete, with smeared and discrete concepts being the most common [16]. The smeared crack approach treats a cracked solid as a continuum, while the discrete-crack approach treats a crack as a geometrical discontinuity. An objection to the discrete-crack approach is the uncertainty concerning the location at which cracks form. However, in DT models the pre-existing crack patterns at the surface are readily available, at least at inspectable surfaces, which partly invalidates this objection. Smeared crack models commonly use a total-strain based approach [17], with either rotating or fixed crack direction. Both have been shown to give reasonable results [18].
Nevertheless, neither the discrete nor the smeared crack approach are sufficiently well developed to describe pre-existing cracks in concrete structures. Discrete cracks are commonly used when assessing dams for the incorporation of critical pre-existing cracks [19]. For such large-scale structures with high self-weight, a simple frictional law without including aggregate interlock in the shear transfer model may be generally acceptable, while it is not for other types of structures. For example, aggregate interlock is known to be an important shear transfer mechanism in beams of low shear reinforcement ratio [20]. Discrete cracks have also been used to represent through-cracks in reinforced concrete slabs tested by bending for fatigue assessment [21]. However, the treatment of aggregate interlock was not generally applicable but tailored to represent the fatigue response of the slab.
Moreover, only a few works were found in the literature concerning the incorporation of pre-existing cracks into the smeared crack approach. In [22], pre-existing cracks in concrete columns were recreated using a special load arrangement which was applied before the loading of interest. This approach represented cracking reasonably well in the studied case but required complicated implementation and was deemed computationally expensive for the treatment of multiple cracks. In the methodology proposed in [23], the constitutive law for the cracked finite elements was derived from the uncracked state, using a scalar damage parameter to reduce stiffness and strength. Variation in damage among cracks of different widths was discussed but not implemented in the analyses, which assumed a similar damage parameter for all cracks. It may be noted that the methodology proposed in [23] shares similarities to a proposal by the authors [24]. However, the derivation of damage parameters is different.
The current work aimed to develop a modelling methodology to incorporate pre-existing cracks in FE analysis and obtain improved structural assessments. The variation in damage depending on the individual crack width and shear retention (aggregate interlock) was explicitly included, to address weak points identified in the literature. The development was based on previous work by the authors [24], where the capability to represent the effect of corrosion-induced splitting cracks on the confinement from the concrete cover in assessment of anchorage capacity was investigated. The proposed methodology includes methods for both smeared and discrete cracks and does not require the load history to be known as, in practice, this is expensive or impossible to map. Rather, it is based on metrics (crack pattern and widths) that are available upon inspection. Furthermore, direct incorporation of pre-existing cracks is more computationally efficient, compared to modelling the full load-history. The proposed methodology allows for capacity estimations of structures with pre-existing cracks and is suitable for use in, say, DT models.

Methods for pre-existing crack incorporation
Two different approaches for considering pre-existing cracks in FE analysis were investigated. These were: a) weakening the continuum elements at the position of a crack and b) introducing discrete crack elements with weakened properties. Both approaches involved choices as to tensile properties and are detailed below.

Weakened-element approach
In the weakened-element approach, finite elements coinciding with pre-existing cracks are assigned weakened tensile properties, as compared to the intact concrete. These weakened tensile properties were determined from an assumed bilinear, mode-I, stress-to-crack width relationship for the undamaged concrete, with the kink point as per [25]: The left-hand graph shows the stress-crack opening relationship and the residual fracture energy. The centre graph shows the relationship converted into a stress-strain relationship, using the crack bandwidth h. The right-hand graph shows the resulting stress-strain relationship, to be used for the weakened elements.
Based on the measured crack width, the tensile strength and remaining crack opening until stress-free cracking was derived from the bilinear relationship, as shown in Fig. 1. Using the measured crack width, denoted by w c in the figure, the tensile stress (f ct,c ) and residual fracture energy (G F,Wc ) were determined from the bilinear, mode-I, stress-to-crack width relationship for sound concrete. The stress-strain relationship corresponding to the stress-to-crack width relationship was derived using the modulus of elasticity (E cm ) and assumed crack bandwidth (h). In determining the crack bandwidth for the weakened elements, localisation within one element row was assumed. This led to a bandwidth of 10 mm. To construct the stress-strain relationship used as input for the weakened elements (to the right in Fig. 1), only the ultimate strain (ε wc,ult ) needed to be determined. This was chosen so that the area under the stress-strain relationship equalled the residual fracture energy, divided by the crack bandwidth (G F,wc /h). This corresponds to: A Poisson's ratio of zero was used for the weakened elements. Furthermore, the cracks with measured widths larger than w ult were assigned low tensile properties, by using a crack width of 0.99w ult in the calculations.
The fixed-crack approach was used for weakened elements. In contrast to the rotating-crack approach, where no shear stresses are transferred in the crack, the fixed-crack approach uses a so-called shearretention factor to specify the remainder of the initial shear modulus. The choice of shear-retention factor is not trivial and may be influenced by several factors, such as the type of aggregates used. A shear-retention factor which decreases with increased normal strain is generally recommended [26] based on, say, the decay of normal stiffness with decay of tensile stress after (damage-based) cracking. This was used for the sound concrete. However, a fixed shear-retention factor of 0.01 was used for the weakened concrete elements. This was done to enable comparison between analyses of weakened elements and discrete cracks. These used the commercial FE software (DIANA 10.3 [27]), which does not allow for variable shear retention.

Discrete-crack approach
The tensile properties of the discrete cracks were derived in a fashion similar to that of the weakened-element approach, but with some slight modifications. The initial elastic properties of the discrete cracks were determined as E cm /t, where t was the assumed thickness of the interface layer (taken as 0.1 mm). The crack stress (f ct,c ) and residual fracture energy (G F,Wc ) were determined in similar fashion to the weakened elements. The stress-crack opening relationship was inputted in the FE software as a linear relationship, based on the stress and residual fracture energy. Thus, the crack opening, corresponding to a stress-free crack, was generally determined as: For pre-existing crack widths larger than the kink point w s, this can be simplified to: As with the weakened elements, the tensile properties for cracks with measured widths larger than w ult were assigned low values by using a crack width of 0.99w ult in the calculations.
The initial shear modulus was determined as: where t was the interface thickness (assumed as 0.1 mm) and ν was the Poisson's ratio for the concrete (taken as zero). It should be noted that the crack widths started at zero in the FE analysis, even though the width of the physical crack was non-zero before mechanical loading. Due to the aforementioned limitations of the FE software, a constant shearretention factor of 0.01 was used. This reduced the shear modulus of the discrete crack elements upon numerical cracking.

Application to cracked beams tested in bending
The experimental campaign is presented first, followed by an account of the FE implementation, as per the approaches described in Section 2.

Summary of experimental campaign
The experimental campaign included three dog-bone-shaped reinforced concrete specimens, see Fig. 2. These were pre-cracked before structural testing and are, therefore, denoted PC. An additional three reference specimens, without their bone-ends (in other words, only the 2 m beam segment), were cast from the same concrete batch. The three specimens within each set (PC and reference) received similar treatment.
The experimental campaign consisted of three main stages: 1) restrained shrinkage, 2) tensile loading and 3) testing in four-point bending. These are described below and augmented with relevant information on the FE implementation, such as material properties, crack information and four-point bending test procedure.

Restrained shrinkage
Each concrete beam segment in the dog-bones was prevented from shortening by two steel struts, placed on each side of the specimens. These struts were connected to the concrete beam through the end blocks, see Fig. 2. As the concrete drying and autogenous shrinkage progressed, the much stiffer steel beams prevented any shortening. Thus, cracks formed in the concrete beam when the stresses reached the tensile strength of the concrete. The degree of initial restraint was estimated at 0.6. This was done using the area of the concrete cross-section, anticipated tensile strength and modulus of elasticity, plus the area of the steel struts and their modulus of elasticity.
The beam segments of the dog-bone specimens and reference specimens had four ϕ10 reinforcement bars, one in each corner with 30 mm of concrete cover in both directions. The longitudinal reinforcement bars extended into the end regions. The ends were also reinforced by three horizontal layers of ϕ12 closed-loop bars and six ϕ10 vertical closedloop stirrups, see Fig. 2.
The dog-bone and reference specimens, plus material test specimens were cast using concrete, mixed as per the proportions specified in Table 1.
After casting, the specimens encased in plywood formwork were wrapped in plastic sheets to reduce water evaporation until de-moulding at seven days. During the 162-day restrained shrinkage period, the relative humidity was lowered to around 30% at an average temperature of 30 • C. During this phase, the top sides of the specimens were monitored using a DIC system, which took pictures every 90 mins.

Material characterisation
Specimens for material characterisation were cast from the same batch as the dog-bone and reference specimens. The compressive and tensile strength, plus fracture energy and modulus of elasticity of the concrete were tested. Based on the fracture energy and the tensile strength and using Eq. (1), w ult was determined as 0.204 mm. The yield and ultimate strength, plus the modulus of elasticity of the reinforcement were also characterised. The concrete and reinforcement material properties at the time of the bending tests (184 days after casting), plus specimen types and test methods are presented in Table 2. It should be noted that the modulus of elasticity of the concrete was tested at 28 days and extrapolated as per EC2 [28]. Furthermore, the characterizations of compressive strength, tensile strength and tensile fracture energy were all based on samples of three specimens, while the characterizations of the modulus of elasticity and the reinforcement properties were performed on samples of five and six specimens, respectively.

Tensile loading
The crack widths due to concrete shrinkage were small; about 0.08 mm when still under restraint. A tensile loading procedure was executed after the shrinkage phase, to increase the crack widths and investigate how they were influenced by an external load. Short segments of the steel struts were cut and manually operated hydraulic jacks inserted, see Fig. 3. Linear variable differential transformers (LVDTs) were placed over some of the cracks to measure the increase in crack width upon force exertion. The force was incrementally increased to around 70 kN per hydraulic jack, a procedure that lasted 1-2 min. This resulted in approximately 140 kN carried through the concrete beam. Thereafter, the load was removed by releasing the pressure and the same procedure was repeated a second time.
The tensile loading increased the crack widths formed in the restrained shrinkage stage, as well as inducing more cracks. The widths of six cracks per beam were monitored using LVDTs under tensile loading. Their average crack width increased from 0.09 mm before loading, to 0.31 mm at maximum load and then diminishing to 0.20 mm after loading. All cracks, not just those monitored by LVDTs, were measured before and after tensile loading. The procedure increased the average crack width, from 0.08 mm under external restraint to 0.15 mm without external restraint.

Crack patterns
The crack patterns resulting from the restrained shrinkage and tensile loading are presented in Fig. 4. The crack patterns are presented in top view, with sides 1 and 2 folded out from the top surface. The blue lines on the specimens' sides represent restraint cracks, labelled with the prefix w r . Red lines represent cracks forming under tensile loading,  *Estimated based on 28-day value. labelled with the prefix w. The top side shows a contour plot of the tensile strains at the end of the restrained shrinkage phase (cracks formed under tensile loading are not included). Note that some noise is visible as a scattered pattern in the contour plots. The width of the cracks shown after tensile loading and unloading in Fig. 4 are presented in Table 3. These widths were measured using a handheld digital microscope at the level of the top reinforcement or, for cracks of widely varying widths, at the point of maximum width. For cracks formed in the restrained shrinkage phase, widths measured before removal of the restraint are also shown in parentheses.

Bending test
After the restrained shrinkage and tensile loading phases, the beams (the middle 2 m segment in Fig. 2) were cut from the dog-bone specimens and tested in four-point bending. The test setup, plus moment and shear diagrams are shown in Fig. 5. The test setup was designed to make crack influence likely. This was done by subjecting the beam to a bending moment/shear ratio that was close to the moment/shear capacity ratio. The two point-loads P 1 and P 2 had different magnitudes, with the ratio P 1 /P 2 chosen as equal to 1.2. This promoted shear failure on one side. Moreover, the point-loads were placed asymmetrically to produce a constant bending-moment region between them. During the test, a DIC system (GOM Aramis DIC 12M) was set up to monitor the region between mid-span and right-hand support (the side with the highest shear force). An LVDT was also placed at mid-span and midheight, with a load cell measuring the total load (P 1 plus P 2 ). The beams were loaded through a load distribution beam, to which a hydraulic jack applied vertical deformation at a rate of 1 mm/min. When tested, the beams were orientated with "side 1" in Fig. 4 as the face seen in Fig. 5.

Implementation of cracks in FE analyses
The FE analyses were conducted using the DIANA 10.3 commercial software [27]. The three-dimensional (3D) geometry was modelled using continuum elements, with an average element size of 10 mm. The beam was placed on support plates, restrained along the centrelines against translation in the vertical and transverse directions, see Fig. 6. One node in the left-hand support was also longitudinally restrained to prevent rigid body movement.
The load was applied to the beam through load plates resting on thin pieces of wood fibreboard. A stiff distribution beam was used to apply displacement-controlled loads of differing magnitude (P 1 /P 2 = 1.2). The ends of the load distribution beam were tied to the centrelines of the load plates, enforcing the same vertical displacement. Moment equilibrium was used to determine the point along the load distribution beam where vertical displacement would be applied. Since the ratio of 1.2 was between P 1 and P 2 , this point was positioned at 387.3 mm from the left-hand support, as indicated in Fig. 6. The element mesh consisted mostly of eight-node brick elements (HX24L, 90% by count) for the concrete. However, five-node pyramids (PY15L, 3%), four-node tetrahedrons (TE12L, 2%) and six-node wedge elements (TP18L, 5%) were also used by the FE software's meshing algorithm. For analysis of the reference case (without pre-existing cracks), strain localisation was promoted by reducing the fracture energy to G F /10 for the light grey elements in Fig. 6.
Interface elements were placed between the fibreboard and concrete and between the support plates and concrete. These were assigned nonlinear elastic properties, with initial normal stiffness calculated from E cm (divided by an assumed interface thickness of 0.1 mm) as 3.565•10 13 N/m 3 . The shear stiffness was assumed to be a small portion (10 − 8 ) of the normal stiffness, resulting in 3.565•10 5 N/m 3 . No tension could be transferred in the interface and the shear stiffness was reduced to zero, at an interface opening of 1.0•10 − 3 mm. This was done to allow for some movement between load and support plates and the concrete and to avoid large stress concentrations. The reinforcement bars were modelled with beam elements. Von Mises plasticity was used, and the stress-strain behaviour was specified using a bilinear relationship based on the properties given in Table 2. The interaction between reinforcement bars and concrete described using a bond stress-slip relationship. The bond stress-slip relationship was specified as per fib Model Code 2010 [33] but with the residual as per [34].
The deformation was imposed in 200 steps of 0.05 mm, followed by 150 steps of 0.1 mm. A line search (as per [27]) was enabled, to promote convergence. The equilibrium iterations were conducted using the Secant (Quasi-Newton) BFGS iteration method. The maximum number of iterations was set to 200, and convergence was considered as fulfilling either an energy norm of 0.0001 or an unbalanced force norm of 0.01.
The sound concrete was modelled using a total-strain-based smeared crack approach, with fixed crack orientation. Damage-based shear retention was used for the sound concrete elements; in other words, the shear stiffness was reduced with increasing normal strain after cracking. The compressive behaviour was modelled as parabolic [35], with strength reduction due to lateral cracking as per [36] and the lateral confinement was incorporated as per [37]. The tensile behaviour was modelled as per Hordijk [38], with a crack bandwidth manually specified to 2 corresponding to the length of two element diagonals. This was verified against the size of the localisation zones in the analyses. Moreover, the Poisson's ratio was reduced based on the cracking (damage) [27].

Weakened-elements approach
The tensile stress-strain relationships for the weakened elements were determined as per the procedure presented in Section 2. The crack widths of each crack are specified in Table 3. A crack bandwidth of 10 mm was used. This corresponds to strain localisation in one element row and was later confirmed for the weakened elements in the analyses. Note that the crack bandwidth differed between the sound concrete and the weakened elements. This since the lower tensile properties caused the strain to localize in the weakened element row while for the sound concrete elements the localization took place over two elements (diagonally). The modelling of the compressive behaviour for the weakened elements excluded strength reduction due to lateral cracking, but was otherwise similar to that of sound concrete.
The surface crack patterns (shown in Fig. 4) were extruded in the width direction of the beam towards the centreline, as a simplified representation of the cracks inside the beam. Side 1 was extruded 55 mm and side 2 was extruded 65 mm. This was done to overlap the cracks by one element at the mid-section, so that cracks on the two sides could join even if one was shifted one row in the longitudinal direction. The elements coinciding with the extruded crack were assigned weakened element properties, using an automated selection procedure based on manually mapped crack patterns. Furthermore, additional elements were manually weakened to allow the path of weakened elements to open without cracking the sound concrete elements. These were selected so that the weakened elements along a crack path were always connected side-to-side, and not only along an edge. Fig. 7 shows the mesh, including weakened (cracked) elements in red for side 1 (diagram a) and the top side of specimen PC1 (diagram b). The manually selected elements in a) and the overlapping elements in b) are indicated in orange.

Discrete-crack approach
In the discrete-crack approach, discrete crack planes were placed at the position of the pre-existing cracks. When surface cracks were present on both sides of the specimen, a through-crack plane was created manually using the geometrical information from the mapped crack pattern. A discrete crack plane was defined by using points on three corners of the beam belonging to the cracks. Two points were chosen on the side with the largest crack width, a third point was chosen on the opposing side and the fourth point was calculated. The width of the discrete crack (assumed to be constant over the entire crack) was taken as the average of the crack width, measured on sides 1 and 2 of the specimen. Cracks only seen on one side were assumed to extend to half the width. The FE implementation of the discrete crack planes was achieved using plane 4 + 4-node quadrilateral interface elements (Q24IF). Furthermore, nodal lumping was used for the discrete crack elements, as suggested by [39]. The discrete crack elements for PC2 appear in Fig. 8.

Table 3
Widths (measured after the tensile tests) of cracks depicted in Fig. 4. The values in parentheses correspond to the width of restrained shrinkage cracks before removal of external restraint. Measurements in mm.

Comment on the implementation methods
The pre-processing of FE models for the weakened-elements approach was straightforward. Using the standard FE mesh as a basis, the properties of certain elements were changed if their position corresponded to that of a crack. The implementation of discrete cracks was more cumbersome as interface elements had to be added to the mesh. In the DIANA 10.3 software [27], this was done by cutting the sections up in pieces prior to meshing and assigning properties to the interfaces between the various parts. The geometry was then meshed into continuum and interface elements. Simple geometry for the discrete cracks was chosen because more complex geometries were deemed too timeconsuming (and had caused numerical problems in early trials). However, it should be noted that using another mesher (whereby interface elements can be placed directly into the mesh) might make preprocessing more straightforward.
Only visual information was used to determine the crack patterns of concrete specimens. This meant that only the external crack pattern was known for the specimens. It is possible that radiographic or electromagnetic methods [40,41], or a combination of both, could provide some information about the internal crack pattern. However, the   equipment necessary for these methods was not available for the experimental campaign and such techniques might prove difficult to apply to real large-scale structures. Thus, it was necessary to make an assumption regarding the specimen's internal crack pattern. As stated earlier, for the weakened elements this assumption meant extruding the surface crack pattern to the centreline of the beam. Meanwhile, the discrete cracks were assumed to penetrate to full or half width of the specimen, depending on the crack pattern on the opposing sides. For the discrete cracks, use of an approach similar to that of weakened elements was investigated but with linearisation of the non-smooth crack paths. However, this led to numerical problems due to ill-shaped elements forming near the beam's centreline.

Results
The load-deflection curves resulting from the four-point bending tests are shown in Fig. 9. A difference in stiffness between the reference and pre-cracked beams may be observed. The reference beams had high initial stiffness before the bending cracks developed, after which the stiffness was reduced. For all reference specimens, the bending moment firstly caused yielding and some hardening of the reinforcement, followed by shear failure. The pre-cracked specimens (PC1-PC3) did not show high initial stiffness. Rather, the load increased in fairly linear fashion until the reinforcement yielded. The stiffness was similar for specimens PC2 and PC3, while the stiffness of PC1 was slightly lower. There was a large difference in ductility among the pre-cracked specimens; for PC1 a shear crack (left-hand shear span) limited the ultimate midspan deflection to around 23 mm. For PC2 crushing between load plates occurred at over 40 mm of deflection. However, for PC3 a shear crack developed in the right-hand shear span at less than 10 mm of midspan deflection.
The results of the FE analyses from the four-point bending test are presented below, alongside the experimental results.

Reference specimens
The tests revealed similar behaviour in all three reference beams tested in four-point bending. After an initial stiff response in the uncracked stage, the beams were able to bear increased load at lower stiffness until yielding. Thereafter, hardening of the reinforcement steel led to a slight increase in load before a shear crack formed in the righthand shear span prevented further deflection. However, there was a difference in ultimate deflection between the three test specimens.
The results in terms of load-deflection and crack pattern are shown in Fig. 10. The FE analysis predicted the ultimate load and ductility well and the stiffness was captured adequately. The failure limiting the deformation capacity obtained in the FE analyses matched the experiments; namely, shear failure in the right-hand shear span after yielding and some hardening of the reinforcement. It is noted that the inclination

Table 4
Crack pattern and failure limiting the deformation capacity in experiments and analyses. Yielding of the reinforcement bars occurred in all cases. Experiments: crack pattern from DIC (red indicates a strain of 0.01 and green 0) and observed failure limiting the deformation capacity. Note that DIC results were unavailable for PC1. Weakened elements: pre-existing cracks implemented as weakened elements, crack pattern at ultimate load (red indicates a strain of 0.01 and blue 0) and predicted failure limiting the deformation capacity. Discrete cracks: position and widths of discrete cracks (red indicates 0.5 mm and blue 0), with crack pattern in continuum elements at ultimate load (red indicates a strain of 0.01 and blue 0) and predicted failurelimiting the deformation capacity mode. Note that the figures of the experimental results only cover half the specimens while the full specimens are shown for the FE results.

Pre-cracked specimens
The load-deflection curves for specimen PC1-PC3 (obtained from experiments and FE analyses) are shown in Fig. 11. Results from the weakened-elements approach are shown in the top row of plots and the discrete-crack approach in the bottom row. The FE analyses of the precracked beams using weakened elements correctly captured the lack of an initial stiff response (corresponding to an uncracked beam), while the discrete-crack approach showed high initial stiffness. Nevertheless, the stiffness in the FE analyses differed when compared to the experimental results. This was observed in both weakened elements and discrete cracks. A discussion of plausible reasons is provided in Section 5.
For specimen PC1, the analyses with weakened elements and discrete cracks failed to capture the reinforcement's yielding and hardening. A shear crack formed on the right-hand side, just as the reinforcement started to yield. Hence, the ultimate deflection was underestimated, as was the ultimate capacity in the weakened-elements approach. The yielding and hardening of specimen PC2 was captured by both approaches, although the ultimate load and deflection at failure were underestimated. The failure in the concrete compression zone occurred in both analyses. The ultimate load and ductility of specimen PC3 was well captured using the weakened-elements approach, while the discrete-crack approach overestimated both metrics. Moreover, the failure limiting the deformation capacity differed between the approaches; the weakened-elements analysis failed in shear on the righthand side, whereas the discrete-cracks analysis failed in compression at the right-hand load plate.
An overview of the experimental crack patterns and failures plus data from the FE analyses are shown for specimens PC1 to PC3 in Table 4. As yielding of the reinforcement occurred in all cases, the failure is referring to the factor limiting the deformation capacity. The table's top row shows a difference between the failure characteristics and crack patterns in the experimental specimens.
For the weakened-elements approach, shear cracks developed in both shear spans. However, the failure ultimately occurred on one side or in the compression zone. It should be noted that the number of cracks at ultimate load was roughly the same between the analyses of weakened elements and discrete cracks. The number of cracks in the latter case includes also cracks formed in continuum elements. The finite element contour plots shown in Table 4 need to be read with care as they include existing and new cracks for the weakened element approach but only new cracks for discrete crack approach.

Discussion
The four-point bending tests were conducted on specimens cast from the same concrete and subjected to the same load-history under the same environmental conditions. Nevertheless, the resulting capacities and failure modes varied between the test specimens. A quantitative comparison of the modelling approaches is given below. The cause of the observed difference in bending stiffness between the experimental results and FE analyses is then investigated. There then follows a sensitivity study of the shear retention used in the FE analyses.

Comparison of modelling approaches
The main results (in terms of ultimate capacity and maximum deflection) of all analyses are compiled in Table 5. Better correlation between FE analysis and experimental results was found for analyses with pre-existing cracks included by weakened elements compared to the discrete crack approach and Standard FEA (cracks omitted) as they had the lowest mean (or average) absolute error, both for ultimate capacity and maximum deflection. Implementation of pre-existing cracks with the discrete approach gave less accurate prediction than the standard FEA. The load at yielding of the reinforcement was similarly compared. All analysis types were able to estimate this load to a margin of error of around 1%.
The prediction of the failure limiting the deformation capacity was also compared, see Table 6. The incorporation of pre-existing cracks led to improved predictions of the failure, as compared to overlooking the cracks. Analyses of weakened elements or discrete cracks correctly predicted the failure for two of the three specimens, while the analyses without pre-existing cracks were correct in one case. It is worth noting that the greatest improvement from adding pre-existing cracks by weakened elements was a more accurate estimate of the maximum Table 5 Ultimate capacities and maximum deflections for specimen PC1-PC3, analysed without pre-existing cracks and with cracks included as weakened elements and discrete cracks. The difference of the FE modelling compared to the experiments was also calculated, as well as the absolute error. The values in parentheses are the relative absolute errors.  deflection and failure limiting the deformation capacity. However, the ultimate capacity was only slightly improved. The computational time for the FE analyses with discrete cracks was around 15% larger compared to analyses with weakened elements. This was based on the average time-consumptions of around 8 and 7 h respectively. Moreover, as described in Section 3.2.3 the discrete cracks were also more time-consuming to implement in the analysis. In conclusion, the weakened-elements approach was most straightforward to implement and less time-consuming. Compared to the discrete-crack approach, it also led to results that were in closer agreement with experimental results for the studied beams. It can be noted that the analyses were stopped when convergence was lost according to the criteria specified in Section 3.2. These criteria constitutes a part of the solution strategy and there is no consensus on the tolerances [26]. It is possible that the use of other convergence criteria could influence the results in terms of estimated ductility and ultimate capacity.

Investigation of difference in stiffness between experimental results and FE analyses
As stated above, a difference in stiffness was observed between the FE analyses and the experimental results for the specimens with pre-existing cracks. In other words, the slope of the load-deflection curves differed (as observed in Fig. 11). Pre-cracked specimens PC1 to PC3 were loaded in tension prior to the four-point bending test (as presented in Section 3). The maximum achieved jacking force of 70 kN meant a stress of around 450 MPa in the reinforcement bars at cracked sections. At such stress levels, it is plausible that the tensile loading damaged the bond between the reinforcement bars and surrounding concrete. Four different bond stress-slip relationships were therefore considered, to investigate whether the difference in stiffness between the experimental results and FE analyses could be explained by a weaker steel-concrete bond. The bond stress-slip relationships considered are shown in Fig. 12. The same relationship was used over the entire length of the reinforcement bars. The first of these, denoted "Original", corresponds to the bond stress-slip relationship given in Model Code 2010 but with the residual as specified in [34]. The weakened bond stress-slip relationships are denoted "Low 1" to "Low 3", with Low 3's initial stiffness and peak stress being the most weakened. The initial stiffness of Low 1, Low 2 and Low 3 stiffnesses were chosen as 8%, 5% and 2% of the initial stiffness of the original relationship. This to qualitatively study a range of possible damage degrees to the steel-concrete bond.
The load-deflection curves obtained from the (weakened-elements approach) analyses with weakened bond relationships are shown in Fig. 13. The original bond relationship gives the stiffest load-deflection response, while the load-deflection behaviour becomes less stiff in the case of weakened bond properties. This is expected since a weaker bond relationship requires larger slip to mobilize the same bond stress compared to a stiffer one; this influences the tension stiffening. Moreover, the influence of the bond stress-slip relationship on ultimate capacity seems relatively small, except for the Low 3 relationship, which leads to premature anchorage failure. It should be noted that no signs of anchorage failure (such as major bar slippage) were noted in the experiments.
The compressive behaviour of the concrete material used for the cracks was identified as another potential factor having a major effect on initial bending stiffness. The FE analysis therefore investigated the compressive properties of the weakened elements. The comparison was made based on the notion that, in the physical test, the cracks in the compressive zone needed to close before any significant compressive force could be transmitted. Thus, the compressive stress-strain curve was adjusted by introducing a lower stiffness of 0.11 GPa, for strains of between 0 and 0.015. Meanwhile, for greater strains, the stiffness was increased to its original value of E cm at 35.7 GPa (see left-hand plot in Fig. 14). Setting the strain value at which the stiffness changed to 0.015 Fig. 12. Bond stress-slip relationships used for sensitivity study of bond stiffness and strength. Fig. 13. Load-deflection curves for PC1 to PC3. These were obtained using four different bond stress-slip relationships, with varying stiffness and peak stress. allowed an average crack width of 0.15 mm to be implied, based on the crack bandwidth of 10 mm. This value corresponds to the average crack width measured. Fig. 14 shows the load-deflection curves for the FE analyses, with original and adjusted compressive stress-strain curves. The bending stiffness was markedly lower for the case in which the adjusted stress-strain relationship was used and shear failure (right-hand side for both analyses) occurred at a lower load level.
It may be concluded that both the bond stress-slip relationship and compressive material model for the pre-existing cracks influence the bending stiffness of the beam. Thus, a combination of both is considered the likely cause of the disparity in stiffness between the FE analyses and experimental results.

On the shear retention of cracked concrete
The shear-retention factor β may be interpreted as the aggregate interlock consideration in the FE analysis [42], when a fixed crack model is used. This factor may be difficult to quantify and may depend on several factors, such as the types of aggregates and size of member analysed. This work used a variable shear-retention factor for the sound concrete elements in all analyses, with the shear modulus decay similar to the normal stiffness decay after cracking. However, for the weakened elements and discrete cracks, the shear retention was constant. This facilitated comparison between both approaches, as explained in Section 2. To investigate the model sensitivity in regard to the shear-retention factor used for pre-existing cracks, a comparison was made of the results obtained using different values. This investigation included shearretention values of 0.1, 0.01, 0.001, 0.0001 and 0 for both the weakened elements and the discrete-crack approach.
Moreover, an individual shear-retention factor for each crack (based on the measured crack width) was also used to investigate both weakened elements and discrete cracks. The weakened-element approach included an additional case in which the individual shear-retention factors were reduced during the analysis, based on the normal crack strain. These individual shear-retention factors were calculated based on the assumption that shear stiffness gradually reduces after cracking and reaches zero at a crack width of half the aggregate size [26]. The shearretention function was derived based on a reduction of the secant stiffnesses in a bi-linear relationship; similar to the stress-strain relationship in Fig. 1 but with the x-axis scaled to reach zero at a crack width of 2.5 mm (corresponding to 50% of the approximate average aggregate size). The resulting shear-retention function (with logarithmic y-scale) appears in Fig. 15, alongside examples of a constant shear-retention factor and decreasing input shear-retention function, both based on a crack width of 0.1 mm. Note that the decreasing shear-retention function is a function of the crack strain. This is calculated by subtracting the measured crack width and dividing by the crack bandwidth (10 mm), see top axis of Fig. 15.
The load-deflection curves for specimens PC1 to PC3 are shown in Fig. 16. These were analysed using the weakened-element approach and with the shear behaviour described above. A similar plot for the discretecrack approach is shown in Fig. 17.
For the weakened-element approach, the three lowest constant values for the shear-retention factors led to similar ultimate capacity and ductility, while the two higher values led to higher capacities. Furthermore, the analyses of individually derived shear-retention factors showed similar results to the analyses of reducing shear-retention functions. In turn, both were similar to analyses conducted with a constant β of 0.001. However, the exception was PC2, for which the individually derived shear-retention factors reached a higher capacity and ductility. Regarding the failure limiting the deformation capacity, the  highest shear-retention factor (of 0.1) led to a compression-type failure in all specimens. This failure also occurred for a β of 0.01 for PC2, with shear-type failure occurring in all other cases.
For the discrete-crack approach, the results were similar for all shearretention factor values greater than zero, albeit with some variation in ductility. PC1 failed in shear and PC2 failed in compression, for all shear retention choices investigated. PC3 failed in compression for all but zero shear retention, for which shear failure occurred. It should be mentioned that since the shear-retention factor specifies the portion of the initial shear modulus being retained, the observed indifference among shearretention factors may be caused by the particular choice of initial shear modulus (see Section 2.2). It may not extend to other initial shear modulus choices.
To conclude, the choice of shear retention was shown to be influential in the analyses involving weakened elements; in these, a high shear retention led to a change of failure mode, plus greater capacity and ductility compared to lower values. In most cases, using individual, constant or decreasing shear retention led to underestimation of the experimental results. The discrete-crack approach was less sensitive to the choice of shear retention. However, this observation may depend on the choice of initial shear stiffness.

Conclusions
This paper has presented a methodology for including pre-existing cracks in FE modelling, with the aim of enhancing structural assessments. Two different approaches were investigated; namely, weakening elements at the crack positions and placing discrete crack elements at the crack positions. The methodology was validated with experiments in which reinforced concrete beams were pre-cracked using restrained shrinkage and tensile loading and thereafter tested in four-point bending. Based on the study, the following main conclusions may be drawn: • influence of pre-existing cracks on failure limiting the deformation capacity, ultimate capacity and ductility was demonstrated; Fig. 16. Load deflection for PC1-PC3, analysed using the weakened-elements approach, with general constant shear retention plus individually assigned constant and decreasing shear retention. Fig. 17. Load deflection for PC1-PC3, analysed using the discrete-crack approach, with general constant shear retention plus individually assigned constant shear retention.
• incorporating pre-existing cracks in weakened elements led to improved estimates compared to traditional FEA, particularly of the ductility and failure characteristics in the present study. However, the ultimate capacity was also slightly improved; • traditional FEA provided better estimates of the ultimate capacity and ductility compared to analyses using discrete cracks in the present study. However, these gave improved predictions of the failure limiting the deformation capacity; • the analyses using weakened elements required less implementation effort and less computational time compared to those with discrete cracks; • modification of the compressive behaviour of the weakened elements (to reflect closure of cracks in the compressive zone), plus reduced reinforcement bond stiffness and strength (to represent damage from previous tensile loading) was shown to influence the bending stiffness in the analyses towards better correlation with the experimental observations; • the choice of shear retention for the weakened elements was shown to influence the analyses results. However, its influence was minor in the analyses using discrete cracks.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.