Size effect in concrete beams under bending – influence of the boundary layer and the numerical description of cracks

In the paper the size effect phenomenon in concrete is analysed. The results of numerical simulations of using FEM on geometrically similar un-notched and notched concrete beams under bending are presented. Concrete beams of four different sizes and five different notch heights under three-point bending test were simulated. In total 18 beams were analysed. Two approaches were used to describe cracks in concrete. First, eXtended Finite Element Method (XFEM) describing cracks as discrete cohesive ones with bilinear softening was chosen. Alternatively, an elasto-plastic constitutive law with Rankine criterion, associated flow rule and bilinear softening was defined. In order to ensure mesh-independent FE results, a non-local theory in an integral format as a regularisation technique was applied in the softening regime. In both approaches the influence of the decrease of the material parameters (mainly fracture energy) in the boundary layer on obtained maximum loads was studied. Additionally the influence of the averaging method in non-local plasticity was also examined. Obtained results were compared with experimental outcomes available in literature.


Introduction
The behaviour of heterogeneous quasi-brittle materials, like concrete, is very complex. In these materials a softening regime with decreasing force with increasing displacement is observed. The response of the material is govern by fracture process zones, where localization of deformation takes place, in the form of cracks (mode I) or shear zones (mode II). As a consequence, a strong size effect (dependence of strength and other mechanical properties on the size of the specimen) is observed. For real concrete structures it is between two extremes: the plastic limit load theory and linear elastic fracture mechanics. The understanding and a proper description of the size effect is very important in extrapolating experimental results obtained on small specimens into bigger ones.
Within continuum mechanics, there exist two main approaches to describe cracks. The first one describes them in a smeared sense as localized zone of microcracks with a non-zero finite width. The material can be described using e.g. elasto-plastic, damage mechanics or coupled constitutive laws. However, classical FEsimulations with material with softening are not able to model localisation properly. The obtained results suffer from the mesh sensitivity and produce unreliable results. The strains concentrate in one element wide. Constitutive laws have to be equipped with a characteristic length of microstructure to preserve the well-posedness of the boundary value problem. It can be done by means of e.g. a micro-polar, non-local or gradient theories. In addition, a size effect can be properly captured. The alternative approach introduces displacement jumps across cracks while keeping the remaining region as a continuous one. The oldest solutions used interface elements defined along element edges. The modern ones allow for considering cracks in the interior of finite elements using embedded discontinuities or XFEM (eXtended Finite Element Method) based on a concept of the partition of unity. A smeared approach is more proper when describing a micro-crack formation phase while a discontinuous one allows for a better simulation of a macro-crack propagation. Usually, only one approach is used to simulate a crack growth in concrete during the entire deformation process.
Recently, some experiment campaigns were performed to investigate fracture and size effect in concrete beams under bending. Hoover and et al. [1] examined 18 un-notched and notched specimens with 4 different sizes and 5 notch heights made from one batch of concrete. Similar experiments were performed by Çağlar and Şener [2]. Grégoire at al. [3] studied the same problem by testing 12 beams with 4 different sizes and 3 notch heights. Next these experiments were numerically simulated with cohesive elements [4] or non-local damage models [3,5].
Paper presents results of numerical simulations of experiments [1] using two different constitutive laws. In both approaches the influence of the decrease of the material parameters (mainly fracture energy) in the boundary layer on obtained maximum loads was studied. Additionally the influence of the averaging method in non-local plasticity was also examined.

Elasto-plasticity
The elasto-plasticity model was defined within a standard plasticity theory. Standard Rankine criterion was used [6][7]. The yield function for 2D case was defined as: where 1 and  2 are the principal stresses, t is the tensile yield stress and  stands for the hardening/ softening parameter (equal to the maximum principal plastic strain). An associated flow rule was assumed. In order to define the tensile yield strength t, a bilinear softening curve was assumed: where f is the value of the herdening parameter when yield stress t reduces to 0 and stress k and hardening parameter define k the knee point at which the postpeak softening slope desreases. The similar curve was assumed by Hoover and Bažant [4]. The bilinear softening yield curve is shown at Fig. 1. Strain softening causes the ill-posedness of the boundary value problem and FE-results become mesh dependent. In order to achieve mesh-independent results, a non-local theory in an integral format was used as a regularization technique. Rates of the softening parameter d were treated non-locally according to Brinkgreve [8]: where quantity d is defined as: In above equations x are the coordinates of considered (actual) point,  are the coordinates of the surrounding points, m is the non-locality parameter (it should be greater than 1) and 0 is the weighting function. Here the Gauss distribution was chosen: where r is the distance between two points and l defines characteristic length of the microstructure. It should be noted that averaging operation affects only to a small area around each material point (the influence of points at the distance of r=3l is only of 0.01%, see Fig. 2). In FE-simulations, an approximated method was used to evaluate non-local quantities. In the given integration point, the influence of its neighbors was determined using the values from the previous iteration. A split method based on the proposal by Strömberg and Ristinmaa [9], called the Newton split, was used (superscripts (i) and (i-1) stand for iteration counters): The terms with the (i-1) iteration counter were taken from the previous global iteration and they were frozen during local iteration at the material point during the global actual iteration (i) while the terms with the superscript (i) were active and they could change their values. It allowed to simplify FE calculations and to preserve the locality of plasticity algorithms.

Extended Finite Element Method
The Extended Finite Element Method (XFEM) was used to describe discrete cracks in continuum. It is based on a local partition of the unity (PUM) concept. The method assumes that a displacement field is discontinuous across the crack just after its initiation. It enables adding 'ad hoc' extra terms to a standard FE displacement field interpolation. These extra functions are responsible for capturing displacement jumps. Cracks do not have to be placed along finite element edges; they may pass through elements. The formulation used in the paper follows (with some improvements and modifications [10][11]) the original concept by Wells and Sluys [12]. In order to describe jumps in the displacement field, the so-called shifted-basis enrichment proposed by Zi and Belytschko [13] was used. The shifted-basis enrichment simplifies the implementation of XFEM (two types of the finite elements exist only). Moreover the total nodal displacements are equal to the standard ones.
Two constitutive laws were defined to describe the behaviour of a solid body with a cohesive crack. In a bulk (un-cracked) continuum, a linear elastic relationship between stresses and strains was assumed.
During active loading the softening of the normal component of the traction vector tn was described by the yield curve n using an bilinear softening law, similarly as Eqn.
(2) in elasto-plasticity. For XFEM this formula was redefined with the aid of initial and total fracture energies Gf and GF (see Fig. 3 for explanation): In order to improve convergence a tension-compression transition of the traction tn the yield stress y was multiplied by a factor Df defined as [14]: where df is a drop factor. With increasing df, the value of Df approaches 1 and the standard bilinear softening law is recovered. In a compressive regime, the penalty stiffness in a normal direction was defined as A new crack could be activated or an existing crack could propagate, if the standard Rankine criterion was fulfilled at least in one integration point at the element at the front of the crack tip: A crack could develop only in a vertical direction (due to the nature of the problem analysed, see next Section).

Problem data
In the experiment by Hoover et all [1] concrete beams of similar geometry and different sizes under three point bending were tested. Figure 4 presents the geometry of a beam. The depths (heights) H were taken as 40 mm, 93 mm, 215 mm and 500 mm. The length-to-depth ratio was equal to 2.4 in all cases. The span length L was taken as 1.176 of the beam height H. The width of all beams was equal to 40 mm. The relative notch depth 0

Elasto-plasticity
First numerical calculations with smeared cracks using elasto-plasticity with non-local softening were executed.
The characteristic length was equal to l=5 mm and the non-locality parameter was taken as m=2. The failure softening parameter was f=7.52·10 -3 . The kink point was defined at 15% of the tensile strength and 22.67% of the f value. All parameters were calibrated to reflect experimental ad numerical data from [4], especially fracture energies (calibrated in 1D problem). Figure 5 presents force-displacement diagrams obtained for standard isotropic averaging scheme. Numerical results overestimated both peak loads and resultant fracture energies (no agreement in the softening regime). Figure 6 presents comparison of peak stresses (calculated as a maximum moment divided by a section modulus without taking into account the notch) obtained from numerical calculations versus averaged experimental values. The average relative error of the peak stress was about 32%.  In order to improve the results, an approach proposed by Giry et al. [15] was used, in which the influence of the stress state in the neighbours is taken into account. As a consequence the averaging operator is not isotropic any longer. Obtained results are presented at Figs. 7 and 8. It can be seen that maximum forces were slightly smaller, but they still overestimated experimental outcomes. The average relative error of the peak stress was about 25%. This method, however, significantly improved the behaviour of the model with to the respect of properly reproducing the softening branch of the force-CMOD curve.
Finally, simulations with Giry et al [15] approach were repeated with adding the influence of the boundary layer along the specimen's edges. In this band the fracture energies linearly decrease with decreasing the distance to outer edges till the specified fraction of the original values. The width of the boundary layer was Fig. 7. Force-CMOD diagrams for all beams from simulations with elasto-plasticity and Giry et el [15] averaging scheme. Fig. 8. Comparison of experimental and numerical maximum peak stresses for elasto-plasticity and Giry et al [15] averaging. k assumed as 1 cm with the reduction fraction equal to 0.5. The average relative error of the peak stress was again smaller, about 18%, but no significant improvement was observed.

Extended Finite Element Method
Next numerical simulations with the discrete cracks using XFEM were performed. To describe the behaviour of the cracks in the softening phase, the initial fracture energy was taken as Gf=49.56 N/m and the total fracture energy was assumed as GF=70 N/m. The kink stress  was defined as 15% of the tensile strength (all values after [4]). The drop factor was equal to df=10 4 . Figure 9 presents obtained force -CMOD curves for all geometries compared to experimental bounds. Figure  10 presents comparison of peak stresses obtained from numerical calculations versus averaged experimental values. In general, curves from numerical simulations are in good agreement with experimental ones during the whole loading (Fig. 9).
The analysis of peak stress values shows some discrepancies. While for the huge beams (H=500 mm, the lowest set of points) a very good agreement was obtained, the same set of parameters was not able to capture experimental values for the small beams (H=40 mm, highest set of points). The relative error of the peak stress was between -1.8% and 19%.
In order to improve the results, simulations with boundary layer with reduced fracture energies were performed. The width of the boundary layer was taken as 0.5 cm and 1.0 cm with the corresponding maximum reduction factors at the specimen's boundary equal to 0.2 and 0.5, respectively. Figures 11 and 12 present maximum peak stresses obtained for the boundary layer equal to 0.5 cm and 1.0 cm, respectively. It can bee seen that for both cases good agreement between numerical and experimental results was obtained. The relative error of the peak stress was between -5.6% and 4.8% for the boundary of 0.5 cm and between -5.6% and 6.3% for the boundary of 1.0 cm. Fig. 11. Comparison of experimental and numerical maximum peak stresses for XFEM with the boundary layer 0.5 cm.

Conclusions
In the paper the simulations of size effect of concrete beams under bending were presented. Both smear and discrete cracks' descriptions were used. Calculations have shown that direct application of the experimental data into elasto-plastic constitutive law led to overestimating maximum experimental forces. A partial solution was the use of non-isotropic averaging schemes. This problem requires further investigations. On the contrary, simulations with XFEM showed very good agreement with experimental outcomes. The activation of the boundary layer with reduces fracture energies resulted in reducing the relative error of the peak stresses.