FE investigations of the effect of fluctuating local tensile strength on coupled energetic–statistical size effect in concrete beams

Engineering Structures 103 (2015) 239–259

Contents lists available at ScienceDirect

Engineering Structures journal homepage: www.elsevier.com/locate/engstruct

FE investigations of the effect of fluctuating local tensile strength on coupled energetic–statistical size effect in concrete beams E. Syroka-Korol a, J. Tejchman a,⇑, Z. Mróz b a b

´ sk University of Technology, Poland Faculty of Civil and Environmental Engineering, Gdan Institute of Fundamental Technological Research, Polish Academy of Sciences, Warsaw, Poland

a r t i c l e

i n f o

Article history: Received 16 January 2015 Revised 23 July 2015 Accepted 8 September 2015 Available online 26 September 2015 Keywords: Autocorrelation length Concrete beam Elasto-plasticity Non-local softening Random field Size effect Strain localization Tensile strength

a b s t r a c t The effect of fluctuating local tensile strength on a coupled energetic–statistical size effect in plain concrete beams under bending was numerically investigated. First, the influence of varying autocorrelation length of the random field describing a spatial variation of local tensile strength was studied. Next, the influence of the coefficient of variation of local tensile strength was analyzed. The numerical FE investigations were performed for unnotched concrete beams of similar geometry under quasi-static three point bending within elasto-plastic model assumptions, accounting for non-local softening. The FE analyses were carried out for four different beam sizes. In the calculations, the tensile strength took the form of spatially correlated random field described by a Gaussian distribution. Ó 2015 Elsevier Ltd. All rights reserved.

1. Introduction A size effect is an inherent property of the behavior of many engineering materials. It means that the nominal strength (or nominal stress at failure) varies with a characteristic size of a structural member. In the case of concrete materials, both the nominal structural strength and structural ductility always decrease with increasing element size under tension [1–3]. Thus, concrete becomes ductile on a small scale and perfectly brittle on a sufficiently large scale. The results from laboratory tests which are scaled versions of the actual structures cannot be directly transferred to them. For simply supported beams under bending, the nominal strength can be expressed in terms of the elastic bending stress at maximal load or the effective stress accounting for coupled tension and shear effects [3]. Two mechanical size effects are of major importance in quasi-brittle and brittle materials: energetic (or deterministic) and statistical (or stochastic) one. The deterministic size effect is caused by the formation of a region of intense strain localization with a certain volume (micro-crack region – called also fracture ⇑ Corresponding author. E-mail addresses: [email protected] (E. Syroka-Korol), [email protected] (J. Tejchman), [email protected] (Z. Mróz). http://dx.doi.org/10.1016/j.engstruct.2015.09.011 0141-0296/Ó 2015 Elsevier Ltd. All rights reserved.

process zone FPZ) which precedes macro-cracks [3,4]. The nominal structural strength which is sensitive to the size of FPZ cannot be appropriately estimated in laboratory tests since it differs for various specimen sizes. Strain localization area is not negligible to the element cross-section dimensions and is large enough to cause a significant stress redistribution in the structure and an associated energy release. The specimen strength increases with increasing ratio lc/D (lc – characteristic length of the microstructure influencing both the size and spacing of localized zones, D – characteristic structure size). In turn, a statistical (stochastic) effect is caused by the spatial variability/randomness of the local material strength. The first statistical theory was introduced by Weibull [5] (called also the weakest link theory) which postulates that a structure is as strong as its weakest component. The structure fails when its strength is reached in the weakest spot, since the stress redistribution is not considered. The Weibull’s size effect model is a power law and is of particular importance for large structures that fail as soon as macroscopic fracture initiates in one small material element. However, it is not able to account for a statistical correlation of local material properties, does not include any characteristic length of micro-structure (i.e. it ignores a deterministic size effect) and it underestimates the effect of small- and intermediate-sizes. Combining the energetic theory with the Weibull statistical theory, a general energetic–statistical

240

E. Syroka-Korol et al. / Engineering Structures 103 (2015) 239–259

theory was developed [6]. The deterministic size effect was obtained for not too large structures and the Weibull statistical size effect was obtained as an asymptotic limit for very large structures. The numerical FE analysis is aimed at investigating the influence of the spatial variability of local tensile strength modeled by a random field on a combined energetic–statistical size effect in unnotched plain concrete beams under three-point bending by taking properly strain localization into account. It is noted that the beams belong to the most important structural elements. Secondly, due to its simplicity, the three-point bending test is the most popular indirect experimental test for measuring concrete flexural tensile strength. Thirdly, our analysis constitutes the continuation of our previous simulations of a coupled energetic– statistical size effect in 3-point bending concrete beams within stochastic elasto-plasticity and non-local softening [7–10]. The spatial fluctuation of local tensile strength described by a random field is under a great influence of an assumed autocorrelation length. Thus, the effect of the varying autocorrelation length was taken into account in our analysis of the strength and postpeak behavior of concrete beams. The influence of varying material coefficient of variation of local tensile strength was also investigated. The plane stress calculations were performed with 4 different concrete beam sizes of a similar geometry scaled in 2 dimensions (height and length). The results were compared with the existing size effect rules proposed by Bazant. The finite element method with a simple elasto-plastic constitutive model using the Rankine criterion coupled with a non-local softening rule was used, which was suitable to capture strain localization under tension [7]. The deterministic calculations were performed assuming the constant tensile strength. The statistical analyses were carried out using solely random fields with a squared exponential autocorrelation function to model the varying local tensile strength. The random fields were generated using the Karhunen–Loeve expansion and wavelet Galerkin approach [11,12]. The approximated results were obtained using a stratified sampling method (SSM) belonging to a group of variance reduced Monte Carlo methods [13]. This approach enables a significant reduction of the number of simulations without losing the accuracy of calculations. It was demonstrated that for quasi-brittle materials such as concrete that the most realistic approach should assume the local material strength distribution with the Gauss core and power law left tail with the zero threshold [14]. The strength of small and medium beam sizes is mainly affected by the core of the Gaussian pdf whereas very large structures which fail in a perfectly brittle manner are governed by the left tail and follow the Weibull size effect. Since we were mainly focused on the intermediate beam sizes in the engineering practice (beam height between 80 mm and 1920 mm), then the material strength variation was assumed to be described by a normal distribution (core and both tails) for all beam sizes. The calculations for the significantly larger beam sizes were limited by the long computation time caused by the necessity of the use of small finite elements to properly describe strain localization. According to Grassl and Bazant [14] the autocorrelated distribution of local material strength is responsible for the statistical size effect observed in concrete. It was also suggested [14] that the autocorrelation length should be several times greater than the FPZ size. A physical justification to use the autocorrelated distribution of local material strength supported by experimental results for multi-filament yearns was given by Vorechovsky and Chudoba [15]. They have demonstrated [15] that the finite autocorrelation length creates un upper limit for the statistical size effect i.e. the mean strength of small specimens is limited, contrary to the Weibull theory. In our previous paper [10], the spatial variation of local axial tensile strength was defined using a homogeneous second order autocorrelation function [16], which produced a weaker

correlation than that specified by a square exponential function. The decay parameters provided a stronger autocorrelation in a horizontal (
80 mm) than in the vertical direction (
30 mm). The mean tensile strength was 3.6 MPa and the standard deviation 0.424. The difference between the deterministic flexural strength and mean stochastic strength increased with increasing specimen size whereas the standard deviation decreased. The deterministic flexural strength approached a horizontal asymptote for all beam sizes. The mean statistical strength corresponded to the deterministic strength for small beam sizes whereas for large beam sizes approached the inclined asymptote with a slope n/m (n – the number of spatial dimensions in which the structure is scaled and m – the dimensionless Weibull modulus responsible for the slope of a large-size asymptote). For the assumed homogenous autocorrelation function [16], the Weibull modulus (calculated based on the nominal strength variation) was m = 48 for the largest beam which differed from the value recommended by Bazant and Novak m = 24, based on experiments [17]. It has been also found out that the Weibull modulus m strongly depends on the beam size. The combined stochastic–deterministic size effects were simulated among others by Carmeliet and Hens [18], Carmeliet and de Borst [19], Gutierrez and de Borst [20], Frantziskonis [21], Vorechovsky [22], Bazant and Novak [17], Bazant et al. [23], Yang and Xu [24], Grassl and Bazant [14] and Elias et al. [25]. The most comprehensive combined calculations were performed by Grassl and Bazant [14] using the random discrete particle model for geometrically similar concrete prisms of 5 sizes (characteristic size D = 40–800 mm) subjected to direct tension. The interaction of the failure zone size and autocorrelation length lcor was investigated with the square exponential autocorrelation function with 3 different values of 40 mm, 20 mm and lcor ? 0. The effect of the varying material coefficient of variation 0.1–0.4 was also analyzed. It was concluded [14] that the main parameter controlling the statistical size effect was the ratio of the autocorrelation length to the fracture process zone size or to the representative volume element (RVE) size (found to be equal to the two–three distances of large aggregate grains). The results demonstrated that the statistical size effect was significant for lcor >> RVE and negligible for lcor << RVE. The higher material coefficient of variation cf, the stronger was the statistical size effect. In the paper by Yang and Xu [24], a heterogeneous cohesive crack model to predict the macroscopic strength based on meso-scale random fields of fracture properties was proposed. A concrete notched beam subjected to mixed-mode fracture was modeled. The effect of the various important parameters on the crack path, peak load, macroscopic ductility and overall reliability (including the variance of random fields, autocorrelation length, and shear fracture resistance) was discussed. It has to be noted that statistical FE analyses are very difficult due the lack of experimental results concerning a autocorrelation function and autocorrelation length in engineering materials of a different size. The innovative points in our numerically oriented paper are the comprehensive statistical calculations of both the failure load and the softening modulus for different concrete beam sizes. The statistical analysis takes into account the varying autocorrelation length of a random field describing the spatial fluctuation of local material tensile strength and the coefficient of variation of local material tensile strength. The effect of the autocorrelation length on the mean nominal strength is analyzed with respect to the size of a localized zone wloc  hloc (width  height) dependent upon a characteristic length of micro-structure. The FE results were compared with the combined energetic–statistical size effect by Bazant [23]. The magnitude of Weibull modulus was also discussed. The outline of the present paper is as follows. First, after Introduction (Section 1) and a description of size effect laws (Section 2), the employed constitutive elasto-plastic model with non-local

241

E. Syroka-Korol et al. / Engineering Structures 103 (2015) 239–259

softening is summarized (Section 3). The description of simulation of discrete random fields, finite element discretization and boundary conditions are provided in Sections 4–6. The numerical results of deterministic and stochastic size effects are described and discussed in Section 7. The conclusions are stated in Section 8. 2. Existing size effect laws There exist currently two different theories of the size effect in quasi-brittle structures: the energetic–statistical theory [3,4,17] and fractal theory [2]. Two size effect laws (SEL) proposed by Bazant [3,4,23] for geometrically similar structures allow for determining their nominal strength. The size effect of Type I (Fig. 1) applies to structures of a positive geometry without notches or pre-existing cracks for which the maximum load occurs as soon as the fracture process zone (FPZ) is fully developed and propagates along the fracture surface, followed by a macroscopic crack. The requirement of a positive geometry enables to incorporate the weakest link theory by Weibull. The structures obeying the size effect of Type I are sensitive to the material randomness. The nominal strength is strongly affected by the material heterogeneity and decreases with increasing structure characteristic dimension. The energetic size effect of Type I assumes that the material strength is bounded for small sizes by a plasticity limit whereas for large sizes the material follows linear elastic fracture mechanics. The following general analytical formulae for energetic size effect predicted by the asymptotic matching was proposed by Bazant [3,23]


1 rDb r rN ðDÞ ¼ f 1 1 þ ; r D þ lp

ð1Þ

where rN is the nominal strength, D is the characteristic structure size while f1 r , Db, lp and r denote unknown parameters to be determined. The parameter f1 represents the elastic-brittle strength r reached for large structures, r controls the strength evolution for varying size, Db is the deterministic scaling length having the meaning of the length of the damaged layer (if Db = 0, the behavior is elastic-brittle) and lp introduces a horizontal asymptote for extremely small sizes (lp is correlated with Db). For large sizes D ? 1, Eq. (1) predicts an asymptotic approach to the horizontal asymptote f1 r . The extended universal formula for a coupled deter ministic–statistical size effect law reads as follows [23]:

Fig. 1. Size effect law of Type I by Bazant for geometrically similar concrete structures failing at crack initiation from smooth surface (D – characteristic size, rN – nominal strength): (a) deterministic size effect and (b) coupled deterministic– statistical size effect [3,6,17].

rN ðDÞ ¼

1 fr




Ls D þ Ls

rnm

rDb þ D þ lp

!1r ;

ð2Þ

where Ls is the statistical scaling length producing an upper bound for fr when D ? 0. For large sizes D ? 1, Eq. (2) asymptotically reaches the dominating Weibull size effect rule with the slope equal to n/m (Fig. 1). 3. Constitutive elasto-plastic model with non-local softening In order to describe the behavior of concrete under tension during three-point bending, the Rankine criterion was used, for which the yield function f with isotropic softening was defined as

f ¼ maxfr1 ; r2 ; r3 g  rt ðjÞ;

ð3Þ

where ri, i = 1, 2, 3, are the principal stresses, rt is the tensile yield stress and j denotes the softening parameter equal to the maximum principal plastic strain e1p [7,26,27]. The associated flow rule was assumed. In order to model concrete softening under tension, the exponential curve by Hordijk [28] was chosen (Fig. 2):

rt ðjÞ ¼ f 0t ½ð1 þ ðA1 jÞ3 Þ expðA2 jÞ  A3 j;

ð4Þ

ft0

where is the uniaxial tensile strength of concrete. The constants A1, A2 and A3 were assumed in the form

A1 ¼

c1

ju

; A2 ¼

c2

ju

; and A3 ¼

1 

ju

 1 þ c31 expðc2 Þ;

ð5Þ

wherein ju denotes the ultimate value of the softening parameter (ju = 0.005). The constants ci were c1 = 3 and c2 = 6.93. A non-local theory has been used as a regularization technique [29–31] wherein the principle of a local action does not take place any more. Polizzotto et al. [32] laid down a thermodynamically consistent formulation of non-local plasticity. In the calculations, the softening parameter was assumed to be non-local [33]

R

xðkx  nkÞjðnÞdn with i ¼ 1; 2; j
ðxÞ ¼ ð1  mnl ÞjðxÞ þ mnl V R xðkx  nkÞdn V ð6Þ
i (x) is the non-local softening parameter, V denotes the where j body volume, x is the coordinate vector of the considered point, n is the coordinate vector of the surrounding points, x denotes the weighting function and mnl is the additional non-local parameter controlling the size of the localized plastic zone. As a weighting function x, the Gauss distribution function was used

Fig. 2. Evolution of softening curve under tension rt = f(j) (rt – tensile yield stress, j – softening parameter) by Hordijk [28].

242

E. Syroka-Korol et al. / Engineering Structures 103 (2015) 239–259

Fig. 3. Stratified sampling scheme (cumulative distribution function cdf against sampling parameter a).

 2 1  lrc xðrÞ ¼ pffiffiffiffi e lc p

ð7Þ

where lc is the characteristic length of micro-structure (related usually to the aggregate size [3]) and the parameter r is the distance between material points. The characteristic length is usually determined with an inverse identification process of experimental data [34]. The averaging in Eq. (7) was restricted to a small representative area around each material point (the influence of points at the distance of r = 3  lc was only of 0.01%). The softening nonlocal parameters near boundaries and at both sides of a localized zone were always calculated on the basis of Eqs. (6) and (7) (which satisfy the normalizing condition). Other approaches were proposed by Polizzotto [35] and Jirásek et al. [36]. In order to simplify the calculations, non-local rates were replaced by their approximations calculated with known total strain increments [36]. The non-local

parameters lc and mnl affect the load–displacement diagram (peak value and softening rate) and the width of localized zones [7]. The non-local model was implemented in the commercial finite element code ABAQUS [37] with the aid of subroutine UMAT (user constitutive law definition) and UEL (user element definition) for efficient computations. The edge and vertex in the Rankine yield function were taken into account by the interpolation of plastic multipliers for 2–3 yield planes according to the Koiter’s rule [38]. To capture the snap-back behavior in large and very large size beams, the so-called arc-length technique may be used. The actual load vector P is defined as kPmax, where k is the multiplier and Pmax is the maximum constant load vector. In general, the determination of the length of the arc in the Pu space (u – displacement vector) involves the displacements of all nodes (as e.g. the Riks method available in [37]). For problems involving strain localization, it is more suitable to use an indirect displacement control method proposed by Jirasek and Bazant [39] where only selected nodal displacements are considered to formulate an additional condition in the Pu space. This method was used in our analyses. The horizontal distance between two nodes lying on the opposite sides of the notch was chosen as the control variable (CMOD, crack mouth open displacement). Our elasto-plastic model with nonlocal exponential softening was verified in size effect calculations of concrete beams [8], based on corresponding experimental results by Le Bellego et al. [34] and size effect calculations of reinforced concrete beams [40], based on corresponding experiments by Syroka-Korol [9]. 4. Random field of tensile strength The spatially correlated random fields of the tensile concrete strength were described by the Gauss distribution function and the homogenous squared exponential autocorrelation function C

Cðx1 ; x2 Þ ¼ exp 

ðx1  x2 Þ2 2

lcor

!

;

Fig. 4. Calculation of function a(x): (a) critical area selection and (b) definition of searching area and scheme of a(x) estimation.

ð8Þ

E. Syroka-Korol et al. / Engineering Structures 103 (2015) 239–259

243

Fig. 5. Prediction of localized zone position for small-size beam (D = 80 mm): (a) exemplary distributions of local tensile strength, (b) calculated a(x)-functions and (c)
in FE simulations. distributions of tensile non-local softening parameter j

Fig. 6. Prediction of localized zone position for very large size beam (D = 1920 mm): (a) exemplary distributions of local tensile strength, (b) calculated a(x)-functions and (c)
in FE simulations. distributions of non-local tensile softening parameter j

244

E. Syroka-Korol et al. / Engineering Structures 103 (2015) 239–259

where |x1  x2| is the distance between two points in a Cartesian space and lcor is the autocorrelation length related to the material randomness at macro-scale due to casting and mixing techniques (the same in the both directions xi). It was assumed that lcor was not linked to lc. The auto-covariance function had the following spectral decomposition

Cðx1 ; x2 Þ ¼

1 X ki f i ðx1 Þf i ðx2 Þ;

ð9Þ

i¼1

where the eigenvalues ki and eigenfunctions fi(x) were the solution of the Fredholm integral equation of the second kind

Z Cðx1 ; x2 Þf i ðx1 Þdx1 ¼ ki f i ðx2 Þ

ð10Þ

D

The analytical solution of Eq. (10) is available with e.g. the exponential autocorrelation function [41]. In other cases, Eq. (14) has to be solved numerically. In our study, the wavelet Galerkin method proposed by Phoon et al. [11,12] was used wherein conventional bases (e.g. trigonometric or polynomials) were replaced by Haar wavelets. The spatially correlated random fields H(x, h) (here with the zero mean and unit variance) were defined according to the Karhunen–Loëve expansion [42,43] by an infinite linear combination of orthogonal functions with random coefficients

Hðx; hÞ ¼

1 pffiffiffiffi X ki ni ðhÞf i ðxÞ;

ð11Þ

i¼1

where ni(h) is the vector of uncorrelated random variables sampled ^ hÞ was from N(0, 1) distribution. The approximated solution Hðx;

Fig. 7. Convergence of flexural tensile strength fr (a) and its coefficient of variation COVfr (b) with increasing number of simulations n compared to results obtained for 124 and 12 simulation with SSM (small-size beam D = 80 mm).

obtained by truncating the series in Eq. (11) after M terms. The resulting error was minimized by sorting eigenvalues in a decreasing order and controlling the sum of ki for i = 1, 2, . . . , M. 5. Stratified sampling method

Fig. 8. Geometry of free-supported unnotched concrete beams subjected to three-point bending.

Initially 2000 samples of a single random field were generated in order to reproduce the assumed Gauss probability distribution function and auto-covariance matrix. The stratified sampling method (SSM) originated by Newman [44] was adopted to reduce the number of samples used in further FE analysis. After the original set of 2000 samples was generated, a sampling parameter a was defined. Next, samples were classified and arranged in an increasing order with reference to the a parameter and the corresponding cumulative probability (frequency) was calculated. Afterwards, the representative set of 12 samples was chosen from 12 equally probable intervals (Pi, i = 1, . . . , 12) using the mid-point method (Fig. 3). The scalar sampling parameter a reflects the features of a two-dimensional random field and its definition is

Fig. 9. FE mesh in small-size concrete beam (s – width of region with finer mesh).

E. Syroka-Korol et al. / Engineering Structures 103 (2015) 239–259

245

Fig. 10. Effect of varying autocorrelation length lcor: (a) evolution of autocorrelation function (Eq. (8)) with lcor = 10, 20, 30 and 40 mm and (b) random distribution of local tensile strength ft0 with lcor = 10 mm and lcor = 40 mm.

Fig. 11. Normalized flexural tensile stress–deflection curves fr/ft0 = f(u/D) from deterministic analyses for: small D = 80 mm (‘a’), medium D = 160 mm (‘b’), large D = 320 mm (‘c’) and very large beam D = 1920 mm (‘d’). (A) line fr/ft0 = 1.5 and (B) line fr/ft0 = 1.0, fr – flexural tensile strength) [10].

Fig. 12. The height of localized zone hloc at peak load versus beam size D from deterministic analysis (in black [10]) and average values from statistical simulations (in red). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

crucial for the selection of the representative set of samples. In our analysis the sampling parameter a was defined as to reflect the probable failure load of the beam. In this way the original set of samples was arranged in an increasing order with reference to the expected failure force and as a result the chosen 12 samples

led to the sufficient representation of the entire range of the probable failure force. The a-parameter was evaluated as follows. After the random field spatial discretization was performed (similar to the FE mesh discretization), the critical area was defined (Fig. 4a). The critical area covered the length of a random field

246

E. Syroka-Korol et al. / Engineering Structures 103 (2015) 239–259

Fig. 13. Evolution of normalized flexural tensile stress fr/ft0 versus normalized localized zone height hloc/D in deterministic simulations [10].

Fig. 14. Calibrated energetic (Eq. (1), solid line) and energetic–statistical (Eq. (2), dashed line) SEL Type I by Bazant compared with results from deterministic (black circles) and statistical FE analyses with lcor = 40 mm (red squares). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

and the height of a localized zone hloc at the peak load (Fig. 4a). This area was selected due to the dominant failure mode in plain concrete beams which always occurred at the crack initiation. Next, the average value f t;i (i = 1, . . . , N, N – the number of FE elements along the critical area, i.e. columns) of the local tensile strength ft 0

was calculated which was further normalized as f N;i ¼ f t;i =f t using the theoretical mean tensile strength ft0 . Afterwards, the value of f N;i was divided by the corresponding normalized elastic stress due to bending ri. Finally, the value a(x) was calculated as the average value in the searching area wloc  hloc (wloc is the width of a localized zone at the peak load) (Fig. 4b). The searching area was shifted along the critical area in order to determine the lowest value of a(x). The value of the minimum of a(x) function defined the sampling parameter a = amin(x). The position of the minimum of a(x) function indicated a probable place of a localized zone. Figs. 5 and 6 present the exemplary distribution of the local tensile strength in the small (D = 80 mm, Fig. 5a) and very large beam (D = 1920 mm, Fig. 6a) compared to the calculated a(x)-profiles in Figs. 5b and 6b. The distributions of the non-local tensile
1 obtained from FE analyses are shown in softening parameter j

Fig. 15. Post-critical behavior in beams with different height D in deterministic analyses: (a) method of calculating softening angle b or softening angle c for normalized stress–deflection curve fr/ft0 = f(u/D) (red line) and (b) values of b and c against D. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Table 1 Results from FE simulations with varying lcor (small-size beam D = 80 mm). lcor (mm)

lF (kN)

sF (kN)

COVF

2 3 4 5 6 7 8 9 10 20 30 40 60 80 100 120 140

3.77 3.68 3.79 3.74 3.69 3.68 3.64 3.74 3.62 3.69 3.71 3.77 3.85 3.82 3.84 3.81 3.83

0.107 0.129 0.178 0.137 0.154 0.171 0.178 0.222 0.107 0.342 0.438 0.379 0.377 0.369 0.370 0.347 0.364

0.03 0.04 0.05 0.04 0.04 0.05 0.05 0.06 0.03 0.09 0.12 0.10 0.10 0.10 0.10 0.09 0.09

Figs. 5c and 6c. The analytical a(x)-function was able to correctly predict the position of a localized zone and the value of amin(x) was related to a probable failure force. The bearing capacity of

E. Syroka-Korol et al. / Engineering Structures 103 (2015) 239–259 Table 2 Results from FE simulations with varying lcor (medium-size beam D = 160 mm). lcor (mm)

lF (kN)

sF (kN)

COVF

5 10 20 30 40 60 80 100 120 140

6.33 6.42 6.24 6.28 6.49 6.67 6.54 6.73 6.64 6.73

0.182 0.381 0.501 0.529 0.653 0.586 0.599 0.684 0.723 0.677

0.03 0.06 0.08 0.08 0.10 0.09 0.09 0.10 0.11 0.10

Table 3 Results from FE simulations with varying lcor (large-size beam D = 320 mm). lcor (mm)

lF (kN)

sF (kN)

COVF

5 10 20 30 40 60 80 100 120 140

11.45 11.31 11.02 10.88 11.33 11.51 11.60 11.80 12.06 12.24

0.489 0.451 0.591 0.413 0.630 1.104 1.157 1.245 1.218 1.496

0.04 0.04 0.05 0.04 0.06 0.10 0.10 0.11 0.10 0.12

Table 4 Results from FE simulations with varying lcor (very large-size beam D = 1920 mm). lcor (mm)

lF (kN)

sF (kN)

COVF

m

10 20 30 40 60 80 100 120 140

61.28 58.88 56.13 55.00 55.60 56.25 56.77 57.21 58.65

2.897 2.439 2.861 1.580 3.452 3.358 5.263 4.847 6.544

0.05 0.04 0.05 0.03 0.06 0.06 0.09 0.09 0.11

26 30 24 44 20 20 13 14 11

the small size beam corresponding to amin(x) = 0.88 (Fig. 5A) was Fmax = 3.17 kN while for amin(x) = 1.12 was Fmax = 4.13 kN (Fig. 5B). In the very large size beam, the failure force was Fmax = 48.4 kN for amin(x) = 0.75 (Fig. 6A) and Fmax = 57.7 kN for amin(x) = 0.97 (Fig. 6B). A similar dependence was observed in many our initial FE simulations for all beam sizes. Hence a general conclusion might be drawn that the lower amin(x), the lower was the failure force Fmax. A preliminary FE stochastic analysis was performed in order to investigate the solution convergence of SSM. Fig. 6 shows the FE results including a convergence process of the mean flexural tensile strength fr (Fig. 7a) and its coefficient of variation COVfr (standard deviation/mean value) (Fig. 7b) for the small-size plain concrete beam D = 80 mm with the autocorrelation length lcor = 80 mm (Eq. (8)). The computed mean flexural tensile strength was fr = 5.306 MPa and fr = 5.318 MPa while the coefficient of variation was COVfr = 0.104 and COVfr = 0.100 from 124 and 12 simulations, respectively. The presented results indicate that SSM provides a convergent solution for 12 FE simulations. 6. FE input data The two-dimensional FE plane stress analysis was performed for unnotched plain concrete beams subjected to three-point

247

bending. Four beam sizes simultaneously scaled in two dimensions i.e. height and length were analyzed. The beam height was D = 80 mm (called small size beam), D = 160 mm (called medium size beam), D = 320 mm (called large size beam) and D = 1920 mm (called very large size beam). The beam span L was proportional to its height L = 3D (Fig. 8) while the thickness was kept constant t = 40 mm independently on the beam size. The geometry of the small (D = 80 mm), medium (D = 160 mm) and large size beam (D = 320 mm) was similar as in the corresponding experiments by Le Bellego et al. [34]. The concrete constants were adopted from our previous FE investigations [8]: the uniaxial tensile strength of concrete ft0 = 3.6 MPa, the modulus of elasticity E = 38.5 GPa and the Poisson ratio t = 0.24. The tensile fracture energy was assumed Gf = gf  wloc = 52 N/m with gf – area under the softening function and wloc  3lc – the width of the localized zone (wloc is independent of the beam size and material randomness [10]). The ultimate softening parameter was ju = 0.005 and the characteristic length of micro-structure was lc = 5 mm (assumed based on our experimental results with the DIC technique [45,46]). In order to properly capture strain localization in concrete, the mesh was very fine in the mid-part of beams (Fig. 9) (the quadrilateral element size was 1.5  2 mm2). The quadrilateral elements were divided into triangular finite elements to avoid volumetric locking. Totally, 13,820 (small-size beam), 39,900 (medium-size beam), 104,780 (large-size beam) and 521,276 (very large-size beam) triangular elements were used, respectively. The width of this region s (Fig. 9) was determined with preliminary calculations: s = 120 mm (D = 80 mm), s = 180 mm (D = 160 mm), s = 240 mm (D = 320 mm) and s = 1920 mm (D = 1920 mm). A quasi-static deformation of the small and medium beam was imposed through a constant vertical displacement increment Du prescribed at the upper mid-point of the beam top. In the case of large and very large beams (to capture the snap-back behavior), a procedure described in Section 3 was used. The deterministic calculations were performed with the uniformly distributed local tensile strength ft0 = 3.6 MPa [10]. The statistical studies were performed with the small, medium, large and very large size beams having the spatially correlated randomly distributed local tensile strength ft only. The mean value of the local tensile strength lft = 3.6 MPa (lft = ft0 ) was kept constant for all beam sizes. Note that due to the fixed value of the ultimate softening parameter ju = 0.005 (Eq. (5)), the tensile fracture energy Gf slightly increased (approximately in a proportional way) with growing local tensile strength ft (the Irwin’s characteristic length was not constant). In the reality, the fracture energy should rather diminish with increasing ft (concrete brittleness increases with its increasing strength). The effect of the varying Gf = f( ft) on the mean flexural beam strength was found to be negligible based on preliminary statistical calculations. First, the effect of strength of the spatial autocorrelation of ft was studied by changing lcor (Eq. (8), Fig. 10a). The material coefficient of variation COVft = 0.12 (COVft = sft/lft, sft – the standard deviation) was kept constant when lcor varied. The spatial autocorrelation became weaker with decreasing lcor, hence the distribution of ft was more heterogeneous (Fig. 10b). Whereas with increasing lcor, the spatial autocorrelation became stronger and the distribution of ft was more homogenous (Fig. 10c). The effect of the varying autocorrelation length was studied with lcor = 2–150 mm, lcor = 5–140 mm, lcor = 5–140 mm and lcor = 10–140 for the small (D = 80 mm), medium (D = 160 mm), large (D = 320 mm) and very large (D = 1920 mm) beam, respectively. Next, the effect of the material coefficient of variation was investigated with 0.12, 0.16 and 0.20. The mean value lft = 3.6 MPa was kept constant while the standard deviation sft varied. The autocorrelation length was

248

E. Syroka-Korol et al. / Engineering Structures 103 (2015) 239–259

lcor = 80 mm. By changing the material coefficient of variation COVft, the fluctuation of ft around its mean value was slightly modified. The spatial variation of ft was not affected by COVft.

7. Deterministic simulation results The deterministic simulation results were described in detail by Syroka-Korol et al. [10]. For the sake of clarity, some of them are also shown in this paper. The evolution of the normalized flexural tensile stress 1.5FL/( ft0 D2t) (F – vertical force, L – beam length, D – beam height, t-beam thickness, ft – tensile strength) versus the normalized deflection u/D during three-point bending for four different beam sizes with the constant tensile strength ft is shown in Fig. 11 (lc = 5 mm and m = 2). The calculated maximum deterministic vertical forces were: Fmax = 3.83 kN (D = 8 cm), Fmax = 6.75 kN

(D = 16 cm), Fmax = 12.57 kN (D = 32 cm) and Fmax = 66.18 kN (D = 192 cm), respectively. The normalized nominal (flexural) strength rN/ft = 1.5FmaxL/(D2t ft0 ) varied between 1.1 (D = 192 cm) and 1.5 (D = 8 cm). The strength and ductility strongly increased with decreasing beam height. For the large and very large-size beam, the snap-back behavior occurred (strength decrease with decreasing deformation), in particular very pronounced for the very large-size beam. The two lines fr/ft0 = 1 and fr/ft0 = 1.5 were also marked in Fig. 11. They constitute the lower and upper bounds to the normalized maximum strength (assuming the classical beam theory with a linear axial strain distribution in the cross-section and the limit analysis design for perfectly plastic symmetric materials). The height of the localized zone hloc at the peak load slightly grew from 24 mm (hloc/D = 0.3) up to 48 mm (hloc/D = 0.024) when the beam height increased from 80 mm up to 1920 mm (Figs. 12 and 13).

Fig. 16. Theoretical Gauss cdf compared to FE results for small-size beam (D = 80 mm) with autocorrelation length lcor: (a) 10 mm, (b) 20 mm, (c) 30 mm, (d) 40 mm, (e) 60 mm, (f) 80 mm, (g) 100 mm, (h) 120 mm and (i) 140 mm.

E. Syroka-Korol et al. / Engineering Structures 103 (2015) 239–259

The outcomes from deterministic simulations indicate a pronounced energetic size effect on the normalized flexural strength fr/ft0 = 1.5FmaxL/(D2tft0 ) (Fig. 14). Our numerical results were compared with Eq. (1). The additional deterministic calculations were performed for the beam sizes D = 2 cm, D = 4 cm and D = 384 cm. Eq. (1) includes 4 unknown parameters ( f1 r , Db, lp and r) which are to be determined. By taking the advantage of a mechanical definition of parameters described in [17], the number of unknown parameters may be reduced to 2. The flexural tensile strength fr(D) approaches for D ? 1 the uniaxial tensile strength ft0 . For rectangular beams subjected to bending, the ratio Mpl/Mel = gp = 3 and lp = rDb/(grp  1). Assuming r = 1.0, the following parameters were found to fit Eq. (1): f1 r = 3.6 MPa, Db = 45 mm, lp = 22.5 mm (Fig. 14). The agreement of FE results with Eq. (1) is very good (Fig. 14). Except of the nominal strength, the other mechanical parameter which strongly depends on the structure size is the post-critical stiffness modulus characterizing the load decrease related to the progressive failure deformation. The global softening modulus on the load–displacement diagram usually decreases with the

249

structure size, leading to the effect of a transition from the snap-through to snap-back, with the displacement control loss in a test. In order to determine the change of the softening rate, we introduced the new parameter – the softening angle b which is the angle measured between the horizontal line for the normalized peak strength fr/ft0 and the line linking the peak with the point on the softening curve with the increment of Dfr/ft0 = 0.2 (Fig. 15a). It increases in a parabolic form with the beam height D: b = 72° (small beam), b = 86° (medium beam), b = 97° (large beam) and b = 102° (very large beam) (Fig. 15b). The softening angle b larger than 90° denotes the snap-back behavior. The evolution of the softening angle b against D (Fig. 15b) resembles the change of the height of the localized zone at the peak strength against D (Fig. 12). The change of the softening rate can be also expressed by the softening angle c, i.e. the angle measured between the elastic part of the curve fr/ft0 and the line linking the peak with the point on the softening curve again with the increment of Dfr/ft0 = 0.2 (Fig. 15a). The softening angle c decreases in a parabolic form from 29.6° (small beam) down to 0.2° (very large beam) (Fig. 15b).

Fig. 17. Theoretical Gauss cdf compared to FE results for medium-size beam (D = 160 mm) with lcor: (a) 10 mm, (b) 20 mm, (c) 30 mm, (d) 40 mm, (e) 60 mm, (f) 80 mm, (g) 100 mm, (h) 120 mm and (i) 140 mm.

250

E. Syroka-Korol et al. / Engineering Structures 103 (2015) 239–259

where l is the expected value and s is the standard deviation,

8. Statistical simulation results


 m  x ; FðxÞ ¼ 1  exp  k

8.1. Effect of varying autocorrelation length lcor Tables 1–4 include results from the FE statistical analyses for the small (D = 80 mm, Table 1), medium (D = 160 mm, Table 2), large (D = 320 mm, Table 3) and very large (D = 1920 mm, Table 4) beam with the varying autocorrelation length (Eq. (8)). Due to the lack experimental data, the wide range of lcor was assumed (lcor = 2–140 mm). The expected values and standard deviations of the failure forces in Tables 1–4 were evaluated by calibrating the theoretical probability distribution functions based on the results from 12 FE simulations in each case (Figs. 16–19 with lcor = 10–140 mm). For the small, medium and large size beam, the Gauss cumulative distribution function (cdf) was adopted (Eq. (12)) whereas for the very large size beam, the Weibull cumulative distribution function (cdf) was used (Eq. (13))

FðxÞ ¼


 1 xl 1 þ erf pffiffiffi ; 2 s 2

ð12Þ

ð13Þ

where k is the scale parameter. The expected value and standard deviation for the Weibull probability distribution were estimated as follows




l ¼ kC 1 þ

 1 ; m


 2 2  l2 : s2 ¼ k C 1 þ m

ð14Þ

ð15Þ

The expected failure force lF for the small size beam (D = 80 mm) varied from 3.62 kN up to 3.85 kN (Fmax = 3.83 kN in the deterministic analysis) while the coefficient of variation COVF changed from 0.03 up to 0.10 (Fig. 20a). For the medium size beam (D = 160 mm), the lowest expected failure force lF was 6.24 kN (Fmax = 6.75 kN in the deterministic analysis) and the highest value reached 6.73 kN,

Fig. 18. Theoretical Gauss cdf compared to FE results for large-size beam (D = 320 mm) with lcor: (a) 10 mm, (b) 20 mm, (c) 30 mm, (d) 40 mm, (e) 60 mm, (f) 80 mm, (g) 100 mm, (h) 120 mm and (i) 140 mm.

E. Syroka-Korol et al. / Engineering Structures 103 (2015) 239–259

251

Fig. 19. Theoretical Weibull cdf compared to FE results for very large-size beam (D = 1920 mm) with lcor: (a) 10 mm, (b) 20 mm, (c) 30 mm, (d) 40 mm, (e) 60 mm, (f) 80 mm, (g) 100 mm, (h) 120 mm and (i) 140 mm.

the coefficient of variation COVF increased from 0.03 up to 0.11 (Fig. 20b). For the large size beam (D = 320 mm), the expected failure force lF varied between 10.88 kN and 12.24 kN (Fmax = 12.57 kN in the deterministic analysis) while the coefficient of variation COVF was in the range of 0.04–0.12 (Fig. 20c). The very large beam achieved the lowest expected failure force lF equal to 55 kN and the highest one equal 61.28 kN (Fmax = 66.18 kN in the deterministic analysis), the coefficient of variation COVF varied between 0.03 and 0.11 (Fig. 20d). The Weibull modulus m calibrated for the very large size beam changed from 11 up to 44. In some cases the theoretical cdf did not fit the numerical results satisfactorily (compare Figs. 16a, 18c and 19b, d, f), the tails of cdfs tended to diverge from the FE results. Thus, the calibrated values of the standard deviation or Weibull modulus m might be affected by some error. The expected failure force lF and its coefficient of variation strongly depended on the autocorrelation length lcor (Fig. 20) and the same general trend could be observed for all beam sizes. First the expected failure force lF (Fig. 20A) decreased and further increased

together with increasing autocorrelation length lcor. In the case of the small size beam, a plateau was observed where the failure force was close to the deterministic value Fmax = 3.83 kN. The lowest expected failure force was obtained with lcor equal to 10, 20, 30 and 40 mm for the small, medium, large and very large-size beam, respectively. The coefficient of variation COVF (Fig. 20B) gradually increased with increasing lcor and next stabilized (small, medium and large size beam) at the level of about COVF = 0.10 which was slightly below of the assumed material coefficient of variation COVft = 0.12. The single localized zone in the beams was always nonsymmetric and curved [10]. The localized zone width wloc was independent of the beam size D and autocorrelation length lcor and was similar as the deterministic value of wloc = 15 mm. The height hloc at the peak load fluctuated for the given autocorrelation length lcor (the relationship between the average height hloc and lcor was not found). The mean localized zone height hloc for the varying lcor was 24 mm, 35 mm, 45 mm and 71 mm for the small, medium,

252

E. Syroka-Korol et al. / Engineering Structures 103 (2015) 239–259

Fig. 20. Calculated expected failure force lF (A) and its coefficient of variation COVF (B) from FE statistical simulations depending on autocorrelation length lcor: (a) small (D = 80 mm), (b) medium (D = 160 mm), (c) large (D = 320 mm) and (d) very large beam (D = 1920 mm).

253

E. Syroka-Korol et al. / Engineering Structures 103 (2015) 239–259

Fig. 23. Calibrated SEL Type I (Eq. (2)) by Bazant compared with results from statistical FE analyses with varying autocorrelation length lcor.

Table 5 Parameters used to calibrated energetic–statistical SEL Type I (Eq. (2)) (m - Weibull modulus, lcor - autocorrelation length, Ls - statistical scaling length, W - geometry parameter).

Fig. 21. FE statistical results: (a) expected flexural tensile strength lfr and (b) its coefficient of variation COVfr versus normalized autocorrelation length lcor/hloc.

m

lcor (mm)

W

Neq (D = 192 m)

Ls (mm)

24 24 10

40 80 40

0.0024 0.0024 0.012

5.53E4 1.382E4 2.856E5

817 1635 359

large and very large beam, respectively (Fig. 12). The expected nominal flexural tensile strength lfr strongly decreased with increasing specimen size D (Fig. 21). The lowest lfr was as follows: 5.09 MPa (reduction by 6% referring to the deterministic fr = 5.39 MPa), 4.39 MPa (reduction by 8% referring to fr = 4.75 MPa), 3.82 MPa (reduction by 8% referring to fr = 4.44 MPa) and 3.22 MPa (reduction by 17% referring to fr = 3.88 MPa) for the small, medium, large and very large beam, respectively (with lcor = 10–40 mm). The strongest reduction of the nominal strength

Fig. 22. Random distribution of local tensile strength ft for small-size beam (D = 80 mm) for lcor (a) 3 mm, (b) 10 mm and (c) 40 mm (with marked area wloc  hloc).

254

E. Syroka-Korol et al. / Engineering Structures 103 (2015) 239–259

with increasing beam size was computed with lcor = 40 mm while the weakest reduction was with lcor = 10 mm. Fig. 21 compares the calculated expected nominal flexural tensile strength lfr and coefficient of variation COVfr depending on the ratio lcor/hloc ratio for all beam sizes. The nominal strength lfr decreases with increasing lcor/hloc ratio up to about 0.5 and further increases. The coefficient of variation COVfr is more or less stable when lcor is higher than hloc. This phenomena can be explained follows. As it was demonstrated in Section 4, the failure load (or the

nominal strength) in statistical analyses depends upon the localized zone size at the peak load wloc  hloc, its position with respect to the support and the average local tensile strength ft inside the localized zone. When the spatial autocorrelation is too weak (lcor  lcor << wloc  hloc) to affect the entire area wloc  hloc, the local tensile strength ft inside the localized zone fluctuates and the average value of ft in the area wloc  hloc is more concentrated around the assumed theoretical mean local value lft = ft0 (i.e. the mean ft has the decreasing variance with decreasing lcor), thus the mean

Fig. 24. Normalized flexural tensile stress–deflection curves fr/ft0 = f(u/D) from statistical simulations for different autocorrelation length lcor: (A) lcor = 40 mm, (B) lcor = 80 mm and (C) lcor = 120 mm corresponding to the lowest, average and highest flexural tensile strength: (a) small (D = 80 mm), (b) medium (D = 160 mm), (c) large (D = 320 mm) and (d) very large size beam (D = 1920 mm).

E. Syroka-Korol et al. / Engineering Structures 103 (2015) 239–259

255

average ft inside the area wloc  hloc is more scattered around the assumed local value lft = ft0 (i.e. the average ft with increasing lcor has the increasing variance) that is reflected by the increasing coefficient of variation of the failure force. The increase of lcor induces the more uniform distribution of ft and consequently the lower number of weak sub-areas (Fig. 22). Hence, the mean failure load increases. In our calculations the width of the localized zone wloc = 5 mm was constant (wloc depends on the characteristic length lc only). Although we assumed one characteristic length lc = 5 mm in calculations, the influence of the varying wloc = f(lc) may be intuitively explained. With increasing localized zone width wloc and constant lcor, the average local material strength inside the area of wloc  hloc becomes closer to the assumed theoretical mean local value lft = ft0 (similar effect to the decrease of lcor for the constant area wloc  hloc). Hence, the decrease of the ratio lcor/wloc would increase the mean failure force and decrease its coefficient of variation. According to the Weibull theory, the modulus m is strictly a material property i.e. it is independent of the beam size and shape. However, our results showed that the modulus m depended on the beam size and autocorrelation length (i.e. by varying COVfr). The varying modulus m indicates the beam statistical strength to be affected by a localized zone size and the local Weibull theory cannot be applied for analyzed intermediate beam sizes. The combined energetic–statistical SEL (Eq. (2)) is supported by the constant Weibull modulus m. Thus, in order to properly calibrate the energetic–statistical SEL formula (Eq. (2)), we used a theoretical approach described by Bazant et al. [23] and Vorechovsky and Sadilek [45,46]. The procedure estimates the unknown parameter Ls in Eq. (2) for the given m. The parameter Ls which corresponds to the intersection of the deterministic horizontal asymptote 1 n/m rN(D) = f1 r and the Weibull statistical asymptote rN(D) = fr D may be determined as (Bazant et al. [23]):


Ls ¼ Dst

Fig. 25. The softening angle b of Fig. 14b from deterministic (red triangles) and statistical analyses (black circles) for different autocorrelation length lcor: (a) lcor = 40 mm, (b) lcor = 80 mm and (c) lcor = 120 mm and beam height D = 80–1920 mm. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

failure force increases. This also explains the decreasing coefficient of variation of the failure load. On the other hand, when the spatial autocorrelation becomes stronger (lcor  lcor >> wloc  hloc), the

rst 1

fr

mn

ð16Þ

where rst is the mean strength of very large structure sizes Dst calculated from the Weibull integral [23]. We used the generally accepted m-value for concrete, namely m = 24. It is close to the average value of m for the largest beam calculated based on the scatter of fr (the average value m  23 for all lcor, Table 4). For comparative purposes we used also the value m = 10 resulting directly from the assumed material coefficient of variation COVft. Fig. 23 compares our statistical results with the calibrated energetic– statistical SEL curves (Eq. (2)) with m = 24 and m = 10 (lcor = 40 mm and lcor = 80 mm) (Table 5). The general tendency of the mean nominal strength with increasing beam size is well reproduced by the calibrated SEL curve with m = 24. The increasing value of lcor causes an increase of the calculated statistical scale length Ls which affects the mean nominal strength predicted by SEL (Eq. (2)). With respect to extremely large sizes, the reduction of the mean nominal strength with increasing beam size is independent of lcor since the slope of the Weibull asymptote expressed in terms of m is constant. On the other hand the value of m = 10 indicates an unrealistically low value of the mean nominal strength for analyzed intermediate beam sizes. Fig. 14 presents the deterministic and mean statistical results with lcor = 40 mm and COVft = 0.12 as compared to the pure energetic (Eq. (1)) and coupled energetic–statistical (Eq. (2)) rules. The mean statistical and deterministic flexural tensile strength predicted by SELs are similar for small sizes and bounded by the plasticity limit. With increasing beam size, the effect of the material randomness becomes stronger and the mean statistical flexural strength diverge from the deterministic one. The statistical outcomes indicate a further reduction of the nominal flexural strength

256

E. Syroka-Korol et al. / Engineering Structures 103 (2015) 239–259

Table 6 The softening angle b (°) from statistical analyses corresponding to lowest, average and highest flexural tensile strength fr for different correlation length lcor and beam height D. D (mm)

lcor = 40 mm

80 160 320 1920

65.94 74.61 84.19 101.77

lcor = 80 mm 67.59 80.60 90.98 102.30

77.49 82.15 93.34 102.37

66.64 73.36 88.50 101.26

with increasing beam size while the deterministic ones reach their lower limit. The mean statistical flexural strength of the largest beam D = 1920 mm was reduced by 17% as compared to the deterministic calculations. The statistical softening angle b was determined for the different curves fr/ft0 = f(u/D) corresponding to the lowest, average and highest flexural tensile strength (lcor = 40–120 mm) (Fig. 24). The statistical angle may be smaller by 12° than the deterministic value for D = 80–320 mm and is almost same for D = 1920 mm (Fig. 25
may be smaller and Table 6). In turn, the mean statistical angle b maximum by 7° than the deterministic value (Fig. 26 and Table 7).
increases with increasing Generally the mean statistical angle b autocorrelation length lcor and approaches a deterministic value of b, except of the very large size beam where the both values are almost the same. The mean statistical angle with lcor = 40–120 mm
= 70.34–72.43° (deterministic for the small beam varies between b
= 79.12– value b = 73.41°), for the medium beam varies between b 85.55° (deterministic value b = 85.51°), for the large beam
= 89.50–95.30° (deterministic value b = 96.53°) and for the very b
= 102.42–102.14° (deterministic value b = 102.42°). large beam b
with respect to the deterministic The smallest reduction of mean b

Fig. 26. The mean softening angle b of Fig. 14b from statistical analyses for different autocorrelation length lcor = 40–120 mm and beam height D = 80–1920 mm as compared to deterministic values.

lcor = 120 mm 74.54 75.13 95.39 102.18

73.41 88.25 97.92 102.35

70.08 82.61 92.88 101.71

72.47 85.17 94.87 102.33

74.74 88.85 98.15 102.38

value takes place for the autocorrelation length of lcor = 120 mm (D = 80–320 mm). 8.2. Effect of varying material coefficient of variation Fig. 27 presents the statistical analysis results with the varying material coefficient of variation COVft from 0.12 up to 0.20 with lcor = 80 mm. The calibrated cdfs (details in Tables 8–11) fit properly the FE results in all cases except of the medium-size beam with COVft = 0.16 where the standard deviation is underestimated. For all beam sizes the expected failure force decreases while its coefficient of variations COVF increases with increasing COVft. This behavior can be simply explained by a growing standard deviation of the average ft inside the area of the localized zone wloc  hloc at the peak load (Section 5). The larger the beam, the stronger is the effect of the material coefficient of variation COVft on lF and COVF. This phenomena is caused by an increasing number of weak sub-areas where a localized zone may initiate. The expected failure force lF for the small size beam (D = 80 mm) decreases from 3.82 kN down to 3.77 kN while the coefficient of variation COVF increases from 0.097 up to 0.180 kN. For the medium-size beam (D = 160 mm), lF decreases from 6.54 kN down to 6.36 kN and COVF increases from 0.092 kN up to 0.162 kN. the large size beam has lF between 11.60 kN and 11.02 kN while COVF changes from 0.10 (the Weibull modulus m = 20.2) up to 0.153 (Weibull modulus m = 6.7). For the very large size beam lF varies between 54.78 kN and 44.13 kN and COVF rises from 0.061 kN up to 0.176. The localized zone width wloc and height hloc are not affected by COVft. The expected flexural tensile strength lfr strongly decreases with increasing specimen size and increasing material coefficient of variation COVft. Fig. 28 compares the calibrated energetic–statistical SEL Type I curves (Eq. (2) with Ls calculated as described in Section 8.1 and m assumed from COVfr for the very large beam D = 1920 mm, Table 11) with the calculated lfr for the varying beam size D. The strongest size effect was obtained with the highest material COVft = 0.20. The calculated expected flexural tensile strength decreased from 5.37 MPa (D = 80 mm) down to 3.20 MPa (D = 1920 mm) and from 5.30 MPa (D = 80 mm) down to 2.58 MPa (D = 1920 mm) for COVft = 0.12 and COVft = 0.20, respectively (Fig. 28). The reduction of the expected lfr with COVft = 0.20 for the very large size beam was 33% as compared to deterministic simulations. The calibrated SEL curve captured well numerical results lfr with COVft = 0.12. With the higher COVft, the difference between numerical lfr and predicted by the calibrated

Table 7
(°) from Table 6 (statistical analyses) for different correlation length lcor and beam height D compared with deterministic value. The mean softening angle b D (mm)

Deterministic value of b (°)

Deterministic value of c (°)


(°) b lcor = 40 mm


(°) b lcor = 80 mm


(°) b lcor = 120 mm

80 160 320 1920

70.34 79.12 89.50 102.14

29.62 16.38 5.59 0.20

71.53 78.92 93.94 101.93

72.43 85.55 95.30 102.14

70.34 79.12 89.50 102.14

257

E. Syroka-Korol et al. / Engineering Structures 103 (2015) 239–259

Fig. 27. Results from statistical FE analyses with varying material coefficient of variation COVft = 0.12, 0.16 and 0.20 compared with calibrated theoretical cdfs: (a) small (D = 80 mm), (b) medium (D = 160 mm), (c) large (D = 320 mm) and very large beam (D = 1920 mm) (cases ‘a’–‘c’) Gauss cdf and case (d) Weibull cdf.

Table 8 Results from FE simulations with varying COVft (small-size beam D = 80 mm). COVft

lF (kN)

sF (kN)

COVF

0.12 0.16 0.20

3.82 3.81 3.77

0.37 0.51 0.68

0.097 0.133 0.180

Table 11 Results from FE simulations with varying COVft (very large size beam D = 1920 mm, m – Weibull modulus). COVft

lF (kN)

sF (kN)

COVF

m

0.12 0.16 0.20

54.78 49.17 44.13

3.358 4.247 7.747

0.061 0.086 0.176

20.22 14.17 6.68

Table 9 Results from FE simulations with varying COVft (medium-size beam D = 160 mm). COVft

lF (kN)

sF (kN)

COVF

0.12 0.16 0.20

6.54 6.62 6.36

0.60 0.69 1.03

0.092 0.104 0.162

Table 10 Results from FE simulations with varying COVft (large-size beam D = 320 mm). COVft

lF (kN)

sF (kN)

COVF

0.12 0.16 0.20

11.60 11.32 11.02

1.16 1.50 1.69

0.100 0.133 0.153

SEL became stronger. However, the size effect, i.e. the reduction of the strength lfr with increasing size D indicated by the slope of the SEL curve, is in satisfactory agreement with numerical results. Our numerical results are in agreement with the similar ones obtained by Grassl and Bazant [14] with respect to the effect of varying lcor on the mean strength (which becomes negligible when lcor is small and pronounced when lcor is large). They are also in agreement with respect to the strong influence of the varying material coefficient of variation (the higher the material coefficient of variation, the lower is the mean nominal strength and higher coefficient of variation of the strength). The highest reduction of the mean nominal strength was 25% between the beam D = 40 mm and D = 800 mm [14].

258

E. Syroka-Korol et al. / Engineering Structures 103 (2015) 239–259

Fig. 28. Calibrated SEL Type I (Eq. (2)) by Bazant compared with results from statistical FE analysis with varying material coefficient of variation COVft.

9. Conclusions  The decrease of the nominal strength, ductility and the ratio of the height of the localized zone at the peak load to the beam height with increasing beam size is significant in deterministic calculations (the nominal strength decreased by 27% when the beam size was increased from 80 mm up to 1920 mm). Their further decline with increasing beam size is caused by a random distribution of the local tensile strength (the mean statistical strength reduction by 17% for the largest beam D = 1920 ). The beam flexural strength is affected by the size of the localized zone wloc  hloc at the peak load, the average local tensile strength inside the area wloc  hloc and the corresponding elastic normal stress due to bending. The position of the localized zone in statistical simulations is obviously non-symmetric and its shape can be a more or less curved. The width of the localized zone wloc is independent of the beam size D and material randomness. The localized zone height hloc at the peak load is affected by both the beam size and material randomness. The average normalized hloc/D from statistical simulations decreases from 0.30 down to 0.037 for the small and very large-size beam.  The expected flexural tensile strength reduction with respect to the deterministic analyses becomes stronger when the beam size increases. The expected flexural tensile strength is affected by both the autocorrelation length and material coefficient of variation. With increasing autocorrelation length, the expected flexural tensile strength first decreases until the autocorrelation length is equal to the half of the localized zone height at the peak load and next increases. The strongest reduction of the flexural tensile strength with increasing beam size (by about 40% for the very large size beam with respect to the small beam) occurs with the autocorrelation length lcor = 40 mm. The coefficient of variation of the flexural tensile strength gradually increases with increasing autocorrelation length and next stabilizes when the autocorrelation length is larger than the localized zone height at the peak load. The varying ratio lcor/wloc has the similar influence on the flexural tensile strength and its variation as the varying ratio lcor/hloc. The highest nominal strength variation for the largest beam COVfr = 0.11 was reached for the highest correlation length lcor = 140 mm which was close to the assumed material coefficient of variation COVft = 0.12. Hence, the Weibull modulus m reflecting the nominal strength variation of beams depended on lcor. The lowest Weibull

modulus obtained for the largest beam was m = 11 (with lcor = 140 mm). The SEL curve calibrated with m = 10 (based on the assumed material COVft = 0.12) strongly underestimated the mean strength of the analyzed beams of intermediate sizes. The nominal strength reduction of beams caused by material randomness may roughly be estimated with m = 24.  The increasing material coefficient of variation causes a further decline of the expected flexural tensile strength. The resulting coefficient of variation of the flexural tensile strength decreases with decreasing material coefficient of variation. The reduction of the expected flexural tensile strength is 20% when comparing the very large size beam with COVft = 0.12 and COVft = 0.20. The energetic–statistical SEL function calibrated with the Weibull modulus m resulting from the nominal strength variation COVfr of the largest beam may be used to approximate the mean strength reduction with increasing beam size. The SEL with the Weibull modulus m calculated from the assumed material coefficient of variation unrealistically provides the low estimation of the mean strength of intermediate beams sizes.  The beam brittleness increases with increasing beam height. For the beam height D P 320 mm the snap-back phenomenon appears. The mean statistical softening angle b may be smaller by 3° (small beam), 6° (medium beam) and 7° (large beam) than the deterministic value. This maximum reduction happens again with the autocorrelation length lcor = 40 mm. The mean statistical softening angle generally increases with increasing autocorrelation length lcor and approaches the deterministic value. The evolution of the softening angle b is similar to the change of the height of the localized zone (versus the beam height).

Acknowledgments The research work has been carried out within the project: ‘‘Experimental and numerical analysis of coupled deterministic–statis tical size effect in brittle materials” financed by the National Research Centre NCN (UMO-2013/09/B/ST8/03598) and the scholarship for young doctors financed from the project: ‘‘The Centre for Advanced Studies – the development of interdisciplinary doctoral studies at the Gdansk University of Technology in the key areas of the Europe 2020 Strategy” (POKL04.03.00-00-238/12). References [1] Bazant ZP. Size effect in blunt fracture: concrete, rock, metal. J Eng Mech ASCE 1984;110:518–35. [2] Carpinteri A. Decrease of apparent tensile and bending strength with specimen size: two different explanations based on fracture mechanics. Int J Solids Struct 1989;25(4):407–29. [3] Bazant ZP, Planas J. Fracture and size effect in concrete and other quasi-brittle materials. CRC Press LLC; 1998. p. 1–640. [4] Bazant ZP. Probability distribution of energetic–statistical size effect in quasibrittle fracture. Prob Eng Mech 2004;19:307–19. [5] Weibull W. A statistical theory of the strength of materials. J Appl Mech 1951;18(9):293–7. [6] Bazant ZP, Novak D. Proposal for standard test of modulus of rupture of concrete with its size dependence. ACI Mater J 2001;98(1):79–87. [7] Tejchman J, Bobin´ski J. Continuous and discontinuous modelling of fracture in concrete using FEM. Berlin-Heidelberg: Springer; 2013. [8] Bobinski J, Tejchman J, Górski J. Notched concrete beams under bendingcalculations of size effects within stochastic elasto-plasticity within non-local softening. Arch Mech 2009;61(3–4):1–25. [9] Syroka-Korol E. Theoretical and experimental study on size effect in concrete beams reinforced with steel and basalt bars. PhD thesis. Gdansk (PL): Gdansk University of Technology, Poland; 2012. [10] Syroka-Korol E, Tejchman J, Mróz Z. FE calculations of a deterministic and statistical size effect in concrete under bending within stochastic elastoplasticity and non-local softening. Eng Struct 2013;48:205–19. [11] Phoon KK, Huang SP, Quek ST. Simulation of second-order processes using Karhunen–Loeve expansion. Comput Struct 2002;80:1049–60.

E. Syroka-Korol et al. / Engineering Structures 103 (2015) 239–259 [12] Phoon KK, Huang SP, Quek ST. Implementation of Karhunen–Loeve expansion for simulation using a wavelet-Galerkin scheme. Prob Eng Mech 2002;17:293–303. [13] Hurtado JE, Barbat AH. Monte Carlo techniques in computational stochastic mechanics. Arch Comput Method Eng 1998;5(1):3–30. [14] Grassl P, Bazant ZP. Random lattice-particle simulation of statistical size effect in quasi-brittle structures failing at crack initiation. J Eng Mech 2009;135 (2):85–92. [15] Vorechovsky M, Chudoba R. Stochastic modeling of multi-filament yarns: II. Random properties over the length and size effect. Int J Solids Struct 2006;43 (3–4):435–58. [16] Bielewicz E, Górski J. Shell with random geometric imperfections – simulationbased approach. Int J Non-linear Mech 2002;37(4–5):777–84. [17] Bazant ZP, Novak D. Probabilistic nonlocal theory for quasibrittle fracture initiation and size effect. II: application. J Eng Mech 2000;126(2):175–85. [18] Carmeliet J, Hens H. Probabilistic nonlocal damage model for continua with random field properties. J Eng Mech 1994;120(10):2013–27. [19] Carmeliet J, de Borst R. Stochastic approaches for damage evolution in quasibrittle materials. Prob Mater 1994;269:491–503. [20] Gutiérrez MA, de Borst R. Deterministic and stochastic analysis of size effects and damage evaluation in quasibrittle materials. Arch Appl Mech 1999;69:655–76. [21] Frantziskonis GN. Stochastic modeling of hetereogeneous materials – a process for the analysis and evaluation of alternative formulations. Mech Mater 1998;27:165–75. [22] Vorechovsky M. Interplay of size effects in concrete specimens under tension studied via computational stochastic fracture mechanics. Int J Solids Struct 2007;44:2715–31. [23] Bazant ZP, Vorechovsky M, Novak D. Asymptotic prediction of energetic– statistical size effect from deterministic finite-element solutions. J Eng Mech ASCE 2007;133(2):153–62. [24] Yang Z, Xu XF. A heterogeneous cohesive model for quasi-brittle materials considering spatially varying random fracture properties. Comput Methods Appl Mech Eng 2008;197(45–48):4027–39. [25] Elias J, Vorechovsky M, Skocek J, Bazant ZP. Stochastic discrete meso-scale simulations of concrete fracture: comparison to experimental data. Eng Fract Mech 2015;135:1–16. [26] Majewski T, Bobin´ski J, Tejchman J. FE-analysis of failure behaviour of reinforced concrete columns under eccentric compression. Eng Struct 2008;30:300–17. [27] Syroka E, Bobin´ski J, Tejchman J. FE analysis of reinforced concrete corbels with enhanced continuum models. Finite Elem Anal Des 2011;47:1066–78.

259

[28] Hordijk DA. Local approach to fatigue of concrete. PhD dissertation. Delft (NL): Delft University of Technology; 1991. [29] Pijaudier-Cabot G, Bazant ZP. Nonlocal damage theory. ASCE J Eng Mech 1987;113:1512–33. [30] Bazant ZP, Jirasek M. Nonlocal integral formulations of plasticity and damage: survey of progress. J Eng Mech 2002;128:1119–49. [31] Bobinski J, Tejchman J. Numerical simulations of localization of deformation in quasi-brittle materials within non-local softening plasticity. Comput Concr 2004;4:433–55. [32] Polizzotto C, Borino B, Fuschi P. A thermodynamic consistent formulation of nonlocal and gradient plasticity. Mech Res Commun 1998;25:75–82. [33] Brinkgreve RBJ. Geomaterial models and numerical analysis of softening (Ph.D. thesis). Delft University of Technology; 1994. [34] Le Bellego C, Dube JF, Pijaudier-Cabot G, Gerard B. Calibration of nonlocal damage model from size effect tests. Eur J Mech A/Solids 2003;22:33–46. [35] Polizzotto C. Remarks on some aspects of non-local theories in solid mechanics. In: Proceedings of the 6th national congress SIMAI. Chia Laguna, Italy, CD-ROM; 2002. [36] Jirásek M, Rolshoven S, Grassl P. Size effect on fracture energy induced by nonlocality. Int J Numer Anal Methods Geomech 2004;28(7–8):653–70. [37] ABAQUS. Theory manual, version 5.8. Hibbit, Karlsson & Sorensen Inc.; 1998. [38] Pramono E. Numerical simulations of distributed and localized failure in concrete (PhD thesis). University of Colorado-Boulder; 1988. [39] Jirasek M, Bazant ZP. Inelastic analysis of structures. John Wiley & Sons; 2001. p. 758. [40] Syroka-Korol E, Tejchman J, Mróz Z. FE analysis of size effects in reinforced concrete beams without shear reinforcement based on stochastic elastoplasticity with non-local softening. Finite Elem Anal Des 2014;88:25–41. [41] Ghanem RG, Spanos PD. Stochastic finite elements – a spectral approach. New York, US: Springer Verlag; 1991. [42] Karhunen K. Über lineare Methoden in der Wahrscheinlichkeitsrechnung. Ann Acad Sci Fenn Ser A I Math Phys 1947;37:1–79. [43] Loève M. Functions aleatoires du second ordre. Supplement to P. Levy, Processus Stochastic et Mouvment Brownien, Paris, Gauthier Villares; 1948. [44] Neyman J. On the two different aspects of the representative method: the method of stratified sampling and the method of purposive selection. J R Stat Soc 1934;97:558–625. [45] Skarzynski L, Syroka E, Tejchman J. Measurements and calculations of the width of the fracture process zones on the surface of notched concrete beams. Strain 2011;7:319–22. [46] Vorechovsky M, Sadilek V. Computational modeling of size effects in concrete specimens under uniaxial tension. Int J Fract 2008;154(1–2):27–49.