Cracking analysis of reinforced concrete structures

The study of a reinforced concrete beam tested under four-point bending is proposed in this article. An analysis of the behaviour of this beam from the global response down to local information such as cracking is performed. In order to describe the progressive degradation of the beam, a damage model is used, associated to a stress-based non-local regularisation method. A post-treatment of the finite element analysis is then performed to characterise the cracking pattern (crack spacing and crack opening). In this article, two different methods of post-treatment are compared: the topological search and continuous/discontinuous crack opening and the global/local analysis. Results show the capability of both methods to give a good estimation of cracking. Finally, the main advantages and drawbacks of both methods are underlined.


Introduction
The durability analysis of reinforced concrete structures requires the quantification of cracking.To deal with this problematic, several approaches are available.
On one hand, the description of cracking can be explicitly introduced in the modelling, as shown in X-FEM (Moës, Dolbow, & Belytschko, 1999), G-FEM (Strouboulis, Babuška, & Copps, 2000) or E-FEM (Oliver, Huespe, & Sánchez, 2006) for example.In these approaches, the location of the crack and the crack opening are directly quantified.However, the modelling of crack initiation is still under discussion.
On the other hand, the initial micro-cracking process and the macro-cracking failure can be modelled in a continuous framework using regularised damage models.The description of cracking is then performed by means of a post-processing procedure.In this paper, the continuous approach is explored with two particular post-treatment methods used to quantify cracking in reinforced concrete elements.Both methods are compared on a case study of a reinforced concrete beam subjected to a four-point bending test.
In the first part of this contribution, the mechanical models used to describe the degradation of the structure in a continuous framework are presented, including recent developments regarding the description of the strain localisation.Then, two original methods designed to extract information on cracking are described.Finally, these methods are compared within the framework of the modelling of a reinforced concrete beam subjected to a four-point bending test.

Description of degradation in concrete
The micro-cracking process observed in concrete is described at the macro-scale in a continuous framework with a non-local damage model.The main details of the model are described in this part.

Damage model
The progressive degradation of concrete is characterised by a scalar damage variable D. This internal variable links the stress tensor r to the strain tensor e, following Equation (1).r ¼ð1 À DÞC:e (1) where C represents the elasticity tensor.This damage scalar variable D ranges from 0 (virgin material) to 1 (failed material).According to the model proposed by Mazars (1986), this variable is driven by an internal variable κ that corresponds to the maximum value reached by the equivalent strain e eq during the loading and is initially equal to e D 0 (strain at first crack in tension).The formula of the equivalent strain e eq is given in Equation ( 2).
where e i are the principal strain values and hÁi þ are the Macaulay brackets.Equation (2) expresses that damage is driven by strain extensions.The damage loading function f D (ε eq , κ) associated to the damage variable D corresponds to the maximum value reached by the equivalent strain e eq during the loading and is initially equal to e D 0 (strain at first crack in tension).The damage loading function f D ðe eq ; jÞ follows the Kuhn-Tucker conditions (Equation (3)).
e eq À j 0; _ j !0; _ jðe eq À jÞ¼0 In order to account for the non-symmetrical behaviour of concrete under uniaxial loading, the damage variable D is expressed as a linear combination of two variables (Equation ( 4)): D t for the degradation in tension and D c for the degradation in compression.
where a c and a t represent the recombination factors for tri-axial loading and b is a parameter accounting for shear behaviour.Equation ( 5) gives the evolution of the variables D t and D c .
This model offers a robust and accurate representation of concrete behaviour and is therefore largely used for industrial studies.Unfortunately, like every model presenting a softening behaviour, it suffers from localisation process during failure.By describing this process with a constitutive model expressed in a local formulation, it leads to non-objective results in numerical applications.In order to address this problem and to introduce an explicit description of the fracture process zone (FPZ) (i.e.introduction of characteristics length of the FPZ in the modelling), several regularisation methods have been proposed (gradient enhanced media, non-local model on internal variable, etc).

Non-local model
A non-local regularisation method on the internal variable is used in this study in order to maintain the objectivity of the results.This model, originally proposed by Pijaudier-Cabot and Bažant (1987), replaces a local variable X by its non-local counterpart X following Equation ( 6).
X ðxÞ¼ where /ðx; sÞ is a weighting function, x is the position of the point where the non-local variable is calculated and s is the position of a point in its neighbourhood.Classically, the weight function can be chosen as the Gaussian function Equation (7).
Different choices can be made regarding the variable to regularise.In the framework of isotropic damage model as the one presented previously, the internal variable j is a classical choice (Saouridis & Mazars, 1992) that does not lead to locking problem (Jirásek, 1998).
This regularisation method allows the objectivity of the results.Nevertheless, it fails to properly describe both the strain field and the damage profile at complete failure and cracking initiation close to boundaries.In order to improve the description of the continuous fields at failure and thus the estimation of cracking by post-treatment of these fields, an evolution of non-local interactions should be introduced during computation.Among the recent developments relative to this pathology inherent to the original nonlocal model, one can quote the works of Desmorat and Gatuingt (2007), Pijaudier-Cabot and Dufour (2010) and Gregoire, Rojas-Solano, and Pijaudier-Cabot (2013).However, in the actual contribution, the method proposed by Giry, Dufour, and Mazars (2011)i s used.The influence of the stress state is introduced in the description of the non-local interactions through a stress factor ρ defined in Equation (8).qðx; rðsÞÞ ¼ 1 where h is the angle between u 1 and the projection of ðx À sÞ over the plane defined by u 1 and u 2 and u is the angle between u 2 and ðx À sÞ.u 1 , u 2 and u 3 are the eigenvectors of the stress tensor r and r 1 , r 2 and r 3 are the associated eigenvalues.The internal length of the non-local model is thus modified following Equation (9).
This stress-based non-local model allows us to describe the progressive decrease of non-local interactions across the FPZ during failure.As a consequence, Giry et al. (2011) have shown that the strain field across the damaged zone is much closer to a strong discontinuity than the original version by Pijaudier-Cabot and Bažant (1987).Furthermore, the interactions close to free boundaries are better described by avoiding the attraction phenomena.
3. From continuous calculation to discrete information: a characterisation of cracking In this part, the two post-treatment methods to extract cracking from the continuous calculation are briefly presented.The basic ideas to locate and quantify the cracksthrough their openingsare given.More details can be found in Bottoni and Dufour (2010), Dufour, Pijaudier-Cabot, Choinska, and Huerta (2008) and Dufour, Giry, and Mazars (2013) for the topological search and continuous/discontinuous crack opening (CDCO) and in Oliver-Leblond, Delaplace, Ragueneau, and Richard (2013) for the global/local analysis.

Location of the crack path
The method to locate the crack path is based on a topological search approach.From a 2D scalar field Y representative of micro-cracking, the locus of the maximum values corresponds to crack path.For the continuous model considered here, j is considered as the representative field of cracking Y.
To perform the crack path search, three main steps need to be defined: the initiation step, the current step and the ending step.The numerical parameters introduced in the procedure are the searching step a and the length of regularization l smooth used in a Gaussian function / introduced for the convolution product.

Initiation step
The Gauss point P ini bearing the maximum value of the post-treated field is used to initiate the searching procedure.As no direction of search is initially defined, a circle centred on P ini with a radius equal to a is used to project the field of interest.The maximum value of the convoluted profile Y ðsÞÃ/ðs; l smooth Þ, where Ã is the convolution product, defines the point P 2 of the crack.The first point P ini is re-evaluated, since its position depends on the location of Gauss points, by taking the maximum of the convoluted profile along the profile which is perpendicular to P 2 P ini ! at P ini .Then, current steps can be performed.As two directions are defined by points P ini and P 0 , each is considered in turn.

Current step
From a point P i and a searching direction P iÀ1 P i !, the point P iþ1 is defined as the locus of the maximum value of the convoluted profile Y ðsÞÃ/ðs; l smooth Þ along the line which is perpendicular to the search direction at point P 0 defined as It is worth noting that the crack path does not depend on the space discretisation of the initial simulation since the crack path is found out from the convoluted field.

Ending step
Different criteria are considered to stop the searching procedure: point P 0 is out of boundaries or, the maximum of the convoluted profile used to identify P iþ1 is smaller than e D 0 (strain at first crack in tension).
As this procedure can only be performed for 2D fields, an intermediate step must be performed by projecting the field of interest on several cut planes along the depth of the considered structure.

Crack opening
In order to estimate crack opening, a post-treatment of the strain field obtained from the FE calculation is performed.After identifying the crack path, the crack opening is calculated by comparing the computed strain field e FE and the analytical strain profile e FE ¼½U dðxÞ corresponding to a strong discontinuity (displacement jump) along a line which is perpendicular to the crack path.As the analytical strain profile is a Dirac-like function, a direct comparison cannot be performed with the numerical strain profile.Thus, a convolution product is applied on both profiles using a Gaussian function /.I n order to get an evaluation of crack opening, one should make a hypothesis between the analytical and the finite element strain profiles in order to compare them.In the following work two different hypotheses have been considered: the "strong approach" and the "weak approach".The "strong approach" considered an equality of the convoluted profiles at the crack location x 0 (Dufour et al., 2008), whereas the "weak approach" considered an equality in average for the convoluted profiles (Dufour, Pijaudier-Cabot, & Legrain, 2009).One can then obtain an estimation of the crack opening ½U (see Equation (10.a) for the "strong approach" and (10.b) for the "weak approach"). (10.a) This method also gives an error indicator Dðx 0 ; ½U Þ along the crack path for each perpendicular profile (Equation ( 11)).It quantifies the error made by comparing the numerical profile to a strong discontinuity.
Equation ( 11) gives only a model error and should be seen by the user as an indicator of the quality of the hypothesis made by comparing the finite element strain profile and the strain profile of a strong discontinuity.
3.2.Location and quantification of crack: global/local analysis 3.2.1.Process of analysis After the continuous computation of the damage pattern at the global scale, a non-intrusive reanalysis at the local scale can be performed with a discrete model in order to extract fine information about cracking.The successive steps of this strategy are summarised in Figure 1.
The obtained damage pattern is studied and several regions of interest (ROIs) are defined corresponding to distinct areas of damage.Then, the loading steps for the extraction of crack features are determined and will correspond to the steps of reanalysis.For those loading steps, boundary conditions are extracted from the continuous displacement field and applied on the non-free surfaces of the ROIsnamely the ones which cut the whole domain.The natural way to transfer the displacement field from the global scale to the local scale is to use the shape functions of the finite elements used for the global computation.Then, the displacement u L ðx L Þ applied at a local node x L of a non-free surface of a ROI is directly obtained with Equation ( 12).
where N j are the shape functions of the finite element model, u j is the displacement vector computed at the global scale and N FE is the number of finite element nodes.
The cracking pattern obtained at the previous step will be retrieved at the current step of the local reanalysis in order to follow the crack propagation accurately.The discrete computation of the chosen ROIs can be parallelised.

The local model
The discrete model used at the local scale has been proposed by Delaplace and Desmorat (2007) and offers a reliable description of concrete behaviour for tensile loadings.
The material is described as an assembly of polyhedral particles linked by Euler-Bernoulli beams.The quasi-brittle behaviour of the material is obtained through a brittle behaviour for the beams.The breaking threshold P ab of a − b, the beam linking the particles a and b, not only depends on the beam extension e ab but also on the rotations of the two particles θ a and θ b (see Equation ( 13)).
The six parameters of the beam a − b need to be calibrated.First, the length and the area are imposed by the geometry.Then, the inertia and the elastic modulus are identified in order to retrieve the elastic behaviour of the global computation.calibration of the breaking thresholds of the beam e cr and h cr allows us to fit the peak and post-peak behaviour of the global model.The calibration is performed on an independent case study (Oliver-Leblond et al., 2013).
Our study focuses on a fine description of crack pattern and on the measurement of the crack opening.The crack pattern is defined as the common side of the particles initially linked by the breaking beams.The opening of the crack is computed by considering the relative displacement u b À u a of the unlinked particles a and b.The measure of the opening between those particles e ab is projected on the normal n ab of the local discontinuity (Equation ( 14)).

Analysis of a reinforced concrete structure
The structure studied in this section is a reinforced concrete beam subjected to a fourpoint bending test.This test has been proposed as a numerical experiment within the French project CEOS.fr.

Description of the test
A reinforced concrete beam -5.1 m long, 0.8 m high and 0.47 m deephas been loaded in a four-point bending test.Three layers of three rebars of 32 mm in diameter are placed at the bottom of the beam (area under tension) and three rebars of 25 mm in diameter are placed at the top of the beam (area under compression).The vertical concrete cover is equal to 50 mm and the horizontal rebar spacing is equal to 35 mm.The two loading points are 1.6 m away from each other.
The various symmetries of the structure and assuming a symmetrical crack pattern allow us to only study a quarter of the beam.A displacement control has been considered for the loading (Figure 2).
The concrete and the steel reinforcement used for this beam have the following characteristics: E cm = 40.8GPa, f cm = 64 MPa, f tm = 4.4 MPa, E s = 200 GPa and f yk = 500 MPa.The model parameters considered for the computation of the continuous analysis (Section 2) are summarised in Table 1.
In order to reduce the computational time, the non-linear behaviour of concrete is only considered for the central part of the beam (1.7 m long) which includes the loading points.Cubic elements with linear interpolation functions have been used to describe concrete and lower rebars.Bar elements have been used to describe upper rebars.The average mesh size is equal to 15 mm in the centre part of the beam and 40 mm outside.The mesh used for the computation and the loading scheme is given in Figure 2. 4.2.Behaviour of the structureforce/displacement analysis One can observe in Figure 3 a classical behaviour of a reinforced concrete structure: the progressive degradation of the flexural modulus with the appearance of the different cracks in the specimen.

Behaviour at the material scalecracking analysis
One can observe during the degradation of the beam, the development of three main areas of cracking (Figure 4).The left crack is initiated at first on the lower rebars under the loading plate due to the perturbation of the strain field above it.A damage evolution along the lower rebars is then observed until the initiation of the second cracks.
In Figure 4, the damage field at failure stays relatively sharp and no spurious diffusion of damage is observed at crack levels and close to rebars.This precise definition of the different damaged areas is obtained thanks to the stress-based non-local model that leads to a progressive decrease of non-local interactions in the crack vicinity.
The cracking profilescrack path and openingobtained with the cut planes method (Figure 5   good agreement with each other.The appearance of a secondary crack in the neighbourhood of the central crack is also caught by both methods.Table 2 summarises the averaged location from the loading plates of the three main cracks for both methods on the surface of the beam.It shows again a good agreement in the determination by both methods of the crack location. These results are also in good agreement with the location of the maximum stress state observed along the lower rebars (Figure 6).The small discrepancies observed for the left and the central cracks are due to the change in the crack path at the level of the rebars.Indeed, when only the upper part of the macro-crack (i.e. the area above the rebars) is considered, the same crack location is obtained for both methods.Figure 7 gives a comparison between both methods for the quantification of crack opening at two different loading steps for a crack at the surface of the beam (average stress in the lower rebars in the central part of the beam: 150 and 219 MPa).
One can observe a good agreement between both methods at the bottom of the beam.At the location of the rebars, a small discrepancy is seen.This difference can be explained by the damage area observed.Indeed, the rebars tends to distribute the damage in their vicinity.At the surface, the influence is still felt for this specimen, leading to a small level damaged area inside which we find a higher damaged more localised area corresponding to the macro-crack.The "strong approach" of the CDCO method supposed an equality of the convoluted profile at the crack location whereas the "weak approach" considered an equality of the profile in average.These two hypotheses can be seen as limit conditions and the estimation of the crack opening of a macro-crack in a micro-cracked field is probably a combination of these two approaches.For the GLA method, the area of broken elements corresponds to the micro-cracked area and the macro-crack.This last is identified among the broken elements as the path with the main crack opening.It is interesting to note that the crack opening obtained with the "strong" and the "weak" approaches of the CDCO method surround the values obtained with the GLA method.At failure, the strain field is fully localised, leading to a decrease of the difference observed between both methods.Furthermore, the difference observed close to rebars can be also explained by the presence of a partial mode II failure that is not caught by the GLA method.In Figure 7, the difference at the level of the rebars between both methods tends to increase.Over this area, the FPZ is no longer influenced by the rebars and a good agreement is recovered between both methods.
The CDCO method identified the crack path from the last step of the computation.One can see in Figure 7 (Average stress in the rebars: 219 MPa) a good agreement between both methods in the identification of the crack tip.For the previous step, the CDCO method tends to give a crack tip location higher than for the GLA method as the diffused micro-cracks area observed in front of the FPZ tends to give a small crack opening along the crack path identified from the last step of the calculation.For the GLA method, the identification of the crack path has been made at each computed step of the discrete reanalysis using graph theory combined to Bellman-Ford algorithm.
In Figure 8, one can see the model error for the CDCO method.At the front of the FPZ, a high error is observed due to the fact that the finite element strain profile is still diffused compared to analytical Dirac-like profile.At the bottom of the beam, the macro-crack is well established and the model error is lower.The error limit value observed is relative to the mesh discretisation considered.Indeed, the mesh used for the finite element computation is coarse as a consequence one gets at failure a finite element strain profile closer to a step function than a Dirac-like function.Nevertheless, this error should be seen as an error indicator for the model only and not be directly taken as a difference with a global error indicator as the different methods used tends to give similar values for crack opening at the bottom of the beam.

Conclusion
Continuum damage mechanics, largely developed these last two decades, succeed in predicting the strongly non-linear behaviour of reinforced concrete structures.The steel rebars employed for large scale structures lead to the occurrence of a large number of micro-and macro-cracks.This phenomenon is the key point to describe so as to account for the structural elements ductility.Some interesting mechanisms may be handled within this framework such as: Young's modulus degradation, permanent strain, crack closure, induced anisotropy, strain rate effects and frictional sliding (Richard & Ragueneau, 2013).Unfortunately, due to the post-peak softening branch, the uniqueness properties of the mechanical solution are lost and regularised models have to be introduced.Moreover, if one is concerned with more precise features of cracks like spacing, openings or tortuosity, the macroscopic model such as damage mechanics fails.This contribution aimed to give some insights regarding these two drawbacks.An original non-local procedure has been presented allowing the introduction of the anisotropic interaction in the FPZ along the crack tip.The energy dissipation when the crack propagates is more physically motivated.Regarding crack opening determination, the use of damage mechanics at the structural scale is kept, for robustness evidence and engineering point of view.Some post-treatments are then needed.This paper exhibited two different approaches for post-treating a Finite Element Analysis performed at the structural scale.
On the one hand, a topological search method is used to locate cracks, analysing a damage strain field.The crack openings are computed by comparing the regularised strain field and the analytical strong discontinuity displacement field induced by a discrete crack.On the other one hand, the global/local method exposed in this paper introduces a discrete mechanical model at the finer scale to obtain crack properties of a Region of Interest using the finite element displacement field as Dirichlet boundary conditions.
Using the same macroscopic damage model, the two approaches have been compared regarding the crack opening prediction for a reinforced concrete beam subjected to a four-point bending test.The results in terms of crack location and opening at several loading steps are very close to each other.This means that regarding cracks refined features, the use of continuum damage mechanics at the higher scale makes sense, since two completely different post-treatment methods give the same cracks pattern and crack opening profiles.
The main advantage of the first method lies in the fact that a unique mechanical model is used and identified.One major difficulty remains in the crack opening computation, implying the post-treatment of only mode I propagations.The use of cut plans only gives discrete values of openings instead of a continuous line in 3D.Concerning the second method, two mechanical models (one finite, the other discrete) have to be handled and identified generating more difficulties.The main advantage in introducing the discrete element methods is the kinematic of such elements: the cracks are naturally described in 3D.Some important phenomena such as tortuosity can be accounted for in case of mechanicaldiffusive coupled problems.

Figure 2 .
Figure 2. Mechanical scheme and mesh of a quarter of the reinforced concrete beam.
(a)) and with the global/local analysis (Figure 5(b)) are shown to be in

ForceFigure 3 .
Figure3.Evolution of the force applied on a loading plate vs. the mid-span deflection.

Figure 4 .
Figure 4. Damage field for the first macro-crack and at the end of the computation.

Figure 5 .
Figure 5. Crack field at the last step of the computation: (a) topological search and CDCO estimation and (b) global/local method.

Figure 6 .
Figure 6.Evolution during the loading of the stress profile along the bottom rebars.

Figure 7 .
Figure 7. Evolution of crack opening along the height for a crack (at right in Figure 7) at the surface of the beam for two average stress levels in the lower rebars: 150 and 219 MPa during the loading of the stress profile along the bottom rebars.

Figure 8 .
Figure 8. Evolution of the model error for both approaches of the CDCO method along the height for a crack (at right in Figure 7) at the surface of the beam for two average stress levels in the lower rebars: 150 MPa (a) and 219 MPa (b) during the loading of the stress profile along the bottom rebars.

Table 1 .
Model parameters for concrete and steel.