Modelling localised fracture of reinforced concrete structures

Recently, localised fracture of reinforced concrete members has been of interest to many researchers and engineers. In the case of utilisation of reinforced concrete structures in offshore structures, it is very important to evaluate the effect of individual crack opening on the corrosion of reinforced steel bars, which could cause a significant strength loss of the structural members. When structures are subjected to extreme loading conditions such as fire, the reinforced concrete beams or slabs are forced into high deformation, and large individual cracks can be formed within the structural members. This phenomenon was observed in some previous experimental tests (Bailey and Toh, 2007; Foster, 2006; Foster et al., 2004). Those large individual cracks influence the exposure condition of reinforcing steel bars to fire and in some cases the steel reinforcements are directly exposed to fire which significantly reduces the fire resistance of the structures. In the case of the fibre-reinforced plastic (FRP) strength reinforced concrete beam, the localised concrete cracking could cause possible debonding failure of externally strengthening FRP layer (Chen et al., 2011). In addition to the global responses like ultimate strength and deflection the propagation of localised cracks in the concrete surrounding the reinforcing steel also affects the bond–slip behaviour between the steel bar and concrete (Cervenka et al., 2003; Gao et al., 2013; Rots, 1985). In the past, plenty of researches have been devoted to develop the numerical models for modelling of the reinforced concrete structure based on a continuum approach in which smeared cracking was adopted to simulate cracks (Huang et al., 2003; Suidan and Schnobrich, 1973; Vecchio and DeRoo, 1995). However, the smeared cracking model is not capable of capturing the formations and propagations of individual cracks within a reinforced concrete member. At present, very little research has been done on the modelling of localised fractures for reinforced concrete structural members.

In the past, a discrete-cracking model has been used successfully for modelling cracks in structural members when the crack path was known in advance. In these cases, the finite-element (FE) mesh was constructed in such a way that the crack path coincided with the element boundaries. However, this approach has to limit cracks to interelement boundaries, or requires performing remeshing during the analysis process in order to capture crack propagation. To overcome this problem, extended FE method (XFEM) was developed and introduced (Belytschko and Black, 1999; Möes et al., 1999) based on the partition of unity theory (Mellenk and Babuŝka, 1996). Recently, a number of XFEM approaches in conjunction with cohesive-zone models have been developed to analyse localisation and fracture in engineering materials (Möes and Belytschko, 2002; Verhoose et al., 2009; Wells and Sluys, 2001; Zi and Belytschko, 2003). In the last decade, XFEM approach has been successfully extended to many applications in modelling two-dimensional (2D) and three-dimensional (3D) problems (Duan et al., 2009; Gravouil et al., 2002; Hansbo and Hansbo, 2004; Möes et al., 2002; Sukumar et al., 2000). However, so far the application of XFEM for modelling of reinforced concrete members is still very limited.

The main objective of this paper is to present a robust FE procedure for modelling the localised fracture of reinforced concrete members. The four-node quadrilateral elements are used for modelling of plain concrete part of reinforced concrete members, in which the XFEM is adopted to simulate the formation and growth of cracks. The reinforcing steel bars are modelled by using the three-node beam–column elements and two-node bond-link elements are employed for modelling the interaction between plain concrete and reinforcing steel bar elements. If XFEM is applied to the plain concrete element its nodal displacement will be divided into two parts, continuous part and discontinuous part, both of them cannot be coincident with the nodal displacement of steel bars. To resolve this problem, a shifted enhancement function proposed by Zi and Belytschko (2003) is adopted so that the total nodal displacement of the plain concrete element can be obtained in the procedure rather than only the continuous part. In this way, the nodal displacement of concrete and reinforcing steel bar can coincide with each other in the proposed model. In this paper, the influence of bond characteristics between the concrete and reinforcing steel bar on the localised crack initiation and propagation within a reinforced concrete member is also examined.

As shown in Figure 1, a reinforced concrete beam is modelled as an assembly of finite plain concrete, reinforcing steel bar and bond-link elements. In the present model, a general 3D isoparametric three-node beam–column element developed by Huang et al. (2009) is adopted to represent the main reinforcing steel bars of reinforced concrete members. The material non-linearities of concrete and steel follows the models specified in Eurocode 2 (BSI, 2004).

To model the interaction between the reinforcing steel and concrete within the reinforced concrete structural members, a two-node bond-link element with zero length developed previously by Huang (2010) was used. As shown in Figure 1(b), the bond-link elements are used to link four-node quadrilateral plain concrete elements with three-node reinforcing steel bar elements. The bond-link element is a specialised two-node element of zero length, which has three translational degrees of freedom u, v, w and three rotational degrees of freedom θ_x, θ_y, θ_z at each node, where x, y, z are local coordinates of the reinforcing steel element in which x is the direction of the longitudinal axis of the reinforcing steel bar (Figure 1(c)). It is assumed that the slip between the reinforcing steel and concrete is only related to the longitudinal axis direction (x-direction).

For reinforcing steel bars, apart from the relative slip along the longitudinal axis direction (x-direction) between the concrete and steel bars the concrete prevents relative movement of reinforcing bars in other directions. It is, therefore, reasonable to assume that common nodes of the concrete and reinforcing bar elements have identical rotations and movements in y and z directions. Hence, in this bond-link element the stiffness coefficients of the element stiffness matrix k₂, k₃, k₄, k₅, k₆, which are related to the degrees of freedom for v, w and θ_x, θ_y, θ_z, respectively, are assumed to have infinite magnitude (10¹⁵). The stiffness coefficient k₁, which is related to the degree of freedom for u, can be calculated using an empirical bond stress–slip relationship given in the CEB-FIP Model Code (CEB-FIP, 1993).

The bond-link element is capable of modelling perfect, partial and zero bond between the concrete and reinforcing steel bar within the reinforced concrete structures. For partial bond condition (ribbed steel bar and smooth steel bar), the average bond stress (τ_b) and stiffness coefficient k₁ of the bond element are calculated using the bond stress–slip curves proposed in the CEB-FIP Model Code (CEB-FIP, 1993). For perfect bond condition k₁ = 10¹⁵(N/mm) and no failure of the bond are assumed. For the zero bond condition, it is assumed that there is no interaction between the concrete and reinforcing steel bars (k₁ = 0, τ_b = 0). The details of the bond-link element can be found in Huang (2010).

As shown in Figure 1(b), the quadrilateral plain concrete elements and the steel bar elements are connected to each other using two-node bond-link elements. The displacement interpolation in a four-node quadrilateral element is linear (see the following sections) while the displacement interpolation in the three-node beam element, which represents reinforcing steel bar, is quadratic. Therefore, there would be a displacement incompatibility between the plain concrete and reinforcing bar elements. However, in the current model, the two different elements are not sharing the same nodes (Figure 1(b)); this will reduce the impacts of displacement incompatibility on the model. Moreover, it can be seen from Figure 1(b) that two plain concrete elements are connected to one steel bar element by three bond-link elements. Hence, the mesh of the quadrilateral concrete element has been doubly refined compared with the steel bar element, which may help to reduce the effect of displacement incompatibility between the plain concrete and the reinforcing bar elements on the modelling results. Generally, for FE modelling, if different elements with displacement incompatibility are used the fine mesh is needed to improve the accuracy of the results. It will be demonstrated in the validation section that the current model is not very mesh sensitive. Hence, the influence of the displacement incompatibility of two different elements on the current model is limited.

To transform the 3D problem into a 2D model, all reinforcing steel bars at the same height within the cross-section of the beam are represented by an equivalent steel bar element at the reference plane, as shown in Figure 1(a). The cross-sectional area of the equivalent steel bar element equals the total cross-sectional area of all real steel bars represented. Besides, the contact area of the equivalent steel bar element with concrete is the summation of the total contact areas of all the real steel bars represented.

A general isoparametric four-node quadrilateral element, each node of which contains two translational degrees of freedom, is used herein for representing the plain 2D concrete element (Bathe, 1996). For modelling the concrete cracking, the XFEM is incorporated into the plain concrete element. The key idea of the XFEM is to use the partition of unity to describe the discontinuous displacement and the displacement field is approximated by the sum of the regular displacement field and the enhancement displacement field (Möes et al., 1999). Considering a four-node quadrilateral element which is crossed by a crack (Γ_d), as shown in Figure 2, the domain is divided into two distinct domains represented as Ω⁺ and Ω⁻. Then, the total displacement field u consists of a continuous regular displacement field u_cont and a discontinuous displacement field u_dis (Wells and Sluys, 2001)

where N_i is the shape function, u_i the regular node displacement, a_i the additional node displacement to describe the discontinuity and Ψ_i (x) the enrichment function to realise the displacement jump over discontinuity and is expressed as

in which sign is the sign function and is defined as

Note that the sign function enrichment is equivalent to the Heaviside function enrichment which has been used in many previous works of XFEM, but the sign function appears to be more symmetrical than the Heaviside function H (H = 1 if x ∈ Ω⁺, H = 0 if x ∈ Ω⁻) (Zi and Belytschko, 2003). sign (x_i) is the sign function of the ith node in a quadrilateral element. Taking the quadrilateral element in Figure 2 as an example: sign (x₂) = sign (x₃) = − 1 for nodes 2 and 3 and sign (x₁) = sign (x₄) = + 1 for nodes 1 and 4, respectively. Compared with the conventional XFEM models, the sign function given in Equation 2 is shifted by sign (x_i). According to Zi and Belytschko (2003), using the shifted sign function can make the enrichment displacement field vanish outside the enhanced element but not alter the approximating basis. In this way, only the element cut by the crack needs to be enhanced rather than its adjacent elements which also contain the enhancement nodes. This method could greatly simplify the implementation of extended FE model in modelling the multiple cracks in a reinforced concrete structure. This is because the interaction between the steel bar and the concrete by way of bond action tends to make cracks to distribute more evenly within the structural members, unlike most of previous XFEM simulations for modelling single crack propagation within brittle materials, such as plain concrete. In addition, a more crucial advantage of using shifted enhancement function is capable of obtaining the total nodal displacement directly in the procedure rather than only in the regular part of XFEM displacement (Ahmed, 2009). This makes the compatibility of total nodal displacements of the plain concrete element and the steel bar element feasible. Hence, the bond-link element can be used to link plain concrete element and steel bar element in a conventional way, such as continuous approach, for modelling of localised cracking within a reinforced concrete member.

In many applications of the XFEM, crack tip enrichment functions such as branch functions were used to enrich the element that contains the crack tip (Xiao et al., 2007). In this paper, for simplicity, it is assumed that the crack tip is always located at the edge of an element; thus, the cracked element can be successfully enriched by the sign function only without other enrichment functions. Thus, the crack branching is not included in the proposed model. It is assumed that a particular element contains only one crack. The main purpose of the model developed in this paper is to capture the major localised cracks within the concrete beam. Therefore, in order to enhance the computational efficiency of the proposed model, precise modelling of the crack tip and crack branching is not considered in this paper. It will be demonstrated in the validation section that the current model is not very mesh sensitive.

In the case of 2D four-node quadrilateral element as shown in Figure 2, the element nodal displacement vector û can be represented as

where u_i and ν_i are the regular nodal displacements related to x and y coordinates, respectively, and a_i and b_i are enhancement nodal displacements related to x and y coordinates, respectively.

The total strain within the element in which the enhancement degrees of freedom are included can be expressed as

in which ε_cont is the continuous strain and ε_dis the discontinuous strain; B_sta^u the standard strain–displacement transformation matrix corresponding to the regular degrees of freedom u_i and B_enr^a the enrichment strain–displacement transformation matrix corresponding to the additional enhancement degrees of freedom a_i.

The shape functions of a general quadrilateral element are calculated as (Bathe, 1996)

Hence, the strain–displacement transformation matrix B including the regular part and the enhancement part can be obtained as B = [B_sta^uB_enr^a], in which

where Ψ_i(x) (i = 1 − 4) is the enrichment function given in Equation 2, and the matrix L contains the differential operators. If strains are reasonably small the stresses σ can be written as

in which D is the constitutive matrix of concrete.

In a FE model, the equilibrium conditions between the internal and the external ‘force’ has to be satisfied. To form the element stiffness matrix and internal force vector, the virtual work equation without body forces reads

where f^int is the internal force vector and f^ext is the external force vector. As shown in Figure 2, for the element with crack Equation 10 can be written as

The internal force f^int contains the regular part f_u^int, enhancement part f_a^int and the traction part f_Γ^int in which the regular internal force f_u^int balances the external force f^ext and the enhancement part f_a^int is related to traction of crack f_Γ^int only, that is

where t_a is the traction acting on the discontinuity Γ_d (Figure 3) and can be written as

where t_an and t_as are the traction normal and tangential to a crack, respectively; w_n (⁠$=n⋅(u+−u−)=n⋅2∑14Niai$⁠) is the crack opening of the normal direction n of the discontinuity (Γ_d), and w_s is the slide of two crack faces tangential to the discontinuity (Γ_d), respectively and T_an is the tangent stiffness of the traction-separation law. Using the principle of virtual work on Equations 12 and 13, the element stiffness matrix in terms of incremental displacements can be obtained

Equation 15 can be rewritten as

where K_uu is the element stiffness matrix with reference to the regular degrees of freedom; K_aa is the element stiffness matrix with reference to the enhancement degrees of freedom; K_ua = K_au^T is related to both; and K_Γ is the element stiffness matrix of traction

where t is the thickness of the element, $N¯=N+−N−=2(N)$ and N is the shape function matrix which is defined in Equations 6 and 7. The orthogonal transformation matrix O performs the transformation of the local orientation of the discontinuity to the global coordinate system, as shown in Figure 3 

where cos(x, n) is the cosine of the angle between the x-axis and the normal orientation of discontinuity, cos(x, s) is the cosine of the angle between the x-axis and the tangential orientation of discontinuity.

Using the principle of virtual work, the internal force vectors in Equation 16 can be written as

Before cracking or crushing occurs, the integral concrete is assumed to be isotropic, homogeneous and linearly elastic. A biaxial concrete failure envelope proposed by Barzegar-Jamshidi (1987), which was based on a slight modification of the Küpfer and Gerstle (1973) expressions, was adopted to determine the cracking and crushing of concrete. This is illustrated in Figure 4. Within this model, the initiation of cracking or crushing at any location occurs when the concrete principal stresses reach one of the failure surfaces. Before cracking or crushing, a linear elastic constitutive model is adopted to describe the material property of concrete. The constitutive matrix D in Equation 9 is the usual set of elastic constant. After crushing, concrete is assumed to lose all stiffness. After cracking, the constitutive model of cracked concrete in this research is based on the cohesive crack concept. The linear elastic material properties are still assumed in the continuous solid, but the enhancement internal force which is related to the traction over the crack would decrease with the increase of the crack opening. A concrete bilinear softening curve is used in the current model, as shown in Figure 5, where f_t is the tensile strength of concrete, t_dn the traction normal to crack and w the crack opening. When the crack opening exceeds the traction-free open width (0·68w_ch) the tangent stiffness is set to be zero. In the cohesive interface, the softening curve is governed by the fracture energy, and the fracture energy formulation proposed by Bazant and Becq-Giraudon (2002) is adopted in this research.

In the current model, after a certain concrete element is cracked, there are two kinds of element nodes: regular nodes and enhancement nodes. The regular node contains only regular degrees of freedom, but the enhancement function should be added to the enhanced nodes that support the element crossed by a crack. This is illustrated in Figure 6, where enhanced nodes are indicated by solid circles and regular nodes by hollow circles.

As mentioned in the previous sections, since the enhancement function (sign function) referenced to enhancement nodes is shifted by sign (x_i) the enhanced shape function of the enhancement nodes only influence the displacement field inside the element crossed by the crack. Therefore, only the element enclosing the crack needs to be enhanced rather than all the elements that contain enhanced nodes. The enhanced elements shown in Figure 6 are filled in grey colour, and for other regular elements the element stiffness matrix K contains the regular part K_uu (Equation 17) only.

In the past, for the majority of the extended FE simulations for modelling structural members of the brittle material only one main crack was allowed to form and propagate. This is reasonable to model plain concrete. However, as far as the reinforced concrete is concerned, the reinforcing steel bar can undertake the tensile force transferred from the cracked concrete to other parts of uncracked concrete within the member through bond interaction between the steel bar and concrete. Therefore, several cracks may exist within the structure member. The model developed in this paper allows two or more cracks initiating and propagating within a reinforced concrete member with the following assumptions: only one crack may exist within a particular element and the nodes related to the cracked element are enhanced once only.

In this research, the crack initiation and growth depend on the principal tensile stresses of the concrete element. At each load increment every concrete element is examined one by one and the average principal tensile stress of all Gauss points in an element is checked to determine whether the concrete is cracking. Once the average concrete principal tensile stress reaches one of the failure surfaces, either in the biaxial tension region (segment AB) or in the combined tension–compression region (segment BC), as shown in Figure 4, a straight crack is inserted through the entire element and the orientation of crack is normal to the maximum tensile principal stress. The initial crack is set to go through the centroid point of a quadrilateral element, and then when the average principal stress of the next element reaches one of the tension failure surfaces, the crack will propagate from the tip of the existing crack into the next element by the orientation normal to the corresponding maximum tensile principal stress. Figure 7 illustrates how a crack initiates and propagates. As can be seen, there are two different possible ways that an initial crack cut a quadrilateral element, initial crack 1 in Figure 7(a) and initial crack 2 in Figure 7(b), each of which has possibly three crack propagation paths within the next element when the initial crack extends from element 1 to element 2. The different possibilities of crack propagation in Figure 7 are all included in the current FE model so that various crack patterns can be simulated in a reinforced concrete member (like bending cracks and shear cracks). A crack propagating depends on the location of the existing crack tip and the orientation of the propagated crack.

In this study, Gauss quadrature is employed to calculate the stiffness matrix and the internal force vector of the concrete element. All stresses, strains and the constitutive matrix of the material discussed are related to Gauss integration points. For the regular four-node element, four Gauss integration points are used as recommended by Bathe (1996). For those enhanced elements containing a crack, integration should be performed separately on both sides (Ω⁺ and Ω⁻) of the crack, respectively, that means sign function sign (x) needs to be applied for each Gauss point within the element.

To properly integrate the field (Ω⁺ and Ω⁻) on both sides of the crack, the enhanced elements are partitioned into sub-triangles where usual Gauss quadrature could be used. Figure 8 shows two possible integration schemes. It can be seen that within each sub-triangle, three Gauss integration points are used, and correspondingly each enhanced element employs 24 Gauss integration points to perform integration over the whole element. As shown in Figure 8(a), a quadrilateral element is cut by a crack into two sub-quadrilaterals; four sub-triangles with 12 Gauss points are applied within each sub-quadrilateral. Figure 8(b) shows that a quadrilateral element is cut by a crack into a pentagon and a triangle. In this case, 15 Gauss points distributed over five sub-triangles are applied within the pentagon and 9 Gauss points distributed over three sub-triangles are applied within the triangle.

After the element is partitioned, the integration of the enhanced element is performed at 24 Gauss points distributed over eight sub-triangles, each of which 3 Gauss integration points are applied, that is

where (r_i^j, s_i^j) are the integration point coordinates of the ith Gauss point within jth sub-triangle, and w^j is the corresponding weighting factor. Bathe (1996) provided the integration point coordinates (r′_i, s′_i)_i=1,2,3 and weight factor w′ of a three-point integration over a triangular domain, where r′₁ = r′₃ = 1/6, r′₂ = 2/3, s′₁ = s′₂ = 1/6, s′₃ = 2/3 and w′ = 1/3. However, because the matrices B_sta^u and B_enr^a are represented by natural coordinates (r, s) of the quadrilateral element, the integration of the stiffness matrix and the internal force has to be performed with reference to the natural coordinate system of the quadrilateral element. Therefore, a coordinate transformation should be conducted for changing (r′_i, s′_i) to (r_i^j, s_i^j)

where (r₁^sub,j, s₁^sub,j), (r₂^sub,j, s₂^sub,j) and (r₃^sub,j, s₃^sub,j) are the coordinates of three vertices of the jth sub-triangle in the natural coordinate system of the quadrilateral element; N_i^tri-1(= 1 − r′_i − s′_i), N_i^tri-2(=r′_i) and N_i^tri-2(=s′_i) are the shape functions of the three-node triangular element recommended by Bathe (1996), which are represented by (r′_i, s′_i). The weight factor w′ should also be transformed to w^j with reference to the nature coordinate system of the quadrilateral element as w^j=2w′a^sub,j, in which a^sub,j is the area of the jth sub-triangle and can be expressed as

In the current FE model, a crack is represented by a straight line within the enhancement element, thus two Gauss points are employed to integrate the discontinuity terms (K_Γ and f_Γ^int) using a one-dimensional integration scheme, shown in Figure 3.

Due to the high non-linearity of the current model, a full Newton–Raphon solution procedure is adopted. In the current model, a reinforced concrete member is modelled as an assembly of finite plain concrete, reinforcing steel bar and bond-link elements (Figure 1). Hence, after multiple cracks formed within the plain concrete elements the forces originally resisted by the plain concrete elements can be transferred into the reinforcing steel bar element through bond-link elements. Hence, the modelling behaviour of the reinforced concrete member is much smooth compared with modelling of plain concrete elements only. The full Newton–Raphon solution procedure used in the current model is robust enough to deal with these convergence problems. The analysis can be performed until the fracture of the reinforcing steel bars or the failure of bond-link elements.

To demonstrate the capability of the model described above for modelling the individual crack within a reinforced concrete beam, a simply supported reinforced concrete beam subjected to transverse point load at the centre was modelled. Figure 9 shows the details of the modelled beam. The compressive strength of concrete is 40 MPa and the yield strength of longitudinal reinforcing steel bar is 460 MPa. To simplify the analysis, the stirrups were ignored in this numerical example. Four real reinforcing steel bars were represented by two equivalent steel bar elements at the reference plane. In this example, the compressive failure of concrete was not taken into account in order to generate a large deformation, so that the localised cracking within the beam can be demonstrated more obviously. In this case, a perfect bond between steel bars and concrete was assumed. The analysis was performed using full Newton–Raphon solution procedure under load control. The total point load was 150 kN and was applied at the mid-span of the beam. The analysis was completed by 94 load steps.

Figure 10 demonstrates the process of crack propagations within the beam. It can be seen that the predicted flexural cracks distributed near the mid-span and the cracks were inclined near the supports. Generally, the predicted cracking pattern is reasonable, and it is evident that the model is capable of modelling the different kinds of cracks within a reinforced concrete beam. The cracks were able to propagate freely within the beam depending on the states of principal stresses within the element. As shown in Figures 10(b) and 10(c), after steel bar yielding the crack opening of the element at mid-span is significantly larger than the others. As the deformed mesh shown in Figure 11, the mid-span element is significantly distorted due to a big crack within the mid-span zone. It is evident that the numerical model developed in this paper is capable of modelling the big individual cracks and capture the localised cracking within a reinforced concrete beam. Figure 12 illustrates the predicted history of opening of mid-span crack in which the crack opening increased significantly after the steel bar was yielded at the same location.

To validate the model proposed in this paper, a reinforced concrete beam tested by Esfahania et al. (2007) was modelled herein. The reinforced concrete beam named B5-16D-0L, 2000 mm long, was tested under four-point loads. The beam having a cross-sectional dimension of 150 mm width and 200 mm height was reinforced by two ribbed 16 mm dia. tensile steel bars and two ribbed 10 mm dia. compressive steel bars. The compressive strength of concrete at testing is 23·8 MPa. The yield strengths of 16 and 10 mm dia. bars are 406 and 365 MPa, respectively. The details of the tested reinforced concrete beam are given in Figure 13. Again in this case, the stirrups have not been included in the modelling in order to simplify the analysis. To examine the effects of different bond actions on the cracking pattern and crack opening of the beam, four cases were modelled under perfect bond, partial bond (ribbed steel bar and smooth steel bar) and zero bond conditions, respectively. Same as in the previous examples, the analyses were performed using full Newton–Raphon solution procedure under load control. A total load of P = 100 kN was applied with the load increment of 1 kN for each loading step. The analyses were completed with 91 load steps for perfect bond case; 86 load steps for the partial bond-ribbed bar case; 58 load steps for partial bond-smooth bar case; 10 load steps for zero bond case; and 80 load steps for the partial bond-ribbed bar case considering the compressive failure.

To investigate mesh sensitivity of the current model, the beam was modelled using different meshes – that is 200 concrete elements plus 40 steel bar elements and 400 concrete elements plus 80 steel bar elements. The comparison of predicted loads against mid-span deflection using different FE meshes is given in Figure 14. It can be seen that the results are almost identical to each other. Hence, the current model is not very mesh sensitive.

The comparisons of the load against mid-span deflection relations predicted by the current model for different conditions together with the tested results are shown in Figure 15. Apart from ‘prediction-compressive failure (partial bond-ribbed bar)’ case, all the predicted curves were modelled without considering the compressive failure of concrete. It was found that the predicted ultimate strength with consideration of compressive failure of concrete is around 6% lower than that without consideration of compressive failure, due to the fact that the extreme compressive elements of concrete were crushed after the tensile reinforcing steel bar yielded. The curves of the perfect bond case and partial bond-ribbed bar case are generally smooth before the yield of reinforcing bars. This is due to the fact that the reinforcing bars could pick up the released tensile stresses of the cracked concrete. In the case of zero bonded condition, the load against mid-span deflection relation shows significantly discontinuous characteristic after the concrete elements cracked. Generally, the predicted ultimate strengths of the beam with consideration of compressive failure and partial bond (ribbed steel bar) agree well with the tested results. As can be seen in Figure 15, if the concrete is assumed to be not cracking, the predicted result greatly overestimates the ultimate load and rigidity of the beam. It is obvious that the effects of bond characteristics on the structural behaviour of the beam are significant.

Other than the global response of the beam, it is also interesting to see how the bond–slip characteristics influence the cracking pattern and crack opening of the reinforced concrete beam. Figure 16 illustrates the crack pattern and crack opening of the predicted beam with different bond–slip characteristics. The stronger the applied bond strength between the steel bar and concrete, the wider the range of cracks distributed along the length of the beam. In contrast, with the decrease of bond strength the number of cracks within the beam becomes fewer and correspondingly the crack opening becomes bigger. For the zero bond case, there are only four parallel flexural cracks within the beam and only the two middle cracks showed remarkable openings. It is worth to note that the reinforcing steel bar is more possible to fracture locally when a stronger bond action is provided, and on the contrary, if little bond force is applied between the steel bar and concrete the stresses of the steel bar tends to be uniformly distributed along the length of the bar, so that the steel bar can bear a very large deformation before fracture. The comparison of the predicted opening history is demonstrated in Figure 17. The load against maximum crack opening curves of the perfect bond case and the partial bond-ribbed bar case are generally close to each other. As shown in Figure 18, the partial bond (ribbed bar) case generally shows flexural deformation along the longitudinal direction of the beam, but apparent localised fracture at mid-span can be observed in the zero bond case.

To validate the model in terms of modelling the propagation of shear cracks, another case presented herein is the shear panel tests conducted at the University of Toronto (Xie, 2009). Xie (2009) performed six shear panel tests under different combinations of axial stress and shear stress. The dimensions of the test panels were 890 mm × 890 mm × 70 mm. The panel specimens were reinforced by ten D8 cold-drawn deformed wires in the X-direction and by five D4 cold-drawn deformed wires in the Y-direction. The real-sectional areas of D8 wire and D4 wire are 51·6 and 24·2 mm², respectively, and their yield strengths are 641 and 581 MPa, respectively. The arrangements of reinforcing bars are shown in Figure 19(a). The panel specimens were tested under a combination of shear stress and different levels of axial stresses, where f_x/v varied from −2·8 to 3·0, in which f_x is the axial stress in the X-direction (Figure 19) and v is the shear stress, and a negative value means compressive axial stress and a positive value means tensile axial stress. Two specimens – that is the panel PL4 with f_x/v = −2·8 and the panel PL6 with f_x/v = 3·0 were modelled in this paper. The compressive strengths of concrete at testing are 43·1 MPa for PL4 and 43·5 MPa for PL6, respectively. The loads were applied through 20 shear keys along the perimeter of the panels, as shown in Figure 19(b), each of which was connected to two hydraulic jacks. Either tensile force or compressive force could be applied by the hydraulic jacks, so that different combinations of shear force and normal force were provided by varying the relative magnitude of the applied forces. In the modelling, two D8 bars at the same height were represented by a steel bar element at the reference plane. The model contained 256 quadrilateral concrete elements, 120 steel bar elements and 255 bond-link elements.

Figures 20 and 21 give comparisons of predicted cracking pattern with tested result for the panel specimens PL4 and PL6, respectively. It can be seen that the model can reasonably predict the formations and propagations of shear cracks within the panels. Figure 22 shows the comparison of the predicted shear stress against shear strain relations together with the tested results. It is evident that reasonable agreements were achieved between the tested results and predictions for such complex shear tests. It seems that the predicted results underestimate the shear stress while overestimate the shear strain.

In this paper, a robust FE procedure for modelling the localised fracture of reinforced concrete members has been developed. In this new model, the reinforced concrete member is modelled as an assembly of plain concrete, reinforcing steel bar and bond-link elements. The XFEM is incorporated into the plain concrete elements for modelling the formation and growth of individual cracks. The displacement jump over a crack is represented by multiplying the enhancement function by additional nodal displacement. The crack initiation and growth depend on the maximum principal tensile stress of concrete element. Different possibilities of crack propagation are considered in the model in order to simulate various crack patterns in a reinforced concrete member. A special integration scheme is applied to the enhanced elements so that the integration can be properly performed on both sides of the crack.

The numerical example and the validations show that the model can reasonably predict the localised flexural cracks and shear cracks within the reinforced concrete structures. The model is robust to identify individual crack initiation and propagation within the structural member. Hence, the current model provides an excellent numerical approach for assessing both structural stability (global behaviour) and integrity (localised fracture) of reinforced concrete members under extreme loading conditions. The predicted results show that the bond characteristic between concrete and reinforcing bars has a significant influence on the crack pattern of reinforced concrete members. The crack opening of the beam increases significantly with the decrease in bond strength.

The authors acknowledge the support of the Engineering and Physical Sciences Research Council of Great Britain under grant no. EP/I031553/1.

B_enr^a

enhanced strain–displacement transformation matrix

B_sta^u

regular strain–displacement transformation matrix

D

material constitutive matrix of plain concrete

f^int

element internal force vector

f_a^int

enhanced element internal force vector

f_u^int

regular element internal force vector

f_Γ^int

element internal force vector corresponding to traction

G_f

fracture energy of concrete

K_aa

enhanced element stiffness matrix

K_uu

regular element stiffness matrix

K_Γ

element stiffness matrix corresponding to traction

sign (x)

sign function

T_a

tangent stiffness of traction–separation relation

t_a

traction within the cracks

u_cont

vector of continuous displacement field

u_dis

vector of discontinuous displacement field

Ψ_i (x)

enhancement function