Computational aspects of the discontinuous deformation analysis framework for modelling concrete fracture

Engineering Fracture Mechanics 65 (2000) 283±298

www.elsevier.com/locate/engfracmech

Computational aspects of the discontinuous deformation analysis framework for modelling concrete fracture C.J. Pearce*, A. Thavalingam, Z. Liao, N. BicÂanic Department of Civil Engineering, University of Glasgow, Glasgow G12 8LT, Scotland, UK

Abstract Some computational aspects of the Discontinuous Deformation Analysis (DDA) modelling framework as applied to the modelling of concrete fracturing are discussed. In particular, the ability of DDA to model concrete behaviour through the entire range of fracture evolution from continuum to discontinuum is considered. Related topics including displacement controlled DDA analysis and the monitoring of progressive failure via eigenvalue analysis are also presented. # 2000 Elsevier Science Ltd. All rights reserved. Keywords: DDA; Concrete; Fracturing; Prescribed displacement; Failure indicator

1. Introduction Computational modelling of concrete structures has for some time moved well beyond providing only an estimate of the ultimate load behaviour Ð the importance of tracing the structural behaviour beyond the peak and an ability to follow the development of the actual failure mechanism including the transition to a discontinuum state in a robust and stable manner has been recognised as vital in any rational assessment of structural integrity of complex structures. A comprehensive modelling approach capable of covering both continuum and discontinuum states of concrete structures is still a long way o€. Such a process would comprise damage and fracturing of an originally continuous medium, the transition from a continuum to a discontinuum in a rational manner, followed by total separation and possible interaction of fragmented structural parts. However, there are a number of noticeable developments where the best practices of continuum based modelling and discontinuum based modelling are rapidly converging closer together. * Corresponding author. Fax: +44-141-330-4557. E-mail address: [email protected] (C.J. Pearce). 0013-7944/00/$ - see front matter # 2000 Elsevier Science Ltd. All rights reserved. PII: S 0 0 1 3 - 7 9 4 4 ( 9 9 ) 0 0 1 2 1 - 6

284

C.J. Pearce et al. / Engineering Fracture Mechanics 65 (2000) 283±298

Modelling of conditions for gradual fracturing (strain localisation prior to separation by cracking or shear slip) presents a complex physical problem in its own right. Nonlinear continuum based material models adopted are usually based on concepts of damage mechanics, strain softening plasticity formulations or some higher order continuum theory; however, often no discontinuities are admitted in the displacement ®eld. In general, continuum based approaches have been stretching their applicability by either enhancing the mathematical framework used for material models (higher order continua, gradient models, softening plasticity), or by enrichment of the trial function space in the actual discretisation process (mesh adaptivity and alignment, discontinuous shape functions, internal bands). In all such formulations, the identi®cation of the actual failure mode implies some form of continuity, i.e. separation of parts of the structure is not considered. On the other hand, the discrete crack framework implies a fracture mechanics based criteria for advancement of the crack front, coupled with extensive local and global remeshing. Such frameworks introduce a displacement discontinuity across the crack and the ®nal failure mode is readily identi®ed. Single crack front propagation is typically controlled by an energy based advancement criterion [1] and tracing of multiple cracks leads to both topological and algorithmic complexities. Usual formulations assume that cracks, once open, do not close, i.e. there is no consideration of a possible re-establishment of contact between the crack surfaces. In a similar way, modelling of concrete via the lattice based approaches [2] has also been very successful as a discontinuum framework, especially in interpreting meso mechanics of concrete fracture. Discontinuous modelling frameworks have become increasingly researched in concrete fracture simulations; such frameworks assume a discontinuum state a priori and come under the broad name of discontinuous or discrete element methods. The distinct feature of these formulations is a treatment of a structure comprising a large number of interacting parts with a constantly changing problem geometry and contact conditions. The Discontinuous Deformation Analysis (DDA) method, the subject of this paper, falls under this umbrella. DDA originated in the ®eld of rock mechanics [3] as an ecient framework for modelling jointed rock as a deformable blocky medium. It has more recently been reformulated as a subset of a more generalised framework which is seen as an alternative approach, along with a number of approximation procedures, for modelling discontinuous media. The most recent generalisation is the Manifold Method [4,5] which advocates similar ideas as the ones used in the Element Free Galerkin Method [6]. Here, the original simply deformable DDA formulation is adopted [3], which can be seen as an alternative to discrete ®nite elements [7] by introducing solid deformability into the discrete element framework. The original formulation is restricted to arbitrary shaped blocks (either convex or nonconvex) in two-dimensions with a constant strain state over the entire area of each block. Improved model deformability within DDA may be achieved by either introducing higher order strain ®elds for each block or by the so-called sub-block concept [8], in which an individual block is subdivided into a set of simply deformable sub-blocks. The latter approach introduces a piecewise constant strain ®eld into each parent block thereby enhancing its deformability. Block to block contact is controlled, whereby a no penetration condition and a friction controlled sliding condition are both satis®ed with varying degrees of accuracy depending on the method adopted Ð Penalty Method [3], Lagrange Multiplier Method or Augmented Lagrangian Method [8,9]. Di€erent ways of imposing these contact constraints lead to di€erent algorithmic sti€ness matrices in DDA and some of the resulting implications on monitoring the evolution of failure will be addressed later. In simulating progressive fracturing of concrete structures it is necessary to apply not only force controlled but also displacement controlled loading. Displacement controlled loading enables the simulation of the entire response spectrum of progressively fracturing media, including post peak softening; several issues related to the displacement controlled approach for DDA are outlined here.

C.J. Pearce et al. / Engineering Fracture Mechanics 65 (2000) 283±298

285

Finally, the use of linearised eigenvalue analyses within the DDA framework is investigated in order to provide a means for rational assessment of structural integrity.

2. Sti€ness formulation in DDA Within the basic DDA framework individual blocks, of arbitrary geometry, are modelled as a simply deformable medium with six degrees of freedom per block. These six variables comprise the centroidal translations, u0 , v0 , centroidal rotation r0 and the three strain components, Ex , Ey , Exy : The block deformation vector, di , is therefore ÿ d i ˆ u0

v0

r0

Ex

Ey

Exy


T

…1†

The displacement components at any position within the block are expressed as ui …x, y † ˆ … u

v† T ˆ T i d i

…2†

where the transformation matrix T i is de®ned as  1 0 ÿ…y ÿ y0 † …x ÿ x 0 † 0 Ti ˆ 0 0 1 …x ÿ x 0 † …y ÿ y0 †

…y ÿ y0 †=2 …x ÿ x 0 †=2

 …3†

The equilibrium equations are derived from the minimisation of potential energy where the total potential energy is composed of contributions from block strain energy, point and body loads, block to block contacts, displacement constraints etc. For a full derivation of each of these contributions the reader is directed to the seminal work of Shi [3]. The di€erentiation of the individual energy contributions results in sti€ness terms corresponding to each of these contributions. All sti€ness contributions can then be collected to form a global system of linear equations of the form 2

K 11 6 K 21 6 6 6 .. 4. K n1

K K 12 K 22 .. . K n2







.. . ...

3

K 1n K 2n 7 7 7 .. 7 . 5 K nn

0

d

1

f

1 d1 f1 B d2 C B C B C ˆ B f2 C B C B C B .. C B .. C @. A @. A dn fn 0

…4†

For an n block system, where each block has six degrees of freedom, the global sti€ness matrix will be 6n  6n in size. In Eq. (4) the o€-diagonal contributions to the system matrix are present only if and when blocks are in contact. In other words, the existence of sub-matrices K ij and K ji (6  6 in size) is the result of contact between blocks i and j. More generally, the diagonal sub-matrices will be a€ected by any source of potential energy which is restricted to the block itself, whereas both the diagonal and the o€-diagonal sub-matrices comprise contributions from the potential energy associated with any interaction between blocks.

286

C.J. Pearce et al. / Engineering Fracture Mechanics 65 (2000) 283±298

3. Active contact constraints and algorithmic representation Imposition of contact constraints represents a source of potential energy and is consequently responsible for a contribution to the system sti€ness matrix. Within the DDA framework these conditions are satis®ed only in an approximate manner. Clearly, di€erent algorithmic treatment of contact constraints results in di€erent contributions to the resulting DDA system sti€ness matrix, which in turn has implications on the system conditioning and convergence. For example, when the Mohr±Coulomb with tension cut-o€ is considered as an interface failure condition, depending on the inter-block contact stress state, the contact conditions can be either . No penetration, limited tension. Sliding allowed, or . No penetration, limited tension. Sliding not allowed. These conditions can be incorporated as block displacement constraints which are algorithmically reduced to an interaction problem between the vertex of one block and the side of another. Considering only the zero normal penetration condition, D ˆ 0, the simplest approach, initially utilised by Shi [3], is the Penalty Method, where a contact spring, with sti€ness pc , is introduced between the vertex of one block normal to the side of the other, implying a contact force equal to f ˆ pc D

…5†

The advantage of this method is its simplicity, however, the problem often leads to an ill-conditioned system matrix when a suitably large penalty sti€ness is adopted to ensure a near zero penetration; too low a sti€ness and the contact condition will not be satis®ed to any degree of accuracy. Moreover, a non-zero penetration is necessary for a ®nite contact force to exist. The Lagrange Multiplier Method introduces contact forces, l, wherever contact between blocks occur, as an additional variable. This approach ensures the contact condition is satis®ed exactly, but it increases the number of equations which need to be solved and indeed can change as the solution progresses. For a single contact between two blocks, Eq. (4) is therefore augmented as follows ! # ! " f d K kTl ˆ …6† l kl 0 fl It should be noted that this method introduces a zero onto the leading diagonal of the system matrix for every constraint introduced and therefore causes the system matrix to be non-positive de®nite. The third approach, the Augmented Lagrangian Method, solves for the contact forces, l, via an iterative combination of a Lagrange Multiplier and a contact penalty spring. This method, advocated by Lin, [8], is probably the most advanced treatment of contact constraints in the context of DDA. Unlike the Lagrange Multiplier Method, Eq. (6), the Augmented Lagrangian Method does not increase the size of the system; instead l is obtained via consecutive iterations li‡1 ˆ li ‡ pc Di

…7†

where pc can theoretically be set to any value. For each iteration, Eq. (4) now becomes K d ˆ f ‡ f li

…8†

where  T f li ˆ 0, 0, . . . , f lki , . . . , f ll i , . . . ,0

…9†

C.J. Pearce et al. / Engineering Fracture Mechanics 65 (2000) 283±298

287

In Eq. (9) f lki and f ll i are equivalent load vectors corresponding to the current iterate of li resulting from contact between blocks k and l. The iterative procedure continues until the zero penetration condition for Di is satis®ed within some tolerance. As the size of the system matrix does not increase, the method retains the simplicity of the penalty method; however, the contact force can now be found without problems associated with an ill-conditioned system matrix, although the choice of the penalty sti€ness a€ects the convergence rate.

4. Displacement constraints and displacement controlled loading In order to simulate a wider spectrum of structural problems and to enable tracing of the complete structural response past the ultimate load carrying capacity, it is necessary to control structural loading via prescribed displacements. To constrain the movement of a point within a block i in a speci®ed direction, a very sti€ spring, in the same direction, is assumed to be connected to a fully ®xed imaginary boundary. In order to ensure a  is utilised. The directionality of the sti€ near zero displacement a relatively large spring sti€ness, p, spring is de®ned by direction cosines lx and ly : Any deformation of the spring d is therefore the projection of the point's displacement in the direction of the spring
ÿ …10† d ˆ lx u ‡ ly v where u and v are the displacement components of the constrained point. Therefore, the strain energy of this sti€ spring is
2 1 2 1 ÿ  ˆ p lx u ‡ ly v Pc ˆ pd 2 2 or

 1 ÿ Pc ˆ p u 2

 
lx ÿ lx v ly

 
u ly v

Recalling u ˆ Td, Eq. (2), and de®ning l ˆ … lx matrix form as  1  Pc ˆ p dTi T Ti l lT T i di 2

…11†

…12† ly †T , the strain energy of the spring may be written in …13†

Minimising the spring strain energy with respect to the unknown deformation vector di results in a set of six linear equations @ Pc  Ti l lT T i di ˆ 0 ˆ pT @ di

…14†

These equations de®ne the equilibrium condition corresponding to the displacement constraint from which the contribution to the global sti€ness matrix is can be written as  Ti l lT T i ˆ)K ii pT

…15†

Fully ®xed displacement boundary conditions for a given point can easily be established with two noncollinear directional springs and complex boundary conditions for DDA blocks by various combinations

288

C.J. Pearce et al. / Engineering Fracture Mechanics 65 (2000) 283±298

of directional sti€ springs. High values for these boundary springs have similar implications on the system conditioning as the active constraint penalty terms, emanating from inter-block contacts. In a similar way, in the standard DDA framework an additional set of three weak springs (two translational and one rotational for 2D) is assumed to exist at every block centroid, connected to a fully ®xed imaginary boundary. Thereby the block sub-matrix for any block is rendered non-singular a priori, even for the case where an individual block is free from the contact with any other block or boundary. The DDA system sti€ness matrix is therefore always positive de®nite, even for the fragmented state. In order to allow for displacement controlled analysis, the above formulation is augmented by prescribing the displacement in a speci®ed direction. Again a sti€ spring connected to a fully ®xed imaginary boundary is utilised but the displacement of the point is now prescribed to be d in the direction de®ned by lx and ly : The force due to the prescribed displacement, resulting from the connection to the sti€ spring, is f ˆ  The corresponding potential energy due to this force is pd: ÿ
 lx u ‡ ly v Pd ˆ ÿpd

…16†

Therefore, the total potential energy due to the prescribed displacement is
2 ÿ
1 ÿ  lx u ‡ ly v Pc ˆ p lx u ‡ ly v ÿpd 2

…17†

where the ®rst term is the potential energy of the sti€ spring and the second term is the potential energy Pd of the imposed force. Eq. (17) can be rewritten as  1 ÿ Pc ˆ p u 2

v

 
lx ÿ lx ly

ly

   
u ÿ
lx  u v ÿ pd ly v

…18†

or in matrix form as
1 ÿ  Ti T Ti l Pc ˆ p dTi T Ti l lT T i di ÿ pdd 2

…19†

Minimising the spring strain energy results in a set of equilibrium equations @ Pc  Ti l lT T i di ÿ pdT  Ti l ˆ 0 ˆ pT @ di

…20†

Thus, for a prescribed displacement, the contributions to the global sti€ness matrix and the global force vector are as follows  Ti l lT T i ˆ) K ii pT

…21†

 Ti l ˆ) f i pdT

…22†

A simulation of the pull-out of an anchor bolt from mass concrete under displacement controlled loading is presented in Section 5 where it can be seen that the post peak response is captured to some extent.

C.J. Pearce et al. / Engineering Fracture Mechanics 65 (2000) 283±298

289

5. Fracturing in DDA Analysis of progressively fracturing media, like concrete, requires both the de®nition of the fracture initiation criterion, as well as some control of the direction of fracture evolution and propagation. In most modelling frameworks, the stress state is checked against the fracture initiation criteria, whereas the propagation direction is either assumed to be along predetermined discontinuities or it is de®ned by some energy based directionality criterion. Topologically, DDA simulations of fracturing can be modelled by releasing predetermined block interfaces or by through block fracturing, thereby creating new blocks. The DDA formulation deals naturally with discontinuities along block boundaries, typically in the sense of a Mohr±Coulomb criterion (with or without a tension cut-o€) and such a fracturing formulation resembles early concrete fracturing simulations by Ngo and Scordelis [10] in the ®nite element context. The principal di€erence and advantage here is that the continuous rede®nition of contact constraints, as the fracturing transforms the initially solid continuum structure into an assembly of partially or fully discontinuous parts, is easily implemented. The traditional algorithmic argument in DDA is simple; the normal and tangential contact springs, which are added to the system whenever a vertex is in the vicinity of a side of another block, are released instantaneously if a pointwise Mohr-Coulomb or tension cut-o€ criterion (in terms of stress resultants) is violated. Similar possibilities exist with the combined FEM/DEM, where the Rankine criterion was utilised to invoke local fracturing and remeshing [11], or with joint/interface ®nite elements [12]. A block fracturing algorithm through block centroids has been proposed recently by Lin [8] in the context of rock fracture, comprising again the Mohr±Coulomb fracturing criterion with a tension cut-o€, where newly formed discontinuities are introduced and are treated in the same way as the original discontinuity planes. The limited deformability associated with the basic DDA framework can lead to solution locking (for example, as a result of too many contact constraints that need to be satis®ed all at the same time) and thus various ways of enhancing the description of the block deformability have been suggested. Similar to the combined FEM/DEM formulations, the DDA sub-block framework [8] assumes that each sub-block is again restricted to a constant strain state. The additional constraint condition between sub-block boundaries remain valid at all times and their sole purpose is to enhance the description of the block deformability. Through block fracturing is also valid for sub-blocks. Concrete is increasingly regarded as a quasi brittle, rather than brittle, material and the fracture energy controlled reduction of material strength has been successfully adopted in continuum based softening simulations of localised failure of concrete. Implications for discontinuous modelling frameworks are that the contacts need to be released gradually and should be associated with some energy dissipation. Within the discontinuous RBSM (Rigid Body-Spring Model) framework, Kitoh et al. [13] have introduced fracture energy control into a block separation criterion by utilising the relative displacement between the two rigid blocks, which is in turn responsible for a reduced interface strength. A similiar technique to allow for energy dissipation in block seperation has been utilised by Meguro and Tagel-Din in the context of the Extended Distinct Element Method. [14]. Equivalent ideas can be applied to the DDA framework, which di€ers from the RBSM framework principally in the way that the medium deformability is accounted for Ð in RBSM it is concentrated on block interfaces, whereas in DDA it is located in blocks themselves. An equivalent cracking strain across block boundaries in the DDA could be used to determine a non-zero crack opening across a given interface (i.e. a ®nite gap, rather than zero) which the DDA algorithm needs to satisfy, and which is at the same time associated with a reduced interface strength. Full separation of DDA blocks would only follow after the crack opening

290

C.J. Pearce et al. / Engineering Fracture Mechanics 65 (2000) 283±298

Fig. 1. Kitoh plain concrete beam …h ˆ 100 mm). (a) Dimensions, loading and constraints. (b) Failure mode for DDA coarse mesh (130 blocks) (c) Failure mode for DDA ®ne mesh (455 blocks).

C.J. Pearce et al. / Engineering Fracture Mechanics 65 (2000) 283±298

291

reaches a critical value, which is in turn controlled by the fracture energy release rate at the interface in question. Several discontinuous modelling attempts in the simulation of fragmentation of concrete structures have been reported. Recently Kitoh [13] conducted a series of block size sensitivity tests using RBSM, exploring both the geometric size e€ect on fracture for a series of "normalised' concrete beams, as well as the in¯uence of the discretisation e€ect, resulting from di€erent block sizes. The Kitoh plain concrete four point beam bending problem …h ˆ 100 mm), Fig. 1(a), is here modelled with two DDA discretisations based on Voronoi tessellations, where no attempt has been made to study the detailed sensitivity with respect to block shape or size. The failure load predictions from a coarse discretisation using 130 simply deformable blocks, Fig. 1(b), is compared to the failure load prediction from a model comprising 455 simply deformable blocks, Fig. 1(c), and both results are set against the failure load prediction reported by Kitoh [13] using RBSM, see Fig. 2. All analyses were conducted using the dynamic relaxation technique. The block material characteristics and concrete/concrete interface material data are identical in both cases: Young's modulus E ˆ 27:5 MPa Poisson's ratio n ˆ 0:20 Cohesion c ˆ 0:0047 MPa

Fig. 2. Kitoh plain concrete beam …h ˆ 100 mm). Comparison of load-displacement diagrams for coarse and ®ne DDA discretisations, together with RBSM failure load.

292

C.J. Pearce et al. / Engineering Fracture Mechanics 65 (2000) 283±298

Tensile strength ft ˆ 0:0029 MPa Friction angle f ˆ 37:08 The analysis has been conducted using an enhanced version of Shi's original DDA source code, under force controlled loading, with adaptive load incrementation approaching failure. An additional simulation of a hypothetical reinforced concrete beam based on the Kitoh beam geometry, Fig. 3, comprises a steel/concrete interface law of Cohesion and c ˆ 0:00188 MPa

Fig. 3. Hypothetical reinforced concrete beam Ð identical geometry as Kitoh beam …h ˆ 100 mm). (a) Failure mode. (b) Detail of crack zone.

C.J. Pearce et al. / Engineering Fracture Mechanics 65 (2000) 283±298

293

Tensile strength ft ˆ 0:00145 MPa Friction angle f ˆ 22:58 whereas the steel/steel interface is fully bonded and the concrete/concrete interface is identical to the one utilised in the plain concrete beam analysis. The ®nal failure mode indicates an arrest of the initial crack followed by a slip failure of the concrete/steel interface resulting in a shift of the major crack location along the beam. The second example reported is the prescribed displacement controlled loading of the pull-out of an anchor-head from plain concrete. The experiments, similar to the RILEM pull-out round-robin analysis, were carried out in the Stevin Laboratory, Delft University of Technology and reported by Vervuurt et al. [2] and Feenstra [15]. In the simulations reported here, only the bolt head was modelled to which a prescribed displacement was applied. Fig. 4(a) shows the initial DDA block assembly and Fig. 4(b) shows a close up of the anchor bolt head when the analysis terminated. Fig. 5 shows the resulting load-displacement diagram from the DDA analysis. Comparison with the experimental results show that the peak load reported here is approximately 28% less than those of the experiments and that the post peak response, which terminates as a result of non-convergence of the algorithm, is overly brittle. However, this example does highlight the ability of the displacement controlled loading to allow the post-peak response to be followed to some degree but also illustrates the need to model interface fracturing in a progressive manner as opposed to the instantaneous release approach adopted here. Furthermore, it is clear that more work must be done on the development of

Fig. 4. DDA simulation of anchor pull-out. [2]. (a) Geometry (dimensions in mm) and DDA block assembly. (b) Close-up of anchor head at the end of analysis.

294

C.J. Pearce et al. / Engineering Fracture Mechanics 65 (2000) 283±298

Fig. 5. DDA simulation of anchor pull-out. [2]. Load vs. displacement of anchor head (mean experimental peak load 22 kN corresponding to mean experimental displacement of anchor head of 0.008 mm).

algorithms which can enable phenomena such as softening, snap-back and snap-through to be captured. In the context of the Finite Element Method a lot of research has been done into such problems (such as the arc-length method) and similar strategies could be adopted for DDA. 6. Monitoring of progressive fracture and failure indicators Various computational frameworks deal di€erently with an indication of failure or a limit state, either on a local level (singularity of the material operator and the acoustic tensor) or on a global level (singularity of the tangent sti€ness matrix). In concrete fracturing, the transition criteria from a continuum to discontinuum are more or less rigorously de®ned in the context of di€used or localised failure, smeared or hard discontinuities. Here a failure indicator strategy in DDA, similar to the current sti€ness concept, [16], is suggested. The strategy is based on an eigenvalue analysis of the linearised DDA system sti€ness matrix, where speci®c DDA dependent features and sti€ness contributions will be accounted for. Various potential energy contributions to the system sti€ness matrix have been considered earlier and Fig. 6 illustrates all contributions, associated with the Augmented Lagrangian DDA formulation, for a simple three block DDA assembly (Blocks A, B and C ). The algorithmic sti€ness matrix in this case comprises sti€ness terms emanating from block deformabilities K 0 , terms associated with an approximate imposition of displacement contact constraints K c (penalty spring pc ), as well as sti€ness terms emanating from displacement boundary conditions K s (penalty spring ps ). These constraints, together with the sti€ness contributions from the weak link or `anti-singularity' springs, K w (penalty

C.J. Pearce et al. / Engineering Fracture Mechanics 65 (2000) 283±298

295

Fig. 6. Example of DDA block assembly illustrating springs associated with weak links, boundary constraints and interface constraints.

spring pw † appear in the algorithmic DDA system sti€ness matrix. Eq. (23) shows symbolically the algorithmic matrices associated with the Augmented Lagrangian DDA formulation for the three block problem at the ith iteration, 3 2 …AÿB † c…AÿB † ‡ K sAA ‡ K w Š ‰K Š ‰0Š ‰K 0AA ‡ K cAA AA AB 7 6 7 6 c…AÿB † …BÿC † …AÿB † c…BÿC † 7 6 ‰K AB Š ‰K 0BB ‡ K cBB ‡ K cBB ‡ Kw Š ‰K Š BB BC 5 4 … † c…BÿC † c BÿC 0 w ‰K CC ‡ K CC ‡ K CC Š ‰0Š ‰K BC Š …23† 1 0 1 0 f A ‡ f A …pc , li † dA  @ dB A ˆ @ f B ‡ f B …pc ,li † A dC f C ‡ f C …pc , li † Here f A is the force vector associated with external loads acting on block A and f A …pc , li † is the force vector associated with contact with other blocks which includes the contact force from the previous iteration. In the context of the Augmented Lagrangian format, it is interesting to consider the di€erence between the algorithmic and physical importance of the K c terms, associated with the contact spring pc : It has been argued earlier that the Augmented Lagrangian procedure advocates the iterative use of a low penalty term, which ensures a correct contact force l in an iterative way, and that this leads to a better conditioned system matrix, as opposed to the penalty method. However, the presence of the low penalty term does not correspond to the block contact conditions, and therefore does not allow for the linearisation to be made using this algorithmic system sti€ness matrix. Consequently, the algorithmic DDA system sti€ness matrix used for advancing the nonlinear solution within the Augmented Lagrangian context cannot be used to assess the current eigenvalue of the

296

C.J. Pearce et al. / Engineering Fracture Mechanics 65 (2000) 283±298

physical system sti€ness matrix. It needs to be linearised in such a way that it re¯ects the current block topology as well as block contact constraints in terms of sti€ness contributions only. This leads to the requirement to replace the low penalty terms pc (desirable for the Augmented Lagrangian) with the very high values of contact penalty terms pc for the purpose of the current sti€ness eigenvalue analysis. Eq. (24) shows the DDA system sti€ness matrix for a linearised eigenvalue analysis for the three block problem, where K c has replaced K c : 2

…AÿB † c…AÿB † ‰K 0AA ‡ K cAA ‡ K sAA ‡ K w Š AA Š ‰K AB

6 6 …AÿB † K  ˆ 6 ‰K cAB Š 4 ‰0Š

3

‰0Š

…AÿB † …BÿC † ‰K 0BB ‡ K cBB ‡ K cBB ‡ Kw BB Š

…BÿC † ‰K cBC Š

…BÿC † Š ‰K cBC

…BÿC † ‰K 0CC ‡ K cCC ‡ Kw CC Š

7 7 7 5

…24†

Fig. 7. Cantilever problem, illustrating DDA block assembly and deformed shape depending on extent of propagating fracture.

C.J. Pearce et al. / Engineering Fracture Mechanics 65 (2000) 283±298

297

Fig. 8. Eigenmodes and eigenvalues corresponding to di€erent stages of fracturing and the evolution of the current sti€ness ratio l=l0 :

A benchmark problem involving the progressive fracturing of a cantilever beam subjected to a point load, Fig. 7, is given to illustrate the changes in the eigenmodes associated with a linearised eigenvalue analysis at various stages. Deformed shapes associated with each stage of fracturing are also shown in Fig. 7. Fig. 8 illustrates the changes that occur in the four lowest eigenvalues of the cantilever beam during progressive fracturing and illustrate the evolution of the eventual discontinuous failure mode, where it is interesting to note that some deformation modes are una€ected by the presence of progressive discontinuities. The current sti€ness ratios l=l0 are also given for various stages and various deformation modes.

7. Conclusions Several computational aspects associated with the use of the Discontinuous Deformation Analysis in modelling progressive fracturing of concrete have been discussed. The displacement controlled loading, required for tracing an overall structural response beyond the limit point, has been considered in detail and illustrated on plain concrete four point bending and anchor pull-out problems. It is noted that the algorithmic DDA sti€ness matrix in the Augmented Lagrangian approach needs to be replaced by an appropriate DDA sti€ness matrix, to properly re¯ect the block contact constraints for the linearised eigenvalue analysis. Current eigenvalues and current sti€ness ratios can be used as failure indicators in simulations of progressive fracturing of concrete.

298

C.J. Pearce et al. / Engineering Fracture Mechanics 65 (2000) 283±298

References [1] Ingra€ea A. FE methods for discrete crack modelling of concrete structures. Lecture notes, CISM Course, Udine, Italy, 1997. [2] Vervuurt A, Schlangen E, Van Mier JGM. A numerical and experimental analysis of the pull-out behaviour of steel anchors embedded in concrete. Technical Report 25.5-93-1/VFL, Faculty of Civil Engineering, Delft University of Technology, 1993. [3] Shi G-H. Discontinuous deformation analysis: a new numerical model for the statics and dynamics of block systems. Ph.D., University of California, Berkeley, 1988. [4] Shi, G-H. Numerical manifold method. In: Ohnishi Y, editor. Proceedings of ICADD-2, The Second International Conference on Analysis of Discontinuous Deformation, 1997, pp. 1±35. [5] Chen G, Ohnishi Y, Ito T. Development of high-order manifold method. International Journal for Numerical Methods in Engineering 1998;43:685±712. [6] Belytschko T, Lu Y, Gu L. Element Free Galerkin method. International Journal for Numerical Methods in Engineering 1994;37:229±56. [7] Ghaboussi J. Fully deformable discrete element analysis using a ®nite element approach. International Journal of Computers and Geotechnics 1997;5:175±95. [8] Lin CT. Extensions to the discontinuous deformation analysis for jointed rock masses and other blocky systems. Ph.D., University of Colorado, 1995. [9] Amadei B, Lin C, Dwyer, J. Recent extensions to the DDA method. In: Salami MR, Banks D, editors. DDA and simulations of discontinuous media. TSI Press, 1996, pp.1±30. [10] Ngo D, Scordelis AC. Finite element analysis of reinforced concrete beams. ACI Journal 1967;64(3):152±63. [11] Munjiza A, Owen DRJ, BicÂanic N. A combined ®nite/discrete element method in transient dynamics of fracturing solids. Engineering Computations 1995;12:145±74. [12] Rots JG, Schellekens JCJ. Interface elements in concrete mechanics. In: BicÂanic N, Mang H, editors. Computer aided analysis and design of concrete structures. Swansea: Pineridge, 1990, pp. 909±18. [13] Kitoh H, Takeuchi N, Ueda M, Higuchi H, Kambayashi A, Tomida M. Size e€ect analysis of plain concrete beams by using RBSM. In: Ohnishi Y, editor, Proceedings of ICADD-2, The Second International Conference on Analysis of Discontinuous Deformation, 1997, pp. 373±82. [14] Meguro K, Tagel-Din T. A new simpli®ed and ecient technique for fracture behaviour analysis of concrete. In: Mihashi H, Rokugo K, editors, Fracture mechanics of concrete structures (FRAMCOS-3), Aedi®catio, 1999, pp. 911±20. [15] Feenstra P. Computational aspects of biaxial stress in plain and reinforced concrete. Ph.D., Civil Engineering Department, Delft University of Technology, The Netherlands, 1993. [16] Ramm E. Strategies for tracing nonlinear responses near limit points. In: Wunderlich W, Stein E, Bathe K, editors. Nonlinear ®nite element analysis in structural mechanics. Berlin: Springer-Verlag, 1981.