Computational modeling of fracture in concrete using a meshfree meso-macro-multiscale method

Computational Materials Science 69 (2013) 204–215

Contents lists available at SciVerse ScienceDirect

Computational Materials Science journal homepage: www.elsevier.com/locate/commatsci

Computational modeling of fracture in concrete using a meshfree meso-macro-multiscale method A. Ghosh ⇑, P. Chaudhuri IIT Madras, Department of Engineering, India

a r t i c l e

i n f o

Article history: Received 18 May 2012 Received in revised form 4 November 2012 Accepted 14 November 2012 Available online 8 January 2013 Keywords: Applied mechanics Fracture mechanics Damage Concrete

a b s t r a c t The focus of this manuscript is on multiscale modeling of material failure of concrete. The concrete material far away from the fracture process zone is modeled by a macroscopic homogenized material model while a non-linear mesoscopic constitutive model is used in the fracture process zone. On the mesoscopic scale, the constituents of the composites material (concrete), i.e. the cement matrix and the aggregates, are modeled. Fracture occurs only in the cement matrix and is modeled by a cohesive zone model (CZM) in the context of an enriched partition of unity meshfree method. We also test the influence of the interface model between the cement and the aggregates. Numerical results demonstrate the importance of the interface model. The method is validated by comparison to experimental data for different failure modes. Ó 2012 Elsevier B.V. All rights reserved.

1. Introduction Due to its low cost and recyclability, concrete is one of the most popular materials in Civil Engineering. Concrete can be considered as composite consisting of aggregates and a porous cement matrix with fluid-filled cavities. The macroscopic behavior of concrete depends on its structure at the mesoscopic scale. For example, the interface behavior between the aggregates and the cement matrix on the mesoscopic scale determine the macroscopic material response. Several numerical methods were developed for the complex composite material concrete. Most of the models are based on homogenized constitutive models on macroscopic scale neglecting the complex mesoscopic structure of concrete materials [1–11]. Those models are well suited to capture the global response of a structure. However, they are not suitable to predict material properties. One of the key phenomena in concrete materials is cracking. Another characteristics of concrete is its high compressive strength compared to its low tensile strength. Modeling the cracking process in concrete is particularly challenging for several reasons. The big size of the process zone [12] prohibits the use of LEFM (Linear Elastic Fracture Mechanics). Moreover, a huge number of cracks commonly occur even under static loading. Smeared crack models [4,3,13,14] were one of the first attempts to model fracture in

⇑ Corresponding author. Tel.: +91 0361 258 2609 2621; fax: +91 0361 258 2609 2633. E-mail address: [email protected] (A. Ghosh). 0927-0256/$ - see front matter Ó 2012 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.commatsci.2012.11.025

concrete. The boom of PU1-methods [15–26] in recent years appear particularly promising for modeling discrete fracture, see also [27– 32]. While those methods have been applied to model fracture on the macroscopic scale, they never have been used for fracture on the mesoscopic scale. Modeling concrete on the mesoscopic scale is computationally demanding [33,34] and restricts the size of the model. On the other hand, a mesoscopic model is needed for the prediction of material properties. In this manuscript, we present a concurrent multiscale methods for fracture in concrete. The constituents of concrete are modeled around the fracture process zone while a macroscopic description is used elsewhere. We employ the XEFGM (extended element-free Galerkin method) [35,36] in the cement matrix. XEFGM is a meshless Galerkin method that has been successfully applied to complex macroscopic material failure in concrete and reinforced concrete under static and dynamic loading [37–42]. It was also extended to other applications in [43–45]. Besides cracking in the cement matrix, we use two interface models: (1) Rigid model and (2) A CZM taking into account the slip between the cement paste and the aggregates. We validate our method through comparison to experimental results. Moreover, we use our multiscale method to extract material parameters such as fracture energy and tensile strength and show the influence of the interface model. In the next section, we describe the XEFGM. The constitutive model and CZM for cracks is explained subsequently. Section 4 explains the interface model used. The discrete equations are derived

1

Partition of unity.

A. Ghosh, P. Chaudhuri / Computational Materials Science 69 (2013) 204–215

205

in Section 5 before we present numerical results in Section 7. The manuscript ends with conclusions in Section 8. 2. Extended element-free galerkin method The XEFGM [35,36] is based on decomposition of the displacement field into a standard part and into an enriched part:

uh ðXÞ ¼ uhS ðXÞ þ uhE ðXÞ

ð1Þ

The first term on the RHS of Eq. (1) is the standard part of the XEFG-formulation while the second term on the RHS, the enriched part, accounts for the discontinuous displacement field. The standard part is the EFG [15] approximation based, that is based on moving least squares (MLSs) shape functions [46]. The MLS approximation is given by

uhS ðXÞ ¼

m X

pI ðXÞ aI ðXÞ ¼ pT ðXÞ aðXÞ

ð2Þ

I¼1

where the vector pT(X) commonly contains polynomial functions and m indicates the number of functions in p; the vector a(X) contains the unknown coefficients that depend on the material coordinates X. We use linear functions in pT(x) = [1 X Y]. The MLS shape functions are obtained by minimizing a discrete weighted L2 norm with respected to the unknowns a: n X J¼ ðpT ðXI Þ aðXÞ  uI Þ2 wðX  XI Þ

ð3Þ

I¼1

yielding the final EFG approximation

uðXÞ ¼ NðXÞ D

ð4Þ

with the shape function matrix N(X):

NðXÞ ¼ pT ðXÞ A1 ðXÞ CðXÞ

ð5Þ

and

AðXÞ ¼

n X wðX  XI Þ pðXI Þ pT ðXI Þ

Fig. 1. Modeling of crack by XEFG.

introduced as additional DOFs (degrees of freedoms) that need to be solved for. The enrichment functions w(fX) and /(X(h, r)) are defined as:

/ðXðh; rÞÞ ¼ r sin  wðfX Þ ¼

h 2

1 f ðXÞ > 0 0 f ðXÞ < 0

ð9Þ

Commonly, a local level set function fX is used to describe the topology of the crack surface; r in Eq. (9) is the distance of a material point to the crack tip and the angle h is explained in Fig. 1. More detailed description can be found in [36,39,18]. Several approaches were developed omitting the crack tip enrichment [50,19,51]. They are based on modifying the dilation parameter or Lagrange multiplier approach. The first approach is restricted to 2D problems while the second approach introduces other difficulties related to the Lagrange multiplier approach. On the other hand, integration the non-polynomial terms of the tip enrichment is avoided. We have decided to use a crack tip enrichment.

I¼1

3. Constitutive and cohesive crack model

CðXÞ ¼ ½wðX  X1 Þpðx1 Þ wðX  X2 Þpðx2 Þ; . . . ; wðX  Xn Þpðxn Þ D ¼ ½u1 u2 ; . . . ; un 

3.1. Constitutive model

ð6Þ where w(X  XI) denotes the kernel function centered at XI, uI and XI are the nodal values of the displacement field at node I and its position, respectively; n represents the number when the shape functions are unequal to zero. We employ the Gaussian kernel function:

( wðdI Þ ¼

2 2 eðdI =cÞ eðdmI =cÞ 2 1eðdmI =cÞ

dI 6 dmI

0

dI > dmI

ð7Þ

We consider linear elastic material behavior of the constituents. Concrete is modeled using the Lemaitre isotropic elastic damage model [52]. The model was developed for cement materials. It has the classical damage format:

r ¼ ð1  DÞC : 

ð10Þ

where C denotes the fourth-order elasticity tensor, r is the stress tensor,  the strain tensor and D (0 6 D 6 1) is a scalar damage var~: iable that depends on an effective strain 

where dmI is the domain of influence of node I, and c is a parameter used to control the dilation of the weight function. Other common weighting functions used in the EFG or related methods are the cubic B-spline, see e.g. [47,6], that are also used in a different context (not only meshfree methods), see e.g. [48,48,49]. The enriched or discontinuous approximation of the displacement field can be further decomposed into a step-enrichment uhES and a tip-enrichment uhET :

where I denotes the principal strain and w the step function, see Eq. (9).

uD Eh ðXÞ ¼ uhES þ uhET

3.2. Cohesive crack model

¼

nC X I

nT X NI ðXÞ wðfX Þ aI þ NI ðXÞ /ðXðh; rÞÞ bJ

ð8Þ

J

where nC denotes the number of step-enriched nodes and nT is the number of tip-enriched nodes; the parameters a and b are

~Þ ¼ 1  ð1  AÞ0 ~1  AeBð~0 Þ Dð

ð11Þ

with the effective strain

~ ¼

rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi X3 2 wðI Þ I¼1 I

ð12Þ

The CZM is derived from a potential U [53] defined as

U¼

tc0 expðaDb Þði4 þ a½i1  i4 Þ 2D

ð13Þ

206

A. Ghosh, P. Chaudhuri / Computational Materials Science 69 (2013) 204–215

gbeing the ‘tensile strength of the interface’, G is the interface

where

i1 ¼ I  ð½½u  ½½uÞ

ð14Þ

i4 ¼ ð½½u  ½½uÞ : n  n

and n is the vector orthogonal to the crack surface; [[u]] is the jump in the displacement field, t0c, a, b and a are material parameters. Taking derivatives of U with respect to [[u]] gives

energy (related to the fracture energy above), Du0 ¼ KgC with penalty stiffness KC and Du is the effective relative displacement between cement and aggregate material at the interface:

Du ¼

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi hDun i2 þ ðaDut Þ2

ð21Þ

un and ut being the relative discontinuous displacement field normal and tangential to the crack surface:

where un and ut are the normal relative displacements and tangential relative displacements; the parameter a controls the relationship between normal and tangential relative displacements and is set to 1 in our computations. The normal and tangential components are finally obtained by:

un ¼ ½½u  n n |fflfflfflffl{zfflfflfflffl}

t n ¼ rðDuÞ

tc ¼

tc0 expðaDb Þðun þ aut Þ D

ð15Þ

un

ð16Þ

½½u  un ut ¼ ½½u  et |fflfflfflffl{zfflfflfflffl} j½½u  un j ut

et denoting the unit vector of the associated crack segment/surface. The state variable D is obtained through the damage yield surface:

W ¼ j½½uj  D

Dun Du a2 Dut t n ¼ rðDuÞ Du

ð22Þ

In compression, the normal interface traction is computed by the penalty stiffness tn = KCDun. More details can be found in the excellent manuscript of [54].

ð17Þ

The following expressions are needed for the cohesive stiffness tensor:

@tc t c0 expðaDb ÞðaI þ ½1  an  nÞ ¼ @½½u D @tc t c0 ¼ ð1  aÞ expðaDb Þðun I þ n  ½½uÞ DI ¼ @un D @tc 1  abDb ¼ tc DII ¼ @D D2 D¼

ð18Þ

For efficiency reasons, the cohesive model is implemented in incremental form.

5. Weak form and discretization S T Consider the domain X0 = Xcs Xfs, Xfs Xcs = with Lipschitz boundary C as shown in Fig. 2. The indices fs and cs indicate the fine scale and the coarse scale. We use a continuum model at both scale. In the fine scale, the concrete material is modeled as heterogeneous material consisting of aggregates and cement paste while in the coarse scale, the concrete is considered as homogeneous material [55–57,11] ensuring computational efficiency. To ensure consistent material behavior [58], we use computational homogenization to obtain the material behavior of the homogenized concrete [59]. The weak form of the problem is written as

Z

3.3. Cracking criterion

dðÞ : r dX0 

X0

As suggested by [38], we use the loss-of-material-stability criterion. The material fails when the minimum eigenvalue of the acoustic tensor Q is smaller or equal to zero:

min eig v alðQ Þ 6 0; |fflfflffl{zfflfflffl}

Q ¼ n  ðC

tan

þ d  rÞ  n

ð19Þ

n

Ctan denoting the consistent material tangent tensor of the Lemaitre damage model. In this case, the vector n is the normal to the crack surface. We opt to control the crack length and introduce a crack segment of length DL = a DX, DX being the minimum distance between two particles in the vicinity of a propagating crack and a was set to 1.2. The distance DX was evaluated in a circle of radius r = bh where h is the size of the domain of influence of the particle closest to the crack tip and b = 1.5. We tested different values of b and a for different discretizations and the results are insensitive to those parameters.



Z

Z

dðuÞ  t0 dC0 þ

C0t

 Þ dC  dðkÞ  ðu  u

Cu

þ

nc Z X i¼1

Ci c

~ Þ  ~tic dC þ dðu

dðuÞ  b dX0

dðuÞ  k dC

Cu

~ Þ  ~ti ðu Ci

¼0

ð23Þ

where b denotes the body force. The terms in first two lines in Eq. (23) are the well-known terms found in the weak form for the EFG method which requires the use of Lagrange multipliers to enforce essential boundary conditions [15]. Here, k = [kx ky]T denotes the vector containing the Lagrange multiplier increments. The terms in the last line of Eq. (23) are associated to the cohesive crack

4. The interface model The interface model between the aggregates and the cement can be considered also as cohesive zone model as described above. We adopt a model developed originally by [54], particularly design for (mesoscopic) concrete materials. As will be seen subsequently, it is similar to the cohesive crack model described above. It is also R Du ~ Þ dDu ~ with based on a potential P ¼ 0 rðDu

rðDuÞ ¼

8 < K C Du : gexp

g½DuDu0 

eG



Du < Du0 else



X0

Z

Z

Z

ð20Þ Fig. 2. Coarse scale and fine scale domain.

A. Ghosh, P. Chaudhuri / Computational Materials Science 69 (2013) 204–215

and the energy dissipation when the cohesive crack opens (or slides) and to the cohesive tractions at the interface, respectively. The EFG approximation explained in the previous sections can be written more compactly in a vector–matrix form. The displacement u and the strain  are given by

uðXÞ ¼ NðXÞ U

ð24Þ

 ¼ BðXÞ U

where U = [Djajb]T is the vector containing the nodal parameters and N and B represent the shape function matrix and its derivative matrix, respectively:

N ¼ ½NS jNES jNET 

ð25Þ

B ¼ ½BS jBES jBET 

where the indices C, DS and DT indicate the continuous, discontinuous step-enriched and discontinuous tip-enriched part of the EFG-approximation. The Lagrange multiplier k to impose essential boundary conditions is interpolated from its nodal values by using Lagrange interpolation as

kðXÞ ¼ Nk ðr u ÞK

ð26Þ

where Nk is a Lagrange-interpolant matrix (based on Dirac delta functions) and ru denotes the arc length along the boundary Cu; K is a vector containing the nodal Lagrange multipliers of all the nodes on the boundary Cu. After substituting the test and trial functions into the weak form and linearization we obtain the equilibrium equation in incremental form after some algebra:

"

K G G

0

#(

DU DK

(

) ¼

DR b DR

) ð27Þ

50mm

100mm

F

255mm

255mm

Fig. 3. Three-point bending beam [64].

207

where

K¼

Z

BT CB dX þ

XZ 0

Z

ðNÞT ðD þ DI þ DII ÞT N dC

ð28Þ

Cc

G¼ NT Nk dC Cu Z Z DR ¼ NT Db dX þ NT Dt0 dC X0 C0 Z e ¼  dC DR NTk Du

ð29Þ ð30Þ ð31Þ

Cu

The integrals are evaluated by integration over back ground cells as explained e.g. in [60,61]. 6. Validation In this section, two validation examples are presented. Therefore, we apply our method in the context of single-macroscopic scale simulations. The results obtained are compared to experimental data. Validation is an important step for a numerical method [62,63]. The first example is a concrete specimen under mode I-fracture while the second example deals with (more complicated) mixed-mode fracture. 6.1. Three-point bending test The most commonly used configuration to investigate the mode I crack propagation in concrete is a notched beam subjected to three-point bending. In this study, the beam shown in Fig. 3 is used. The depth, length and thickness of the beam are 100 mm, 510 mm and 100 mm, respectively. The ratio between the notched depth and the beam depth is equal to 0.5. The crack mouth opening displacement (CMOD) is used as a controlling parameter of loading. We use discretizations from 386 nodes up to 18,423 nodes. At 2756 nodes, convergent results are obtained in load–deflection curve while crack paths agrees well for only 386 nodes. The material parameters are: Young’s modulus 20 GPa, Poisson ratio 0.2, tensile strength 2.4 MPa and fracture energy 113 N/m. The results are compared with the experimental results by [64]. Similar numerical studies were done e.g. by [65]. The obtained response is in good agreement with the experimental result. No discretization dependence occurs, Fig. 5. The displaced configuration of three-point bending beam for 386 nodes and 2756 nodes is shown in Fig. 4. (See Fig. 5).

Fig. 4. Displaced configuration and contour plot of principal tensile stress of three-point bending beam.

208

A. Ghosh, P. Chaudhuri / Computational Materials Science 69 (2013) 204–215

1.5

6.2. Four-point double-notched shear test

1.4

The second and last validation example is the four-point double-notched shear test. The geometry of the tested beam is shown in Fig. 6. The depth, length and thickness of the beam are 200 mm, 800 mm and 100 mm, respectively. In this problem, the crack mouth sliding displacement (CMSD) is used as a controlling parameter of loading. The material parameters are: Young’s modulus 27 GPa, Poisson ratio 0.1, tensile strength 2 MPa and fracture energy 100 N/m. Results with 1687 and 2967 nodes are presented. The obtained results are compared with experimental and numerical results by [66] as shown in Fig. 8. Good agreement between the obtained results can be observed. The crack path agrees very well, too. The displaced beam is illustrated in Fig. 7.

1.3 1.2

Kormeling 1983 8232 nodes 2756 nodes 4121 nodes

1.1

Load (kN)

1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2

7. Results

0.1 0

0

0.25

0.5

0.75

In this section, we apply our method to multi-scale modeling of concrete. The method can be regarded as concurrent multi-scale method where areas of interest are modeled on the mesoscopic scale while macroscopic scale is used elsewhere. The area of interest is around the expected crack path. On mesoscopic scale,

1

Deflection (mm) Fig. 5. Load–deflection response of the three-point bending beam [64].

F Rigid bar A 40mm 200mm

40mm 2x80mm 400mm

400mm

Fig. 6. Four-point double-notched shear beam [66].

Fig. 7. Displaced configuration of four-point double-notched shear beam at different load steps.

A. Ghosh, P. Chaudhuri / Computational Materials Science 69 (2013) 204–215

50 45

Boccaetal. (1991) 1687 nodes 2967 nodes 5434 nodes

40

LoadF (kN)

35 30 25 20 15 10 5 0

0

0.02

0.04

0.06

0.08

0.1

Load point displacement A (mm) Fig. 8. Load–displacement response of the four-point double-notched shear beam.

constituents of concrete are modeled, i.e. aggregates, cement paste and pores or micro-voids. We carried out random studies and varied the positions of aggregates and pores and show that the dissipated energy due to cohesive cracks is bounded (see Fig. 8). We use the interface models between the aggregates and cement paste described in Section 4. We assume that the crack can only propagate through the cement but not through the aggregates. This assumption seems valid for static applications while fractured aggregates were observed for high load rates only. 7.1. L-shaped panel Consider L-shaped concrete panel shown in Fig. 9. This experiment was done by [67] and the experimental crack patterns observed in a series of experiments are shown on the right hand side in Fig. 9. We first show results for a pure macroscopic simulation: On macroscopic scale, the material parameters are: Young’s modulus E = 20 GPa, Poisson’s ratio m = 0.2, tensile strength ft = 2.8 MPa and fracture energy Gf = 140 N/m. We used three

different numbers of nodes: 4500 nodes, 12,500 nodes and 34,000 nodes. The load deflection curve of these three macroscopic simulations in shown in Fig. 13a and the crack pattern for finest discretization is shown in Fig. 10. For other discretizations, crack pattern looks similar. We note that we get one sharp crack pattern and can capture the global behavior fairly well. Next, we model the area around the expected crack path on mesoscopic scale, consisting of cement, aggregates and voids. The material parameters of aggregates are [67]: Young’s modulus: E = 38 GPa and Poisson’s ration m = 0.2. The cement’s material parameters are given in accordance with the experiments: Young’s modulus E = 19 GPa, Poisson’s ratio m = 0.2, tensile strength ft = 2.6 MPa, and fracture energy Gf = 120 N/m. We note that material parameters of the constituents of concrete are more reliable due to the less heterogeneous structure. We tested different numbers of nodes starting from 80,000 nodes up to 420,000 nodes. Since the results were within the experimental scatter, we present the results of the finest discretization. We found that the spatial distribution of the aggregates has a more pronounced influence of the number of nodes for the models we tested. Of course, a certain refinement is needed to resolve the shape of the finest aggregates increasing computational expense. Nevertheless, we believe that mesoscopic modeling has the potential gaining knowledge on failure process. [67] specified the decomposition of concrete used in this experiment. The volume fraction of cement was 0.3 while 0.7 was contribution of aggregates. We first assume absence of voids that is idealization and later on test the effect of voids on the results. The decomposition of different aggregate sizes for this concrete is given by volume fraction: 0–2 mm: 0.5; 2–4 mm:0.17 and 4– 8 mm:0.03. First, we tested the influence of spatial distribution of aggregates on results. Typical crack patterns for different geometric distributions of the aggregates are shown in Figs. 11 and 12. Fig. 11 shows the crack paths and some close-ups in the deformed configuration for different spatial distributions of the aggregates. We note that the deformed configuration is shown exaggerated with a factor of 250 and overlaps of the aggregates and cement paste are due to post-processing. The crack paths of all simulations is similar but somehow different and depends on the position of the aggregates. However, the different crack paths all are within experimental observations. Exemplary, we show final crack path of one simulation in Fig. 12. While crack path is smooth in

0.25m thickness 0.1m area of experimental crack path

0.25m

0.25m

209

0.25m Fig. 9. L-shaped panel experiment.

210

A. Ghosh, P. Chaudhuri / Computational Materials Science 69 (2013) 204–215

Fig. 10. Crack path of L-shaped panel experiment using macroscopic model.

Fig. 11. Displaced configuration (factor 250), close ups and contour plot of maximum principal stress of L-shaped panel for different spatial distribution of aggregates; note the slight difference of the crack path.

A. Ghosh, P. Chaudhuri / Computational Materials Science 69 (2013) 204–215

211

Fig. 12. Displaced configuration (factor 250), close ups and contour plot of maximum principal stress of L-shaped panel.

macroscopic model, crack path is less smooth in mesoscopic model and moves around the stiffer aggregates. This creates an increased energy dissipation that needs to be compensated by an increase fracture energy in the continuum model. While this effect is less pronounced in statics, it might play a significant role in dynamics where many microcracks occur. This was recently pointed out by [68] and experimentally by [69]. Therefore, we believe mesoscopic modeling can bring more insight into failure mechanism of concrete. We note that the fracture energy in macroscopic analysis is around 15% higher compared to the fracture energy of the cement paste. Without this ’increase’, post-peak behavior is more brittle in macroscopic model and post-localization is captured less accurately. We believe that material parameters such as fracture energy are more reliable and have less scatter for homogeneous material such as cement. This is advantage for mesoscopic modeling over macroscopic modeling. The upper and lower bound of load deflection curve is shown in Fig. 13b and matches well experimental results. We also studied the influence of voids and concrete material and assume volume fraction of voids is 0.04 and deduct 0.01 from every constituent. However, the influence on results are minor and maximum deviation is less than 3% in global behavior, i.e. load– deflection curve (e.g. peak load) and dissipated energy. The influence of pores might be more pronounced in dynamic simulations when several cracks and micro-branches are expected.

7.2. Tensile test The last example is a uniaxial tensile test of a concrete specimen. The concrete consists of 29% cement, approximately (estimated) 2% micro-voids and 69% aggregates of different sizes: 0–2 mm 21%, 2–8 mm 26% and 8–16 mm 22% according to the experiment performed by [70]. In the experiment, crack nucleation occurs in the center of the specimen with the slimmest cross section. Therefore, we use mesoscopic model in the vicinity where cracking is expected while the rest of the specimen is represented on macroscopic scale. The thickness of the specimen is 100 mm and we assume plane stress conditions and model the specimen in two dimensions. The material parameters for concrete on macroscopic level are: Young’s modulus 38 GPa and Poisson’s ratio 0.2. The material parameters of the constituents are given after experimental testing by [70]. We assume fracture occurs only in cement paste that is suitable for static applications. The material properties of aggregates are: Young’s modulus 78 GPa and Poisson’s ratio 0.2 while the material for the cement paste are: Young’s modulus 36GPa, Poisson’s ratio 0.2, tensile strength 4.5 MPa and fracture energy 400 N/m. We note that not more material parameters are needed for the mesoscopic model compared to the macroscopic model besides the fact that parameters are needed for every ingredient that is more reliable since the material is more homogeneous.

A. Ghosh, P. Chaudhuri / Computational Materials Science 69 (2013) 204–215

8

8

7.5

7.5

7

7

6.5

6.5

6

6

5.5

5.5

Vertica lLoad [kN]

Vertical Load [kN]

212

5 4.5 4 3.5 3 Experiment 4,500 nodes 12,500 nodes 34,000 nodes

2.5 2 1.5

5 4.5 4 3.5 3 2 1.5

1

1

0.5

0.5

0

0

0.1

0.2

0.3

Experiment lower bound upper bound

2.5

0.4

0.5

0

0

Vertical Deflection [mm]

0.1

0.2

0.3

0.4

0.5

Vertical Deflection [mm]

Fig. 13. Vertical Load-vertical deflection curve for L-shaped panel problem for (a) pure macroscopic model with different number of nodes and (b) mesoscopic model with different spatial distribution of aggregates; only the upper and lower bound of two models are illustrated, the other 13 models tested lie within this scatter.

Fig. 16 shows load deflection curves of simulations with and without bond model. It can be seen that the bond model has a significant influence on the overall behavior. In contrast to macroscopic modeling, mesoscopic modeling (or multi-scale modeling) allows the extraction of material parameters. Simulations on macroscopic scale only requires for example the knowledge about the fracture energy. It is well known that size effects occur when the specimen is modeled only on macroscopic level. Assuming that a more fundamental knowledge about the material parameters of the constituents of the composites exists, it is possible to extract macroscopic material parameters such as the ’macroscopic’ fracture energy of the composite. It is defined by:

4

Experiment Simulation

Nomial stress [MPa]

3.5 3 2.5 2 1.5

R umax F du R Gf ¼ u¼0 dA A

1 0.5 0

0

0.001

0.002

0.003

0.004

0.005

Nominal strain [ ]

ð32Þ

meaning by dividing the area under the load–deflection curve by the cross-sectional area. For the above example, the fracture energy of the composite is 0.22 N mm/mm2. We note that it is difficult to extract the fracture energy from experiments due to difficulties in reliably measuring the load–deflection curve at post-localization.

Fig. 14. Nominal stress-nominal strain for uniaxial tensile test.

8. Discussions We used a rigid bond and an interface model as described in Section 4. The parameters of the interface model are [70]: e ¼ 0:1 N mm=mm2 . K C ¼ 0:5 MN=mm3 ; g ¼ 3:5 MPa; G We studied different refined discretizations and obtained convergent results when reaching around 20000 nodes. Note that the high number of nodes is needed due to the mesoscopic model. We also modified position of the aggregates to test influence on the results. However, the effect in this example is marginal, probably due to short crack length. This will change for the problem studied in the next section. The nominal stress–strain curve is illustrated in Fig. 14 and matches well with experimental result. The final crack paths is shown in Fig. 15. More realistic crack paths is observed as compared to macroscopic model where the crack is exactly perpendicular to the load direction. This is due to mesoscopic modeling and since crack has to propagate around aggregates.

We have presented a multi-scale approach for fracture in concrete materials. The XEFG-method has been used to model the concrete material on macroscopic scale and mesoscopic scale. On the mesoscopic scale, the constituents of concrete, i.e. cement paste and aggregates, are modeled explicitly while concrete is considered as homogenized material on the macrscopic scale. For efficiency reasons only areas around the fracture process zone are modeled on mesoscopic scale. On the mesoscopic scale, the aggregates are assumed to behave linear elastic which is a reasonable assumption under static loading conditions. A non-linear isotropic damage model is used to describe the tensile behavior of the cement paste. A cohesive crack model is introduced that can capture complex mixed-mode fracture. The interface behavior between the cement and the aggregates is modeled by a phenomenological model. One of the key features in the future will be to improve the interface behavior

A. Ghosh, P. Chaudhuri / Computational Materials Science 69 (2013) 204–215

213

Fig. 15. Displaced configuration, principal tensile stress contour plot for tensile test; displacements are shown exaggerated.

as it was demonstrated that it can have a significant influence on the overall material behavior. We have applied the method to three problems. The first (benchmark) problem was used to show the validity of our XEFG-method. The simulations were based on a pure macroscopic model. In the second problem, an L-shape specimen was loaded until failure.

We showed that a homogenized pure macroscopic description was sufficient to capture the overall behavior. The key purpose of considering finer scales was demonstrated in the last example. A specimen was loaded under uni-axial tension. Based on the mesostructure of the material we were able to capture the experimentally measured load–deflection curve. Moreover, we were able to

214

A. Ghosh, P. Chaudhuri / Computational Materials Science 69 (2013) 204–215

4 Rigid bond Cohesive Interface

Nomial stress [MPa]

3.5 3 2.5 2 1.5 1 0.5 0

0

0.001

0.002

0.003

0.004

0.005

0.006

Nominal strain [ ] Fig. 16. Nominal stress-nominal strain for uniaxial tensile test: Comparison between rigid bond model and interface model.

extract the fracture energy. It is not possible to obtain the fracture energy with a macroscopic model. Therefore, the approach can be used to extract material parameters. In the future, we intend to extend the method to fiber-reinforced concrete. As already mentioned, the interface behavior is significant. Therefore, we intend to consider even finer scales in order to gain a better understanding of those attractive materials. Moreover, we intend to improve the efficiency due to adaptivity that can be easily implemented in meshfree methods [71,72]. References [1] R. Kunigahalli, Journal of Computing in Civil Engineering 11 (1997) 92–101. [2] A. Hillerborg, M. Modeer, P.E. Peterson, Cement and Concrete Research 6 (1976) 773–782. [3] Z.P. Bazant, B.H. Oh, Materials and Structures 16 (1983) 155–177. [4] M. Jirasek, T. Zimmermann, Journal of Engineering Mechanics 124 (1998) 842– 851. [5] M.V.K.V. Prasad, C.S. Krishnamoorthy, Computer Methods in Applied Mechanics and Engineering 191 (2002) 2699–2725. [6] T. Rabczuk, J. Eibl, International Journal for Numerical Methods in Engineering 56 (10) (2003) 1421–1444. [7] T. Rabczuk, J. Eibl, L. Stempniewski, Engineering Fracture Mechanics 71 (4-6) (2004) 547–556. [8] T. Rabczuk, J. Eibl, International Journal of Solids and Structures 41 (3-4) (2004) 1061–1080. [9] T. Rabczuk, J. Eibl, International Journal of Impact Engineering 32 (11) (2006) 1878–1897. [10] J. Oliver, A.E. Huespe, P.J. Sanchez, Computer Methods in Applied Mechanics and Engineering 195 (2006) 4732–4752. [11] K.T. Danielson, S.A. Akers, J.L. O’Daniel, Journal of Computing in Civil Engineering 22 (2008) 140–146. [12] M. Slowik, Computational Materials Science 50 (2011) 1347–1352. [13] T. Rabczuk, J. Akkermann, J. Eibl, International Journal of Solids and Structures 42 (5-6) (2005) 1327–1354. [14] L.J. Malvar, M.E. Fourney, Engineering Fracture Mechanics 35 (1-3) (1990) 251–260. [15] T. Belytschko, Y.Y. Lu, L. Gu, International Journal for Numerical Methods in Engineering 37 (1994) 229–256. [16] J.M. Melenk, I. Babuska, Computer Methods in Applied Mechanics and Engineering 139 (1996) 289–314. [17] T. Belytschko, T. Black, International Journal for Numerical Methods in Engineering 45 (1999) 601–620. [18] T. Rabczuk, S. Bordas, G. Zi, Computational Mechanics 40 (3) (2007) 473–495. [19] T. Rabczuk, S. Bordas, G. Zi, Computers and Structures 88 (2010) 1391–1411. [20] T. Chau-Dinh, G. Zi, P.S. Lee, J. Song, T. Rabczuk, Computers & Structures 92-93 (2012) 242–256. [21] P.M.A. Areias, T. Rabczuk, International Journal for Numerical Methods in Engineering 74 (3) (2008) 475–505. [22] S. Bordas, T. Rabczuk, H. Nguyen-Xuan, S. Natarajan, T. Bog, Nguyen Vinh Phu, Quan Do Minh, Hiep Nguyen Vinh, Computers and Structures 88 (23-24) (2010) 1419–1443.

[23] S. Natarajan, P.M. Baiz, S. Bordas, T. Rabczuk, P. Kerfriden, Composite Structures 93 (2011) 3082–3092. [24] S. Bordas, S. Natarajan, S. Dal Pont, T. Rabczuk, P. Kerfriden, D.R. Mahapatra, D. Noel, Z. Gao, International Journal for Numerical Methods in Engineering 86 (4-5) (2011) 637–666. [25] N. Vu Bac, H. Nguyen-Xuan, L. Chen, S. Bordas, P. Kerfriden, R.N. Simpson, G.R. Liu, T. Rabczuk, CMES-Computer Modeling in Engineering & Sciences 1898 (1) (2011) 1–25. [26] L. Chen, T. Rabczuk, G.R. Liu, K.Y. Zeng, P. Kerfriden, S. Bordas, Computer Methods in Applied Mechanics and Engineering 212 (4) (2012) 250–265. [27] M. Arrea, A.R. Ingraffea, Cornell University, Department of Structural Engineering (81–13) (1982). [28] S. Hao, W.K. Liu, P.A. Klein, A.J. Rosakis, International Journal of Solids and Structures 41 (7) (2004) 1773–1799. [29] X.-P. Xu, A. Needleman, Journal of the Mechanics and Physics of Solids 42 (1994) 1397–1434. [30] D. Organ, M. Fleming, T. Terry, T. Belytschko, Computational Mechanics 18 (1996) 225–235. [31] T. Belytschko, M. Tabbara, International Journal for Numerical Methods in Engineering 39 (6) (1996) 923–938. [32] G. Ventura, J. Xu, T. Belytschko, International Journal for Numerical Methods in Engineering 54 (6) (2002) 923–944. [33] A. Lachihab, K. Sab, Computational Materials Science 33 (2005) 467–490. [34] T. Sadowski, G. Golewski, Computational Materials Science 43 (2008) 119– 126. [35] T. Rabczuk, T. Belytschko, International Journal for Numerical Methods in Engineering 61 (13) (2004) 2316–2343. [36] T. Rabczuk, G. Zi, Computational Mechanics 39 (6) (2007) 743–760. [37] T. Rabczuk, T. Belytschko, International Journal of Fracture 137 (1-4) (2006) 19–49. [38] T. Rabczuk, T. Belytschko, Computer Methods in Applied Mechanics and Engineering 196 (2007) 2777–2799. [39] T. Rabczuk, P.M.A. Areias, T. Belytschko, International Journal for Numerical Methods in Engineering 72 (5) (2007) 524–548. [40] T. Rabczuk, G. Zi, S. Bordas, H. Nguyen-Xuan, Engineering Fracture Mechanics 75 (2008) 4740–4758. [41] T. Rabczuk, G. Zi, S. Bordas, H. Nguyen-Xuan, Computer Methods in Applied Mechanics and Engineering 199 (2010) 2437–2455. [42] H. Talebi, C. Samaniego, E. Samaniego, T. Rabczuk, International Journal for Numerical Methods in Engineering 89 (9) (2012) 1009–1027. [43] T. Rabczuk, P.M.A. Areias, International Journal for Numerical and Analytical Methods in Geomechanics 30 (11) (2006) 1159–1172. [44] T. Rabczuk, P.M.A. Areias, T. Belytschko, International Journal for Numerical Methods in Engineering 69 (5) (2007) 993–1021. [45] T. Rabczuk, R. Gracie, J.-H. Song, T. Belytschko, International Journal for Numerical Methods in Engineering 81 (1) (2010) 48–71. [46] P. Lancaster, K. Salkauskas, Mathematics and Computation 37 (1981) 141–158. [47] T. Belytschko, Y. Krongauz, D. Organ, M. Fleming, P. Krysl, Computer Methods in Applied Mechanics and Engineering 139 (1996) 3–47. [48] N. Nguyen-Thanh, J. Kiendl, H. Nguyen-Xuan, R. Wüchner, K.U. Bletzinger, Y. Bazilevs, T. Rabczuk, Computer Methods in Applied Mechanics and Engineering 200 (47-48) (2011) 3410–3424. [49] R. Simpson, S. Bordas, J. Trevelyan, P. Kerfriden, T. Rabczuk, Computer Methods in Applied Mechanics and Engineering 209–212 (2012) 87–200. [50] G. Zi, T. Rabczuk, W. Wall, Computational Mechanics 40 (2) (2007) 367–382. [51] S. Bordas, T. Rabczuk, G. Zi, Engineering Fracture Mechanics 75 (2008) 943– 960. [52] J. Lemaitre, Evaluation of Dissipation and Damage in Metal Submitted to Dynamic Loading, 1971. [53] T. Gasser, G. Holzapfel, Computer Methods in Applied Mechanics and Engineering 194 (2005) 1873–1905. [54] V. Tvergaard, Engineering Fracture Mechanics 70 (2003) 1859–1868. [55] T. Zimmermann, Nuclear Engineering and Design 92 (3) (1986) 389–410. [56] T. Rabczuk, S.P. Xiao, A. Sauer, Communications in Numerical Methods in Engineering 22 (10) (2006) 1031–1065. [57] J.M. Sancho, J. Planas, A.M. Fathy, International Journal for Numerical and Analytical Methods in Geomechanics 31 (2007) 173–187. [58] H. Talebi, G. Zi, M. Silani, E. Samaniego, T. Rabczuk, A simple circular cell method for multi-level finite element analysis, Journal of Applied Mathematics, pages Article ID 526846, 2012, doi: 10.1063/1.4735246. [59] T. Rabczuk, J.Y. Kim, E. Samaniego, T. Belytschko, International Journal for Numerical Methods in Engineering 61 (7) (2004) 1009–1027. [60] T. Rabczuk, T. Belytschko, S.P. Xiao, Computer Methods in Applied Mechanics and Engineering 193 (12-14) (2004) 1035–1063. [61] V.P. Nguyen, T. Rabczuk, S. Bordas, M. Duflot, Mathematics and Computers in Simulation 79 (2008) 763–813. [62] L. Wang, N. Zhao, S. Ma, Computational Materials Science 46 (2009) 672– 676. [63] I. Babuska, T.J. Oden, International Journal of Numerical Analysis and Modeling 3 (2006) 255–272. [64] H.A. Kormeling, H.W. Reinhard, Determination of the Fracture Energy of Normal Concrete and Epoxy Modified Concrete. Report No. 5-83-15, Delft University of Technology, Delft, 1983. [65] J. Alfaiate, E.b. Pires, J.A.C. Martins, Computers and Structures 63 (1997) 17–26. [66] P. Bocca, A. Carpinteri, S. Valente, International Journal of Solids and Structures 27 (1991) 1139–1153.

A. Ghosh, P. Chaudhuri / Computational Materials Science 69 (2013) 204–215 [67] B Winkler, Traglastuntersuchungen von unbewehrten und bewehrten betonstrukturen auf der grundlage eines objektiven werkstoffgesetzes fur beton. PhD thesis, University of Innsbruck, Austria, 2001. [68] T. Rabczuk, J.-H. Song, T. Belytschko, Engineering Fracture Mechanics 76 (2009) 730–741. [69] S. Sharon, J. Fineberg, Physical Review B 54 (1996) 7128–7139.

215

[70] M. van Vliet, J. van Mier, Engineering Fracture Mechanics 65 (2000) 165–188. [71] T. Rabczuk, T. Belytschko, International Journal for Numerical Methods in Engineering 63 (11) (2005) 1559–1582. [72] T. Rabczuk, E. Samaniego, Discontinuous modelling of shear bands using adaptive meshfree methods, Computer Methods in Applied Mechanics and Engineering 197 (6-8) (2008) 641–658.