- Email: [email protected]
Finite element modelling of multiple cohesive discrete crack propagation in reinforced concrete beams
Engineering Fracture Mechanics 72 (2005) 2280–2297 www.elsevier.com/locate/engfracmech
Finite element modelling of multiple cohesive discrete crack propagation in reinforced concrete beams Z.J. Yang
a,*
, Jianfei Chen
b
a
b
Institute of Hydraulic Structures and Water Environment, College of Civil Engineering and Architecture, Zhejiang University, Hangzhou 310027, China Institute for Infrastructure and Environment, School of Engineering and Electronics, Edinburgh University, The King’s Buildings, Edinburgh EH9 3JN, United Kingdom Received 22 August 2003; received in revised form 5 January 2005; accepted 8 February 2005 Available online 25 April 2005
Abstract This paper presents a finite element (FE) model for fully automatic simulation of multiple discrete crack propagation in reinforced concrete (RC) beams. The discrete cracks are modelled based on the cohesive/fictitious crack concept using nonlinear interface elements with a bilinear tensile softening constitutive law. The model comprises an energybased crack propagation criterion, a simple remeshing procedure to accommodate crack propagations, two state variable mapping methods to transfer structural responses from one FE mesh to another, and a local arc-length algorithm to solve system equations characterised by material softening. The bond-slip behaviour between reinforcing bars and surrounding concrete is modelled by a tension-softening element. An example RC beam with well-documented test data is simulated. The model is found capable of automatically modelling multiple crack propagation. The predicted cracking process and distributed crack pattern are in close agreement with experimental observations. The load–deflection relations are accurately predicted up to a point when compressive cracking becomes dominant. The effects of bond-slip modelling and the efficiency and effectiveness of the numerical algorithms, together with the limitations of the current model, are also discussed. 2005 Published by Elsevier Ltd. Keywords: Finite element analysis; Local arc-length method; Cohesive crack model; Multiple crack propagation; Reinforced concrete beams
*
Corresponding author. Tel.: +86 571 87978025; fax: +86 571 87951356. E-mail addresses: [email protected] (Z.J. Yang), [email protected] (J. Chen).
0013-7944/$ - see front matter 2005 Published by Elsevier Ltd. doi:10.1016/j.engfracmech.2005.02.004
Z.J. Yang, J. Chen / Engineering Fracture Mechanics 72 (2005) 2280–2297
2281
Nomenclature A fc0 ft G Gf l Nd Nint w k b f fe K P r u D d $( )
crack surface area concrete cylindrical compressive strength concrete tensile strength total strain energy release rate (SERR) concrete fracture energy arc length desired optimum iteration number number of interface elements crack opening displacement (COD) loading factor tolerance of convergence cohesive force vector reference loading vector elastic stiffness matrix equivalent nodal force vector out-of-balance force vector displacement vector incremental change of a scalar or vector iterative change of a scalar or vector relative displacement vector (RDV)
1. Introduction Finite element (FE) modelling of RC structures such as beams, plates and shells presents significant challenges to the engineering community because of their complex nonlinear behaviour. The nonlinearities arise mainly from two major material factors: plasticity of reinforcements and compressive concrete, and cracking of concrete. Other nonlinearities are due to bond-slip interaction between reinforcements and concrete, aggregate interlock of cracked concrete, dowel action of reinforcements, concrete creep and shrinkage etc. One of the most challenging issues has been accurate prediction of concrete cracking. Although great effort has been made in developing numerical models in recent years, a recent review [1] indicates that only trivial structural cracking problems with simple geometries have so far been solved. It is now well known that in order to accurately model cracking behaviour of normal-sized structures such as RC beams, the fracture process zone (FPZ) must be properly modelled to consider the gradual energy dissipation during cracking [1]. This means models based on nonlinear fracture mechanics rather than linear elastic fracture mechanics should be used. Two types of crack models, i.e., smeared crack model and discrete crack model, have been most widely used for modelling FPZ. In modelling single tensile crack propagation in plain concrete beams, both discrete and smeared models are able to satisfactorily predict crack trajectories and load–displacement responses [2]. However, both models have only achieved limited success in modelling multiple distributed crack propagation in structures such as RC beams, with the smeared crack model being far more popular because of its computational convenience [1,3,4]. It may be noted that discrete crack models have distinct advantages over smeared models where information such as detailed cracking path, crack spacing and crack width is of interests.
2282
Z.J. Yang, J. Chen / Engineering Fracture Mechanics 72 (2005) 2280–2297
Due to the necessity of a remeshing procedure to accommodate crack propagation, the discrete crack model has been rarely used to model multiple concrete cracking. Its first application to multiple crack modelling seemed to be a study of the bond-slip behaviour of cylindrical ‘‘tension-pull’’ specimens by Ingraffea et al. [5]. Crack patterns of similar specimens were predicted by Yao and Murray [6]. Prasad and Krishnamoorthy [7] recently modelled a tension specimen with known cracking paths. Shi et al. [8] modelled a plain concrete beam using artificial notches to represent initial cracks. Feenstra et al. [9] carried out multiple discrete crack simulations of RC beams by inserting cohesive cracks at predetermined locations. The authorÕs previous study [10] presented an automatic modelling of multiple discrete crack propagation in RC beams strengthened with externally bonded composite beams based on linear elastic fracture mechanics. In general, most of these very limited studies of multiple discrete crack modelling were focused on bond-slip behaviour and based on predetermined crack paths or pre-specified initial crack locations. An automatic finite element modelling of multiple discrete crack propagation in RC structures exhibiting complex fracture behaviour is still not available. This study presents a nonlinear FE model for automatically modelling multiple discrete crack propagation in RC beams. It first describes the key aspects of the model and its computer implementation. A RC beam is then modelled. Numerical predictions of both cracking process and load–displacement relations are compared with experimental records. The effect of the bond-slip model, the efficiency and effectiveness of the numerical algorithms and the limitations of the current model are also discussed.
2. FE model for multiple cohesive discrete crack propagation The discrete crack FE model for RC beams consists of six key components: a cohesive crack model simulating the softening behaviour of tensile concrete, a bond-slip model considering interaction between steel reinforcement and concrete, an energy based crack propagation criterion, a remeshing procedure, meshmapping techniques to transfer structural responses from one FE mesh to another, and a numerical algorithm to solve the nonlinear equation systems. These key aspects are briefly described in the following sections. 2.1. Cohesive crack model (CCM) of concrete The CCM or fictitious crack model first proposed by Hillerborg et al. [11] can accurately model the energy dissipation process in quasi-brittle materials such as concrete. It assumes that a fictitious crack or a FPZ exists ahead of a real crack tip. The FPZ has the capability of transferring stresses through mechanisms such as aggregate interlock and material bonding until the crack opening displacement (COD) reaches a critical value. The CCM has become the basis of nonlinear discrete crack modelling and has been incorporated into some finite element codes in the form of 2-dimensional four-node, six-node and 3-dimensional eight-node interface elements to model mode-I and mixed-mode crack propagation (e.g., [2,5,12– 14]). The four-node interface element proposed by Gerstle and Xie [12] is adopted to represent the cohesive cracks in this study. Fig. 1 schematically shows the FPZ in concrete structures modelled with two interface elements. PeterssonÕs bi-linear curve [15] is used here to model the softening behaviour of the cohesive interface elements. Fig. 2 shows the bi-linear COD–traction curve with unloading path indicated. The initial stiffness before the concrete tensile strength is reached should be high enough to represent the uncracked material but not too high to cause numerical ill-conditioning. The model assumes an ‘‘irreversible unloading path’’, i.e., after reaching a value w*, for decreasing value of w an elastic unloading occurs with the corresponding scant stiffness. In the model, when COD is negative (compression) in the iterations on unloading, a high compressive stiffness will be assigned to the interface element to prevent the crack surfaces from kinking into each other.
Z.J. Yang, J. Chen / Engineering Fracture Mechanics 72 (2005) 2280–2297
2283
Aggregate Concrete bulk Real crack
Fracture process zone (FPZ)
Traction distribution along FPZ
ft
4 3
wc 2 1 Interface element 2
Fictitious crack tip node y x Interface element 1
Fig. 1. Modelling of FPZ with four-node interface elements (after Xie [13]).
σ σ0=ft
Loading Unloading
σ1=ft/3 w w0
w*
2 w1 = wc 9
wc
Fig. 2. COD–traction curve of nonlinear interface elements.
2.2. Bond-slip model The reinforcing bars are modelled using elastic–perfectly plastic truss elements. A simple yet rational bond-slip model for the bond behaviour between reinforcement bars and concrete proposed by Ingraffea et al. [5] is adopted. Assuming that the bond-slip behaviour between a reinforcing bar and its surrounding concrete is mainly governed by localised secondary cracking of concrete near a primary crack, the model lumps all the nonlinearities caused by bond-slip due to secondary cracking onto a special two-node ‘‘tension-softening’’ element. Fig. 3 illustrates a tension-softening element connected with two four-noded cohesive interface elements modelling the primary crack and two truss elements modelling a reinforcing bar. The tension-softening
2284
Z.J. Yang, J. Chen / Engineering Fracture Mechanics 72 (2005) 2280–2297
Secondary crack
concrete
2 Cohesive interface element Continuum element
Steel bar
∆
Primary crack
Truss element
tension-softening element Stress σ
Stress σy
σy ∆y
Crack opening ∆
Strain ε
Fig. 3. Bond-slip modelling with tension softening element (after Xie [13]).
element is characterised by a stress-opening relation, which is dependent upon concrete strength, the diameter and form of the reinforcing bar and loading conditions [13]. 2.3. Crack initiation and propagation criteria When the principal tensile stress at a node exceeds the tensile strength of the concrete, a crack perpendicular to the principal tensile stress is assumed to initiate. A proper crack propagation criterion is needed to determine when and in which direction a crack will propagate. The original CCM [11] assumes that the crack propagates when the maximum principal stress of the crack tip node reaches the concrete tensile strength. This stress-based assumption has been used in most previous studies (e.g., [16–18]), but it requires very fine crack-tip meshes to predict accurate nodal stresses. In order to reduce the required mesh density near the crack tip so that the remeshing procedure can be simplified, Xie [13] developed an energy-based cohesive crack propagation criterion G uT
of ¼0 oA
ð1Þ
in which u is the global displacement vector of the whole structure, f is the global cohesive force vector and A is the crack surface area and G is the total strain energy release rate (SERR) which can be calculated using [14]: 1 oK oP G ¼ uT u þ uT 2 oA oA
ð2Þ
where K is the total stiffness matrix of the elastic bulk and P the total equivalent nodal force vector due to external and body forces. The value of G in Eq. (2) can be calculated using a simple coplanar virtual crack extension (VCE) technique [13]. The second term in Eq. (1) can be explicitly derived when the four-node interface elements in Fig. 1 are used [13]. Because the convergence rate for energy entities is faster than that for stresses in FE analysis, crack propagation can be modelled more accurately even when coarse meshes are used [2,13,14]. As regards the crack propagation direction, Xie [13] found that using the direction of the maximum energy release rate as the crack propagation direction did not lead to reasonable results for nonlinear mixed-
Z.J. Yang, J. Chen / Engineering Fracture Mechanics 72 (2005) 2280–2297
2285
mode crack propagation problems. Instead, using the direction of maximum principal stress proved reasonable (e.g., [2,12–14,17]) and was adopted in this study. Therefore, the current crack propagation criterion is not a pure energy-based criterion. 2.4. Remeshing procedure An efficient remeshing procedure is paramount in discrete crack modelling to accommodate crack propagation. Available procedures can be generally classified into two categories. One may be termed ‘‘removerebuild’’ algorithm. In this algorithm, a new crack-tip node is determined by extending a specified crack growth increment in the calculated propagation direction. The original mesh within a certain range around the new crack-tip node is then completely removed. A complex procedure is followed to form the new crack and regenerate the mesh within this range where a regular rosette is added at the crack tip. Examples of this category include those developed by Wawrzynek and Ingraffea [19] and Bocca et al. [16,17]. The other category of algorithms may be termed ‘‘insert-separate’’ algorithms (e.g. [13,14]). In this procedure, a new edge from the old crack-tip node is first inserted into the local mesh in the propagation direction. The intersection point of this edge with the original mesh is the new crack-tip node. The new crack is then formed by separating those nodes along the line through the new and old crack-tip nodes. A rosette can be finally added to refine the tip-node mesh. Because this procedure neither completely removes nor rebuilds the new crack-tip mesh as in the ‘‘remove-rebuild’’ procedure, fewer elements are affected and the procedure is much simpler. The ‘‘insert-separate’’ [13] is used in this study. 2.5. Mesh mapping techniques After remeshing, the structural state variables from the old mesh need to be accurately mapped to the new mesh as their initial values to be used in the next loading step to ensure accuracy and numerical convergence. The most widely used mapping methods are inverse isoparametric mapping (e.g. [20]) and direct interpolation (e.g. [21]). Both methods are implemented in the model. 2.6. Local arc-length method The material softening represented by the nonlinear traction–COD relation (Fig. 2), together with continuous change of geometry due to discrete crack propagation, poses significant challenges in finding a robust and efficient solution procedure, especially when post-peak responses are desired. A recent comparative study [2] shows that a variety of numerical procedures have been used to solve nonlinear equation systems associated with material softening. Among them, arc-length controlled procedures are found to be more robust in tracing equilibrium paths characterised by snap-back which often occurs in mixedmode fracture modelling. It was also shown that the local arc-length algorithms are superior to global ones in terms of numerical robustness and efficiency [2]. The term ‘‘global’’ here means that the arc-length constraint equations include all the degrees of freedom whereas only selected degrees of freedom of dominant nodes are included in a ‘‘local’’ arc-length method. DuanÕs local updated normal plane constraint equation [22,23] is used in this study. Within one load increment in the local arc-length method, the constraint equations, the loading factor at the beginning of the load increment dk1 and the iterative loading factors dkk can be expressed as X e
rðDuT Þ rðDuk Þ ¼ l2 ;
ðk ¼ 2; 3; 4. . . . .Þ
ð3aÞ
2286
Z.J. Yang, J. Chen / Engineering Fracture Mechanics 72 (2005) 2280–2297
l dk1 ¼ rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi P rðd uT Þ rðd uÞ
ð3bÞ
e
P k
1
dk ¼ dk
e
rðd uT Þ ðrðDuk1 Þ þ rðdu ÞÞ P ; rðd uT Þ rðd uÞ
ðk ¼ 2; 3; 4. . . . .Þ
ð3cÞ
e
where l is the specified arc-length, u is the displacement vector and k is the loading factor; the symbols D and d represent incremental and iterative change, respectively; k is the iterative number; du and du* are iterative displacement vectors [24]. The summation in Eq. (3) is calculated in an element-by-element manner. Only the elements contributing to structural nonlinearity are included in the constraints and they are termed dominant elements. In the present study, these are the nonlinear interface elements. The symbol $ denotes the relative displacement vector (RDV) of dominant elements, rðaÞ ¼ ½ a1 an
a2 a1
a3 a2
an an1
T
ð4Þ
where vector a is any displacement vector (e.g. du) in Eq. (3) and n is the degrees of freedom of an interface element. E.g., for a four-node interface element, Eq. (4) results in a vector with 8 members. When the nonlinear interface elements are treated as the dominant elements in the present discrete crack modelling, the RDV may include only two crack opening displacements and two crack sliding displacements of the two pairs of nodes for each four-node interface element (Fig. 1), i.e. rðaÞ ¼ ½ a2 a8
a4 a6
a1 a7
T
a3 a5
ð5Þ
This local arc-length method represented by Eqs. (3)–(5) have three advantages [22,23]: (i) Eq. (3a) limits the iteration trajectory on a spherical surface, which always has an intersection with the equilibrium path; (ii) it does not need to choose a proper root as spherical/cylindrical constraints do since there is only one root for the loading factor; and (iii) the sign of the loading factor changes automatically if necessary within the iterations if the total secant stiffness is used. These advantages have also been demonstrated in [2]. At the beginning of each loading step, the arc-length l (Eq. (3)) must be determined. Bellini and Chulya [25] found that the definition of l had a direct effect on the performance of cylindrical/spherical arc-length algorithms applied to geometrically nonlinear problems. The arc-length of the ith loading increment li is often determined from [24] m Nd li ¼ li1 ; i ¼ 2; 3; 4 . . . . . . ð6Þ N i1 where Ni1 and li1 are the iteration number and arc-length in the (i 1)th loading step, respectively, Nd is a desired optimum iteration number which is problem-dependent and coefficient m ranges from 0.25–0.5 [25]. Eq. (6) tends to adjust l to keep the iteration number around Nd. The arc-length in the first loading increment l1 may be determined from Eq. (3b) by specifying an initial loading factor based on a proper reference loading condition, e.g., dk1 = 0.1. In the present model, the structure is linear elastic initially. The global arc-length method is used before any crack is initiated (thus there is no interface element yet) for the reason of compatibility with the subsequent analyses. The local arc-length method is used thereafter. Special considerations need to be given at this transition. Assume that interface elements are first added before the loading increment j, a relatively large l1 is calculated from the global arc-length method as all nodal displacements are used in Eq. (3b). If a comparable lj calculated from Eq. (6) is used in Eq. (3b), dk1 for the jth loading increment can be very large because only RDVs of the added interface elements are used in Eq. (3b) in the local arc-length
Z.J. Yang, J. Chen / Engineering Fracture Mechanics 72 (2005) 2280–2297
2287
method. This can result in a large number of iterations or even numerical divergence. Eq. (6) is modified in this research to ensure a smooth transition from global to local arc-length method 8 1 > i ¼ 1; . . . . . . ; j 1 < dki ¼ kini m ð7Þ Nd > li1 i ¼ j; j þ 1; . . . . . . : li ¼ N i1 where kini is the initial loading factor for the first j increments. The following default values are used: m = 0.5 and kini = 0.1. The value of Nd appears to have significant effect on numerical efficiency and stability in solving problems involving nonlinear interface elements. As cracks propagate, more interface elements are added and thus more RDVs are included in Eq. (3b). Using a constant Nd therefore tends to reduce the arc-length l calculated from Eq. (7) and leads to smaller incremental loading factors and more increments. Based on simulation experience, the following relationship is proposed N d ¼ N gint þ N d0
ð8aÞ
where Nint is the number of interface elements varying with loading level and Nd0 is the initial reference iteration number. Nd0 = 10 is used in this study. The coefficient g is dependent on Nint and the type of stiffness matrix used. For secant-stiffness based algorithms 1.5 if N int 6 10 g¼ ð8bÞ 1.2 if N int > 10 For tangential stiffness based algorithms 1.1 if N int 6 10 g¼ 1.0 if N int > 10
ð8cÞ
The higher g for secant stiffness based algorithms is used to compensate its slow convergence. The model adopts the modified Newton–Raphson iterative procedure with the convergence criterion based on the out-of-balance force factor, i.e., krðuÞk 6b k0 kf e k
ð9Þ
where r is the out-of-balance force vector, f e is the reference loading vector, b is the tolerance, k0 is the converged total loading factor at the last loading increment and the norms are Euclidean. The whole loading procedure can be automatically simulated using Eqs. (3)–(9) if the following four parameters are specified: the reference loading vector f e, the initial load factor kini, the initial reference iteration number Nd0 and the convergence tolerance b.
3. Computer implementation of the model The cohesive crack model in forms of nonlinear interface elements (Section 2.1), the energy based crack propagation criterion (Section 2.3) and the remeshing procedure (Section 2.4) proposed by Xie [13] have been implemented in his program AUTOFRAP. This program, designed originally for workstations, has been transplanted to PCs and integrated with the other above-mentioned aspects into a FE analysis program with pre-processing and post-processing functions. Fig. 4 shows a simplified flowchart of the program. During remeshing, the FE information changes constantly. This requires a robust method to execute a large number of manipulations on the mesh topology. Six arrays of type TYPE in Fortran 90, representing
2288
Z.J. Yang, J. Chen / Engineering Fracture Mechanics 72 (2005) 2280–2297
Pre-processing
Nonlinear Analysis Nonlinear Fracture Propagation Yes
end
Collapsed? No Remeshing Mapping of State Variables
Next loading step
Fig. 4. Key steps of discrete crack FEM.
the six basic finite element entities, i.e., elements, nodes, edges, cracks, boundaries and materials, have been designed to carry out this task. Due to the flexibility of dynamic memory allocation technique in FORTRAN90, the sizes of these six arrays can be decided by multiplying their sizes of the initial FE mesh by an estimated size-expanding ratio to accommodate remeshing. The operations on the six arrays include enquiry, modification, deletion and addition of their sizes and properties. The most complex yet important manipulation is probably the deletion and addition of elements, which leads to the change of all relative topology information. In this model, all the topology information and response information, together with their histories caused by continuous remeshing are stored in and dealt with by only manipulating the six arrays. This greatly facilitated the program development.
4. Numerical example 4.1. Example beam and material properties A series of tests on RC beams carried out Bresler and Scordelis [26] have been widely used for validating nonlinear finite element models to simulate the behaviour of RC structures. For example, Cedolin and Dei Poli [27] modelled two of the test beams OA-1 and A-1. Barzegar and Maddipudi [28] also modelled the beam OA-1. Both studies used the smeared crack models. Test beam OA-2 was modelled in this study because its crack pattern at failure is available for comparison [26]. The beam OA-2 was simply supported and tested under a concentrated load at the mid-span. It had a length of 5029.2 mm with a shear span of 4572.0 mm and a cross-section of 305.8 · 561.3 mm (12 · 22 ft) (Fig. 5). It was reinforced with five 9# high strength deformed steel bars with total nominal cross-sectional area of 3290 mm2 at an effective depth of 466.2 mm. The beam was not reinforced with stirrups. The cylindrical compressive strength of concrete was fc0 ¼ 23.7 MPa. Its YoungÕs modulus was evaluated pffiffiffiffi from Ec ¼ 4730 fc0 to be 24 GPa according to ACI 318-95 [29]. The tensile strength was evaluated from ft ¼ 0.1fc0 ¼ 2.37 MPa. For ordinary concrete, the fracture energy Gf has values between 65 and 200 N/m [30]. Gf = 100 N/m was used in this study. Because Gf equals the area under the PeterssonÕs bilinear
Z.J. Yang, J. Chen / Engineering Fracture Mechanics 72 (2005) 2280–2297 F
9# steel bars 466.2
2514.6 561.3
2514.6
2289
4572.0
305.8
Fig. 5. RC beam OA-2 tested by Bresler and Scordelis [26] (unit: mm).
softening law (Fig. 2), the following parameters were obtained from Gf. r0 = ft = 2.37 MPa, wc = 0.152 mm, r1 = 0.79 MPa and w1 = 0.034 mm. The value of w0 is assumed as 0.0001 mm with a high initial stiffness k0 = 2400 MPa/mm. The steel bars had a YoungÕs modulus of 200 GPa and an average yield strength 552 MPa. The PoissonÕs ratios for the concrete and the steel bars were assumed here to be 0.18 and 0.3, respectively. 4.2. Finite element modelling The reinforcing bars were located at different levels in the cross section (Fig. 5). For convenience and simplicity, they were modelled by a single two-node truss element at the centroid of the bars in a given section. The bulk of concrete was assumed to behave linear-elastically and modelled with four-node isoparametric and three-node constant strain elements. The cracks were modelled with four-node nonlinear interface elements. A plane stress condition was assumed for the analysis. Only half of the beam was modelled considering symmetry. Fig. 6 shows an initial FE mesh modelled for the example beam. The tension-softening elements described in Section 2.2 were used to model the bond-slip behaviour. The force-opening relation reported by Ingraffea et al. [5] was slightly modified to consider the different steel yield strength (Fig. 7). In order to investigate the effects of bond-slip modelling on the simulation, two other relations were also modelled. One assumed a strong bond condition with a high constant stiffness (5400 MPa/mm) whereas another assumed a weak bond with a much lower stiffness (1700 MPa/mm) after the same initial stiffness. The three bond-slip data are termed ‘‘Mild Bond’’, ‘‘Strong Bond’’ and ‘‘Weak Bond’’ hereafter for convenience of discussion. The analysis was carried out according to Fig. 4. Firstly the nonlinear equation system at the current configuration was solved using the local arc-length algorithm. The crack initiation criterion and energy crack propagation criterion were then checked based on the converged stress and displacement solutions. Nonlinear interface elements were then inserted at the nodes satisfying the initiation criterion and the remeshing procedure was triggered at the existing crack tips where the energy crack propagation criterion is matched. The mesh mapping techniques were applied to transfer the stress and displacement fields to the new mesh, which was followed by another nonlinear analysis on the newly formed configuration. The iteration continued until the structure collapsed or some FE elements were too ill-shaped to result in accurate predictions.
Fig. 6. Initial finite element mesh: 410 elements and 420 nodes.
2290
Z.J. Yang, J. Chen / Engineering Fracture Mechanics 72 (2005) 2280–2297 600
Stress (MPa)
500 400 300 200
Strong Bond Mild Bond (Ingraffea et al, 1984) Weak bond
100 0 0
0.05
0.1
0.15
0.2
0.25
0.3
0.35
Opening (mm)
Fig. 7. Bond-slip assumptions for tension-softening element.
5. Results and discussion 5.1. Cracking process Fig. 8a–e illustrate the cracking process of the beam using the Strong Bond data, which was fully automatically modelled using the presented discrete crack model. A scale of 50 has been applied to the deformation for clarity of presentation. The dotted lines are truss elements representing all the five steel bars. The FE simulation shows that a few flexural cracks first appear near the mid-span at a load about F = 60 kN (Fig. 8a). As the load increases, these initial flexural cracks begin to propagate upwards with increasing crack widths. More cracks away from the mid-span are also initiated. At a load of about 110 kN (Fig. 8b), the major flexural cracks have propagated into the upper half of the beam and the cracks gradually become curved towards the loading point. When the load is about 200 kN (Fig. 8c), the major cracks have extended to about one-third of the beam depth below the compression face. The widths of most cracks continue to increase steadily and more cracks initiate away from the mid-span. The width of major cracks increases rapidly thereafter (Fig. 8d), but the number of cracks has been stabilised by now. It is interesting to note that, as the load further increases (Fig. 8e), the crack width increases much faster in the middepth of the beam than at the steel reinforcement level because of the restraints of the reinforcing bars on the crack. The analysis was terminated because some FE elements became too ill-conditioned due to remeshing. Although this remains a challenge to avoid such ill-conditioned meshes at late stage of analysis, this may become less a problem if the plastic behaviour of concrete under compression is considered as discussed in the next section. It can be seen that the FE mesh has significantly changed in Fig. 8a–e from the initial mesh (Fig. 6). The numbers of elements, nodes and interface elements increase from 410, 420 and 0 to 926, 696 and 115, respectively. It is remarkable that the adopted remeshing procedure can adapt such dramatic changes. The predicted cracking process in general agrees well with experimental observations [2]. Fig. 9a shows the test crack pattern [26]. This may be compared with the predicted one as shown in Fig. 9b which is the same as Fig. 4d except that the un-deformed mesh is shown here. Within the shear span, the model predicted 17 cracks (a number of which may be still too small to be seen in real test) compared with 13 observed in test. It may be noted that the shear crack near the support is not modelled yet because in the current model only cracks on the tensile face of the beam are allowed to initiate and they are not allowed
Z.J. Yang, J. Chen / Engineering Fracture Mechanics 72 (2005) 2280–2297
Fig. 8. Cracking process.
2291
2292
Z.J. Yang, J. Chen / Engineering Fracture Mechanics 72 (2005) 2280–2297
Fig. 9. Crack pattern at failure.
to propagate downwards to avoid crack intersection and remeshing difficulties. Further research is being carried out to alleviate these restrictions. 5.2. Load versus deflection and loading capacity Fig. 10 shows the predicted load versus deflection relationship at the mid-span using the Strong Bond assumption. Test results from [26] are also shown for comparison. It can be seen that the prediction is in very close agreement with the test data for loading up to about 285 kN. Beyond this point the stiffness of the beam is slightly over-predicted. This is inevitable because the plasticity of concrete under compression and compression cracking have not been considered in the current model. Fig. 11a and b show contours of the principal compressive stress. At F = 200 kN (Fig. 11a), the compressive stress near the loading point has exceeded the concrete compressive strength (fc0 ¼ 23.7 MPa) already. Because the plastic zone is confined in a very small zone at this stage even if it is considered, neglecting the plasticity has insignificant effect on the deflection. As the load is increased to F = 285 kN, the area where the 400 350
Load F (kN)
300 250 200 150 100
Strong Bond assumption Test data (Bresler and Scordelis, 1963)
50 0 0
2
4
6
8
10
12
Deflection at midspan (mm)
Fig. 10. Load–deflection curve using the strong bond assumption.
14
Z.J. Yang, J. Chen / Engineering Fracture Mechanics 72 (2005) 2280–2297
2293
Fig. 11. Principal compressive stress.
principal compressive stresses exceeds fc0 has expanded considerably. The actual plastic zone would be slightly larger than this area because of stress redistribution due to plastic deformation. The effect of plastic deformation of concrete on the stiffness of the beam thus starts to be visible. According to [26], a longitudinal splitting (compression cracks) started to occur in the compression zone near the loading point at a load about 80% (285 kN) of the ultimate failure load (358 kN) of the beam. ‘‘The failures occurred as a result of longitudinal splitting in the compression zone near the load point, and also by horizontal splitting along the tensile reinforcement near the end of the beam’’ [26]. The ultimate failure of the beam was due to the formation of a critical diagonal tension crack, joining a compressive splitting crack with a splitting crack at the reinforcement bar level. The current discrete crack model cannot yet simulate this final failure process because neither the plasticity of the concrete under compression nor compression cracking is considered. Therefore, the predicted crack pattern (Fig. 9e) and load–deflection curve above F = 285 kN (Fig. 10) do not accurately reflect the actual behaviour of the beam. Further research is being undertaken to properly consider the plasticity of concrete under compression and concrete compression cracking, aiming at simulating the entire loading process of RC structures. 5.3. Effect of bond-slip modelling All the three bond-slip assumptions as shown in Fig. 7 were used to simulate the behaviour of the beam. The resulting cracking process and crack patterns are very similar in all three cases, but at a same loading level the crack widths predicted from the weak bond model are slightly larger than those from the other two bond-slip assumptions. Fig. 12 shows the predicted load–deflection curves using the three bond-slip assumptions. All three curves are virtually identical up to a load about 110 kN. This is understandable because at the initial
2294
Z.J. Yang, J. Chen / Engineering Fracture Mechanics 72 (2005) 2280–2297 400 350
Load F (kN)
300 250 200 150
Strong Bond Bond-slip data (Ingraffea et al, 1984) Weak bond Test data (Bresler and Scordelis, 1963)
100 50 0 0
2
4
8 10 6 Deflection at midspan (mm)
12
14
Fig. 12. Effect of bond-slip assumption.
cracking stage the cracks are still within the concrete cover or have just crossed the reinforcement bars. The effect of the bond-slip assumption does not take place yet if a crack is still within the concrete cover, or is still insignificant if the crack has just crossed the reinforcing bar. Discrepancies between the weak bond model and the other two become significant when the load is above 160 kN when most of the cracks have penetrated a significant part of the beam depth. Because the significant lower stiffness of the tension-softening elements in the weak bond model allows the cracks to open more freely than the stronger bond-slip models, it considerably under-predicts the stiffness of the beam. The Mild Bond assumption slightly under predicts the stiffness. From this analysis, it may be concluded that the bond-slip between reinforcements and concrete has a significant effect on the overall structural behaviour of the beam. It also demonstrates that the two-node tension-stiffening element is able to model the bond-slip behaviour in an efficient and simple way, but further research is necessary to determine the accurate bond-slip relationship for the model. 5.4. Effectiveness and efficiency of the local arc-length algorithm Both tangential and secant iterative stiffness matrices for interface elements modelling the FPZ were tested with the local arc-length method. The results showed that although the tangential stiffness matrix based method was far more efficient than that based on secant stiffness matrix, it was less reliable than the latter and often led to numerical instability and divergence, especially when a number of cracks initiate or propagate simultaneously. This is perhaps due to its incompatibility with the physical unloading path assumed in the model (Fig. 2), i.e., in one loading increment, the sign of tangential stiffness at a Gauss point (Fig. 1) may alter from negative on loading to positive on unloading or the vice versa. In contrast, the secant stiffness is always positive during loading-unloading iterations, leading to better numerical convergence and stability. Another reason perhaps lies in the complex of the multiple cracking problems such as the modelled RC beam. Some cracks in such a complex problem may experience opening-closure loops from one iteration to another. A previous study [2] showed that the tangential stiffness based arc-length algorithms succeeded in modelling a shear beam with a single crack which keeps opening during the whole loading process. Therefore, the use of secant iterative stiffness matrices is recommended. The secant iterative stiffness based local arc-length methods using RDVs in form of both Eqs. (4) and (5) have resulted in successful simulations. The predicted load–deflection curves and crack patterns are
Z.J. Yang, J. Chen / Engineering Fracture Mechanics 72 (2005) 2280–2297
2295
virtually identical. For the strong bond assumption, the analysis using Eq. (4) required a total number of load increments of 33 with a total of 507 iterations for a load up to F = 285 kN. The corresponding numbers for an analysis using Eq. (5) are 66 and 369, respectively. Although using Eq. (5) doubles the total load increment number, the total number of iterations is reduced by nearly 30%. Therefore Eq. (5) is more efficient than Eq. (4), which has also been demonstrated in [2]. 5.5. Limitations of the model and further research needs Although the remeshing procedure adopted in the model is simple and efficient, it may generate illshaped elements, especially in dealing with multiple closely distributed cracks. Fig. 13a–c show three FE meshes during the cracking process. It can be seen that the initially regular mesh gradually becomes irregular, especially in areas where two cracks are close. Considerable challenge remains to devise a robust algorithm which can deal efficiently with multiple cracking and crack intersection problems. Although mesh objectivity is not a problem in single discrete crack modelling as demonstrated by the previous study [2], it becomes a major concern in modelling problems with multiple distributed cracks. The initial element size may influence the prediction of the crack propagation process and the final crack pattern because cracks are only allowed to initiate at nodes and cracks cannot initiate at two nodes in one element (or otherwise mesh refining is needed before cracks initiate). More investigation is being undertaken to clarify this issue.
Fig. 13. Progressive FE meshes.
2296
Z.J. Yang, J. Chen / Engineering Fracture Mechanics 72 (2005) 2280–2297
6. Conclusion This paper has presented a cohesive discrete crack model based finite element model for fully automatic simulation of multi-crack propagation in reinforced concrete beams. It comprises an energy-based crack propagation criterion, a simple remeshing procedure, two mesh mapping techniques and a local arc-length numerical solution technique. Four-node interface elements are used to model the cohesive cracks with bilinear concrete tensile softening constitutive laws. The bond-slip behaviour between reinforcing bars and their surrounding concrete is modelled using a two-node tension-softening element. The finite element model has been applied to simulate an RC beam with well-documented test data. The model was found capable of automatically modelling multiple crack propagation with good robustness and efficiency. The predicted cracking process and crack pattern are in close agreement with experimental observations. The load–deflection relations are accurately predicted up to a point when compressive cracking becomes dominant. It is also shown that the bond-slip modelling has significant effect on the predicted behaviour of the beam. The limitations of the model are also briefly discussed.
Acknowledgement This work was partly supported by the National Science Foundation of China (No. 50279046). The authors would like to thank Dr. Ming Xie in Weidlinger Associate INC and Prof. Walter H. Gerstle in the Department of Civil Engineering of New Mexico University in the USA for their kind provision of program AUTOFRAP, a modified version of which has been incorporated into the program developed and used in this research.
References [1] ACI Report 446.3 R-97. Finite element analysis of fracture in concrete structures: state-of-the-art, reported by ACI Committee 446; 1997. [2] Yang ZJ, Proverbs D. A comparative study of numerical solutions to nonlinear discrete crack modelling of concrete beams involving sharp snap-back. Engng Fract Mech 2003;71(1):81–105. [3] Abdollahi A. Numerical strategies in the application of the FEM to RC structures—I. Comput Struct 1996;58(6):1173–82. [4] Abdollahi A. Investigation of objectivity in the application of the FEM to RC structures—II. Comput Struct 1996;58(6):1183–211. [5] Ingraffea AR, Gerstle WH, Gergely P, Saouma V. Fracture mechanics of bond in reinforced concrete. ASCE J Struct Engng 1984;110(4):871–90. [6] Yao BD, Murray DW. Prediction of distributed concrete cracking in RC analysis. ASCE J Struct Engng 1993;119(10):2813–34. [7] Prasad MVKV, Krishnamoorthy CS. Computational model for discrete crack growth in plain and reinforced concrete. Comput Methods Appl Mech Engng 2002;191(25–26):699–2725. [8] Shi ZH, Ohtsu M, Susuki M, Hibino Y. Numerical analysis of multiple cracks in concrete using the discrete approach. ASCE J Struct Engng 2001;127(9):1085–91. [9] Feenstra PH, De Borst R, Rots JG. Numerical study on crack dilatancy. II: Applications. ASCE J Engng Mech 1991;117(4):754–69. [10] Yang ZJ, Chen JF, Proverbs D. Finite element modelling of concrete cover separation failure in FRP plated RC beams. Construct Build Mater 2003;17(1):3–13. [11] Hillerborg A, Modeer M, Petersson PE. Analysis of crack formation and crack growth in concrete by means of fracture mechanics and finite elements. Cem Concr Res 1976;6:773–82. [12] Gerstle WH, Xie M. FEM modelling of fictitious crack propagation in concrete. ASCE J Engng Mech 1992;118(2):416–34. [13] Xie M. Finite element modelling of discrete crack propagation, PhD thesis, University of New Mexico, USA; 1995. [14] Xie M, Gerstle WH. Energy-based cohesive crack propagation modelling. J Engng Mech, ASCE 1995;121(12):1349–458. [15] Petersson PE. Crack growth and development of fracture zone in plain concrete and similar materials. Report TVBM-1006, Lund Institute of Technology, Lund, Sweden; 1981.
Z.J. Yang, J. Chen / Engineering Fracture Mechanics 72 (2005) 2280–2297
2297
[16] Bocca P, Carpinteri A, Valente S. Size effects in the mixed-mode crack propagation: softening and snap-back analysis. Engng Fract Mech 1990;35(9):159–70. [17] Bocca P, Carpinteri A, Valente S. Mixed-mode fracture of concrete. Int J Solid Struct 1991;27(9):1139–53. [18] Cendon DA, Galvez JC, Elices M, Planas J. Modelling the fracture of concrete under mixed loading. Int J Fract 2000;103:293–310. [19] Wawrzynek PA, Ingraffea AR. An interactive approach to local remeshing around a propagation crack. Finite Elem Anal Des 1989;5:87–96. [20] James MA. A plane stress finite element model for elastic–plastic mode I/II crack growth. PhD thesis, Kansas State University; 1998. [21] Harbaken AM, Cescotto S. An automatic remeshing technique for finite element simulation of forming process. Int J Numer Methods Engng 1990;30:1503–25. [22] Duan Y. Non-linear material finite element analysis for concrete and masonry structures. PhD thesis, University of Bradford, UK; 1994. [23] May IM, Duan Y. A local arc-length procedure for strain softening. Comput Struct 1997;64(1–4):297–303. [24] Crisfield MA. Non-linear finite element analysis of solids and structures. Advanced topics, vol. II. New York: John Wiley and Sons; 1997. [25] Bellini PX, Chulya A. An improved automatic incremental algorithm for the efficient solution of nonlinear finite element equations. Comput Struct 1987;26:99–110. [26] Bresler B, Scordelis AC. Shear strength of reinforced concrete beams. J Am Concr Inst 1963;60-4:51–72. [27] Cedolin L, Del Poli S. Finite element studies of shear-critical R/C beams. ASCE J Engng Mech Div 1977;103(EM3):395–410. [28] Barzegar F, Maddipudi S. Three-dimensional modelling of concrete structures. II: Reinforced concrete. ASCE J Struct Engng 1995;123(10):1347–56. [29] ACI 318-95. Building code requirements for structural concrete (318-95) and commentary (318R-95), American Concrete Institute (ACI); 1999. [30] Elfgren L, editor. Report of the technical committee 90-FMA fracture mechanics to concrete-applications, RILEM. London: Chapman and Hall; 1989.
Finite element modelling of multiple cohesive discrete crack propagation in reinforced concrete beams Z.J. Yang
a,*
, Jianfei Chen
b
a
b
Institute of Hydraulic Structures and Water Environment, College of Civil Engineering and Architecture, Zhejiang University, Hangzhou 310027, China Institute for Infrastructure and Environment, School of Engineering and Electronics, Edinburgh University, The King’s Buildings, Edinburgh EH9 3JN, United Kingdom Received 22 August 2003; received in revised form 5 January 2005; accepted 8 February 2005 Available online 25 April 2005
Abstract This paper presents a finite element (FE) model for fully automatic simulation of multiple discrete crack propagation in reinforced concrete (RC) beams. The discrete cracks are modelled based on the cohesive/fictitious crack concept using nonlinear interface elements with a bilinear tensile softening constitutive law. The model comprises an energybased crack propagation criterion, a simple remeshing procedure to accommodate crack propagations, two state variable mapping methods to transfer structural responses from one FE mesh to another, and a local arc-length algorithm to solve system equations characterised by material softening. The bond-slip behaviour between reinforcing bars and surrounding concrete is modelled by a tension-softening element. An example RC beam with well-documented test data is simulated. The model is found capable of automatically modelling multiple crack propagation. The predicted cracking process and distributed crack pattern are in close agreement with experimental observations. The load–deflection relations are accurately predicted up to a point when compressive cracking becomes dominant. The effects of bond-slip modelling and the efficiency and effectiveness of the numerical algorithms, together with the limitations of the current model, are also discussed. 2005 Published by Elsevier Ltd. Keywords: Finite element analysis; Local arc-length method; Cohesive crack model; Multiple crack propagation; Reinforced concrete beams
*
Corresponding author. Tel.: +86 571 87978025; fax: +86 571 87951356. E-mail addresses: [email protected] (Z.J. Yang), [email protected] (J. Chen).
0013-7944/$ - see front matter 2005 Published by Elsevier Ltd. doi:10.1016/j.engfracmech.2005.02.004
Z.J. Yang, J. Chen / Engineering Fracture Mechanics 72 (2005) 2280–2297
2281
Nomenclature A fc0 ft G Gf l Nd Nint w k b f fe K P r u D d $( )
crack surface area concrete cylindrical compressive strength concrete tensile strength total strain energy release rate (SERR) concrete fracture energy arc length desired optimum iteration number number of interface elements crack opening displacement (COD) loading factor tolerance of convergence cohesive force vector reference loading vector elastic stiffness matrix equivalent nodal force vector out-of-balance force vector displacement vector incremental change of a scalar or vector iterative change of a scalar or vector relative displacement vector (RDV)
1. Introduction Finite element (FE) modelling of RC structures such as beams, plates and shells presents significant challenges to the engineering community because of their complex nonlinear behaviour. The nonlinearities arise mainly from two major material factors: plasticity of reinforcements and compressive concrete, and cracking of concrete. Other nonlinearities are due to bond-slip interaction between reinforcements and concrete, aggregate interlock of cracked concrete, dowel action of reinforcements, concrete creep and shrinkage etc. One of the most challenging issues has been accurate prediction of concrete cracking. Although great effort has been made in developing numerical models in recent years, a recent review [1] indicates that only trivial structural cracking problems with simple geometries have so far been solved. It is now well known that in order to accurately model cracking behaviour of normal-sized structures such as RC beams, the fracture process zone (FPZ) must be properly modelled to consider the gradual energy dissipation during cracking [1]. This means models based on nonlinear fracture mechanics rather than linear elastic fracture mechanics should be used. Two types of crack models, i.e., smeared crack model and discrete crack model, have been most widely used for modelling FPZ. In modelling single tensile crack propagation in plain concrete beams, both discrete and smeared models are able to satisfactorily predict crack trajectories and load–displacement responses [2]. However, both models have only achieved limited success in modelling multiple distributed crack propagation in structures such as RC beams, with the smeared crack model being far more popular because of its computational convenience [1,3,4]. It may be noted that discrete crack models have distinct advantages over smeared models where information such as detailed cracking path, crack spacing and crack width is of interests.
2282
Z.J. Yang, J. Chen / Engineering Fracture Mechanics 72 (2005) 2280–2297
Due to the necessity of a remeshing procedure to accommodate crack propagation, the discrete crack model has been rarely used to model multiple concrete cracking. Its first application to multiple crack modelling seemed to be a study of the bond-slip behaviour of cylindrical ‘‘tension-pull’’ specimens by Ingraffea et al. [5]. Crack patterns of similar specimens were predicted by Yao and Murray [6]. Prasad and Krishnamoorthy [7] recently modelled a tension specimen with known cracking paths. Shi et al. [8] modelled a plain concrete beam using artificial notches to represent initial cracks. Feenstra et al. [9] carried out multiple discrete crack simulations of RC beams by inserting cohesive cracks at predetermined locations. The authorÕs previous study [10] presented an automatic modelling of multiple discrete crack propagation in RC beams strengthened with externally bonded composite beams based on linear elastic fracture mechanics. In general, most of these very limited studies of multiple discrete crack modelling were focused on bond-slip behaviour and based on predetermined crack paths or pre-specified initial crack locations. An automatic finite element modelling of multiple discrete crack propagation in RC structures exhibiting complex fracture behaviour is still not available. This study presents a nonlinear FE model for automatically modelling multiple discrete crack propagation in RC beams. It first describes the key aspects of the model and its computer implementation. A RC beam is then modelled. Numerical predictions of both cracking process and load–displacement relations are compared with experimental records. The effect of the bond-slip model, the efficiency and effectiveness of the numerical algorithms and the limitations of the current model are also discussed.
2. FE model for multiple cohesive discrete crack propagation The discrete crack FE model for RC beams consists of six key components: a cohesive crack model simulating the softening behaviour of tensile concrete, a bond-slip model considering interaction between steel reinforcement and concrete, an energy based crack propagation criterion, a remeshing procedure, meshmapping techniques to transfer structural responses from one FE mesh to another, and a numerical algorithm to solve the nonlinear equation systems. These key aspects are briefly described in the following sections. 2.1. Cohesive crack model (CCM) of concrete The CCM or fictitious crack model first proposed by Hillerborg et al. [11] can accurately model the energy dissipation process in quasi-brittle materials such as concrete. It assumes that a fictitious crack or a FPZ exists ahead of a real crack tip. The FPZ has the capability of transferring stresses through mechanisms such as aggregate interlock and material bonding until the crack opening displacement (COD) reaches a critical value. The CCM has become the basis of nonlinear discrete crack modelling and has been incorporated into some finite element codes in the form of 2-dimensional four-node, six-node and 3-dimensional eight-node interface elements to model mode-I and mixed-mode crack propagation (e.g., [2,5,12– 14]). The four-node interface element proposed by Gerstle and Xie [12] is adopted to represent the cohesive cracks in this study. Fig. 1 schematically shows the FPZ in concrete structures modelled with two interface elements. PeterssonÕs bi-linear curve [15] is used here to model the softening behaviour of the cohesive interface elements. Fig. 2 shows the bi-linear COD–traction curve with unloading path indicated. The initial stiffness before the concrete tensile strength is reached should be high enough to represent the uncracked material but not too high to cause numerical ill-conditioning. The model assumes an ‘‘irreversible unloading path’’, i.e., after reaching a value w*, for decreasing value of w an elastic unloading occurs with the corresponding scant stiffness. In the model, when COD is negative (compression) in the iterations on unloading, a high compressive stiffness will be assigned to the interface element to prevent the crack surfaces from kinking into each other.
Z.J. Yang, J. Chen / Engineering Fracture Mechanics 72 (2005) 2280–2297
2283
Aggregate Concrete bulk Real crack
Fracture process zone (FPZ)
Traction distribution along FPZ
ft
4 3
wc 2 1 Interface element 2
Fictitious crack tip node y x Interface element 1
Fig. 1. Modelling of FPZ with four-node interface elements (after Xie [13]).
σ σ0=ft
Loading Unloading
σ1=ft/3 w w0
w*
2 w1 = wc 9
wc
Fig. 2. COD–traction curve of nonlinear interface elements.
2.2. Bond-slip model The reinforcing bars are modelled using elastic–perfectly plastic truss elements. A simple yet rational bond-slip model for the bond behaviour between reinforcement bars and concrete proposed by Ingraffea et al. [5] is adopted. Assuming that the bond-slip behaviour between a reinforcing bar and its surrounding concrete is mainly governed by localised secondary cracking of concrete near a primary crack, the model lumps all the nonlinearities caused by bond-slip due to secondary cracking onto a special two-node ‘‘tension-softening’’ element. Fig. 3 illustrates a tension-softening element connected with two four-noded cohesive interface elements modelling the primary crack and two truss elements modelling a reinforcing bar. The tension-softening
2284
Z.J. Yang, J. Chen / Engineering Fracture Mechanics 72 (2005) 2280–2297
Secondary crack
concrete
2 Cohesive interface element Continuum element
Steel bar
∆
Primary crack
Truss element
tension-softening element Stress σ
Stress σy
σy ∆y
Crack opening ∆
Strain ε
Fig. 3. Bond-slip modelling with tension softening element (after Xie [13]).
element is characterised by a stress-opening relation, which is dependent upon concrete strength, the diameter and form of the reinforcing bar and loading conditions [13]. 2.3. Crack initiation and propagation criteria When the principal tensile stress at a node exceeds the tensile strength of the concrete, a crack perpendicular to the principal tensile stress is assumed to initiate. A proper crack propagation criterion is needed to determine when and in which direction a crack will propagate. The original CCM [11] assumes that the crack propagates when the maximum principal stress of the crack tip node reaches the concrete tensile strength. This stress-based assumption has been used in most previous studies (e.g., [16–18]), but it requires very fine crack-tip meshes to predict accurate nodal stresses. In order to reduce the required mesh density near the crack tip so that the remeshing procedure can be simplified, Xie [13] developed an energy-based cohesive crack propagation criterion G uT
of ¼0 oA
ð1Þ
in which u is the global displacement vector of the whole structure, f is the global cohesive force vector and A is the crack surface area and G is the total strain energy release rate (SERR) which can be calculated using [14]: 1 oK oP G ¼ uT u þ uT 2 oA oA
ð2Þ
where K is the total stiffness matrix of the elastic bulk and P the total equivalent nodal force vector due to external and body forces. The value of G in Eq. (2) can be calculated using a simple coplanar virtual crack extension (VCE) technique [13]. The second term in Eq. (1) can be explicitly derived when the four-node interface elements in Fig. 1 are used [13]. Because the convergence rate for energy entities is faster than that for stresses in FE analysis, crack propagation can be modelled more accurately even when coarse meshes are used [2,13,14]. As regards the crack propagation direction, Xie [13] found that using the direction of the maximum energy release rate as the crack propagation direction did not lead to reasonable results for nonlinear mixed-
Z.J. Yang, J. Chen / Engineering Fracture Mechanics 72 (2005) 2280–2297
2285
mode crack propagation problems. Instead, using the direction of maximum principal stress proved reasonable (e.g., [2,12–14,17]) and was adopted in this study. Therefore, the current crack propagation criterion is not a pure energy-based criterion. 2.4. Remeshing procedure An efficient remeshing procedure is paramount in discrete crack modelling to accommodate crack propagation. Available procedures can be generally classified into two categories. One may be termed ‘‘removerebuild’’ algorithm. In this algorithm, a new crack-tip node is determined by extending a specified crack growth increment in the calculated propagation direction. The original mesh within a certain range around the new crack-tip node is then completely removed. A complex procedure is followed to form the new crack and regenerate the mesh within this range where a regular rosette is added at the crack tip. Examples of this category include those developed by Wawrzynek and Ingraffea [19] and Bocca et al. [16,17]. The other category of algorithms may be termed ‘‘insert-separate’’ algorithms (e.g. [13,14]). In this procedure, a new edge from the old crack-tip node is first inserted into the local mesh in the propagation direction. The intersection point of this edge with the original mesh is the new crack-tip node. The new crack is then formed by separating those nodes along the line through the new and old crack-tip nodes. A rosette can be finally added to refine the tip-node mesh. Because this procedure neither completely removes nor rebuilds the new crack-tip mesh as in the ‘‘remove-rebuild’’ procedure, fewer elements are affected and the procedure is much simpler. The ‘‘insert-separate’’ [13] is used in this study. 2.5. Mesh mapping techniques After remeshing, the structural state variables from the old mesh need to be accurately mapped to the new mesh as their initial values to be used in the next loading step to ensure accuracy and numerical convergence. The most widely used mapping methods are inverse isoparametric mapping (e.g. [20]) and direct interpolation (e.g. [21]). Both methods are implemented in the model. 2.6. Local arc-length method The material softening represented by the nonlinear traction–COD relation (Fig. 2), together with continuous change of geometry due to discrete crack propagation, poses significant challenges in finding a robust and efficient solution procedure, especially when post-peak responses are desired. A recent comparative study [2] shows that a variety of numerical procedures have been used to solve nonlinear equation systems associated with material softening. Among them, arc-length controlled procedures are found to be more robust in tracing equilibrium paths characterised by snap-back which often occurs in mixedmode fracture modelling. It was also shown that the local arc-length algorithms are superior to global ones in terms of numerical robustness and efficiency [2]. The term ‘‘global’’ here means that the arc-length constraint equations include all the degrees of freedom whereas only selected degrees of freedom of dominant nodes are included in a ‘‘local’’ arc-length method. DuanÕs local updated normal plane constraint equation [22,23] is used in this study. Within one load increment in the local arc-length method, the constraint equations, the loading factor at the beginning of the load increment dk1 and the iterative loading factors dkk can be expressed as X e
rðDuT Þ rðDuk Þ ¼ l2 ;
ðk ¼ 2; 3; 4. . . . .Þ
ð3aÞ
2286
Z.J. Yang, J. Chen / Engineering Fracture Mechanics 72 (2005) 2280–2297
l dk1 ¼ rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi P rðd uT Þ rðd uÞ
ð3bÞ
e
P k
1
dk ¼ dk
e
rðd uT Þ ðrðDuk1 Þ þ rðdu ÞÞ P ; rðd uT Þ rðd uÞ
ðk ¼ 2; 3; 4. . . . .Þ
ð3cÞ
e
where l is the specified arc-length, u is the displacement vector and k is the loading factor; the symbols D and d represent incremental and iterative change, respectively; k is the iterative number; du and du* are iterative displacement vectors [24]. The summation in Eq. (3) is calculated in an element-by-element manner. Only the elements contributing to structural nonlinearity are included in the constraints and they are termed dominant elements. In the present study, these are the nonlinear interface elements. The symbol $ denotes the relative displacement vector (RDV) of dominant elements, rðaÞ ¼ ½ a1 an
a2 a1
a3 a2
an an1
T
ð4Þ
where vector a is any displacement vector (e.g. du) in Eq. (3) and n is the degrees of freedom of an interface element. E.g., for a four-node interface element, Eq. (4) results in a vector with 8 members. When the nonlinear interface elements are treated as the dominant elements in the present discrete crack modelling, the RDV may include only two crack opening displacements and two crack sliding displacements of the two pairs of nodes for each four-node interface element (Fig. 1), i.e. rðaÞ ¼ ½ a2 a8
a4 a6
a1 a7
T
a3 a5
ð5Þ
This local arc-length method represented by Eqs. (3)–(5) have three advantages [22,23]: (i) Eq. (3a) limits the iteration trajectory on a spherical surface, which always has an intersection with the equilibrium path; (ii) it does not need to choose a proper root as spherical/cylindrical constraints do since there is only one root for the loading factor; and (iii) the sign of the loading factor changes automatically if necessary within the iterations if the total secant stiffness is used. These advantages have also been demonstrated in [2]. At the beginning of each loading step, the arc-length l (Eq. (3)) must be determined. Bellini and Chulya [25] found that the definition of l had a direct effect on the performance of cylindrical/spherical arc-length algorithms applied to geometrically nonlinear problems. The arc-length of the ith loading increment li is often determined from [24] m Nd li ¼ li1 ; i ¼ 2; 3; 4 . . . . . . ð6Þ N i1 where Ni1 and li1 are the iteration number and arc-length in the (i 1)th loading step, respectively, Nd is a desired optimum iteration number which is problem-dependent and coefficient m ranges from 0.25–0.5 [25]. Eq. (6) tends to adjust l to keep the iteration number around Nd. The arc-length in the first loading increment l1 may be determined from Eq. (3b) by specifying an initial loading factor based on a proper reference loading condition, e.g., dk1 = 0.1. In the present model, the structure is linear elastic initially. The global arc-length method is used before any crack is initiated (thus there is no interface element yet) for the reason of compatibility with the subsequent analyses. The local arc-length method is used thereafter. Special considerations need to be given at this transition. Assume that interface elements are first added before the loading increment j, a relatively large l1 is calculated from the global arc-length method as all nodal displacements are used in Eq. (3b). If a comparable lj calculated from Eq. (6) is used in Eq. (3b), dk1 for the jth loading increment can be very large because only RDVs of the added interface elements are used in Eq. (3b) in the local arc-length
Z.J. Yang, J. Chen / Engineering Fracture Mechanics 72 (2005) 2280–2297
2287
method. This can result in a large number of iterations or even numerical divergence. Eq. (6) is modified in this research to ensure a smooth transition from global to local arc-length method 8 1 > i ¼ 1; . . . . . . ; j 1 < dki ¼ kini m ð7Þ Nd > li1 i ¼ j; j þ 1; . . . . . . : li ¼ N i1 where kini is the initial loading factor for the first j increments. The following default values are used: m = 0.5 and kini = 0.1. The value of Nd appears to have significant effect on numerical efficiency and stability in solving problems involving nonlinear interface elements. As cracks propagate, more interface elements are added and thus more RDVs are included in Eq. (3b). Using a constant Nd therefore tends to reduce the arc-length l calculated from Eq. (7) and leads to smaller incremental loading factors and more increments. Based on simulation experience, the following relationship is proposed N d ¼ N gint þ N d0
ð8aÞ
where Nint is the number of interface elements varying with loading level and Nd0 is the initial reference iteration number. Nd0 = 10 is used in this study. The coefficient g is dependent on Nint and the type of stiffness matrix used. For secant-stiffness based algorithms 1.5 if N int 6 10 g¼ ð8bÞ 1.2 if N int > 10 For tangential stiffness based algorithms 1.1 if N int 6 10 g¼ 1.0 if N int > 10
ð8cÞ
The higher g for secant stiffness based algorithms is used to compensate its slow convergence. The model adopts the modified Newton–Raphson iterative procedure with the convergence criterion based on the out-of-balance force factor, i.e., krðuÞk 6b k0 kf e k
ð9Þ
where r is the out-of-balance force vector, f e is the reference loading vector, b is the tolerance, k0 is the converged total loading factor at the last loading increment and the norms are Euclidean. The whole loading procedure can be automatically simulated using Eqs. (3)–(9) if the following four parameters are specified: the reference loading vector f e, the initial load factor kini, the initial reference iteration number Nd0 and the convergence tolerance b.
3. Computer implementation of the model The cohesive crack model in forms of nonlinear interface elements (Section 2.1), the energy based crack propagation criterion (Section 2.3) and the remeshing procedure (Section 2.4) proposed by Xie [13] have been implemented in his program AUTOFRAP. This program, designed originally for workstations, has been transplanted to PCs and integrated with the other above-mentioned aspects into a FE analysis program with pre-processing and post-processing functions. Fig. 4 shows a simplified flowchart of the program. During remeshing, the FE information changes constantly. This requires a robust method to execute a large number of manipulations on the mesh topology. Six arrays of type TYPE in Fortran 90, representing
2288
Z.J. Yang, J. Chen / Engineering Fracture Mechanics 72 (2005) 2280–2297
Pre-processing
Nonlinear Analysis Nonlinear Fracture Propagation Yes
end
Collapsed? No Remeshing Mapping of State Variables
Next loading step
Fig. 4. Key steps of discrete crack FEM.
the six basic finite element entities, i.e., elements, nodes, edges, cracks, boundaries and materials, have been designed to carry out this task. Due to the flexibility of dynamic memory allocation technique in FORTRAN90, the sizes of these six arrays can be decided by multiplying their sizes of the initial FE mesh by an estimated size-expanding ratio to accommodate remeshing. The operations on the six arrays include enquiry, modification, deletion and addition of their sizes and properties. The most complex yet important manipulation is probably the deletion and addition of elements, which leads to the change of all relative topology information. In this model, all the topology information and response information, together with their histories caused by continuous remeshing are stored in and dealt with by only manipulating the six arrays. This greatly facilitated the program development.
4. Numerical example 4.1. Example beam and material properties A series of tests on RC beams carried out Bresler and Scordelis [26] have been widely used for validating nonlinear finite element models to simulate the behaviour of RC structures. For example, Cedolin and Dei Poli [27] modelled two of the test beams OA-1 and A-1. Barzegar and Maddipudi [28] also modelled the beam OA-1. Both studies used the smeared crack models. Test beam OA-2 was modelled in this study because its crack pattern at failure is available for comparison [26]. The beam OA-2 was simply supported and tested under a concentrated load at the mid-span. It had a length of 5029.2 mm with a shear span of 4572.0 mm and a cross-section of 305.8 · 561.3 mm (12 · 22 ft) (Fig. 5). It was reinforced with five 9# high strength deformed steel bars with total nominal cross-sectional area of 3290 mm2 at an effective depth of 466.2 mm. The beam was not reinforced with stirrups. The cylindrical compressive strength of concrete was fc0 ¼ 23.7 MPa. Its YoungÕs modulus was evaluated pffiffiffiffi from Ec ¼ 4730 fc0 to be 24 GPa according to ACI 318-95 [29]. The tensile strength was evaluated from ft ¼ 0.1fc0 ¼ 2.37 MPa. For ordinary concrete, the fracture energy Gf has values between 65 and 200 N/m [30]. Gf = 100 N/m was used in this study. Because Gf equals the area under the PeterssonÕs bilinear
Z.J. Yang, J. Chen / Engineering Fracture Mechanics 72 (2005) 2280–2297 F
9# steel bars 466.2
2514.6 561.3
2514.6
2289
4572.0
305.8
Fig. 5. RC beam OA-2 tested by Bresler and Scordelis [26] (unit: mm).
softening law (Fig. 2), the following parameters were obtained from Gf. r0 = ft = 2.37 MPa, wc = 0.152 mm, r1 = 0.79 MPa and w1 = 0.034 mm. The value of w0 is assumed as 0.0001 mm with a high initial stiffness k0 = 2400 MPa/mm. The steel bars had a YoungÕs modulus of 200 GPa and an average yield strength 552 MPa. The PoissonÕs ratios for the concrete and the steel bars were assumed here to be 0.18 and 0.3, respectively. 4.2. Finite element modelling The reinforcing bars were located at different levels in the cross section (Fig. 5). For convenience and simplicity, they were modelled by a single two-node truss element at the centroid of the bars in a given section. The bulk of concrete was assumed to behave linear-elastically and modelled with four-node isoparametric and three-node constant strain elements. The cracks were modelled with four-node nonlinear interface elements. A plane stress condition was assumed for the analysis. Only half of the beam was modelled considering symmetry. Fig. 6 shows an initial FE mesh modelled for the example beam. The tension-softening elements described in Section 2.2 were used to model the bond-slip behaviour. The force-opening relation reported by Ingraffea et al. [5] was slightly modified to consider the different steel yield strength (Fig. 7). In order to investigate the effects of bond-slip modelling on the simulation, two other relations were also modelled. One assumed a strong bond condition with a high constant stiffness (5400 MPa/mm) whereas another assumed a weak bond with a much lower stiffness (1700 MPa/mm) after the same initial stiffness. The three bond-slip data are termed ‘‘Mild Bond’’, ‘‘Strong Bond’’ and ‘‘Weak Bond’’ hereafter for convenience of discussion. The analysis was carried out according to Fig. 4. Firstly the nonlinear equation system at the current configuration was solved using the local arc-length algorithm. The crack initiation criterion and energy crack propagation criterion were then checked based on the converged stress and displacement solutions. Nonlinear interface elements were then inserted at the nodes satisfying the initiation criterion and the remeshing procedure was triggered at the existing crack tips where the energy crack propagation criterion is matched. The mesh mapping techniques were applied to transfer the stress and displacement fields to the new mesh, which was followed by another nonlinear analysis on the newly formed configuration. The iteration continued until the structure collapsed or some FE elements were too ill-shaped to result in accurate predictions.
Fig. 6. Initial finite element mesh: 410 elements and 420 nodes.
2290
Z.J. Yang, J. Chen / Engineering Fracture Mechanics 72 (2005) 2280–2297 600
Stress (MPa)
500 400 300 200
Strong Bond Mild Bond (Ingraffea et al, 1984) Weak bond
100 0 0
0.05
0.1
0.15
0.2
0.25
0.3
0.35
Opening (mm)
Fig. 7. Bond-slip assumptions for tension-softening element.
5. Results and discussion 5.1. Cracking process Fig. 8a–e illustrate the cracking process of the beam using the Strong Bond data, which was fully automatically modelled using the presented discrete crack model. A scale of 50 has been applied to the deformation for clarity of presentation. The dotted lines are truss elements representing all the five steel bars. The FE simulation shows that a few flexural cracks first appear near the mid-span at a load about F = 60 kN (Fig. 8a). As the load increases, these initial flexural cracks begin to propagate upwards with increasing crack widths. More cracks away from the mid-span are also initiated. At a load of about 110 kN (Fig. 8b), the major flexural cracks have propagated into the upper half of the beam and the cracks gradually become curved towards the loading point. When the load is about 200 kN (Fig. 8c), the major cracks have extended to about one-third of the beam depth below the compression face. The widths of most cracks continue to increase steadily and more cracks initiate away from the mid-span. The width of major cracks increases rapidly thereafter (Fig. 8d), but the number of cracks has been stabilised by now. It is interesting to note that, as the load further increases (Fig. 8e), the crack width increases much faster in the middepth of the beam than at the steel reinforcement level because of the restraints of the reinforcing bars on the crack. The analysis was terminated because some FE elements became too ill-conditioned due to remeshing. Although this remains a challenge to avoid such ill-conditioned meshes at late stage of analysis, this may become less a problem if the plastic behaviour of concrete under compression is considered as discussed in the next section. It can be seen that the FE mesh has significantly changed in Fig. 8a–e from the initial mesh (Fig. 6). The numbers of elements, nodes and interface elements increase from 410, 420 and 0 to 926, 696 and 115, respectively. It is remarkable that the adopted remeshing procedure can adapt such dramatic changes. The predicted cracking process in general agrees well with experimental observations [2]. Fig. 9a shows the test crack pattern [26]. This may be compared with the predicted one as shown in Fig. 9b which is the same as Fig. 4d except that the un-deformed mesh is shown here. Within the shear span, the model predicted 17 cracks (a number of which may be still too small to be seen in real test) compared with 13 observed in test. It may be noted that the shear crack near the support is not modelled yet because in the current model only cracks on the tensile face of the beam are allowed to initiate and they are not allowed
Z.J. Yang, J. Chen / Engineering Fracture Mechanics 72 (2005) 2280–2297
Fig. 8. Cracking process.
2291
2292
Z.J. Yang, J. Chen / Engineering Fracture Mechanics 72 (2005) 2280–2297
Fig. 9. Crack pattern at failure.
to propagate downwards to avoid crack intersection and remeshing difficulties. Further research is being carried out to alleviate these restrictions. 5.2. Load versus deflection and loading capacity Fig. 10 shows the predicted load versus deflection relationship at the mid-span using the Strong Bond assumption. Test results from [26] are also shown for comparison. It can be seen that the prediction is in very close agreement with the test data for loading up to about 285 kN. Beyond this point the stiffness of the beam is slightly over-predicted. This is inevitable because the plasticity of concrete under compression and compression cracking have not been considered in the current model. Fig. 11a and b show contours of the principal compressive stress. At F = 200 kN (Fig. 11a), the compressive stress near the loading point has exceeded the concrete compressive strength (fc0 ¼ 23.7 MPa) already. Because the plastic zone is confined in a very small zone at this stage even if it is considered, neglecting the plasticity has insignificant effect on the deflection. As the load is increased to F = 285 kN, the area where the 400 350
Load F (kN)
300 250 200 150 100
Strong Bond assumption Test data (Bresler and Scordelis, 1963)
50 0 0
2
4
6
8
10
12
Deflection at midspan (mm)
Fig. 10. Load–deflection curve using the strong bond assumption.
14
Z.J. Yang, J. Chen / Engineering Fracture Mechanics 72 (2005) 2280–2297
2293
Fig. 11. Principal compressive stress.
principal compressive stresses exceeds fc0 has expanded considerably. The actual plastic zone would be slightly larger than this area because of stress redistribution due to plastic deformation. The effect of plastic deformation of concrete on the stiffness of the beam thus starts to be visible. According to [26], a longitudinal splitting (compression cracks) started to occur in the compression zone near the loading point at a load about 80% (285 kN) of the ultimate failure load (358 kN) of the beam. ‘‘The failures occurred as a result of longitudinal splitting in the compression zone near the load point, and also by horizontal splitting along the tensile reinforcement near the end of the beam’’ [26]. The ultimate failure of the beam was due to the formation of a critical diagonal tension crack, joining a compressive splitting crack with a splitting crack at the reinforcement bar level. The current discrete crack model cannot yet simulate this final failure process because neither the plasticity of the concrete under compression nor compression cracking is considered. Therefore, the predicted crack pattern (Fig. 9e) and load–deflection curve above F = 285 kN (Fig. 10) do not accurately reflect the actual behaviour of the beam. Further research is being undertaken to properly consider the plasticity of concrete under compression and concrete compression cracking, aiming at simulating the entire loading process of RC structures. 5.3. Effect of bond-slip modelling All the three bond-slip assumptions as shown in Fig. 7 were used to simulate the behaviour of the beam. The resulting cracking process and crack patterns are very similar in all three cases, but at a same loading level the crack widths predicted from the weak bond model are slightly larger than those from the other two bond-slip assumptions. Fig. 12 shows the predicted load–deflection curves using the three bond-slip assumptions. All three curves are virtually identical up to a load about 110 kN. This is understandable because at the initial
2294
Z.J. Yang, J. Chen / Engineering Fracture Mechanics 72 (2005) 2280–2297 400 350
Load F (kN)
300 250 200 150
Strong Bond Bond-slip data (Ingraffea et al, 1984) Weak bond Test data (Bresler and Scordelis, 1963)
100 50 0 0
2
4
8 10 6 Deflection at midspan (mm)
12
14
Fig. 12. Effect of bond-slip assumption.
cracking stage the cracks are still within the concrete cover or have just crossed the reinforcement bars. The effect of the bond-slip assumption does not take place yet if a crack is still within the concrete cover, or is still insignificant if the crack has just crossed the reinforcing bar. Discrepancies between the weak bond model and the other two become significant when the load is above 160 kN when most of the cracks have penetrated a significant part of the beam depth. Because the significant lower stiffness of the tension-softening elements in the weak bond model allows the cracks to open more freely than the stronger bond-slip models, it considerably under-predicts the stiffness of the beam. The Mild Bond assumption slightly under predicts the stiffness. From this analysis, it may be concluded that the bond-slip between reinforcements and concrete has a significant effect on the overall structural behaviour of the beam. It also demonstrates that the two-node tension-stiffening element is able to model the bond-slip behaviour in an efficient and simple way, but further research is necessary to determine the accurate bond-slip relationship for the model. 5.4. Effectiveness and efficiency of the local arc-length algorithm Both tangential and secant iterative stiffness matrices for interface elements modelling the FPZ were tested with the local arc-length method. The results showed that although the tangential stiffness matrix based method was far more efficient than that based on secant stiffness matrix, it was less reliable than the latter and often led to numerical instability and divergence, especially when a number of cracks initiate or propagate simultaneously. This is perhaps due to its incompatibility with the physical unloading path assumed in the model (Fig. 2), i.e., in one loading increment, the sign of tangential stiffness at a Gauss point (Fig. 1) may alter from negative on loading to positive on unloading or the vice versa. In contrast, the secant stiffness is always positive during loading-unloading iterations, leading to better numerical convergence and stability. Another reason perhaps lies in the complex of the multiple cracking problems such as the modelled RC beam. Some cracks in such a complex problem may experience opening-closure loops from one iteration to another. A previous study [2] showed that the tangential stiffness based arc-length algorithms succeeded in modelling a shear beam with a single crack which keeps opening during the whole loading process. Therefore, the use of secant iterative stiffness matrices is recommended. The secant iterative stiffness based local arc-length methods using RDVs in form of both Eqs. (4) and (5) have resulted in successful simulations. The predicted load–deflection curves and crack patterns are
Z.J. Yang, J. Chen / Engineering Fracture Mechanics 72 (2005) 2280–2297
2295
virtually identical. For the strong bond assumption, the analysis using Eq. (4) required a total number of load increments of 33 with a total of 507 iterations for a load up to F = 285 kN. The corresponding numbers for an analysis using Eq. (5) are 66 and 369, respectively. Although using Eq. (5) doubles the total load increment number, the total number of iterations is reduced by nearly 30%. Therefore Eq. (5) is more efficient than Eq. (4), which has also been demonstrated in [2]. 5.5. Limitations of the model and further research needs Although the remeshing procedure adopted in the model is simple and efficient, it may generate illshaped elements, especially in dealing with multiple closely distributed cracks. Fig. 13a–c show three FE meshes during the cracking process. It can be seen that the initially regular mesh gradually becomes irregular, especially in areas where two cracks are close. Considerable challenge remains to devise a robust algorithm which can deal efficiently with multiple cracking and crack intersection problems. Although mesh objectivity is not a problem in single discrete crack modelling as demonstrated by the previous study [2], it becomes a major concern in modelling problems with multiple distributed cracks. The initial element size may influence the prediction of the crack propagation process and the final crack pattern because cracks are only allowed to initiate at nodes and cracks cannot initiate at two nodes in one element (or otherwise mesh refining is needed before cracks initiate). More investigation is being undertaken to clarify this issue.
Fig. 13. Progressive FE meshes.
2296
Z.J. Yang, J. Chen / Engineering Fracture Mechanics 72 (2005) 2280–2297
6. Conclusion This paper has presented a cohesive discrete crack model based finite element model for fully automatic simulation of multi-crack propagation in reinforced concrete beams. It comprises an energy-based crack propagation criterion, a simple remeshing procedure, two mesh mapping techniques and a local arc-length numerical solution technique. Four-node interface elements are used to model the cohesive cracks with bilinear concrete tensile softening constitutive laws. The bond-slip behaviour between reinforcing bars and their surrounding concrete is modelled using a two-node tension-softening element. The finite element model has been applied to simulate an RC beam with well-documented test data. The model was found capable of automatically modelling multiple crack propagation with good robustness and efficiency. The predicted cracking process and crack pattern are in close agreement with experimental observations. The load–deflection relations are accurately predicted up to a point when compressive cracking becomes dominant. It is also shown that the bond-slip modelling has significant effect on the predicted behaviour of the beam. The limitations of the model are also briefly discussed.
Acknowledgement This work was partly supported by the National Science Foundation of China (No. 50279046). The authors would like to thank Dr. Ming Xie in Weidlinger Associate INC and Prof. Walter H. Gerstle in the Department of Civil Engineering of New Mexico University in the USA for their kind provision of program AUTOFRAP, a modified version of which has been incorporated into the program developed and used in this research.
References [1] ACI Report 446.3 R-97. Finite element analysis of fracture in concrete structures: state-of-the-art, reported by ACI Committee 446; 1997. [2] Yang ZJ, Proverbs D. A comparative study of numerical solutions to nonlinear discrete crack modelling of concrete beams involving sharp snap-back. Engng Fract Mech 2003;71(1):81–105. [3] Abdollahi A. Numerical strategies in the application of the FEM to RC structures—I. Comput Struct 1996;58(6):1173–82. [4] Abdollahi A. Investigation of objectivity in the application of the FEM to RC structures—II. Comput Struct 1996;58(6):1183–211. [5] Ingraffea AR, Gerstle WH, Gergely P, Saouma V. Fracture mechanics of bond in reinforced concrete. ASCE J Struct Engng 1984;110(4):871–90. [6] Yao BD, Murray DW. Prediction of distributed concrete cracking in RC analysis. ASCE J Struct Engng 1993;119(10):2813–34. [7] Prasad MVKV, Krishnamoorthy CS. Computational model for discrete crack growth in plain and reinforced concrete. Comput Methods Appl Mech Engng 2002;191(25–26):699–2725. [8] Shi ZH, Ohtsu M, Susuki M, Hibino Y. Numerical analysis of multiple cracks in concrete using the discrete approach. ASCE J Struct Engng 2001;127(9):1085–91. [9] Feenstra PH, De Borst R, Rots JG. Numerical study on crack dilatancy. II: Applications. ASCE J Engng Mech 1991;117(4):754–69. [10] Yang ZJ, Chen JF, Proverbs D. Finite element modelling of concrete cover separation failure in FRP plated RC beams. Construct Build Mater 2003;17(1):3–13. [11] Hillerborg A, Modeer M, Petersson PE. Analysis of crack formation and crack growth in concrete by means of fracture mechanics and finite elements. Cem Concr Res 1976;6:773–82. [12] Gerstle WH, Xie M. FEM modelling of fictitious crack propagation in concrete. ASCE J Engng Mech 1992;118(2):416–34. [13] Xie M. Finite element modelling of discrete crack propagation, PhD thesis, University of New Mexico, USA; 1995. [14] Xie M, Gerstle WH. Energy-based cohesive crack propagation modelling. J Engng Mech, ASCE 1995;121(12):1349–458. [15] Petersson PE. Crack growth and development of fracture zone in plain concrete and similar materials. Report TVBM-1006, Lund Institute of Technology, Lund, Sweden; 1981.
Z.J. Yang, J. Chen / Engineering Fracture Mechanics 72 (2005) 2280–2297
2297
[16] Bocca P, Carpinteri A, Valente S. Size effects in the mixed-mode crack propagation: softening and snap-back analysis. Engng Fract Mech 1990;35(9):159–70. [17] Bocca P, Carpinteri A, Valente S. Mixed-mode fracture of concrete. Int J Solid Struct 1991;27(9):1139–53. [18] Cendon DA, Galvez JC, Elices M, Planas J. Modelling the fracture of concrete under mixed loading. Int J Fract 2000;103:293–310. [19] Wawrzynek PA, Ingraffea AR. An interactive approach to local remeshing around a propagation crack. Finite Elem Anal Des 1989;5:87–96. [20] James MA. A plane stress finite element model for elastic–plastic mode I/II crack growth. PhD thesis, Kansas State University; 1998. [21] Harbaken AM, Cescotto S. An automatic remeshing technique for finite element simulation of forming process. Int J Numer Methods Engng 1990;30:1503–25. [22] Duan Y. Non-linear material finite element analysis for concrete and masonry structures. PhD thesis, University of Bradford, UK; 1994. [23] May IM, Duan Y. A local arc-length procedure for strain softening. Comput Struct 1997;64(1–4):297–303. [24] Crisfield MA. Non-linear finite element analysis of solids and structures. Advanced topics, vol. II. New York: John Wiley and Sons; 1997. [25] Bellini PX, Chulya A. An improved automatic incremental algorithm for the efficient solution of nonlinear finite element equations. Comput Struct 1987;26:99–110. [26] Bresler B, Scordelis AC. Shear strength of reinforced concrete beams. J Am Concr Inst 1963;60-4:51–72. [27] Cedolin L, Del Poli S. Finite element studies of shear-critical R/C beams. ASCE J Engng Mech Div 1977;103(EM3):395–410. [28] Barzegar F, Maddipudi S. Three-dimensional modelling of concrete structures. II: Reinforced concrete. ASCE J Struct Engng 1995;123(10):1347–56. [29] ACI 318-95. Building code requirements for structural concrete (318-95) and commentary (318R-95), American Concrete Institute (ACI); 1999. [30] Elfgren L, editor. Report of the technical committee 90-FMA fracture mechanics to concrete-applications, RILEM. London: Chapman and Hall; 1989.