- Email: [email protected]

Contents lists available at ScienceDirect

Comput. Methods Appl. Mech. Engrg. journal homepage: www.elsevier.com/locate/cma

A hybrid ﬁnite element-scaled boundary ﬁnite element method for crack propagation modelling E.T. Ooi *, Z.J. Yang Department of Engineering, The University of Liverpool, Liverpool L69 3GQ, UK

a r t i c l e

i n f o

Article history: Received 17 June 2009 Received in revised form 9 November 2009 Accepted 15 December 2009 Available online 23 December 2009 Keywords: Scaled boundary ﬁnite element method Finite element method Discrete crack model Multiple cohesive crack propagation Remeshing Fracture mechanics

a b s t r a c t This study develops a novel hybrid method that combines the ﬁnite element method (FEM) and the scaled boundary ﬁnite element method (SBFEM) for crack propagation modelling in brittle and quasibrittle materials. A very simple yet ﬂexible local remeshing procedure, solely based on the FE mesh, is used to accommodate crack propagation. The crack-tip FE mesh is then replaced by a SBFEM rosette. This enables direct extraction of accurate stress intensity factors (SIFs) from the semi-analytical displacement or stress solutions of the SBFEM, which are then used to evaluate the crack propagation criterion. The fracture process zones are modelled using nonlinear cohesive interface elements that are automatically inserted into the FE mesh as the cracks propagate. Both the FEM’s ﬂexibility in remeshing multiple cracks and the SBFEM’s high accuracy in calculating SIFs are exploited. The efﬁciency of the hybrid method in calculating SIFs is ﬁrst demonstrated in two problems with stationary cracks. Nonlinear cohesive crack propagation in three notched concrete beams is then modelled. The results compare well with experimental and numerical results available in the literature. Ó 2009 Elsevier B.V. All rights reserved.

1. Introduction Since the smeared [1] and discrete [2] crack models were introduced in the late 1960s, they have been extensively used to model crack propagation in brittle and quasi-brittle materials. In smeared crack models the cracks which are ‘‘smeared” out over the domain are represented by softening stress–strain relationships in the normal and tangential directions of the crack surfaces. Discrete crack models, on the other hand, represent cracks by separating the nodes along the crack path and usually needs remeshing. Although the smeared crack models have been more popular due to their ease in implementation, they suffer from spurious stresses due to the kinematic incompatibility between elements and have to be remedied by introducing a shear retention factor [3] or by rotating the cracks so that they are aligned with the principal stress and strain directions [4]. Recent approaches to model discontinuities caused by cracks enrich either the elements cut by the crack (embedded crack models e.g., [5–7]) or the nodes along the crack path (extended ﬁnite element method (XFEM) e.g., [8–10]) using discontinuous functions. These approaches do not need remeshing and have shown tremendous potential in modelling crack propagation. Discrete crack models use remeshing algorithms to update mesh topology as cracks propagate. They require more implemen* Corresponding author. Tel.: +44 7942 521 333. E-mail address: [email protected] (E.T. Ooi). 0045-7825/$ - see front matter Ó 2009 Elsevier B.V. All rights reserved. doi:10.1016/j.cma.2009.12.005

tation effort. However, they are capable of accurately representing the discontinuities between the crack surfaces, and are useful to study the local behaviour in the vicinity of cracks. Most of the discrete crack models are based on the ﬁnite element method (FEM). As stress singularities at crack tips are not accounted for by conventional FEM formulations, very ﬁne crack tip meshes are often needed to accurately calculate fracture parameters such as the stress intensity factors (SIFs), energy release rates and crack-tip stresses, which are used in crack propagation criteria. Various remedies have been proposed, such as the quarter-point elements [11], hybrid Trefftz elements [12], elements having shape functions with crack tip asymptotic expansions [13], enriched elements [14] and hybrid crack elements (HCE) [15–18]. In general, however, these remedial methods still need ﬁne crack-tip meshes or impose special requirements on the shape and element types of crack-tip meshes. As a result, the remeshing algorithms used in FEM to accommodate crack propagation are usually very complicated. Recently, Yang et al. [19–21] successfully applied the scaled boundary ﬁnite element method (SBFEM) to model discrete crack propagation in concrete structures. The SBFEM is a semi-analytical method developed by Wolf and Song [22] and is very efﬁcient in solving problems with unbounded media, discontinuities and singularities. In modelling crack propagation problems, the SBFEM can extract accurate SIFs from semi-analytical solutions of displacements or stresses using substantially fewer degreesof-freedom (DOF) than the FEM. Moreover, the remeshing and mesh mapping procedures in the SBFEM are much simpler yet

E.T. Ooi, Z.J. Yang / Comput. Methods Appl. Mech. Engrg. 199 (2010) 1178–1192

more accurate. The SBFEM-based methods developed in [19–21] are particularly suitable for problems with single crack or a few cracks. For problems with many cracks that become too close during propagation, the remeshing operation is cumbersome because the subdomains may become so distorted that not all the nodes are visible from their scaling centres. To tackle this problem, this study further develops a hybrid method that combines the FEM and the SBFEM. It can model crack propagation in both linear elastic- and nonlinear-fracture mechanics problems. In the latter, cohesive interface elements (CIEs) are automatically inserted into the FE mesh as the cracks propagate. A very simple yet ﬂexible local remeshing procedure, based on the FE mesh only, is used to accommodate crack propagation. The crack-tip FE mesh is then replaced by a SBFEM rosette to extract accurate SIFs, which are subsequently used to evaluate the crack propagation criterion. Both the FEM’s ﬂexibility in remeshing for multiple crack problems and the SBFEM’s high accuracy in calculating SIFs are thus exploited. 2. The hybrid ﬁnite element-scaled boundary ﬁnite element method

x1 þ x2 2x0 ðx2 x1 Þs x ¼ x0 þ n ¼ x0 þ nxs ðsÞ; þ 2 2 y þ y2 2y0 ðy2 y1 Þs y ¼ y0 þ n 1 ¼ y0 þ nys ðsÞ; þ 2 2

uðn; sÞ ¼ NðsÞuðnÞ;

ð2Þ

ð3Þ

where uðnÞ is the displacement function in the radial direction and N(s) is the shape function matrix. For a two-noded element on the subdomain boundary

NðsÞ ¼

N1

0

N2

0

0

N1

0

N2

ð4Þ

with

1 ð1 sÞ; 2

N2 ¼

1 ð1 þ sÞ: 2

ð5Þ

Substituting Eq. (3) into the principle of virtual work, the governing equations of SBFEM for a subdomain can be obtained as [23]

T p ¼ E0 uðnÞ;n þ E1 uðnÞ

ð6Þ

n¼1

E0 n2 uðnÞ;nn þ ðE0 þ E

2.1. The scaled boundary ﬁnite element method

ð1Þ

ðx1 ; y1 Þ and ðx2 ; y2 Þ are the nodal coordinates of an element on the boundary. In each subdomain, an approximate solution is sought in the form

N1 ¼

The basic idea of this method is to carry out remeshing solely based on the FE mesh, and to calculate the SIFs from the SBFEM rosette which replaces the crack-tip FE mesh. The key elements of the method are described below.

1179

1T

E1 ÞnuðnÞ;n E2 uðnÞ ¼ 0;

ð7Þ

where In the SBFEM, a domain is partitioned into a number of subdomains. Fig. 1a shows a typical two-dimensional SBFEM mesh with three subdomains, in which subdomain 1 contains a crack (cracked subdomain). Only the boundaries of subdomains are discretized using one-dimensional ﬁnite elements. The crack surfaces are not discretized. Each subdomain has a scaling centre ðx0 ; y0 Þ from which all the nodes on the boundary of the subdomain are visible. In the cracked subdomain, the scaling centre is positioned at the crack tip. The coordinates of any point in a subdomain are uniquely deﬁned by a local coordinate system ðn; sÞ. n is the radial coordinate which varies from zero at the scaling centre to one on the subdomain boundary. s is the circumferential coordinate (see Fig. 1b). The local coordinates ðn; sÞ are related to the Cartesian coordinates ðx; yÞ by the scaling equations [23]

(a)

s

Crack edges (side-faces)

p¼

Z

NðsÞT t ds;

ð8Þ

s

E0 ¼

Z

B1 ðsÞT DB1 ðsÞjJjds;

ð9Þ

B2 ðsÞT DB1 ðsÞjJjds;

ð10Þ

B2 ðsÞT DB2 ðsÞjJjds:

ð11Þ

s

E1 ¼ E2 ¼

Z

Zs s

Eqs. (6) and (7) ensure equilibrium is satisﬁed analytically in the radial direction and in the ﬁnite element sense in the circumferential direction. In Eqs. (8)–(11), p is the nodal load vector due to boundary tractions, E0 ; E1 and E2 are coefﬁcient matrices that depend on the geometry and material properties of the subdomain,

(b) s

Nodes

Node 2(x2, y2)

p 1 3

ξ ξ

s=1 ξ=1

2

Radial Similar lines curve Defining Scaling curve centre, ξ = 0

s = -1

Node 1 (x1, y1)

Scaling centre, ξ=0

Internal Boundary Scaling Centre Fig. 1. (a) A typical SBFEM mesh and (b) a two-noded line element on subdomain boundary.

1180

E.T. Ooi, Z.J. Yang / Comput. Methods Appl. Mech. Engrg. 199 (2010) 1178–1192

t is the traction load intensity and D is the material constitutive matrix. jJj is the determinant of the Jacobian matrix at the boundary,

ð12Þ

9 8 > = < wxx ðsÞ > wi ðsÞ ¼ wyy ðsÞ ¼ D½ki B1 ðsÞ þ B2 ðsÞui : > ; : w ðsÞ > xy

7 xs ðsÞ;s 5NðsÞ;

ð13Þ

ys ðsÞ;s 3 0 7 xs ðsÞ 5NðsÞ;s :

It is clear that both the displacement (Eq. (22)) and stress solutions (Eqs. (23) and (24)) are analytical with respect to the radial coordinate n.

ð14Þ

jJj ¼ xs ðsÞys ðsÞ;s ys ðsÞxs ðsÞ;s ; B1 ðsÞ and B2 ðsÞ are the strain–displacement matrices

2 16 B ðsÞ ¼ 4 jJj 2 16 4 jJj

B2 ðsÞ ¼

0 xs ðsÞ;s ys ðsÞ 0 xs ðsÞ

3

0

ys ðsÞ;s

1

ys ðsÞ

uðnÞ ¼ c1 nk1 u1 þ c2 nk2 u2 þ ;

ð15Þ

where ki and ui are eigenvalues and corresponding eigenvectors that closely satisfy equilibrium in the radial direction. The constants ci are the contribution of each mode to the displacement ﬁeld. Omitting the subscript, the displacement of each mode takes the form

uðnÞ ¼ nk u;

ð16Þ

where u is the modal displacement at the boundary nodes and k is the modal scaling factor in the n direction. Substituting Eq. (16) into Eqs. (6) and (7), a linear standard eigenvalue problem can be derived following [23] as 1

E0 E1 1

E E

0 1

E

1T

T

E0 2

þE

1

1

E E0

# u

1

q

u ¼k : q

ð17Þ

For a subdomain having N DOFs, Eq. (17) yields 2N modes of which only the N modes with positive real eigenvalues lead to ﬁnite displacements at the scaling centre. This subset of modal displacements is U ¼ fu1 ; u2 ; . . . ; uN g. The displacement solution of Eq. (7) for a subdomain can then be expressed as

uðnÞ ¼

N X

ð24Þ

2.2. The cohesive interface elements

The solution uðnÞ of Eq. (7) for a subdomain has the form

"

where

ci nki ui

ð18Þ

The fracture process zone (FPZ) is modelled using 4-noded CIEs (shown in Fig. 2) developed by Xie [24]. They are characterised by softening laws relating the normal traction r and shear traction s on the crack surfaces to the crack opening displacement (COD) and crack sliding displacement (CSD), respectively as

r kn 0 ¼ s 0 ks

COD CSD

COD ¼k ; CSD

ð25Þ

where kn and ks are the normal and shear stiffness shown in Fig. 3. The element stiffness matrix Kc in the local coordinate system is derived from the principle of virtual work [24] Ng AX MT ki Mi Hi : 2 i¼1 i

Kc ¼

ð26Þ

A is the crack surface area, Hi the one-dimensional Gaussian weight factor, and N g the number of Gaussian integration points. M is the shape function matrix

Mi ¼

N1

0

N2

0

N2

0

N1

0

0

N1

0

N2

0

N2

0

N1

;

ð27Þ

where the shape functions, N1 and N 2 are deﬁned in Eq. (7). In the present study, shear fracture resistance is neglected because the numerical examples modelled did not report s-CSD relations. This negligence will only slightly affect the post-peak but not the overall load–displacement curves [20].

i¼1

and the stiffness matrix of the subdomain is T

Ks ¼ E0 U½kU1 þ E1 :

Aggregate

ð19Þ

The stiffness matrix and load vector of each subdomain are assembled in a manner similar to the FEM. The global nodal displacement vector of all the subdomains Ub is calculated as

Ub ¼ K1 SBFEM P;

ð20Þ

where KSBFEM and P are the assembled stiffness matrices and load vectors of all the subdomains in the model. The constants c ¼ ðfc1 ; c2 ; . . . ; cN gT Þ in Eq. (18) are computed from the nodal displacements on the subdomain boundary ub which is a subset of Ub by 1

c ¼ U ub :

Concrete bulk Real crack

Traction distribution along FPZ

f

t

ð21Þ

The displacements at any point in the subdomain are recovered from

uðn; sÞ ¼ NðsÞ

Fracture process zone (FPZ)

N X

ci nki ui

ð22Þ

4 3

Fictitious crack tip node

wc

y

i¼1

and the stresses are computed as

rðn; sÞ ¼ NðsÞ

N X i¼1

ci nki 1 wi ðsÞ;

1

ð23Þ

Interface element 2

2

x Interface element 1

Fig. 2. Modelling FPZ with interface elements after [24].

1181

E.T. Ooi, Z.J. Yang / Comput. Methods Appl. Mech. Engrg. 199 (2010) 1178–1192

(a) σ

(b)

ft

τ τu Loading Unloading

kn 0

GfII

ks0

Gf

f1

ks -sc

kn

w0

sc CSD Loading

wc COD

w1

Unloading τu

Fig. 3. (a)

Γ

Γt

r-COD and (b) s-CSD relation for CIEs.

t y

Ω SBFEM

Finite elements

Crack

r θ

x

Ω FEM SBFEM rosette Fig. 4. Discretisation of a domain with a SBFEM rosette surrounding the crack tip and ﬁnite elements in regions away from the crack tip.

2.3. The scaled boundary ﬁnite element rosette Fig. 4 shows a cracked domain that is discretized using ﬁnite elements in regions XFEM away from the crack tip, and a SBFEM rosette, XSBFEM in the vicinity of the crack tip (shaded region). The rosette has two layers of subdomains: a cracked subdomain in the interior with the scaling centre at the crack tip, and an outer layer of normal subdomains with scaling centres at their geometric centres, connecting the cracked subdomain with the FE mesh. Because of the semi-analytical nature and the ﬂexible nodal distributions of the SBFEM, the cracked subdomain can be modelled with more nodes to ensure the accuracy of the SIFs, while the outer layer of subdomains can be modelled with nodes compatible with the surrounding ﬁnite elements. The FEM and SBFEM meshes are coupled through the common nodes on the rosette’s exterior (shown as unﬁlled circles in Fig. 4). Discretising the principle of virtual work in XFEM (following the usual FEM procedures, e.g., [25]) and in XSBFEM (following the procedures in Section 2.1), the coupled system of equations is

KU ¼ F;

ð28Þ

where U and F are the global nodal displacement and load vectors, respectively, and

8 KFEM þ KSBFEM > > > < without CIEs : linear elastic fracture mechanics K¼ > K > FEM þ KSBFEM þ KCIE > : with CIEs : nonlinear-fracture mechanics

ð29Þ

Fig. 5. Cartesian and polar coordinate system with a crack tip.

is the assembled global stiffness matrix with KFEM and KCIE being the assembled stiffness matrices of the ﬁnite elements and CIEs. The assembly process is the same as in the FEM, treating the SBFEM subdomains as super elements. 2.4. Calculation of stress intensity factors The SIFs can be extracted directly from either the semi-analytical displacement solutions [26] or the stress solutions [27]. If the displacement solutions are used, the relations between the SIFs and crack-tip displacement ﬁelds with a coordinate system shown in Fig. 5 are [28]

rﬃﬃﬃﬃﬃﬃﬃ r h 1 2 h cos ðj 1Þ þ sin 2p 2 2 2 rﬃﬃﬃﬃﬃﬃﬃ K II r h 1 h þ sin ðj þ 1Þ þ cos2 ; 2 2 2 G 2p rﬃﬃﬃﬃﬃﬃﬃ KI r h 1 h uy ðr; hÞ ﬃ sin ðj 1Þ þ cos2 2 2 2 G 2p rﬃﬃﬃﬃﬃﬃﬃ K II r h 1 2 h ðj þ 1Þ þ sin : þ cos 2 2 2 G 2p

ux ðr; hÞ ﬃ

KI G

In Eqs. (30) and (31), G is the shear modulus and

j¼

where

3 4m

for plane strain

ð3 mÞ=ð1 þ mÞ for plane stress

m is the Poisson’s ratio.

;

ð30Þ

ð31Þ

j is ð32Þ

1182

E.T. Ooi, Z.J. Yang / Comput. Methods Appl. Mech. Engrg. 199 (2010) 1178–1192

Using h ¼ p to decouple K I and K II in Eqs. (30) and (31), and r ¼ ro where r o is the length from the crack tip to the crack mouth node at h ¼ p, the SIFs are

KI K II

¼

2G jþ1

sﬃﬃﬃﬃﬃﬃﬃ 2p uy ðr o ; pÞ ; r o ux ðr o ; pÞ

Start

Pre-processing

ð33Þ

ux ðr o ; pÞ and uy ðr o ; pÞ can be extracted from the displacement ﬁeld

^ ¼ NðsÞ u

X

ki

ci n ui ;

Save FEM mesh

ð34Þ

i¼I;II

Insert SBFEM rosette

where I and II are the two modes with ki ¼ 0:5: When cohesive cracks are present, a shadow domain procedure [20] is used to calculate the SIFs. In this procedure, the cracked subdomain in Fig. 6a is split into two normal subdomains as shown in Fig. 6b so that the crack surfaces are discretised and CIEs can be inserted into the SBFEM rosette. The SBFEM mesh in Fig. 6b is the shadow domain of the mesh in Fig. 6a. The nodal displacements and cohesive tractions are computed from the shadow domain. The SIFs considering cohesive tractions are then determined from the cracked subdomain after mesh mapping the nodal displacements from the shadow domain as follows. In the SBFEM, if the nodal load distribution on a side-face Ft ðnÞ in a subdomain has the following form of power functions of n

Ft ðnÞ ¼ nt Ft

Generate shadow domain

Yes

Mesh-map state variables

Previously remeshed? No

Insert CIEs

Nonlinear analysis

ð35Þ

the nodal displacement mode for this side-face force ut can be shown to be [23] 2

0

ut ¼ ½ðt þ 1Þ E þ ðt þ 1ÞðE

1T

1

2 1

E Þ E Ft ¼ B1 ðtÞFt

Yes End

Structure failed?

ð36Þ No

Ft is a vector of nodal components of the side-face load and t is a constant depending on the load distribution on the side-face. In order to use Eqs. (35) and (36), the normal traction distribution on the crack faces rðnÞ is assumed to be represented by a summation of M number of power functions of n

rðnÞ ¼ ft

M X

ei nti ;

Calculate stress intensity factors

Any cracks propagate?

ð37Þ

i¼1

Yes

where ft is the tensile strength and ei are coefﬁcients to be determined. The exponents t i are

t i ¼ i 1 þ l ði ¼ 1; 2; . . . ; MÞ;

No

FEM-based remeshing

ð38Þ Fig. 7. Flow chart of the hybrid method.

where l is an arbitrary shift (0.0001 in this study) so that t i is not equal to any of the ki in Eq. (18). In order to solve M number of un-

(b)

(a)

8

8 7

7

2

2 σ(ξ)

6

α

A

ξ

6

1 L0

1

9

3

3

5

5

4

4 CIEs

y x L

Fig. 6. SBFEM rosette for a crack with (a) a crack subdomain in its interior and (b) a shadow domain of the rosette with the cracked subdomain split into two.

1183

E.T. Ooi, Z.J. Yang / Comput. Methods Appl. Mech. Engrg. 199 (2010) 1178–1192

(a)

(b) Case 1

Case 2

Old crack tip

Old crack tip

Drag existing node

(d)

(c)

New crack tip

New crack tip

Old crack tip

Old crack tip

Fig. 8. FEM-based remeshing: (a) mesh from a previous load step, (b) two cases of crack propagation, (c) remeshing in Case 1, and (d) remeshing in Case 2.

(a)

9

8

10

7

1 θ2

6

5

(b)

9

2

7

3

6

10

8 1

2

θ1 4

5

3

4

SBFEM Rosette (c)

9

8

10 2

7 6

(d)

5

4

3

Fig. 9. Generation of SBFEM rosette (a) FE mesh after remeshing, (b) shape optimisation, (c) removal of crack-tip elements, and (d) replacing crack-tip elements with SBFEM rosette.

knowns ei in Eq. (37), M number of equations are built using known values of rðnÞ at the Ng Gaussian points of CIEs in the cracked subdomain, rð0Þ ¼ ft at the crack tip, and r(1) at the crack mouth as

rj ¼ rðnj Þ ¼ ft

M X

t ei nj i

j ¼ 1; 2; . . . ; M:

ð39Þ

and L is the crack length and li is the distance from the ith point to the crack tip. Eq. (39) is then solved to obtain the coefﬁcients ei and the nodal side-face load vector becomes

Ft ðnÞ ¼

M X i¼1

nti Fti ¼

M X

nti Aft ei R;

ð41Þ

i¼1

i¼1

where A is the area of the crack surface and R is a transformation vector

Here,

ni ¼

li L

ði ¼ 1; 2; . . . ; MÞ

ð40Þ

R ¼ f sin a cos a 0 . . . 0 sin a cos a gT

ð42Þ

1184

E.T. Ooi, Z.J. Yang / Comput. Methods Appl. Mech. Engrg. 199 (2010) 1178–1192

σ

(a)

(b)

L a W

σ

(c)

(e)

(d)

Fig. 11. Convergence of K I with increasing DOFs.

rðn; sÞ ¼ NðsÞ

N X

i ðsÞ; ci nki 1 w

ð48Þ

i¼1

where

9 8 > = < wxx ðsÞ > i ðsÞ ¼ w yy ðsÞ ¼ D½ki B1 ðsÞ þ B2 ðsÞu i; w > > ; : wxy ðsÞ Fig. 10. Plate with an edge crack subjected to a far ﬁeld stress, (a) dimensions and loading conditions, (b) SBFEM rosette, (c) structured N M mesh, (d) hybrid mesh and (e) unstructured hybrid mesh.

and a is deﬁned as in Fig. 6a. The complete displacement solution in the cracked subdomain is thus

" uðn; sÞ ¼ NðsÞ

N X

ci nki ui þ

i¼1

M X

# ð43Þ

j¼1

c ¼ U1 ðub Ut eÞ;

ð44Þ

where

Ut ¼ fAft B1 ðt 1 ÞR Aft B1 ðt2 ÞR . . . Aft B1 ðt M ÞRg;

ð45Þ

e ¼ fe1 ; e2 ; . . . ; eM gT :

ð46Þ

The displacement and stress ﬁelds in the cracked subdomain can then be expressed as

uðn; sÞ ¼ NðsÞ

NþM X

ki ¼ k1 ; k2 ; . . . ; kN ; t Nþ1 þ 1; tNþ2 þ 1; . . . ; t NþM þ 1; ¼ u1 ; u2 ; . . . ; uN ; utNþ1 ; utNþ2 ; . . . ; utNþM ; u

ð51Þ

c ¼ fc1 c2 cN e1 e2 eM gT :

ð52Þ

ð50Þ

i replacing Eq. (33) is then used to compute the SIFs with ci ; ki and u ci ; ki and ui in Eq. (34), respectively. 2.5. Crack propagation criterion

nti þ1 utj :

The nodal displacements on the subdomain boundary ub , obtained by mesh mapping the displacements from the shadow domain, can then be used to solve for the coefﬁcients ci in Eq. (43) as

"

ð49Þ

# k i

ci n u i ;

ð47Þ

i¼1

After the SIFs are calculated, the following crack propagation criterion is used to judge whether the crack propagates under the current loading

K I P 0:

ð53Þ

The crack propagation direction is determined using a linear elastic fracture mechanics criterion. In this study, the maximum circumferential stress criterion is used. 3. Aspects of numerical implementation Fig. 7 shows a ﬂowchart illustrating the solution of the hybrid method. The domain is initially discretised using ﬁnite elements in the pre-processing module. This FE mesh is stored in case remeshing is needed. The crack-tip elements are then replaced by SBFEM rosettes and shadow subdomains are generated. If

Table 1 Normalised K I for the edge-cracked plate under mode-I loading. FEM

XFEM

Hybrid

No. Elements

No. DOF

Normalised K I

No. Elements

No. DOF

Normalised K I

No. Elements

No. DOF

Normalised K I

48 8 16 16 32 20 40 32 64

94 314 1138 1742 4322

0.750 0.857 0.924 0.958 0.960

57 9 15 17 31 21 39 33 63

136 368 1216 1832 4448

0.815 0.916 0.963 0.971 0.983

48 8 16 16 32 20 40 32 64

222 442 1266 1870 4450

0.926 0.965 0.985 0.990 0.995

1185

E.T. Ooi, Z.J. Yang / Comput. Methods Appl. Mech. Engrg. 199 (2010) 1178–1192

τ=1

E = 30 × 106 ν = 0.25

L a W

Fig. 15. Clamped plate with an edge crack subjected to unit shear traction. Fig. 12. Convergence of K I with increasing nodes in cracked subdomain.

nonlinear analysis is then carried out to obtain the structure’s responses. To capture the equilibrium path characterised by snapback due to material softening, a local arc-length algorithm [29,30] that uses a local updated normal plane constraint equation is implemented in the nonlinear solution module. After a converged solution is obtained, the SIFs of each crack are calculated using the procedures described in Section 2.4. If any crack satisﬁes the condition in Eq. (53), remeshing is carried out on the previously stored FE mesh. The procedures for FEM-based remeshing, SBFEM rosette generation and mesh mapping are described below. 3.1. FEM-based remeshing

Fig. 13. Effect of relative size of cracked subdomain to rosette on accuracy of K I .

remeshing is performed in the previous loading step, the nodal displacements in the current hybrid mesh are mapped from the old mesh to be used as initial values for the current loading step. A

In this study, a simple local remeshing procedure devised by Xie and Gerstle [31] is used. Fig. 8a shows the FE mesh near a crack tip. First, all the quadrilateral elements that are connected to the old crack tip are split into triangular elements, resulting in the mesh shown in Fig. 8b. There are now two possible cases that a crack can propagate. Case 2 can be simpliﬁed to Case 1 if an existing node is dragged to the location of the new crack tip. To propagate the crack, the old crack-tip node is ﬁrst split into two nodes and the element connectivities that are affected by the addition of the new node at the old crack tip are appropriately updated. Fig. 8c and d show the remeshing procedures in Cases 1 and 2, respectively. For problems with multiple cracks, each crack is remeshed in turn using the same procedure.

(a) KI/KIexact = 0.990

(b) KI/KIexact = 1.003

(c) KI/KIexact = 1.010

(d) KI/KIexact = 0.996

(e) KI/KIexact = 1.008

(f) KI/KIexact = 0.978

(g) KI/KIexact = 1.024

(h) KI/KIexact = 0.972

Fig. 14. Normalised K I from different shapes of cracked subdomains and SBFEM rosettes.

1186

E.T. Ooi, Z.J. Yang / Comput. Methods Appl. Mech. Engrg. 199 (2010) 1178–1192

Table 2 Normalised stress intensity factors for the edge-cracked plate under shear loading. No. Elements in FEM mesh

FEM

48 8 16 16 32 20 40

Hybrid

No. DOF

Normalised K I

Normalised K II

No. DOF

Normalised K I

Normalised K II

94 314 1138 1742

0.807 0.896 0.946 0.957

0.922 0.957 0.977 0.981

222 442 1266 1870

0.908 0.962 0.988 0.992

0.953 0.980 0.989 0.990

3.2. Generation of SBFEM rosette

2H

2H 0.2H

0.1F

F 0.2H

A H

0.2H Fig. 16. A single edge notched concrete beam [37].

After remeshing, the crack-tip elements are replaced by SBFEM rosettes. Fig. 9a illustrates a crack-tip mesh after remeshing where the crack-tip elements are shown in shaded areas. Because the shape of the SBFEM rosettes slightly affects the solution accuracy as will be shown later, the element connectivities of any triangular crack-tip elements and their neighbours are modiﬁed if necessary. For example, if the vertex angles at the crack tip h1 and h2 (Fig. 9a) exceed a critical value hcrit , the crack-tip mesh will be modiﬁed to that shown in Fig. 9b. hcrit ¼ 85 is assumed in this study. The crack-tip node and elements are then removed from the FE mesh (Fig. 9c) and replaced with a SBFEM rosette (Fig. 9d). The nodes in the rosette are generated from the coordinates of nodes 2–9

Fig. 17. Results for SEN beam, (a) load-loading point displacement, (b) load-CSD and (c) load-COD curves.

E.T. Ooi, Z.J. Yang / Comput. Methods Appl. Mech. Engrg. 199 (2010) 1178–1192

1187

Fig. 18. Predicted crack propagation process of SEN beam.

80

F/5

using Eqs. (1) and (2) with the crack tip as the scaling centre and n = 0.5. They are then connected to form the normal and cracked subdomains of the rosette.

F

40 3.3. Mesh mapping

200 40

320

80

400

Fig. 19. A double-edge notched concrete beam [41].

After remeshing, mesh mapping is carried out to transfer the displacements from the old mesh to the new mesh. Since the mesh has both ﬁnite elements and SBFEM subdomains, mesh mapping depends on whether the point in the new mesh belongs to XFEM or XSBFEM in the old mesh. If the point in the new mesh belongs to XSBFEM in the old mesh, Eq. (18) is used to calculate the displacements. Otherwise, an inverse transformation [32] is used to obtain

1188

E.T. Ooi, Z.J. Yang / Comput. Methods Appl. Mech. Engrg. 199 (2010) 1178–1192

Fig. 20. Load-loading point displacement curve of DEN beam.

the coordinates of the point in terms of local coordinates in FEM. The displacements are then computed using the FEM approximation. 4. Numerical examples Five examples are modelled. The ﬁrst two examples of stationary crack problems demonstrate the efﬁciency of the hybrid method in computing SIFs for linear elastic fracture mechanics problems. The other three examples of notched concrete beams illustrate its effectiveness in modelling nonlinear cohesive crack propagation. In the numerical examples involving nonlinear cohesive crack propagation, a local arc length equation solver is used to capture equilibrium paths characterised by snapback and material softening. The local updated normal plane constraint equation [29,30] is used to compute the load increment and iterative loading factors. The controlled variable is the relative displacements of two CODs and two CSDs of the two pair of nodes of the CIEs. 4.1. An edge-cracked plate under mode-I fracture Fig. 10a shows an edge-cracked plate subjected to a far ﬁeld stress r. The dimensions are: W = 7, L = 16 and a = 3.5. The exact solution of the SIF is [33]

pﬃﬃﬃﬃﬃﬃ K I ¼ C r ap;

ð54Þ

where C is a ﬁnite geometry correction factor

C ¼ 1:12 0:231

a

a 2

a 3 21:72 þ 10:55 W W W

a 4 þ 30:39 : W

meshes using the FEM, the XFEM and the hybrid method. Fig. 11 shows the convergence of K I with increasing number of DOFs. The results in Fig. 11 indicate that K I calculated using the hybrid method converges much faster than FEM and XFEM. For example, to obtain about 5% error, the hybrid method needs 442 DOF while the FEM and XFEM needs 1742 DOF and 1216 DOF, respectively. Fig. 12 shows that K I converges monotonically as the number of nodes in the cracked subdomain of the hybrid meshes increase. Fig. 10e shows an unstructured hybrid mesh with 1156 DOF. The normalised K I of this mesh is 0.972. Compared with the structured hybrid mesh that has approximately the same number of DOF (on the third row of Table 1), the difference of K I is only 0.5%. Therefore, the accuracy of the hybrid method appears not to be adversely affected by the arrangement of FE mesh. Fig. 13 shows the effects of the relative size of the cracked subdomain to the rosette, i.e., the ratio n on K I . n P 0:5 has little effect on K I . As n decreases, the accuracy of K I ﬁrst increases, until it reaches a critical value ncrit . This may be because for a low n, the cracked subdomain has nodes that are sufﬁciently close to the crack tip, allowing the singularity to be modelled more accurately. When n < ncrit , the accuracy of K I worsens as n may be too low that it leads to triangular normal subdomains (2–9 in Fig. 10b), which may result in less accurate solutions. For the meshes considered, the optimal ncrit lies between 0.1 and 0.3. Since there is no single optimal value for n; n ¼ 0:5 is used in all the following analyses considering that the increase in accuracy as n ! 0:3 is not signiﬁcant. Fig. 14 examines the shape effects of the cracked subdomain and the SBFEM rosette on the accuracy of the normalised K I using a 20 40 FE mesh. All the ﬁrst four rosettes with different cracked subdomain shapes (Fig. 14a–d) predict less than 1% error. The maximum difference of error from typical shapes of SBFEM rosettes (Fig. 14e–h) is 3%. These results demonstrate that K I is almost independent of the shape of the cracked subdomains and the rosettes. This is important because the remeshing procedure may lead to irregular cracked subdomains or rosettes listed in Fig. 14. 4.2. An edge-cracked plate under mixed-mode fracture A clamped plate, shown in Fig. 15, having W = 7, L = 16 and a = 3.5 subjected to unit shear traction s is considered. Similar meshes as in the ﬁrst example are used. The normalised K I and K II from the FEM and the hybrid method listed in Table 2. The reference solutions are K I ¼ 34:0 and K II ¼ 4:55 [36]. As in the ﬁrst example, the hybrid method converges faster than the FEM for both SIFs. For example, the hybrid method needs 1266 DOF to achieve an error of less than 2% for K I and K II , while the FEM using a mesh with 1741 DOF is still unable to achieve the same level of accuracy. 4.3. Cohesive crack propagation: single edge notched concrete beam

ð55Þ

The plate is meshed using N M 4-noded quadrilateral elements (Fig. 10c) where N and M are the number of elements in the horizontal and vertical dimensions of the plate. The crack-tip elements are replaced with a SBFEM rosette with 67 nodes in the cracked subdomain (Fig. 10b). Thus for a mesh having equal number of elements, the hybrid method has more DOFs than both FEM and XFEM. The hybrid mesh is shown in Fig. 10d, where the rosette is illustrated by thick lines. The J-Integral approach in Abaqus [34] and a domain interaction integral method [35] are used to calculate the SIFs by the FEM and the XFEM, respectively. Table 1 compares the computed K I (normalised with the exact solution of Eq. (54)) of ﬁve

The single edge notched (SEN) beam subjected to a four-point shear test shown in Fig. 16 is modelled. Similar beams were studied in [37–39]. The depth of the beam H is 200 mm and the thickness is 100 mm. The material properties are: Young’s modulus E = 28 GPa, Poisson’s ratio m ¼ 0:1, tensile strength ft ¼ 2:4 MPa and fracture energy Gf ¼ 0:122 N=mm. A bilinear r-COD softening law with wc ¼ 18Gf =5f t ; w1 ¼ 2wc =9 and f1 ¼ ft =3 (see Fig. 3a) according to Petersson [40] is used to model the cohesive crack. The beam is initially meshed using 330 4-noded quadrilateral elements. The crack-tip elements are replaced by a SBFEM rosette with 55 nodes. Fig. 17a–c show the predicted load–displacement curves at point A, load-CSD and load-COD curves, respectively. The experimental data and numerical modelling results [37] are also shown

E.T. Ooi, Z.J. Yang / Comput. Methods Appl. Mech. Engrg. 199 (2010) 1178–1192

1189

Fig. 21. Predicted crack propagation process of DEN beam.

for comparison. The present method predicts a peak load F u ¼ 61:26 kN, which is close to F u ¼ 63:05 kN (averaged from four specimens) as reported in the experiment. The pre-peak response of the load–displacement curve is also close to the experimental data. A sharp snapback is predicted once the peak load is reached. As in [37], the post-peak response almost coincides with the prepeak region, and is different from the experimental data. An explanation for this difference is given in [39]. Both predicted load-CSD and load-COD curves are in good agreement with the experimental data. The simulated crack propagation process is shown in Fig. 18, which is similar to that observed in the experiment.

4.4. Cohesive crack propagation: double-edge notched concrete beam A double-edge notched concrete beam (DEN) subjected to fourpoint shear loading (Fig. 19) is simulated. The material properties are: E = 27 GPa, m ¼ 0:1; ft ¼ 2 MPa and Gf ¼ 0:1 N=mm. A bilinear r-COD softening law: f1 ¼ 0:67 MPa;w1 ¼ 0:04 mm and wc ¼ 0:18 mm is used to model the cohesive cracks. The hybrid mesh has 672 elements and two SBFEM rosettes, each with 55 nodes. Fig. 20 shows that the predicted load–displacement responses using the hybrid method agree well with the experiment [41] and other numerical results [30,41]. The crack propagation process

1190

E.T. Ooi, Z.J. Yang / Comput. Methods Appl. Mech. Engrg. 199 (2010) 1178–1192

CL

P/2

200

10 10 10 A B C 100

200

200

200

CL

P

P/2

100

100

200

200

A 100 160

100 100

10 B 120

100

CL

CL

30

200

BEAM 2

BEAM 1

P/2

200

50 C

10 10 A B

P/2

P/2 200

30

30 A

C 160

160 100

BEAM 3

100 160

100 100

10 B 120

P/2 200

20 C 160

160 100

BEAM 4

Fig. 22. Multiple-notched concrete beams [42–44].

Fig. 23. Load-central deﬂection curves for multiple-notched concrete beams.

is characterised by two anti-symmetric cracks propagating from the notches towards the loading points. The cracks reach approximately half the depth of the beam at the peak load F u ¼ 43:53 kN. This is followed by a snapback and unstable crack propagation. The

predicted crack propagation process is shown in Fig. 21, which resembles the experimental observation. It is noted that unlike the method developed in [21], the present method does not require cumbersome partitioning of subdomains

E.T. Ooi, Z.J. Yang / Comput. Methods Appl. Mech. Engrg. 199 (2010) 1178–1192

1191

Fig. 24. Final crack patterns of the multi-notched concrete beams.

or repositioning of scaling centres during crack propagation as the rosettes are rebuilt after each crack propagation step. 4.5. Cohesive crack propagation: multiple-notched simply-supported beams Fig. 22 show four simply-supported and notched concrete beams that were modelled by Shi et al. [42–44] to validate their FEM-based discrete crack model. The material properties are: E ¼ 27:5 GPa; m ¼ 0:2; ft ¼ 2:8 MPa and Gf ¼ 0:1 N=mm. A bilinear softening r-COD relation with f1 ¼ ft =4; w1 ¼ 0:75Gf =ft and wc ¼ 5Gf =ft is used to model the cohesive cracks. Beam 1 and Beam 2 are meshed with 480 4-noded quadrilateral elements (1394 DOF), and Beam 3 and Beam 4 with 495 elements (1424 DOF). SBFEM rosettes with 55 nodes in the cracked subdomain are used in all the beams. Fig. 23a–d compare the load-central deﬂection curves predicted by the present method with that in [42–44]. The hybrid meshes of the beams at failure are shown in Fig. 24a–d. In Beam 1, the central crack B is dominant. It reaches half the depth of the beam at peak load and then propagates unstably throughout the depth of the beam. The secondary cracks at notches A and C cease to propagate after the peak load is reached. In Beams 2 and 3, the cracks at notches A and C compete with each other during the pre-peak stage. At the peak load, the lengths reach approximately half the depth of the beam. In the case of Beam 2, crack A then becomes dominant and continues to propagate unstably while crack C stops propagating. In Beam 3, Crack A then ceases to propagate while crack C becomes dominant. The central crack B in Beams 2 and 3 stop propagating at the early stages of loading. In Beam 4, crack A dominates the entire cracking process while cracks B and C stop growing after the peak load is reached. The crack propagation processes and ﬁnal crack patterns of all the four beams agree well with the FEM simulations reported in [42–44]. It is noted that the simple crack propagation criterion Eq. (53) is sufﬁcient in predicting the complicated crack propagation processes in all the beams. 5. Conclusions A hybrid method that combines the FEM and the SBFEM has been developed for modelling linear elastic and nonlinear crack propagation problems in quasi-brittle materials. The simple local remeshing procedure is based on the FE mesh only and the crack-tip singularities are modelled using SBFEM rosettes that replace the crack-tip elements after remeshing. The hybrid method thus exploits both the ﬂexibility of FEM in remeshing for multiple cracks and the SBFEM’s high efﬁciency and accuracy in extracting SIFs from the semi-analytical solutions of displacements or stres-

ses. It has been successfully applied to model various concrete fracture problems, demonstrating its unique advantages over conventional FEM procedures. Acknowledgement This research is funded by the EPSRC UK (Project No: EP/ F00656X/1). References [1] Y.R. Rashid, Analysis of prestressed concrete pressure vessels, Nucl. Engrg. Design 7 (1968) 334–344. [2] D. Ngo, A.C. Scordelis, Finite element analysis of reinforced concrete beams, J. Am. Concrete Inst. 64 (1967) 154–163. [3] M. Suidan, W.C. Schnobirch, Finite element analysis of reinforced concrete, ASCE J. Struct. Division 99 (1973) 2109–2122. [4] R.J. Cope, P.V. Rao, L.A. Clark, P. Norris, Modelling of reinforced concrete behaviour for ﬁnite element analysis of bridge slabs, in: C. Taylor, E. Hinton, D.R.J. Owen (Eds.), Numerical methods for non-linear problems, Pineridge Press, Swansea, 1980. [5] H.R. Lotﬁ, P.B. Shing, Analysis of concrete fracture with an embedded crack approach, in: H. Mang, N. Bic´anic´, R. de Borst (Eds.), Computational Modelling of Concrete Structures, Pineridge, Swansea, 1994. [6] J. Oliver, Modelling strong discontinuities in solid mechanics via strain softening constitutive equations. Part 1: Fundamentals, Part 2: numerical simulation, Int. J. Numer. Methods Engrg. 39 (1996) 3575–3624. [7] J. Alfaiate, A. Simnoe, L.J. Sluys, Non-homogeneous displacement jumps in strong embedded discontinuities, Int. J. Solids Struct. 40 (2003) 5799–5817. [8] N. Moës, J. Dolbow, T. Belytschko, A ﬁnite element method for crack growth without remeshing, Int. J. Numer. Methods Engrg. 46 (1999) 131–150. [9] N. Sukumar, N. Moës, B. Moran, T. Belytschko, Extended ﬁnite element method for three-dimensional crack modelling, Int. J. Numer. Methods Engrg. 48 (2000) 1549–1570. [10] T. Hettich, A. Hund, E. Ramm, Modelling of failure in composites by X-FEM and level sets within a multiscale framework, Comput. Methods Appl. Mech. Engrg. 197 (2008) 414–424. [11] R.S. Barsoum, Triangular quarter-point elements as elastic and perfectlyplastic crack tip elements, Int. J. Numer. Methods Engrg. 11 (1977) 86–98. [12] J.A. Teixeira de Freitas, Z.Y. Ji, Hybrid-Trefftz ﬁnite element formulation for simulation of singular stress ﬁelds, Int. J. Numer. Methods Engrg. 39 (1996) 281–308. [13] M. Stern, E.B. Becker, A conforming crack tip element with quadratic variation in singular ﬁelds, Int. J. Numer. Methods Engrg. 12 (1978) 279–288. [14] S.E. Benzley, Representation of singularities with isoparametric ﬁnite elements, Int. J. Numer. Methods Engrg. 8 (1974) 537–545. [15] P. Tong, T.H.H. Pian, S.J. Larsy, A hybrid-element approach to crack problems in plane elasticity, Int. J. Numer. Methods Engrg. 11 (1977) 180–184. [16] B.L. Karihaloo, Q.Z. Xiao, Accurate determination of the coefﬁcients of elastic crack tip asymptotic ﬁeld by a hybrid crack element with p-adaptivity, Engrg. Fract. Mech. 68 (2001) 1609–1630. [17] D. Zeng, N. Katsube, J. Zhang, W. Soboyejo, Hybrid crack-tip element and its applications, Finite Elements Anal. Design 28 (2002) 319–335. [18] Q.Z. Xiao, B.L. Karihaloo, X.Y. Liu, Direct determination of SIF and higher order terms of mixed mode cracks by a hybrid crack element, Int. J. Fract. 125 (2004) 207–225. [19] Z.J. Yang, Fully automatic modelling of mixed-mode crack propagation using scaled boundary ﬁnite element method, Engrg. Fract. Mech. 73 (2006) 1711– 1731.

1192

E.T. Ooi, Z.J. Yang / Comput. Methods Appl. Mech. Engrg. 199 (2010) 1178–1192

[20] Z.J. Yang, A.J. Deeks, Fully automatic modelling of cohesive crack growth using a ﬁnite element-scaled boundary ﬁnite element coupled method, Engrg. Fract. Mech. 74 (2007) 2547–2573. [21] E.T. Ooi, Z.J. Yang, Modelling multiple cohesive crack propagation using a ﬁnite element-scaled boundary ﬁnite element coupled method, Engrg. Anal. Boundary Elements 33 (2009) 915–929. [22] J.P. Wolf, C. Song, Finite-Element Modelling of Unbounded Media, John Wiley and Sons, Chichester, 1996. [23] A.J. Deeks, J.P. Wolf, A virtual work derivation of the scaled boundary ﬁniteelement method for elastostatics, Comput. Mech. 28 (2002) 489–504. [24] M. Xie, Finite Element Modelling of Discrete Crack Propagation. Ph.D. Thesis, University of New Mexico, 1995. [25] K.J. Bathe, Finite Element Procedures, Prentice-Hall, New Jersey, 1996. [26] S.R. Chidgzey, A.J. Deeks, Determination of coefﬁcients of crack tip asymptotic ﬁelds using the scaled boundary ﬁnite element method, Engrg. Fract. Mech. 72 (2005) 2019–2036. [27] C. Song, J.P. Wolf, Semi-analytical representation of stress singularities as occurring in cracks in anisotropic multi-materials with the scaled boundary ﬁnite-element method, Comput. Struct. 80 (2002) 183–197. [28] M.L. Williams, On the stress distribution at the base of a stationary crack, J. Appl. Mech. 24 (1957) 109–114. [29] I. May, Y.A. Duan, A local arc-length procedure for strain softening, Comput. Struct. 64 (1997) 297–303. [30] Z.J. Yang, J.F. Chen, Fully automatic modelling of cohesive discrete crack propagation in concrete beams using local arc-length methods, Int. J. Solids Struct. 41 (2004) 801–826. [31] M. Xie, W.H. Gerstle, Energy-based cohesive crack propagation modelling, ASCE J. Engrg. Mech. 121 (1995) 1349–1358. [32] C. Hua, An inverse transformation for quadrilateral isoparametric elements: Analysis and application, Finite Elements Anal. Design 7 (1990).

[33] H. Ewalds, R. Wanhill, Fracture Mechanics, Edward Arnold, New York, 1989. [34] ABAQUS 6.7, User documentation, Dessault Systems, 2007. [35] J. Yau, S. Wang, H. Corten, A mixed-mode crack analysis of isotropic solids using conservation laws of elasticity, J. Appl. Mech. 47 (1980) 335–341. [36] W.K. Wilson, Combined-mode Fracture Mechanics. Ph.D. Thesis, University of Pittsburgh, 1969. [37] A. Carpinteri, S. Valente, G. Ferrara, G. Melchiorri, Is mode II fracture energy a real material property?, Comput Struct. 48 (1993) 397–413. [38] G. Ferrara, P.A. Morabito, A contribution of the holographic interferometry to studies on concrete failure, in: S.P. Shah, S.E. Swartz, B. Barr (Eds.), Fracture of Concrete and Rock, Elsevier Applied Science, Amsterdam, 1989. [39] F. Barpi, S. Valente, Size-effects induced bifurcation phenomena during multiple cohesive crack propagation, Int. J. Solids Struct. 35 (1998) 1851–1861. [40] P.E. Petersson, Crack growth and development of fracture zone in plain concrete and similar materials. Report TVBM-1006, Lund Institute of Technology, Lund Sweden, 1981. [41] P. Bocca, A. Carpinteri, S. Valente, Size-effects in the mixed mode crack propagation: softening and snap-back analysis, Engrg. Fract. Mech. 35 (1990) 159–170. [42] Z. Shi, M. Ohtsu, M. Suzuki, Y. Hibino, Numerical analysis of multiple cracks in concrete using the discrete approach, ASCE J. Struct. Engrg. 127 (2001) 1085– 1091. [43] Z. Shi, M. Suzuki, M. Nakano, Numerical analysis of multiple discrete cracks in concrete dams using extended ﬁctitious crack model, ASCE J. Struct. Engrg. 129 (2003) 324–336. [44] Z. Shi, M. Suzuki, Numerical studies on load-carrying capacities of notched concrete beams subjected to various concentrated loads, Const. Build. Mater. 18 (2004) 173–180.

Copyright © 2021 COEK.INFO. All rights reserved.