On Some Problems in Determining Tensile Parameters of Concrete Model from Size Effect Tests

Abstract The paper presents results of numerical simulations of size effect phenomenon in concrete specimens. The behaviour of in-plane geometrically similar notched and unnotched beams under three-point bending is investigated. In total 18 beams are analysed. Concrete beams of four different sizes and five different notch to depth ratios are simulated. Two methods are applied to describe cracks. First, an elasto-plastic constitutive law with a Rankine criterion and an associated flow rule is defined. In order to obtain mesh independent results, an integral non-local theory is used as a regularisation method in the softening regime. Alternatively, cracks are described in a discrete way within Extended Finite Element Method (XFEM). Two softening relationships in the softening regime are studied: a bilinear and an exponential curve. Obtained numerical results are compared with experimental outcomes recently reported in literature. Calculated maximum forces (nominal strengths) are quantitatively verified against experimental values, but the force – displacement curves are also examined. It is shown that both approaches give results consistent with experiments. Moreover, both softening curves with different initial fracture energies can produce similar force-displacement curves.


INTRODUCTION
Concrete is one of the most popular materials used in civil engineering to erect buildings, bridges and other structures on land and various port facilities like breakwaters or quays. Its tensile strength is several times smaller than compressive one, therefore it is usually applied with steel reinforcing bars. These bars are responsible to carry (mostly) tensile stresses, while the concrete itself sustains compressive stresses, and the tensile strength of concrete is neglected. The situation is different in plain or weekly reinforced massive concrete structures in civil or coastal engineering, like foundations, dams etc. The tensile strength of concrete cannot be neglected and the proper estimation of the concrete's properties in tension is of great importance.
One of the salient phenomenon observed in concrete structures is the presence of size effect. The strength (and other properties) of material depends on the size of a specimen examined; small elements have a greater nominal strength than bigger ones. It is caused by the existence of the so-called fracture process zone (FPZ), which size is not negligible comparing to the size of the specimen. Deformations localise inside these zones; at the beginning as a set of diffused microcracks and later a discrete macro-crack. The observed size effect results are placed between plastic limit theory and linear elastic fracture mechanics solutions. Therefore, the proper theoretical capture of the size effect phenomenon is crucial when experimental results from tests performed on small specimens are to be extrapolated into bigger ones.
Within continuum mechanics, cracks can be described as a smeared (continuous), discrete (discontinuous) or by coupling these two approaches (Mazars et al. [20], Marzec et al. [18], Unger et al. [26], Bobiński and Tejchman [5]). The first (smeared) approach defines a crack as region (band) with a certain width, while in the second formulation it is presented as a line (surface) with zero thickness and assumed displacement jump across. When using classical constitutive laws with strain softening, a mesh dependency is observed. In order to obtain mesh independent results and to restore the well-posedness of the boundary value problem, a regularisation method is needed. This enrichment introduces a characteristic length of the microstructure into the macroscopic material description. It can be done via a non-local theory, gradient terms or using viscosity in dynamic problems (Brinkgreve [6], Glema et al. [10], Marzec et al. [19], Wang et al. [27], Winnicki et al. [29]). More sophisticated formulations couple continuous and discontinuous descriptions (Simone et al. [22,23]). Non-locality can be also introduced via fractional differential operators (Beda [3], Błaszczyk [4], Lazopoulos i Lazopoulos [17], Sumelka [24], Sumelka et al. [25]).
Despite a huge amount of performed experiments on concrete specimens, there is still no consensus in describing the fracture process in concrete. Due to several factors influencing the results and a small number of specimens tested in a single experimental programme, it is very hard to properly compare the obtained results and to form general conclusions. Therefore, some experimental campaigns were executed recently to overcome these limitations. Hoover et al. [15] examined in plane geometrically scaled unnotched and notched concrete beams under three-point bending. Four different beam sizes and five different notch to depth ratios were analysed. In total 18 different geometries were defined and more than one hundred specimens were tested. Çağlar and Şener [7] examined geometrically identical beams (80 specimens), but cast in a horizontal position. A slightly smaller number of beams (34 specimens) was tested by Grégoire et al. [12]. They analysed beams of four different sizes and three different notch to depth ratios (unnotched and notched).
Different constitutive laws were used later to reproduce obtained experimental results. Hoover and Bažant [16] used the crack band and a bilinear softening law. They stated that an exponential softening curve was not able to give realistic results. Grégoire et al. [12] applied the isotropic damage constitutive law coupled with the integral non-local theory and an exponential curve defined in softening. A similar model was used in numerical simulations presented by Havlásek et al. [13]. In addition, they studied standard and distance based averaging methods in non-locality.
In the paper, numerical simulations of size effect in plain concrete beams under bending are presented. Two alternative crack descriptions are used: a smeared one via an elastoplastic with a Rankine criterion enriched by a non-local theory in the softening regime and a discrete one within XFEM. The influence of defined softening curve and assumed value of the initial fracture energy in both approaches on the obtained results is investigated. Such vast analysis with two different crack descriptions and two different softening curves has not been performed yet. In addition, the deficiencies of the initial fracture energy are pointed out. Numerical results are compared with experimental outcomes.
The paper is organised as follows. The constitutive models are described in Section 2. The problem data are given in Section 3. Obtained results and discussion are presented in Section 4. Conclusions and final remarks are listed in Section 5.

ELASTO-PLASTICITY
As a first option, a smeared description of cracks in concrete is used. An elasto-plastic constitutive law with the classical Rankine criterion is proposed. The yield criterion in 2D is postulated as where: σ 1 , σ 2 -the principal stresses, σ t -the tensile yield stress, κ C -the softening state variable.
The tensile yield stress σ t is defined in Sec. 2.3 and the state variable κ C is equal to the maximum principal plastic strain. Plastic strains are calculated assuming an associated flow rule.
When defining strain softening with the standard constitutive laws, classical finite element calculations fail. Obtained result are mesh dependent and a regularisation method is required to restore the well posedness of the boundary value problem. Here, an integral non-local theory is applied. Non-local rates of state variables dκ C are evaluated as (after Brinkgreve [6]) where: x -point under consideration, m -non-locality parameter, C d -rate of averaged state variable. The non-locality parameter m should be greater than 1 to effectively apply the non-local theory with elasto-plasticity (so called over-non-local formulation). The rate of averaged state variable is calculated as: where: V -volume of the integration domain, ξ -surrounding point coordinates, ω -the weighting function.
The weighting function reflects the influence of the surrounding points on the material's behaviour in the considered point x. Here a Gauss distribution is applied: where: r -the distance between points, l -the characteristic length.
The characteristic length is the link with the microstructure of the material. The characteristic length l (but also the non-locality parameter m, the definition of the weight function ω and the formulation of the constitutive law) influences the width of the localisation zone, which in general depends on concrete mesostructure. The averaging in Eqn. (4) is restricted only to a small neighbourhood around a point considered (the influence of points at the distance of r=3l is only of 0.01%). Therefore, although the weight function ω defined in Eqn (4) has unbounded support, only points at the distance no larger than r=3l from the integration domain V. For points x lying close to the boundary, only points lying within a circle with a radius r=3l and belonging to the specimen are taken into account (see Fig. 1a). In both cases, the denominator in Eqn (3) normalises the averaging operation; the uniform local field remains unchanged after applying the Eqn (3). Near notches so called "shading effect" is taken into account (Fig. 1b).

EXTENDED FINITE ELEMENT METHOD
As an alternative approach, cracks can be described as displacement jumps within continuum. In the paper eXtended Finite Element Method (XFEM) is utilised. It is based on a local partition of the unity (PUM) concept by (Melenk and Babuška [21]). It enables adding 'ad hoc' extra terms to a standard FE displacement field interpolation to capture displacement jumps across a crack. As a consequence, cracks can be placed through elements; they do not have to follow elements' edges.
In the paper the formulation (with minor improvements and modifications) proposed by Wells and Sluys [28] for cohesive cracks is adopted. The only major change is the application of the so-called shifted-basis enrichment proposed by Zi and Belytschko [30]. Theoretically, this improvement is equivalent to the classical model, but it simplifies the implementation (only two types of finite elements exist and total nodal displacements are equal to the standard ones).

Fig. 1. Averaging domain (grey area) for a point near boundary (a) and close to a notch (b)
Two constitutive relationships have to be declared to describe the behaviour of the material. Outside a crack in a solid (bulk) body a linear elasticity law is assumed. Along the crack a constitutive law between displacement jumps [[u]] and tractions t is postulated. The following loading function is assumed: The state variable κ D is calculated as a maximum crack opening [[u n ]] attained during the load history. In active loading (growth of the crack opening) the normal traction t n is equal to the yield traction t y defined as: t f y D t (6) where: D f -correction term, σ t -the yield stress defined in Sec. 2.3.
The correction term D f improves the performance of the model in tension-compression transition cases. It is defined as: where: d f -the drop factor (Cox [8]), f t -the tensile strength, G F -the total fracture energy.
In unloading/reloading a return to the origin (damage formulation) is assumed: In compression a penalty stiffness is taken. It is calculated as: In tangential direction, a linear dependence between shear tractions t s and tangential displacement jump [[u s ]] is assumed: where: T s -the initial shear stiffness.

SOFTENING CURVES
In both smeared and discrete crack descriptions two softening curves are used in the numerical simulations. In the present section, symbol κ should be taken either as the state variable κ C or as the state variable κ D . First, a bilinear relationship is taken (Fig. 2a): where: σ k -the yield stress at the kink point, κ k -the state variable at the kink point, κ u -the state variable at zero stress.
Alternatively, an exponential definition is used (Fig. 2b): u t y f exp (12) where: κ u -controls the slope of the curve.
More advanced exponential relationship was proposed by Hordijk ([14]) based on experimental outcomes.
In eXtended Finite Element Method, state variables κ k (only for bilinear softening) and κ u can be directly related to the initial fracture energy G f and the total fracture energy G F . The initial fracture energy G f (Bažant [1]) is the area under the initial tangent line from the peak at the stress -displacement curve (grey regions in Fig. 2), while the total fracture energy G F is the area under the whole curve. In bilinear softening these parameters can be derived as: and: while in exponential softening:

PROBLEM
In the paper the experiment performed by Hoover et al. [15] is numerically simulated. They examined 128 unnotched and notched plain concrete beams under three-point bending. Figure 3 presents the geometry of a beam. Four different beam sizes were tested with height (depth) D taken as 500, 215, 93 and 40 mm, labelled here as a huge, large, medium and small specimen, respectively. The length to depth ratio was kept constant and it was equal to 2.4 for all beams and the span length L was defined as 2.176D. Five different notch to depth ratios α 0 were assumed: 0.0 (no notch), 0.025, 0.075, 0.15 and 0.30. In total 18 different geometries were defined (small and medium beams with the ratio α 0 =0.025 were not cast). The thickness of all beams was B=40 mm and the width of the notch was equal to 1.5 mm for all specimens.

Fig. 3. Boundary conditions and geometry of a beam
A load was imposed in the middle of the top edge of the specimen. Steel plates were placed at supports and under the load. A hard (rigid) contact was assumed between concrete and steel plates at the supports and under the load. All tests were executed under crack mouth opening displacement control by increasing the distance between two chosen points at the bottom edge symmetrically with the respect to vertical axis of symmetry. The initial distance between these two points depended on the beam size D and the ratio α 0 (it was between 12.7 and 162 mm).

FE-CALCULATIONS INPUT DATA
Numerical simulations are performed in Abaqus Standard programme. The total elongation of the gauge is set to Δ=0.3 mm and an indirect displacement method is used. At least 250 increments are required to complete a simulation. Plane stress conditions are assumed with 3-node constant strain triangular finite elements. In elastoplasticity an approximated method is used in calculations of averaged quantities. In an integration point the influence of the neighbour points is determined with values from the previous iteration. The refined FE mesh along the vertical axis of symmetry is defined in all specimens with the maximum element size not greater than 1 mm for calculations with elasto-plasticity and 2 mm for XFEM simulations. The total number of finite elements varies between 6967 and 55251 and between 4981 and 11660 in calculations with smeared and discrete cracks, respectively. A denser FE-mesh for calculations with elasto-plasticity ensures the effective application of the non-local theory ever for small values of the characteristic lengths and obtaining mesh independent results. Some comparative simulations with XFEM have shown that identical results have been received using coarser and denser meshes. An exemplary FE-mesh used in XFEM calculations in the region around the notch for the medium beam and the notch to depth ratio α 0 =0.15 is shown in Fig. 4.

Fig. 4. Exemplary FE mesh around the notch
In the simulations the Young's modulus is assumed as E=41. 25 GPa and the Poisson's ratio is v=0.172 (taken from experiments by Hoover et al. [15]). The total fracture energy is set to G F =70 N/m (after Hoover and Bažant [16]). In bilinear softening the yield stress at the kink point is always defined as σ k =0.15f t . In XFEM the drop factor is equal to d f =10 4 and the shear stiffness is T s =0.0 N/m 3 . Material for steel plates is taken as a linear elastic with the Young's modulus E s =200 GPa and the Poisson's ratio v s =0.3. where: In order to rate a simulation quality, the following relative error Err 0 is introduced:

RESULT AND ERROR MEASURES
The whole set is evaluated using errors Err 1 and Err 2 : where: N -the number of beams (N=18). In statistics, the errors Err 1 and Err 2 are called the mean percentage error (MPE) and the mean absolute percentage error (MAPE), respectively. The presence of the negative error values allows for distinguish between underestimated and overestimated numerical results.

RESULTS WITH XFEM
First, simulations with bilinear softening curve are executed. In order to find the values of the tensile strength f t and the initial fracture energy G f , seed calculations are performed.
The tensile strength f t changes between 3.6 MPa and 5.6 MPa with an increment 0.4 MPa. The initial fracture energy G f varies between 30 N/m and 54 N/m with an increment 2 N/m. Obtained error Err 1 and Err 2 isolines are depicted at Fig. 5. It can be seen that for both measures the best f t and G f pairs form an inclined line (or region), but they do not indicate a single optimum point. Please note, however, that not only peak values, but also agreement between whole numerical and experimental curves has to be taken into account. Finally, the tensile strength is set to f t =5.0 MPa and the initial fracture energy is taken as G f =48 N/m. With these values the following errors are obtained: Err 1 =0.67% and Err 2 =2.69%. Figure 6a presents a comparison between nominal strength obtained in the experiment and FE-calculations. Generally, a very good agreement can be observed. The largest deviation (the error Err 0 =10.7%) is for the small beam (D=40 mm) and the notch to depth ratio α 0 =0.3. Whole force -displacement curves are presented in Fig. 7. It can be seen that numerical results fit nicely into experimental ranges (grey regions). This fact confirms the proper assumption of the total fracture energy G F =70 N/m. The performance of the exponential softening curve with the following parameters: f t =4.8 MPa and G f =35 N/m is studied also (values based on some initial calculations). In the case of the exponential softening, there is no possibility to control independently both: initial G f and total G F fracture energies. The initial fracture energy G f is always equal to 50% of the total fracture energy G F . FE-calculations with above parameters produce the following errors: Err 1 =0.20% and Err 2 =3.65%, comparable to values obtained with the bilinear softening. Here the worst specimen is the huge beam (D=500 mm) and the notch to depth ratio α 0 =0.3; it gives the error Err 0 =8.9%. Comparison between experimental and numerical nominal strength is done in Fig. 6b. Again, a good agreement is achieved (here some differences occur also for unnotched beams).

RESULTS WITH ELASTO-PLASTICITY
Next the numerical simulations are repeated with the elasto-plastic constitutive law with a Rankine criterion. Two characteristic lengths are assumed: l=5 mm and l=2 mm. The values of κ u (for both softening curve definitions) are chosen in such a way to obtain the identical force-displacement curves with XFEM and elasto-plastic model with any of two characteristic lengths analysed. In calculations with the bilinear softening the ultimate value of the state variable is taken as κ u =4.05·10 -3 and κ u =10.13·10 -3 for the characteristic length l=5 mm and l=2 mm, respectively (to obtain initial fracture energy G f =48 N/m). The tensile strength is f t =5.0 MPa (as in XFEM calculations). Simulations give the following errors: Err 1 =Err 2 =9.76% for the characteristic length l=5 mm and Err 1 =Err 2 =6.80% for the characteristic length l=2 mm. Comparison between numerical and experimental nominal strengths is made in Fig. 8. In general, both parameter sets overestimate experimental outcomes. The worst specimen returns errors Err 0 =28.1% and Err 0 =19.2% (the small beam with α 0 =0.3). Fortunately, the overestimation reduces with decreasing the characteristic length. It suggests the smaller values should be assumed in simulations. Slightly better results are obtained from simulations with the exponential softening law. Here the ultimate value of the softening parameter is assumed as κ u =0.758·10 -3 and κ u =1.894·10 -3 for the characteristic length l=5 mm and l=2 mm, respectively (to obtain initial fracture energy G f =35 N/m). The tensile strength is f t =4.8 MPa (as in XFEM calculations). The following errors are returned: Err 1 =4.91% and Err 2 =6.82% for the characteristic length l=5 mm and Err 1 =2.43% and Err 2 =4.16% for the characteristic length l=2 mm. Again, in general too high peak loads are obtained in calculations, although last errors Err 2 is already comparable with error values produced in XFEM simulations. Figure 9 presents a comparison between numerical and experimental nominal strengths, while force -CMOD curves are shown in Fig. 10. Despite the overestimation of peak loads, numerical curves fit into experimental limits.

CONCLUSIONS
In the paper the size effect phenomenon in plane concrete beams under bending has been investigated. The numerical calculations using smeared and discrete methods describing cracks in concrete have been executed. Obtained results have been compared with experimental outcomes. Both approaches gave results consistent with experiments, although smaller errors have been attained with XFEM. In elasto-plasticity, simulations with smaller values of the characteristic length have returned better results (smaller errors have been obtained). Both softening curves have produced comparable results, assuming different initial fracture energies. This observation points some problems in unique identification of the initial fracture energy, when different curve definitions in softening are used. When assuming a curved softening relationship (e.g. exponential one) purely geometrical definition of the initial fracture energy is not sufficient and the curvature of the curve in the domain after the peak should be taken into account.
The ongoing research is focused on improving the agreement of numerical results obtained with elasto-plastic constitutive law with experimental curves. More advanced, anisotropic averaging schemes in non-locality proposed by Giry et al. [9] and Grassl et al. [11] are investigated. The consequences of decreasing the fracture energy in the boundary layer (Bažant et al. [2]) are also examined.