**TECHNICAL PAPER**

**Modelling of the fracture process zone to improve the crack propagation criterion in concrete**

**S Shahbazpanahi; A A A Ali; F N Aznieta; A Kamgar; N Farzadnia**

**ABSTRACT**

**Keywords:** cohesive, crack, energy, model, stiffness

**INTRODUCTION**

The mechanical behaviour of ductile materials is different from that of quasi-brittle materials. Cracks grow in ductile materials such as metals due to the intersection and coalescence of micro-voids, while in quasi-brittle materials such as concrete, cracks propagate when the aggregates interlock or when micro-crack bridging occurs (Yang & Liu 2008).

In fracture mechanics a crack is assumed to start when there is a notch or a stress concentration in the tension zone. Linear elastic fracture mechanics (LEFM) was first used to study crack propagation during World War II (Esfahani 2007). Later some studies used LEFM to analyse crack propagation in concrete. However, Kaplan (1961) found that LEFM could not be applied to crack problems in normal concrete sizes.

The first model based on nonlinear fracture mechanics in concrete was proposed by Hillerborg *et al* (1976). It was shown that there is a region called the fracture process zone (FPZ) in front of the real crack tip, which is responsible for crack closure (see Figure 1). This significant and relatively large zone contains micro-cracks in matrix-aggregate, gel pores, shrinkage cracks, bridging and branches of cracks that are located ahead of the macro-cracks. Since a significant amount of energy is stored in the FPZ, a crack can have stable growth before the peak load. In addition, the existence of the FPZ accounts for the strain softening behaviour in the stress-crack opening curve that is observed after peak load. In this region interlocking crack surfaces contribute to a gradual decline in stress and prevent sudden failure (Esfahani 2007). The dimension of the FPZ depends on the size of the structure and the length of the initial crack, as well as on the loading and material properties of the concrete. The length of the FPZ is of special interest as compared with its width. The effective modulus of elasticity is reduced when the crack moves from undamaged regions into the FPZ.

]]>

The so-called Griffith energy approach can be used to describe the crack propagation criterion in the fracture process at the crack tip. This approach states that the energy release rate, defined as the amount of the energy stored in the FPZ which is required to form the crack, must be sufficiently larger than the critical fracture energy. Hence to study the crack state, the crack propagation criterion can be defined in terms of the energy release rate.

To simulate the FPZ, Hillerborg *et al* (1976) used cohesive stress, which is a function of crack opening. In the FPZ cohesive stress reaches its maximum at the tip of the crack, which is equal to the tensile strength *f _{t},* and in the critical opening of the crack

*w*it declines to zero. The area under the stress-opening curve is defined as the energy release rate. The cohesive zone model (CZM) was applied by Hillerborg to model the FPZ in normal-size structures using either the nodal force release method or the technique of interface element with zero initial thickness (Ingraffea

_{c}*et al*1984). Since the length of the FPZ,

*l*is of special interest compared to its width, the interface element with zero thickness is preferred (Yang & Liu 2008).

_{p},To model the FPZ, Bazant and Oh (1983) used dummy bands in which micro-cracks were uniformly distributed. This model, called the crack band model (CBM), was simulated using a layer continuum element in the finite element method. As the CBM depends on the width of the element, it was suggested to model only Mode I fracture (Rots & Brost 1987).

The non-local continuum approach is another method that uses the width and the length of the interface element for modelling the FPZ. This approach, however, uses too many degrees of freedom, and for this very reason it is not computationally affordable(Ingraffea *et al* 1984).

From a finite element point of view, to model the FPZ the stiffness of the element should be properly chosen. In practice, compared to undamaged zones, the FPZ has a different stiffness due to micro-cracking, bridging and branching processes, which provide the energy required for the crack growth. In the softening zone, although some resistance is observed, the stress drops dramatically. Parallel use of the constitutive model of the CBM and the cohesive model may lead to a better simulation of the FPZ behaviour in terms of elastic branching and the softening zone. In addition, in order to estimate the nodal force, an accurate constitutive model for normal and shear stresses is needed.

To estimate the crack propagation state in the finite element model, either of two methods are used: the strength-based or the energy-based approach. To achieve higher accuracy, the energy-based approach is often used to simulate the CZM. The energy approach criterion depends on the stiffness matrix, the displacement and the crack geometry (Wua *et al* 2011). Therefore the crack propagation criterion should be modified based on the new stiffness matrix in the FPZ.

In addition, when the FPZ has fully propagated and reaches its maximum length, a stress-free length appears in front of the notch (or macro-crack) behind the FPZ (Wua *et al* 2011). This stress-free length has not been considered in previous research (Bocca *et al* 1991; Xie & Gersle 1995; Prasada & Krishnamoorthy 2002; Yang & Liu 2008; Shi 2009; Ooi & Yang 2011; Guo *et al* 2012).

When modelling cracks, another issue to deal with is the direction of the crack. The initial direction of the propagation is mostly unknown. Numerous reports have proposed the use of approximate re-meshing algorithms as the crack starts to grow (Xie & Gersle 1995). In these algorithms, a significant number of nodes are created for re-meshing and a large stiffness matrix is created. In this approach, the computational complexity is relatively high. An alternative method is the so-called inter-element boundaries technique, which identifies the crack path (Alfaiate *et al* 1997).

*et al*(1997) instead of re-meshing. In this method, the crack propagation direction is identified by following the inter-element boundaries. The proposed approach is capable of modelling the mixed mode of the crack, which is described by the stress versus the crack opening displacement curve. Three examples are modelled to validate the criterion: a plain concrete beam with a notch shear crack, a notched reinforced concrete beam and a reinforced concrete beam with simple support are analysed and compared to recent experimental results and previous modelling.

**NUMERICAL MODEL**

**Stiffness of interface element**

As mentioned earlier, the FPZ has a softening behaviour due to the interlock of aggregates and micro-cracks. Thus different stiffness characteristics are used in the finite element method to model the FPZ. A single four-node interface element with linear variation of the crack opening displacement (COD) is shown in Figure 2. To model the mixed mode, the relationship between the stress and the displacement is represented by matrix *D _{s},* which is given as:

]]> where

*D*and

_{nn}*D*are the normal and tangential stiffness respectively, and

_{tt}*D*represents the additional stiffness due to the interaction between the shear and the normal stiffness in the fracture mechanism.

_{nt}To estimate the normal stiffness, the normal stress versus COD curve was used. Figure 3 illustrates the concrete COD due to normal stress.

The total opening can be separated into two components:

where *dw _{e}* and

*dw*are the opening elastic and the opening softening, respectively.

_{s}The softening parameter is defined as:

]]>This parameter is in fact the slope of stress in the softening portion of the curve and its value is negative. If *E _{e}* and

*E*are the slope of the elastic zone and the slope of the softening zone at the

_{s}^{i}*i*th iteration of the nonlinear solution respectively, we can write:

Eq (4) implies that the softening parameter changes due to changes in *E _{s}^{i}.* On the other hand, the normal stiffness can be expressed as:

(5) where the *L* is the interface element length. From Eq (4) and (5) we obtain the following:

From the study done by Yoshikawa *et al* (1989), shear stiffness and the additional stiffness due to the interaction between the shear and normal stiffnesses are represented as a linear function of the normal stiffness as:

(8) where *µ* and *β* are a frictional coefficient and dilatency factor respectively. Although these two parameters are functions of the normal stiffness, in this paper, for the sake of simplicity, they are assumed to be constant. The stiffness of the interface element *K _{s}* in the softening zone is given by:

where *B* is the strain-displacement matrix used to model the interface element (Desai *et al* 1984), *dA* is the differential element of the crack surface area, and the *T* is transpose.

In material nonlinearity, the matrix *D _{s}* changes with different levels of loading, particularly after the peak load, where the stiffness matrix is of interest. The method of nonlinearity is an incremental-iterative technique. In this method, the stiffness matrix is updated at the beginning of each load increment. Small displacements are assumed to implement the nonlinear solution (Prasada & Krishnamoorthy 2002). Based on small displacements, at each step of the loading the tangential stiffness method is used at the ith iteration

*(i*=

*j*+ 1) of the nonlinear solution by using the

*jth*displacements of the interface element. To estimate the slope of the softening zone, an initial displacement, and the elastic slope,

*E*are used (Gerstle & Xie 1992). Corresponding to the initial displacement, the residual force is calculated linearly, and then the correction of the displacement to the trial value is calculated. The stiffness matrix is updated by accumulating the corrections of the displacements. This iterative procedure is continued until the residual force converges to zero. The stress distribution at all nodes is automatically implemented in the finite element program.

_{e},The vector of the cohesive forces in the nodes is given by:

where σ* _{ij}* = [στ]

*is the stress vector due to the proposed constitutive model for normal and shear stresses in the FPZ. By using Gaussian integration, employing linear shape functions, and setting the determinant of the Jacobian matrix equal to L/2, the stiffness of the interface element and nodal cohesive forces can be obtained.*

^{T}The normal stress, a, is obtained from Petersson's constitutive law (1981) according to the on the interface element. For shear stress, τ, Bazant and Gambarova (1980) assumed that opening displacement takes place prior to the occurrence of the slip. However, in the present study, based on the shear stress versus the crack sliding displacement (CSD) curve obtained by Yoshikawa *et al* (1989), it is proposed that when the normal stiffness is in the elastic zone, the tangential stiffness value is very large and constant. However, when the normal stress in the softening part decreases, the tangential stiffness starts to reduce (Figure 4). It seems more logical to assume that sliding starts when the normal stress begins to decrease.

]]>

To estimate the maximum shear stress, *τ**max,* the initial slip, S0, and the critical slip S_{c}, shown in Figure 4, the following method is used. Previous experimental results have shown that the relationship between the slide and opening and displacements can be approximated as a linear relationship (Paulay & Loeber 1974). The relationship between COD and CSD is assumed to be given by:

where *δ**t* and *δ**n* are the slide and opening displacements respectively, the parameter *a* is greater than one, and *δ**n** is a constant. Based on Eq (11), the values of *S*_{0} and *Sc* (Figure 4) are assumed as a(*w*_{0} - *δ*_{n}*) and a(*w*_{c} - *δ*_{n}*) respectively. In addition, from Eq (7) the maximum shear stress, *τ**max* is given by:

**Crack propagation criterion**

The Griffith criterion for crack propagation is:

where //, *U* , *W* and *E _{s}* are the total potential energy, the strain energy, the work done by applying the load, and the surface dissipated energy respectively. Note that is given by (Xie & Gerstle 1995):

*u*and

*N*are the nodal displacement vector and the shape function matrix respectively, the parameter

*t*is the surface load on the surface

*S*and

_{t}*b*is the body force in the

*SE*

volume *V*. We can also write as:

where matrix *M* is the linear shape function matrix which depends on the crack opening displacement and is related to the geometry of the crack (Xie & Gerstle 1995).

In this study, the effect of the softening strain energy is considered in order to determine the strain energy release rate. The strain energy release rate is assumed to be the same in both mixed mode and Mode I (Bocca *et al* 1991). The rate of change in the strain energy in the interface element with respect to area is expressed as:

where σ and *ε* are normal stress and strain respectively. On the other hand, the strain energy density *W _{d}* is written as:

where *ε** ^{e}* and

*ε*

^{s}are the elastic strain and the softening strain respectively. The rate of change in the strain energy density is:

Using the chain rule in partial differentiation, we obtain:

Thus, using Eq (20), we write Eq (19) as:

Thus:

Substituting Eq (14), Eq (15) and Eq (22) into Eq (13) and rearranging the terms, we obtain the following:

]]> The first bracket is the equilibrium equation and is equal to zero. The body forces and surface load are small and can be ignored (second bracket). Hence, using the fact that (Xie & Gerstle 1995):In the finite difference method, in order to estimate the first part of Eq (25), only crack-tip stiffness elements can be used by applying crack extension (Yang & Liu 2008). The second part can be evaluated by using a Gaussian integration method (Xie & Gerstle 1995). Eq (25) can be applied to mixed-mode and multiple-crack fracture problems.

Furthermore, a stress-free region appears in front of the initial notch or macro-cracks when the FPZ length is fully propagated (Wua *et al* 2011). It has been shown that if the crack-opening displacement reaches 3.6 *Gc/ft,* the stress-free region appears in front of the initial notch. *Gc* is the area under the curve (fracture energy) in Figure 1. Thus the stress-free region length is 0.08 times the ligament length, approximately *h* - *α*_{0} where *h* and *α*_{0} are the depth of the beam and the length of the initial notch respectively. In finite element methods the length of the stress-free region is formulated by:

where *n* is the number of elements that fail behind the crack. When the FPZ is fully propagated, *n* elements will be set to zero behind the crack and the crack will grow along the respective elements upon appropriate identification of the cracking direction in each step, as described in the next section.

**Crack propagation direction**

One of the most important aspects in discrete cracking is the direction of propagation. In this investigation, the crack is assumed to follow the existing inter-element boundaries and no re-meshing algorithm is needed to decrease the computational complexity (Alfaiate *et al* 1997). This method has a simple algorithm. It is assumed that the directions of maximum principal tensile stresses, which were automatically known in the program at each step, were perpendicular to the crack propagation direction (Xie & Gerstle 1995). Thus the angle *(**γ**)* is recognised (Figure 5). Crack propagation follows one of the inter-elements (AB) or (AC). It is assumed that the crack will not stop or intersect the main element and that the direction of propagation is perpendicular to the maximum tensile principal stress (Figure 5). There are two possible cases for the crack path: if the orientation angle (γ) is less than 45°, the path of growth is (AC), otherwise it will be (AB). The stiffness matrix, nodal forces and displacements of the interface element are changed from the local system to the global system by using the transformation matrix. Although the crack paths are non-smooth, the ones found with this method are in good agreement with the correct crack path.

]]>

The FEAPpv program code was developed to analyse cracks in concrete (Taylor 2009). Four-node isoparametric elements are used for bulk concrete, for which the material behaviour is considered to be linear elastic. Figure 6 shows the major steps used to solve fracture in the beam in the present numerical model. To model the post-peak curve of the structure, displacement was assumed to be incremented rather than the load, i.e. displacement-controlled numerical analysis was applied. The displacements were classically those for nodes at the crack mouth on the modelling of crack tip behaviour. The nonlinear dynamic relaxation method is implemented to find the load-displacement curve. This method is preferred to other methods such as Newton (Gerstle & Xie 1992).

**NUMERICAL EXAMPLES**

Three benchmark test specimens were analysed to validate the model. Figure 7 illustrates the previously tested plain beam used to simulate mixed-mode fracture (Arrea & Ingraffea 1982). The boundary conditions and material properties are indicated in Figure 7.

The Young's modulus, Poisson's ratio, and tensile strength of concrete were assumed to be 24 800, 0.18 and 4 MPa respectively. The thickness of the beam was 152 mm and the length of the initial notch was 82 mm. The parameter values of fracture were *G _{c} =* 150 N/m,

*w*0.135 mm and

_{c}=*w*0.0001 mm. Four-node isoparametric elements were used for the bulk concrete with linear elastic behaviour, and plane stress was considered as the analysis condition. The values of parameters

_{0}=*µ*

*,*

*β*

*, a*and

*δ*

*were chosen as 1.16, 1.64, 0.5 and 0.1 respectively. The initial mesh (c) is illustrated in Figure 8.*

_{n}*

Figure 9 shows the result of the load versus the crack mouth sliding displacement (CMSD) curve for the beam, the experiment of Arrea and Ingraffea (1982) and the numerical model of Xie *et al* (1995).

The round dots represent the experimental envelope by Arrea and Ingraffea (1982), the results of Xie and Gerstle (1995) are shown in black dashed lines and the results of the proposed model are shown by coloured lines. As seen in Figure 9, the results of the proposed model show good agreement with those of the experimental method. Mesh (a) has 110 interface elements, mesh (b) has 225 interface elements, and mesh (c) has 306 interface elements. Approximate matching of the three curves demonstrates independence of the model from the mesh size and shows fast convergence of the proposed model. It can be seen from Figure 9 that the peak loads are close to each other, although the mesh size changes.

]]> In the elastic part, the results lay close to the midpoint of the experimental results obtained previously. However, the peak load obtained by the numerical method is slightly shifted below the upper limit of the envelope. The difference between the data of the proposed model and the experimental data is inevitable since the behaviour of concrete is assumed to be linear elastic in fracture mechanics, but in fact it is nonlinear plastic; compression crushing has also been ignored. The peak load obtained by the proposed method differs by almost 7% from that of the experimental method. The peak load in the numerical model is considerably different from that of the experimental peak load.It is seen that after the peak load, the curves in the softening zone (up to 57 kN) are closer to the experimental data than to the numerical model (Xie *et al* 1995), which is slightly more brittle thereafter. In the softening zone after 60 kN, the proposed model shows more agreement with the experimental data in terms of ductility. This may be because the stress-free zone in the tip of the notch was not considered in previous models.

Figure 10 shows the predicted crack path in mesh (c) compared with the experimental data. The FPZ propagation elements are shown in red lines, while the stress-free elements are displayed in black. It should be noted that the crack path is a smooth curve, although in this study the crack path is illustrated by straight lines. It can be seen that the predicted crack path in the mesh (c) is very close to the experimental result.

To further validate the proposed numerical method in mode I cracking, the experimental data reported by Prasada and Krishnamoorthy (2002) were chosen. The test arrangement, the boundary conditions and the geometry of the RC beam are illustrated in Figure 11. The Young's modulus, compressive strength, tensile strength and fracture energy were 29 270 MPa, 30.1 MPa, 4.11 MPa, 100 N/m respectively. The yield strength of the steel was 395 MPa. Other parameters were assumed to be the same as those used in the previous example. The bond between the bars and the concrete was assumed to be perfect, i.e. no bond-slip was considered.

The data from the numerical analysis are compared with experimental records and the numerical model information of Prasada and Krishnamoorthy (2002) (Figure 12). The numerical results are in good agreement with the upper limit of the experimental envelop. This could be due to the assumptions in fracture mechanics such as the linear elastic behaviour of concrete and tension cracks which are not available in practice, especially for low loads. In this study, the bond-slip of the steel bars was ignored, unlike in the experimental test and the numerical model by Prasada and Krishnamoorthy (2002). This is why there is a slight over-estimation in the numerical model for a load of about 13.2 kN.

]]>

The FPZ appeared in front of the notch tip at a loading of about 2.5 kN, while deflection at mid-span was 0.012 mm. The FPZ began to grow with increasing deflection at mid-span as the load increased. Since the crack-opening displacement was smaller than 3.6 *GJf _{t},* the FPZ did not fully propagate. In the initial stages, the load was sustained by the FPZ and reinforcement bars were not involved yet. When the load was 13.2 kN, the reinforcement bars arrested the crack. As shown in Figure 12, at a load of 19 kN the stiffness of the beam was slightly reduced. This may be because the stress in the reinforcement bars reached the yield stress. At a load of 23.4 kN, the crack-opening displacement was equal to 3.6

*G*Thus the FPZ completely propagated and a stress-free crack was created in front of the notch tip.

_{c}lf_{t}.Figure 13 shows the FPZ length, the length of the stress-free region and the COMD at the final load. The CMOD is 0.468 mm, and the deflection at mid-span is 0.576 mm at a loading of 28.1 kN. The FPZ propagation reaches to almost three-quarters of the beam depth after the thirteenth loading step. At the tenth loading step the FPZ fully propagates and the stress-free region length appears.

In addition, a reinforced concrete beam with simple supports (Figure 14), which was tested by Bresler and Scordelis (1963), was analysed using the proposed model. The RC beam was 4 572 mm long and 305.8 mm thick. The modulus of elasticity and the Poisson's ratio of concrete are 24 000 MPa and 0.18 MPa respectively. The modulus of elasticity, the Poisson's ratio, the cross-sectional area and the yield strength of steel are 200 GPa, 0.3 MPa, 3 290 mm^{2} and 552 MPa respectively. The tensile strength of concrete is 2.8 MPa and the critical crack opening displacement is 0.152 mm. A two-node truss element with elastic-perfect plastic behaviour and a four-node isoparametric element with linear elastic behaviour were used to model the steel bar and the concrete respectively. To model the symmetry condition, only half of the beam was simulated. The bond-slip between the bars and the concrete was assumed to be perfect. Also, the crack in the RC beam was simulated by a primary crack, first introduced by Ingraffea *et al* (1984), where the reinforcing bar crosses a primary crack.

]]> In this study, load versus deflection at the middle of the beam was compared to experimental results in Figure 15. It can be seen that the shapes of the curves are similar to those of the experimental data. It can be seen that the stiffness of the beam obtained by the proposed method is slightly greater than that of the experimental observation (almost 6%) and the proposed model underestimates the deflection of the beam. However, as the load is increased, the results no longer agree. The reason could be that the proposed model ignores the compression crushing, the nonlinear behaviour of concrete, plastic deformation and the bond-slip of bar concrete.

Initially a few flexural cracks appeared near the mid-span perpendicular to the longitudinal axis when the loading reached about 55 kN. The length of the biggest crack for this load was 190 mm. The width of the flexural cracks was greater than the shear crack width, which occured near the supports. The first flexural crack reached 245 mm in length and 0.168 mm in width at a loading of 100 kN. In one-quarter of the span to mid-span the cracks tended to grow towards the loading point. When the load was about 200 kN, the first crack propagated with a length of about 320 mm and a width of 0.187 mm. A shear crack was observed in the vicinity of the support and its width was greater than that of the flexural crack. Figure 16(a) shows the crack patterns in the experimental study of Bresler and Scordelis (1963) and Figure 16(b) illustrates the crack paths in this study at a loading of around 285 kN. As can be seen, both shear and flexural crack formation resemble the experimental data. Initially, the cracks grow straight and then slowly propagate towards the loading point. The initiation and location of some of these cracks may change due to the mesh size.

]]>

**CONCLUSIONS**

In this study an alternative stiffness matrix was applied to model the FPZ. The relationship between the normal stress and the shear transfer was considered in the implementation of the finite element method. A new constitutive model for shear stress was proposed in a four-node interface element. The energy-based crack propagation criterion was then improved by using a new stiffness matrix. The load-deflection curve in this numerical model and the curve in a previous experimental study were in good agreement. The crack directions at the tensile face in three recent experimental data sets and the present study were close to one another. Several case studies were considered and global load-deflection responses computed with the proposed model being in reasonable agreement with the results found in the literature. Therefore it can be concluded that the model is applicable, as it has been verified computationally that it is sufficiently able to predict the crack pattern.

**REFERENCES**

Alfaiate, J, Pires, E & Martins, J 1997. A finite element analysis of non-prescribed crack propagation in concrete. *Computers and Structures,* 63(1): 17-26. [ Links ]

Arrea, M & Ingraffea, A 1982. *Mixed-mode crack propagation in mortar and concrete.* Department of Structural Engineering, Cornell University, US, Report No. 81-13. [ Links ]

Bazant, Z & Gambarova, P 1980. Rough cracks in reinforced concrete. *Journal of the Structural Division - ASCE,* 106(4): 819-842. [ Links ]

Bazant, Z & Oh, B 1983. Crack band theory for fracture of concrete. *Material and Structures, RILEM,* 16(3): 155-177. [ Links ]

Bocca, P, Carpinteri, A & Valente, S 1991. Mixed-mode fracture of concrete. *International Journal of Solids and Structures,* 27(9): 1139-1153. [ Links ]

Bresler, B & Scordelis, A 1963. Shear strength of reinforced concrete beams. *Journal of the American Concrete Institute,* 60(4): 51-72. [ Links ]

Desai, C S, Zaman, M M, Lightner, J G & Siriwardane, H J 1984. Thin-layer element for interfaces and joints. *International Journal of Numerical Analytical Methods in Geomechanics,* 8(1): 19-43. [ Links ]

Esfahani, M 2007. *Fracture mechanics of concrete.* Tehran, Iran: Tehran Polytechnic Press. [ Links ]

Gerstle, W & Xie, M 1992. FEM modelling of fictitious crack propagation in concrete. *Journal of Engineering Mechanics - ASCE,* 118(2), 416-434. [ Links ]

Guo, X, Su, R K & Young, B 2012. Numerical investigation of the bilinear softening law in the cohesive crack model for normal-strength and high-strength concrete. *Advances in Structural Engineering,* 15(3): 373-388. [ Links ]

Hillerborg, A, Modeer, M & Petersson, P 1976. Analysis of crack formation and crack growth in concrete by means of mechanics and finite elements. *Cement and Concrete Reserch,* 6: 773-782. [ Links ]

Ingraffea, A, Gerstle, W, Gergely, P & Saouma, V 1984. Fracture mechanics of bonds in reinforced concrete. *Journal of Structural Engineering - ASCE,* 110(4): 871-890. [ Links ]

Kaplan, M 1961. Crack propagation and the fracture of concrete. *ACIJournal,* 58(5): 591-610. [ Links ]

Ooi, E T & Yang, Z J 2011. Modelling crack propagation in reinforced concrete using a hybrid finite element-scaled boundary finite element method. *Engineering Fracture Mechanics,* 78(2): 252-273. [ Links ]

Paulay, T & Loeber, P 1974. Shear transfer by aggregate interlock. Shear in reinforced concrete. *Special Publication - ACI,* 42(1): 1-15. [ Links ]

Petersson, P 1981. *Crack growth and development of fracture zone in plain concrete and similar materials.* Division of Building Materials, Lund Institute of Technology, Lund, Sweden: Report No. TVBM-1006. [ Links ]

Prasada, M & Krishnamoorthy, C 2002. Computational model for discrete crack growth in plain and reinforced concrete. *Computer Methods in Applied Mechanics and Engineering,* 191: 2699-2725. [ Links ]

Rots, J & De Borst, R 1987. Analysis of mixed-mode fracture in concrete. *Journal of Engineering Mechanics,* 113(11): 1739-1758. [ Links ]

Shi, Z 2009. *Crack analysis in stuctural concrete: Theory and application.* Burlington, US: Butterworth-Heinemann. [ Links ]

Taylor, R 2009. *FEAPpv source, A finite element analysis program, personal version.* University of California, Berkeley, US. [ Links ]

Wua, Z, Rong, H, Zheng, J & Xu, F 2011. An experimental investigation on the FPZ properties in concrete using digital image correlation technique. *Engineering Fracture Mechanics,* 78(17): 2978-2990. [ Links ]

Xie, M & Gerstle, W 1995. Energy-based cohesive crack propagation modeling. *Journal of Engineering Mechanics - ASCE,* 121(12): 1349-1458. [ Links ] Xie, M, Gerstle, W H & Rahulkumar, P 1995. Energy- based automatic mixed-mode crack-propagation modeling. *Journal of Engineering Mechanics - ASCE,* 121(8): 914-923. [ Links ]

Yang, Z & Liu, G 2008. Towards fully automatic modelling of the fracture process in quasi-brittle and ductile materials: A unified crack growth criterion. *Journal of Zhejiang University Science,* 9(7): 1862-1775. [ Links ]

Yoshikawa, H, Wu, Z & Ta, T 1989. Analytical model for shear slip of cracked concrete. *Journal of Structural Engineering - ASCE,* 115(4): 771-788. [ Links ]

**Contact details: ** Shahriar Shahbazpanahi

Department of Civil Engineering, Universiti Putra Malaysia

43400, Serdang, Malaysia

T: +60 1 2343 1539 ]]> E: sh.shahbazpanahi@gmail.com

**Contact details: **Dato' Abang Abdullah Abang Ali

Housing Research Centre, Universiti Putra Malaysia

43400, Serdang, Malaysia

T: +60 3 8946 6380

E: aaaa@eng.upm.edu.my

**Contact details: **Farah Nora Aznieta Bt Abd Aziz

Department of Civil Engineering, Universiti Putra Malaysia ]]> 43400, Serdang, Malaysia

T: +60 3 8946 4406

E: farah@eng.upm.edu.my

**Contact details: **Alaleh Kamgar

Department of Civil Engineering, Universiti Putra Malaysia

43400, Serdang, Malaysia

T: +6 01 2334 1538

E: halaleh.kam@gmail.com

**Contact details: ]]>
**Nima Farzadnia

Housing Research Centre, Universiti Putra Malaysia

43400, Serdang, Malaysia

T: +60 1 7634 1221

E: nimafarzadnia@gmail.com

SHAHRIAR SHAHBAZPANAHI is a PhD candidate in structural engineering at Universiti Putra Malaysia. He is a senior lecturer in the civil engineering department of the Islamic Azad University (IAU), Sanandaj Branch, in Iran. His areas of interest are fracture mechanics In concrete and the modelling of FRP strengthening.

]]>

PROF DATO' ABANG ABDULLAH ABANG ALI is the current president of the Malaysian Society for Engineering & Technology (mSET) and Director of the Housing Research Centre, Universiti Putra Malaysia. His areas of interest are structural and materials engineering, affordable quality housing, industrialised building systems and engineering education. More than 200 of his papers have been published in various journals and presented at conferences.

DR FARAH NORA AZNIETA BT ABD AZIZ Is Head of the Structural Department, Universiti Putra Malaysia. Her areas of interest are fibre-reinforced concrete, analysis and design of concrete structures and strengthening of reinforced concrete structures. She is Head of the Structural Department, Universiti Putra Malaysia.

]]>

ALALEH KAMGAR (MSc) is a structural engineer with more than eight years of experience in consultation and design. Her area of interest is the analysis and design of concrete frame structures.

NIMA FARZADNIA Is a PhD candidate in structural engineering and a researcher at the Housing Research Centre, Universiti Putra Malaysia. His areas of interest are construction and building materials, concrete properties and strengthening methods.

]]>