CONSTITUTIVE MODELING OF FIBER REINFORCED CEMENT COMPOSITES by MOHAMED BOULFIZA Ingeniorat d'etat, ENTP Algiers, 1990 M.A.Sc, Universite Laval, 1993 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENT FOR THE DEGREE OF DOCTOR OF PHILOSOPHY in THE FACULTY OF GRADUATE STUDIES (Department of Civil Engineering) We accept this thesis as conforming to the required standard THE UNIVERSITY OF BRITISH COLUMBIA April 1998 © Mohamed Boulfiza, 1998 In presenting this thesis in partial fulfilment of the requirements for an advanced degree at the University of British Columbia, I agree that the Library shall make it freely available for reference and study. I further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the head of my department or by his or her representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission. Department The University of British Columbia Vancouver, Canada DE-6 (2/88) A B S T R A C T The role of fibers in the enhancement of the inherently low tensile stress and strain capacities of fiber reinforced cementitious composites ( F R C ) has been addressed through both the phenomenological, using concepts of continuum damage mechanics, and micromechanical approaches leading to the development of a closing pressure that could be used in a cohesive crack analysis. The observed enhancements in the matrix behavior is assumed to be related to the ability of the material to transfer stress across cracks. In the micromechanics approach, this is modeled by the introduction of a nonlinear closing pressure at the crack lips. Due to the different nature of cracking in the pre-peak and post peak regimes, two different micro-mechanical models of the cohesive pressure have been proposed, one for the strain hardening stage and another for the strain softening regime. This cohesive pressure is subsequently incorporated into a finite element code so that a nonlinear fracture analysis can be carried out. On top of the fact that a direct fracture analysis has been performed to predict the response of some F R C structural elements, a numerical procedure for the homogenization of F R C materials has been proposed. In this latter approach, a link is established between the cracking taking place at the meso-scale and its mechanical characteristics as represented by the Young's modulus. A parametric study has been carried out to investigate the effect of crack patterning and fiber volume fractions on the overall Young's modulus and the thermodynamic force associated with the tensorial damage variable. After showing the usefulness and power of phenomenological continuum damage mechanics ( P C D M ) in the prediction of F R C materials' response to a stimuli (loading), a combined P C D M - N L F M 1 approach is proposed to model (predict, forecast) the complete response of the composite up to failure. Based on experimental observations, this approach assumes that damage mechanics which predicts a diffused damage is more appropriate in the pre-peak regime whereas, N L F M is more suitable in the post-peak NLFM: Non Linear Fracture Mechanics ii stage where the opening and propagation of a major crack w i l l control the response of the material and not a deformation in a continuum sense as opposed to the pre-cracking zone. Tensile and compressive tests have been carried out for the sole purpose of calibrating the constitutive models proposed and/or developed in this thesis for F R C materials. The suitability of the models in predicting the response of different structural members has been performed by comparing the models' forecasts with experimental results carried out by the author, as well as experimental results from the literature. The different models proposed in this thesis have the possibility to account for the presence of fibers in the matrix, and give fairly good results for both high fiber volume fractions (v f low fiber volume fractions (v f > 2%) and < 2%). Use of interface elements in a finite element code has been shown to be a powerful tool in analyzing the behavior of concrete substrateF R C repair materials by the introduction of a zero thickness layer of interface elements to account for the interface properties which usually control the effectiveness of the repair material. iii TABLE OF CONTENTS Abstract ii Table of Contents iv List of Tables vii List of Figures viii Acknowledgement xiii L Introduction and Objectives II. Thermodynamics of Continuous Media - Formulation of Constitutive Equations - Thermodynamics of Continua - The First Principle of Thermodynamics - The Second Principle of Thermodynamics - Complimentary Laws - Application I: Continuum Damage Mechanics - Example: Brittle Elastic Materials - Application JJ: Fracture Mechanics 7 7 11 11 12 15 15 21 23 in. Homogenization of Heterogeneous Media - Introduction - Microscopic and Macroscopic Scales - Representative Volume Element - Description of Few Homogenization Techniques - Variational Methods - Method of Effective M o d u l i - Self-Consistent Method - Periodic Homogenization - Comparative Example 32 33 35 39 42 42 45 48 50 Continuum Damage Mechanics - Introduction - Hypothesis of the effective stress - Damage Models for Cementitious Materials - Damage Model with Anisotropy and Permanent Deformation - Unilateral M o d e l - Scalar Damage M o d e l 57 57 64 66 68 70 72 IV. 1 iv T A B L E OF CONTENTS - V. Application of P C D M to F R C Materials - Beam Bending - Concluding Remarks - Plate Bending - Results and Discussion - Concluding Remarks 78 78 83 84 89 92 Fracture Mechanics - Introduction - Linear Elastic Fracture Mechanics - M i x e d M o d e Interaction Theories - Nonlinear Fracture Mechanics - Fictitious Crack Modeling - Two-Parameter Fracture M o d e l - R-Curve Based Models 95 97 103 107 108 110 112 VI. Numerical Modeling of Fracture Growth - Background and History - Smeared Modeling of Crack Propagation - Discrete Modeling of Crack Propagation - Interface Elements - Dynamic Relaxation Solver - Cohesive Crack Analysis - The Influence Method - Application to F R C Materials - Smeared Approach - Fracture-Based Approach - Concluding Remarks 118 118 120 125 131 137 141 146 153 159 160 162 VII. Micro-Mechanical Modeling of F R C Materials - Introduction - M o d e l Derivation - Strain-Hardening Zone - Strain-Softening Zone - Summary - Experimental Program - Direct Tensile Test - Bending Test - Compressive Test 165 165 169 170 173 178 179 180 184 187 VIH. Numerical Implementation and Validation - N L F M - B a s e d Predictions for F R C Structures - A Combined C D M - N L F M Approach - O n the Overall Behavior of F R C Materials 191 191 204 211 TABLE OF CONTENTS LX. Conclusions and Recommendations - Summary - Conclusions - Recommendations for Future Research 220 220 222 223 Bibliography 229 vi List of Tables Table 2.1 Number of field equations vs. number of unknown variables Table 2.2 Classification of the different field variables 16 Table 4.1 M i x proportions and fiber characteristics 78 Table 4.2 Comparison between model predictions and experimental results. Beam Bending 80 Comparison between model predictions and experimental results Plate Bending 92 Table 4.3 7 Table 7.1 M i x proportions and fiber characteristics 179 Table 7.2 Effect of fiber volume fraction on the tensile behavior 184 Table 7.3 Effect of fiber volume fraction on the bending behavior 185 Table 7.4 Effect of fiber volume fraction on the compressive behavior 188 Table 8.1 M o d e l predictions vs experimental results. Beam bending N L F M approach 193 M o d e l predictions vs experimental results. D C B N L F M approach (Steel fiber-Cement paste) 198 M o d e l predictions vs experimental results. D C B N L F M approach (Steel fiber - Mortar) 198 M o d e l predictions vs experimental results. D C B N L F M approach (Carbon fiber-Cement paste) 201 M o d e l predictions vs experimental results. D C B N L F M approach (Carbon - Mortar) 201 Table 8.2 Table 8.3 Table 8.4 Table 8.5 vii List of Figures Figure 1.1 Different stages in the response of a high V F R C Figure 2.1 A continuum body with a crack 24 Figure 2.2 Line crack in an infinite plate 29 Figure 2.3 Normalized critical stress versus the ratio of the semi-axes of an elliptical crack 30 Figure 3.1 Variation of density with the size of the rve 37 Figure.3.2 Macroscale vs micrscale representation of an rve 39 Figure 3.3 Example of an rve containing a cavity-type inclusion 40 Figure 3.4 Schematic representation of the porous medium and choice of an rve 51 Figure 3.5 f 2 Study of the rve associated to the porous medium subjected to Z n loading 52 Figure 3.6 Unit Cell for the isotropic medium 55 Figure 3.7 Comparative chart for the predictions of three homogenization techniques for a porous medium 56 Figure 4.1 A damaged element 57 Figure 4.2 Schematic Representation of the failure process of materials 59 Figure 4.3 Typical stress-strain curves for low fiber volume and high fiber volume F R C 60 M o d e l ' s capabilities to represent A (tensile behavior, and b) compressive behavior 75 Figure 4.5 Experimental setting of the Three point bending test 79 Figure 4.6 Bending behavior of C F R C materials with L =3 mm 81 Figure 4.7 Bending behavior of C F R C materials with L = 10 mm 82 Figure 4.8 Typical damage distribution as predicted by the model and experimental failure patterns 83 Figure 4.4 f f viii Figure 4.9 F l o w chart of the main steps followed by the finite element code 90 Figure 4.10 Schematic representation of the test setup 91 Figure 4.11 Comparison of the model's predictions to the experimental results for V = 2% 93 f Figure 4.12 Comparison of the model's predictions to the experimental results for V = 3% 93 Failure pattern as predicted by the model and observed experimental failure patterns 94 Figure 5.1 Stress distribution close to an elliptical hole in an infinite plate 98 Figure 5.2 Fracture process in concrete 99 Figure 5.3 The three basic modes of crack growth 99 Figure 5.4 Figure 5.5 Material element near a crack tip Elliptical crack at an angle to a uniform applied stress in an infinite medium 102 Figure 5.6 Principle of the fictitious crack model 109 Figure 5.7 Test for the determination of K Figure 5.8 Typical R- curve representation 114 Figure 5.9 Superposition of stress intensity factors at crack tip 116 Figure 6.1 Example of crack problems that can be represented using quarter symmetry 118 A discrete crack simulation using the technique of nodal release in reinforced concrete beam 119 Random microstructure, scatter of microstress, and crack band or sharp crack model s 121 Test configuration analyzed by Rots together with a representation of the deformed mesh and principal stresses 124 F l o w chart of the iterative scheme 127 f Figure 4.13 Figure 6.2 Figure 6.3 Figure 6.4 Figure 6.5 s ix Ic and CTOD c 101 111 Figure 6.6 Three different crack configurations using the splitting element approach 130 Figure 6.7 Interface idealization and notations 132 Figure 6.8 2 D nodal interface elements 133 Figure 6.9 Linear and quadratic 2 D interface elements 134 Figure 6.10 Schematic representation of stress-strain behavior associated with fracture behavior of strain softening materials 142 Figure 6.11 Cohesive crack problem 145 Figure 6.12 Discretized cohesive crack problem 146 Figure 6.13 Example of a softening model of interface elements 149 Figure 6.14 Typical force displacement responses 149 Figure 6.15 Zero-slope condition for the opening profile at the fictitious crack tip 152 Figure 6.16 Spalled concrete beam 154 Figure 6.17 Direct Tensile behavior of a C F R C Reinforced Repair composite (No Surface Preparation) 161 Direct Tensile behavior of a C F R C Reinforced Repair composite (With Surface Preparation) 161 A n example of the finite element mesh used in the fracture-based approach 163 Numerical Prediction of Crack Evolution in the Tensile test (With Surface preparation) 164 Numerical Prediction of Crack Evolution in the Tensile test (No Surface preparation) 164 Typical stress-strain curves for low fiber volume and high fiber volume F R C 166 Single fiber bridging a matrix crack 171 Figure 6.18 Figure 6.19 Figure 6.20 Figure 6.21 Figure 7.1 Figure 7.2 Figure 7.3 Figure 7.4 Schematic representation of a fiber pull out Pull out analysis of a single fiber bridging a crack 173 174 Figure 7.5 A fiber bridging a matrix crack at an angle 177 Figure 7.6 Shape and dimensions of the tensile test specimen 180 Figure 7.7 Schematic representation of the direct tensile test setup 182 Figure 7.8 Stress-strain curves in uniaxial tension For Lf =3mm and L = 10mm f 183 Figure 7.9 Schematic representation of the bending test setup Figure 7.10 Load deflection Curves for 3 point Bending Test For Lf =3 mm 184 and L =10 mm 186 f Figure 7.11 Schematic representation of the compression test 187 Figure 7.13 Compressive stress-strain curves for L =3 mm 189 Figure 7.14 Compressive stress-strain curves for L = 10 mm 190 Figure 8.1 Geometry and boundary conditions for the three point bending test 192 Figure 8.2 Numerical prediction at an advanced stage of loading 193 Figure 8.3 Experimental results versus numerical predictions for L =3 mm f f f Beam bending Figure 8.4 194 Experimental results versus numerical predictions for L = 10 mm f Beam bending 195 Figure 8.5 Specimen geometry and a schematic of the loading procedure 196 Figure 8.6 Numerical prediction of the crack growth at an intermediate stage 197 Figure 8.7 Numerical prediction vs experimental results. D C B (Steel fiber-Cement paste) 199 Numerical Prediction V S Experimental results. D C B (Steel fiber-Mortar) 200 Figure 8.8 xi Figure 8.9 Figure 8.10 Figure 8.11 Numerical Prediction V S Experimental results. D C B (Carbon fiber-Cement paste) 202 Numerical Prediction V S Experimental results. D C B (Carbon fiber- Mortar) 203 Schematic representation of the C D M - N L F M equivalence 209 Figure 8.11.bis Numerical prediction vs experimental results for carbon fiber- Mortar and steel fiber-Mortar mixes Figure 8.12 210 Numerical visualization of the main difference between dilute and interacting cracks in an R V E 213 Figure 8.13 A typical R V E prior to strain localization 215 Figure 8.14 Uniform crack distributions 217 Figure 8.15. Effect of uniform crack distribution on the effective moduli 218 Figure 8.16 Effect of fiber volume fraction on the effective moduli (uniform crack distribution) 219 Figure 8.17 Effect of uniform crack distribution on the Thermodynamic force associated with damage 220 Figure 8.18 N o n uniformCrack patterns 221 Figure 8.19 Effect of non uniform crack distribution on the effective moduli 222 xii Acknowledgement This research could not be performed without the help and assistance of many people. First and foremost, I would like to thank professor Nemkumar P. Banthia for giving me the opportunity to work on this exciting research project, support, trust and most importantly the freedom he has given me so that my research skills could develop and flourish. I am also grateful to professor S. Mindess and Dr. C. Yan for their constructive comments and help throughout these last few years. I am also particularly grateful to all personnel of the department of civil engineering who had in one way or another contributed to the success of this study. I have been particularly privileged to be surrounded by a very long list of fellow students and friends, the enumeration of which is impossible in this short note. They made my stay at UBC a very enjoyable one. I can never thank them enough for that. xiii The mathematicians and the physics men have their mythology; they work alongside the truth never touching it; their equations are false but the things work. Robinson Jejfers xiv Chapter 1 Introduction and Objectives INTRODUCTION AND OBJECTIVES INTRODUCTION Cement-based materials require reinforcement because of their low tensile strength and strain capacities. The tensile strain capacity of cement-based materials can be dramatically improved by the addition of fibers. In conventional FRCs {V < 1%) fibers are not added to improve the strength. Rather f their role is to affect the behavior of the material once the matrix has cracked and control the cracking of the FRC [3]. The fibers that bridge the cracks formed in the matrix contribute to the energy dissipation through the processes of debonding and pull-out. As a result, the post peak response continues to be of strain softening, but with a slope less steep compared to the unreinforced cementitious material. Thus the fibers improve the energy absorption capacity (or the toughness). It is generally accepted that for FRCs with V < 2%, the major contribution of the fibers is after the matrix strain localization, which f occurs around the peak of the tensile stress strain curve. However new processing techniques have helped in the manufacture of thin-sheet products with fiber volume fractions (V ) as high as 15%. In these composites type, fibers are added not only to f improve the ductility of the material, but also to increase the strength of the composite. This nonlinear increase in strength, termed strain hardening, is associated with the appearance of multiple cracking in the material which requires a higher energy input to open the microcracks. See Figure. 1.1. for a typical stress strain response of high V with f discontinuous fibers. 1 Chapter 1 Introduction and Objectives 00 Multiple Cracking Zone Strain Softening Zone Strain Strain Crack Opening Figure 1.1. Different stages in the response of a high V FRC f The response consists of 3 regions: I - elastic region (up to matrix cracking represented by point 1) II - multiple cracking region (up to max. post cracking point 2) III- Failure region (crack opening localization beyond point 2) The fracture process zone in high fiber volume fraction cementitious composites is not concentrated in a narrow band around the localized crack, but spreads over a large volume of the material. Thus, this kind of composites exhibit a high tensile strain and energy absorption capacity due to the extensive cracking off and parallel to the main crack plane. A n important feature of composite materials lies in the fact that their mechanical and fracture properties can be adjusted or tailored to the needs of specific applications. Indeed, a composite material can be designed for either high strength or 2 Chapter 1 Introduction and Objectives high ductility. As one would expect, the tailoring is done mainly through varying microstructural parameters associated with the three fundamental phases in a fiber composite (the fiber, the matrix, and their interface). Fiber volume fraction, fiber aspect ratio, matrix fracture toughness, fiber-matrix bond strength, etc. are examples of such parameters. Due to their enhanced fracture properties, cementitious composites can be useful in a wide range of engineering applications [6,30]. These include cases of structures under intensive shear loading, where the high ductility of the composite could help avoid a brittle catastrophic failure [107], and earthquake resistant structures, where the high fracture energy could increase the energy dissipation of the structure in a significant way. A fiber reinforced cementitious material would be very beneficial if applied at locations where plastic hinges are likely to occur, such as, near the beam-column connections. Another potential area of application for such materials, is their use as repair or patch materials [7,10] where their high tensile strain and stress capacities combined with their compatibility with the substrate concrete are used to create a high performance repair. Partly due to the difficulty of carrying out structural tests as compared to the material tests, there is a lack of knowledge regarding the behavior of fiber reinforced structural elements. This provides a rational for development of analytical models that provide the missing link between the laboratory test and the field performance. OBJECTIVES Powerful techniques and methods for the modeling of heterogeneous media have mushroomed over the last few decades. These methods which are very well established in the elastic case are developing when dissipative phenomena are considered. Unfortunately, most efforts along these lines deal with media other than cement-based materials. This is particularly the case of the aerospace industry where the need for a precise modeling of the material's behavior is a crucial issue. Despite the fact that 3 Chapter 1 Introduction and Objectives concrete is one of the oldest materials known to mankind, a precise modeling of its behavior has always lagged behind, at least partly, by the complicated nature and poor understanding of phenomena underlying its behavior. The quest for a rational constitutive law for FRC that reproduces its most characteristic features while linking them to what is happening at a mesoscopic level (level of heterogeneities) is still under way. In this thesis the term FRC refers to any kind of fiber reinforced cementitious composite. When modeling composite materials in general, and FRC in particular, two different aspects, which are not necessarily mutually exclusive, are usually considered: 1- Material aspect, related to the optimization of the composite by tailoring its constituents in such a way that the desired benefits are maximized, 2- Structural aspect, related the development of a constitutive law that could be used to characterize the material a given structure is made of, so that its response to a given stimuli can be computed. Depending on the final objective, two main modeling approaches are usually available. First, we have the phenomenological models, dealing with macroscopic quantities where some a priori defined field variables called internal variables (in a thermodynamical sense) are supposed to account for the local mechanisms without any explicit consideration of these latter ones. Secondly, we have the micromechanical-based models which provide a direct link between the structural performance of the composite and the local mechanisms underlying that behavior. This is usually achieved by using the socalled homogenization techniques, which are averaging methods consisting in lumping what is happening around a given material point in the heterogeneous material (microcracking, fiber slippage/breakage/debonding, etc.) into the constitutive law of a material point in an equivalent homogeneous material. The model's accuracy then depends on the adequacy and relevance of the mechanisms we choose to account for. In this thesis we give a description of both approaches along with advantages and inconveniences of each of them while considering their application to FRCs. Among the major benefits of analytical models, one can mention: 4 Chapter 1 Introduction and Objectives Interaction with experimental studies: numerical predictions can be used to explain experimental results or to optimize the actual experimental setup. The analyses may reveal, the existence of phenomena which should be further investigated. - Numerical simulations represent a good alternative to experimental testing which are too expensive or too difficult to carry out. A n example of such applications, is found in the size effect investigations where large specimens are usually needed. Another example is related to material's tailoring and the related parametric studies which would necessitate a very large number of tests. Here, a model assuming a discrete representation of cracking of the material is proposed using the generalized Dugdale-Barenblat cohesive model. This representation includes in a single stress-COD relationship, the complete behavior of the material from the first departure from linearity at the onset of microcracking in the pseudostrain hardening regime up to ultimate failure. This assumes that even in the multiple cracking zone, the tensile behavior of the material could be represented by the introduction of an "equivalent" cohesive crack that would exhibit the same material response as the one related to a representative volume element of the given material. The proposed model has been used to analyze a three point bending beam, and a double cantiliver beam for different fiber volume fractions and two lengths of carbon fibers. These analyses serve to check the validity and limitations of the proposed model and to further clarify the fracture behavior of cementitious composites. A combined continuum damage mechanics-nonlinear fracture mechanics approach has been proposed as an alternative procedure for modeling the nonlinear behavior of F R C composites. A finite element-based numerical procedure is used to compute the effective moduli of F R C composites together with the macroscopic damage variable and its associated thermodynamic force. A parametric study has been performed to investigate the effects of 5 Chapter 1 Introduction and Objectives fiber volume fraction and crack patterns in a cracked R V E on the evolution of the mechanical characteristics and the damage variable. Chapter 2 Thermodynamics of Continuous Media THERMODYNAMICS OF CONTINUOUS MEDIA FORMULATION OF CONSTITUTIVE EQUATIONS Everyday experience shows that constitutive elements of a given system play an essential role in its evolution: in spite of the fact of having the same geometry and subjected to identical loads, a metallic rod does not respond (deform) the same as a timber or an elastomeric rod, water does not flow like oil, . . . . The necessity of this characteristic behavior of materials becomes clear when we consider the movement problem. A s can be seen on Table 1, an inventory of the available field equations and the unknown fields shows that the equilibrium equations together with the equation of continuity form a set of 4 scalar equations in a given coordinate system whereas the unknown fields a,u and p provide 10 scalar unknowns. The 6 missing equations are none other than the constitutive equations that must relate tr to the other unknowns. In a general framework, such constitutive equations allow us to represent in a phenomenological manner the diversity of "continuous" materials encountered in Nature. Indeed, the field equations are valid for all types of media irrespective of their internal constitution. However, material bodies having the same mass and geometry, when subjected to identical external stimuli, respond differently. The internal constitution of matter is responsible for these different responses. Since we are not concerned with the atomic structure of matter in continuum mechanics, the constitutive equations needed must reflect the gross differences in the structural performance. Unknowns Displacement Vector, u Stress Tensor, tr Density, p Balance Equations 1 scalar variable Equations p—-diver = f dt Continuity Equation y Field 3 scalar variables 6 scalar variables 3 scalar equations J —+div(pu) dt =0 H 1 scalar equation Table 2.1 Number offield equations vs. number of unknown variables 7 Chapter 2 Thermodynamics of Continuous Media Despite the fact that atomic concepts in general do not intervene in the formulation of constitutive equations, certain physical and mathematical considerations must be satisfied in the selection of constitutive variables and in the formulation of constitutive equations. Constitutive theory follows along a deductive scheme once a certain number of postulates and basic principles have been adopted. These principles are but the formalization and a priori acceptance of hypotheses which have proven to be efficient and useful both in peculiar applications and in everyday experience. Such principles serve as guidelines and constraints in the construction of the relevant constitutive theory. They thus have a heuristic value. Their foundations cannot be demonstrated and, as such, they cannot be related to any logical frame. Symbolically one can write a(t) = SR e(r) or e(t) = 3 <-r<t (T(T) Tit This equation states that the stress field at instant t depends on the entire history of deformation and vice versa. This functional representation, initially introduced by Volterra in 1912, must however be particularized for each material. We will see in the sequel how certain invariance and thermodynamic principles limit the class of allowable functional. Due to a lack of precise knowledge of bonding etc... at the atomic level, practical identification of constitutive laws at this stage is necessarily experimental. principle of determinism: Past and present causes (independent constitutive determine the present "effects" (dependent constitutive variables), a time functional variables) suggesting the use of for the description of general constitutive equations. However, if no hereditary phenomena are being accounted for, like in this thesis, the principle is usually reduced to the simpler form: present "causes" determine present "effects". principal of contiguity or local action: This principle emphasizes the influence of what is nearby, thus suggesting that only the "causes" defined in an arbitrarily small neighborhood of the point X affects appreciably the dependent constitutive variables at X (so-called axiom of neighborhood). The values taken by independent constitutive variables at distant material points from one material point X do not affect appreciably 8 Chapter 2 Thermodynamics of Continuous Media the values of the dependent variables at X. This leads to the notion of materials of the material-gradient type or theory of n-th gradient of a material. For example, let A and B be an "effect" and a "cause", respectively. Then, accordingly, the constitutive equation x,cr) = A(y B(X,a) R J, with X=X(x, a), corresponds to a first-order gradient theory with respect to the field B. A l l classical theories of continua (classical elasticity, hydrodynamics) are of this type. Higher-order gradient theories are defined by induction, e.g., A = A^7^B(X,cr),V^V^B(X,a)j corresponds to a second-order gradient theory. A counter example of the principle of contiguity is provided by the so-called nonlocal theories (see chapter IV). principle of objectivity: This principle singles out the objective and physical properties of a material from those which are subjective and depend on the observer. The constitutive equations of an ideal "continuous" material must be form-invariant with respect to superimposition of an arbitrary rigid-body motion. In other words, they must be forminvariant with respect to space-time changes of frame. Physically, the principle of objectivity signifies that: "If two observers consider the same motion of a given material the same response of the material, e.g., the same stress and/or body, then they strain determine fields". Although such a principle is applied unconsciously in everyday life, it bears a profound operational significance; the internal forces in a continuum depend on the relative deformations of the continuum with respect to itself and not on the parameters of a rigid body motion. thermodynamical admissibility: The expression of constitutive equations must not contradict the thermodynamical evolution represented by the second principle of thermodynamics. 9 Chapter 2 Thermodynamics of Continuous Media material symmetry: mathematical The constitutive equations of a material must reflect in their very structure the degree of symmetry (isotropy, anisotropy and/or crystallographic group to which it belongs) of the material in its reference configuration or a privileged configuration considered as such. principle of virtual work: Although it is a principle that deals with the concept of energy, the celebrated principle of virtual work, also referred to as d'Alembert's principle has a very intimate relation with the principle of objectivity. Let us consider the case of a continuum subjected to surface tractions T on the external bounding surface dD and to body forces F throughout the region D occupied by matter. The principle of virtual work can be stated as follows SW.+SW.+5W, =0 where SW/is the rate of work performed by the internal forces, SW the rate of work performed by the surface tractions T on dD and body forces F in D, and SWj the rate of work performed by inertial forces. e with SW =- SW t def = \a:SedV SW = ^T.Su*dS e + = [ (a.n).Su * dS - jdiva-Su * dV D j pF.Su*dV D SW =-\ pF.Su*dV j D Hence one can write £ [diva + p (F - y)\. Su * dV + | [T - <j.n]. du * ds = 0 d using Haar's lemma in variational calculus, localization of the above equation for an arbitrary virtual displacement field Su*, yields both the local field equations diva + p(F - -y) the associated boundary conditions T-a.n = 0 10 =0 Chapter 2 Thermodynamics of Continuous Media It is interesting to notice that the virtual displacement field does not have to respect the physical constraints on the real body. For instance, if a degree of freedom is blocked, a virtual displacement which does follow it will yield the corresponding reaction at that point. THERMODYNAMICS OF CONTINUA Before stating the principles of thermodynamics in the form they take for continuous media, we introduce or recall some fundamental notions. We call system, S, a part of the material universe (an open region of R ). The 3 complement of S in R , S , is the exterior of S. 3 ex A system is said to be closed if it does not exchange matter with the exterior. We call thermodynamical system a closed system of which the energy exchanged with its exterior consists only in heat exchange and works performed by body and surface forces acting on S. A thermodynamical system that does not exchange any energy with its exterior is said to be isolated. T H E FIRST PRINCIPLE O F T H E R M O D Y N A M I C S This principle expresses the conservation of energy and could be written as —(U + K)=W dt +Q (2.1) where U is the internal energy of the system under consideration, K its kinetic energy, W the work rate done by external loads, and Q is the rate of thermal energy supplied to the system. These quantities are in turn given by, 11 Chapter 2 Thermodynamics of Continuous Media U= \ pedQ D K=] ^pv'dfi b (2.2) Ivy = [ fv dn+ [rvdr Q= [rd£2- [ qndr JD JdD where p is the mass density, e the internal energy density, v the velocity field, / the body forces, T surface forces applied on the bounding surface 3D of the body, r the non-mechanical heat source per unit volume, q the heat flux supplied through 3D, and n the external normal to 3D .By recalling that T =a n (2.3) and by using Stokes' formula, it is possible to derive the local form of equation (1) p— + div(q -arv) = f v + r dt (2.4) where E represents the total energy density and is given by 1 2 E = e + —v 2 (2.5) THE SECOND PRINCIPLE OF THERMODYNAMICS The first principle introduces the notions of internal energy and heat, whereas the second principle introduces the one of entropy s . More specifically, we assume that the internal energy e is function of the entropy s , the total strain e and the internal variables a (i = \,...,n) i (scalars or tensors) which represent the internal state of matter (crack or dislocation density, residual stresses , etc.). Therefore we have, 12 Chapter 2 Thermodynamics of Continuous Media de = T ds + — x de-—F P P da. (2.6) ' with ds de (2.7) de K=-P de da. where T is interpreted as the absolute temperature of the medium, x represents the reversible stress, and F , the thermodynamic force associated to the internal variable a a. t E q . (2.7) represents the state laws. Introduction of the free energy i// = e-Ts (2.8) allows one to put E q . (2.7) under the following form dy/ s=— dT dy/ \X = P de (2.9) dy/ da. Eqs. (2.5) and (2.6) lead to ds _ p dE v dv x de F, da dt TdT T dt T dt T dt a t (2.10) Combining Eqs. (2.10) and (2.2) and the equilibrium equations given by p— -diva • dt = f J one can deduce 13 (2.11) Chapter 2 Thermodynamics of Continuous Media ds (o ~\ r dt \l ) 1 with 1 1 1 dt where cr* is the entropy production rate, ^ - i s interpreted as the entropy flux, and V represents the gradient symbol. The second principle states that the production of entropy is always positive a*>0 (2.14) Assuming that the absolute temperature is always positive (T > 0 ) , E q . (2.14) can be put under the form d +d >0 i (2.15) T with d = (cr - %): e(v) + F t d = -j-.VT T a (intrinsic dissipation) dt = -q. V(LogT) (thermal dissipation) Traditionally, it is assumed that both the intrinsic dissipation d t dissipation d T and the thermal are positive, not only their sum. A s a general rule the second principle is not sufficient for the complete construction of constitutive laws of materials, since it can only allow the determination of the admissible ones. T o obtain these latter ones, one needs to define the complimentary laws. 14 Chapter 2 Thermodynamics of Continuous Media COMPLIMENTARY LAWS In the sequel we consider only the case of isothermal evolutions (T = Cste). The second principle is then given by d, =(<T- ): % dm s(y) + F —±>0 ' dt (2.16) a The complimentary laws allow the definition of a relationship between the variable dot a = cr - x ir a n d the thermodynamic forces F to the velocities e(v) and — , L a dt respectively. To construct such laws, the normality rule is usually assumed to hold. This rule states that there exists a positive function (p(cr ,F ) ir a assuming its minimum at zero such that (i(u),d) e d(p(o~i , F ) r a (2.17) where u is the displacement and d is the sub-gradient symbol. This assumption, however, must be verified experimentally for each material. According to HalphenNguyen [68], a material which obeys the normality rule and has a free energy \\i and a function <p that are convex of their arguments is said to be a generalized standard material, (p is then called dissipation pseudo-potential. Application of the principles of the thermodynamics of continuous media are given in the next two sections in the context of continuum damage mechanics and fracture mechanics. A P P L I C A T I O N I: Continuum Damage Mechanics Important efforts have been made over the last few decades in developing the area of thermodynamics by establishing the necessary bases required in the construction of physically admissible constitutive laws. We adopt in this section the theory of internal variables. Therefore, the state of a continuum at instant t is fully determined by the knowledge of the value of the observable variables together with the values taken by a set 15 Chapter 2 Thermodynamics of Continuous Media of hidden variables which are assumed to characterize the internal state of matter. The reader interested in a detailed presentation of the content of this section should consult references [41,48,101,103,117] W e limit the presentation in this section to the case of a damage-elasticity coupled behavior, with no permanent deformation (e = e ), e in the framework of isothermal transformations (T = Cste ). The usual classification of the different variables in this case is then as follows Observable variables Hidden Variables Associated variables e CT T s D Y Table 2.2. Classification of the different field variables where s is the entropy of the system and Y the thermodynamic variable associated with D. Given the expression of the specific free energy y/ = y/(e,D,T) (2.18) The state laws are derived as follows a = p s= Y=p dy/ de dy/ (2.19) dy/ 3D Since we are dealing with isothermal evolutions, the thermal dissipation vanishes and the intrinsic dissipation becomes, <j), = -Y : D > 0 16 (2.20) Chapter 2 Thermodynamics of Continuous Media with Y :D = Tr[Y.D] the Clausius-Duhem inequality having its roots i n the second principle of thermodynamics imposes -Y The complimentary evolution :D>0 laws (2.21) which govern the internal variables, are thermodynamically acceptable i f the intrinsic dissipation is always positive or equal to zero. The introduction of a dissipation pseudo-potential cp leads to the adoption of the generalized normality rule and the respect of the positive aspect of the dissipation. Let q> * be a function of the associated variables Y, dual of cp , obtained through a legendre Frenchel Transformation. In the case of generalized standard materials [68], we can state (p(<j,Y;e,b) = where A D is the damage factor, and / X f((T,Y;e,b) (2.22) D is a convex function, centered at the origin representing the damage threshold. This results in an evolution law of the form with or and = 0 and 17 f <0 Chapter 2 Thermodynamics of Continuous Media Two elements govern the implementation, or particularization to a given material, of the constitutive laws • knowledge of the damage threshold function • Knowledge of the expression of the specific free energy The expression of specific free energy y/ is chosen so that it is compatible with the expression of the stress tensor tr through the constitutive law. Therefore, in the case of an elastic-damageable material we have a = A(D): ( s v} i r d = p (2.24) If we choose y/ to be given by py/ = ^A{B):e\e (2.25) then, the thermodynamic variable associated with D can be put under the form Y=P dyi dD IdAiPf = 2 dD . £: £ (2.26) ' Cordebois [48] proposed a damage threshold of the form f(Y) where Y = yjTrY 2 u = Y -K(D) u =0 (2.27) and K(D) represents the damage threshold for a given value of D. If we assume that K(D) = K(D ) (2.28) D.^VTTD (2.29) n with 7 18 Chapter 2 Thermodynamics of Continuous Media Then A,=™ (2.30, A. Consequently, the damage evolution law (2.23) becomes • df Y D = -X ^= -X — dY Y (2.31) Du D u or JTW =A D ^ ^ =Z D (2.32) If we express the fact that a loading case, other than fatigue, corresponds to / = 0 and / = 0 we deduce from E q . (2.27) Y --^-D =0 n (2.33) n or ^x mzn in B dD =0 (2 .34) D Y„ n U Eqs. (2.32) and (2.33) lead to a new expression of the damage factor 4 = d d *"/ K D K n — NlW (2-35) or Y- i iY- uDn _ t t z v (2.36) -Tr[D.Y] dD dK n 19 Chapter 2 Thermodynamics of Continuous Media Further particularization is achieved by the introduction of a damage tensor of the same type as the one proposed by Chaboche [41], and which is governed by a differential equation of the type b = Q(a)D (2.37) where Q(<r) is a second order tensor which depends on the material and the state of stress. In the particular case of a proportional equivalent to a 'fixed' loading this simplifying assumption is anisotropic situation irrespective of the loading (crstaying proportional to a fixed stress cr ). Under these conditions, we have 0 D = Q((T)D (2.38) Noting that VrrZ> = D ^TrQ 2 (2.39) 2 and D = yJTrt) 2 n = D ^TrQ (2.40) 2 dK Let = Mix) be the so-called 'modulus of deconsolidation' 'ii dD, M plays the role of the hardening parameter in the theory of plasticity and is determined experimentally. E q . (2.39) leads to a new expression of the damage factor given by (2.36) 4 = ^ 7 M (2.41) 20 Chapter 2 Thermodynamics of Continuous Media and the evolution law becomes D = Y Y UMY (2.42) U The problem is completely determined, since we know the expression of Y , which is given by Y = ^ : s : s 2 8D K (2.43) J It is worth noticing that in the particular case treated here the intrinsic dissipation is given (p,=-Y Since <j = A:e :Q(<r)b it can be shown [41] that -Y (2.44) :Q(<r) is a quadratic form of the srains. Under these conditions the Clsusius-Duhem inequality (2.21) imposes D>0 Which (2.45) states that damage is an irreversible process and corresponds to a real deterioration phenomenon. E X A M P L E : Brittle Elastic Materials The free energy for these materials is usually given under this form py/ = ^a{D):s:s (2.46) where D is a tensorial or scalar damage parameter representing the state of deterioration or internal damage of the material (cavities, or cracks). The dissipation pseudo-potential is assumed to correspond to the "indicatrice " of the convex C defined by 21 Chapter 2 Thermodynamics of Continuous Media C = {F , D where F D f(F )<\} (2.47) D is the thermodynamic force associated with D. The complimentary and constitutive laws lead to j-i F = 1 da n D 2dD Z = D =A df dF :£•£ a(D).e (2.48) with n \A = 0 if f<l U>0 if f =l When damage is isotropic, the stiffness matrix of the material a(D) is expressed as a function of the initial stiffness matrix a(0) in the following way a(£>) = (l-Z>) a(0) (2.49) This equation is equivalent to the effective stress hypothesis which states that the constitutive law of a damaged medium can be obtained from the constitutive law of the plain material simply by replacing the real stress a by the effective stress a eJf defined by (2.50) It is seen that for these materials, the free energy is no longer convex with respect to all the variables considered together, but is only convex with respect to each variable taken separately. This model applicable only in the case of fragile materials such as plain concrete, is not of great interest as far as F R C materials are concerned. Another model, better suited for this kind of materials, w i l l be presented in more details and applied to predict the response of F R C specimens in chapter 4. 22 Chapter 2 Thermodynamics of Continuous Media A P P L I C A T I O N II: Fracture Mechanics In this sequel section, the general equations of the Griffith theory w i l l be deduced from the coming derivations of the energy balance for a cracked body and w i l l be applied to two examples. W e start with the global energy balance in a continuum during crack growth from which the Griffith criterion is deduced as a special case. In an attempt to extend the principles of linear elastic analysis to situations of highly localized yielding at the crack front, the various irreversibilities associated with fracture are lumped together to define the fracture toughness of the material. This approach allows the applicability of Griffith's theory to metals and other inelastic engineering materials. A general energy balance of a deformable continuum subjected to arbitrary loading and containing a crack is presented. N o particular assumptions regarding the constitutive equations relating stresses and strains are made. Consider the case of a cracked continuum subjected to surface tractions T external bounding surface 8D T on the and to body forces F throughout the region D occupied by matter, as we can see from Figure 4.1. Furthermore, we assume 1. the total bounding surface of the body 3D to be given by dD(t) = dD + dD (t) T c where dD (t) represents the surface of the crack c 2. the cracked surfaces to be stress free 3. the total volume to be unaffected by crack growth 4. The crack is not necessarily stationary but may be propagating and it is assumed that the crack growth is described by the crack area as a simple parameter. 23 Chapter 2 Thermodynamics of Continuous Media Figure 2.1 A continuum body with a crack Application of the first law of thermodynamics to this cracked body leads to W +Q =E +K +t (2.51) Here, W is the work performed per unit time by the surface tractions T on dD T and body forces F in D, Q is the thermal energy supplied to the body per unit time, E is the rate of change of the internal energy, K is the kinetic energy of the continuum, and f the energy per unit time spent in increasing the crack surface dD . c quantities are in turn given by 24 These different Chapter 2 Thermodynamics of Continuous Media ^=[ T D T k^dZ+lpF u k k dV Q=[ qn dZ+[ pkdV D k ) E = — \pedV dt *> = [pedV (2.52) where, p is the mass density, u is the displacement vector, n the unit vector normal to the bounding surface, e the internal energy density per unit mass, q the heat flux vector per unit surface, h the non-mechanical heat source per unit mass, and y the energy required to form a unit of new surface. For the scope of the present study, we assume that a- the applied loads are time independent b- the crack grows slowly in a stable manner c- there is no heat flux into the body d- absence of any non-mechanical heat source in the body Under conditions a and b, the velocity field ii developed in the continuum is small and hence the kinetic energy which is proportional to u 2 may be ignored. In other words, the realm of dynamic fracture mechanics where the kinetic energy term, K, cannot be omitted is not addressed in this study. According to conditions c and d, crack growth takes place under isothermal or adiabatic conditions. Thus, the thermal energy term in equation (2.52) Q can be ignored. In E q . (2.52), the rate of change of internal energy can be rewritten as, (2.53) 25 Chapter 2 Thermodynamics of Continuous Media where cr. and s tj denote the components of the stress tensor and strain tensor respectively. The strain increment de tj may be split into a recoverable elastic part ds.j and a permanent inelastic part ds™ , with de = ds.j + de\" Hence, E q . (2.53) takes the {j form E=lcT ^dV l cr e-dV u where D' n + n (2.54) u denotes the part of the body that has been undergoing an inelastic deformation. The elastic strain energy rate may be given by I)' = f (2.55) £ . dV e CT.. whereas, the rate of inelastic strain work is given by . y« i] (2.56) y According to the above mentioned assumptions, all changes with respect to time are solely caused by changes in crack size. Let A be the crack surface, therefore dt ~ dA 8t~ W = ™A dA with dt = [ T ^dZ+\ F -^dV d JBD k t k pi A y d£ e E=—(u +u )A e dA X = in ' A >0 f u k P k QA , ds" I cr -±dV+ f o-,—?-dV + b '} b'" ii 26 'J P)A (2.57) (2.58) (2.59) Chapter 2 Thermodynamics of Continuous Media Consequently, the energy balance equation E q . (2.51) during quasi-static stable crack growth, where dynamic effects and thermal energy are ignored, takes the form 8W dU 8A dU' • + dA dA e ar + •dA (2.60) or a n _ du ar m ~dA~ dA (2.61) ~dA where i l = U -W is the potential energy of the system. E q . (2.60) shows that the work e rate supplied to the body by the external forces is equal the rate of the elastic strain energy and the inelastic strain work plus the energy dissipated during crack growth. The other form of the energy balance, E q . (2.61), indicates that the rate of decrease in potential energy during crack propagation is equal to the rate of energy dissipated in inelastic deformation and crack growth. Assuming that the energy required to form a unit of new material surface is constant for a given material and environmental conditions, Griffith [65,66] used the minimum potential theorem of elastostatics to approach the problem of fracture of an ideally brittle material. If y represents the energy required to form a unit of new surfaces in such materials where the inelastic energy dissipation is usually negligible, E q . (2.60) takes the form dA dA representing the energy available for crack growth, is given the symbol G in honor of Griffith. 2y represents the resistance the material opposes to crack propagation and is a material constant. Therefore, E q . (2.62) represents the fracture criterion for crack growth. There are two important cases that are usually encountered in the literature. One 27 Chapter 2 Thermodynamics of Continuous Media called 'fixed grips' loading and corresponds to the particular case where the surface of the body on which the loads are applied remains stationary when the crack grows. The other, called 'fixed load' is characterized by the fact that the applied loads on the surface of the continuum are kept constant during crack growth. It is interesting to note that for both cases, the magnitude of the elastic strain energy release rate necessary for crack growth is the same. However, the elastic strain energy of the system decreases for 'fixed grips' and increases for 'fixed load' conditions. In either case, E q . (2.62) can be written as G =- ^ = 2 y dA (2.63) —(n+r)=o (2.64) or dA where IT is the potential energy defined by E q . (2.61). This form of the energy balance states that, the total potential energy of the system IT + T has a stationary value. Example 1 (Line crack in an infinite plate) The change in elastic strain due to the presence of a line crack of length 2a in an infinite plate of unit thickness subjected to a uniform stress ka- perpendicular, and a parallel to the crack as calculated by Sih and Liebowitz [153], is given by 2 2 U =— 8/i (*: + !) 28 (2.65) Chapter 2 Thermodynamics of Continuous Media Figure 2.2 line crack in an infinite plate The critical stress required for unstable crack growth is given by E q . (2.62) 2Ey (2.66) 7ta For generalized plane stress 2Ey (2.67) na(\-v-) } For plane strain A s expected, it is seen that the stress ka parallel to the crack plane has no effect on the critical fracture load. Furthermore, Eqs. (2.66) and (2.67) Confirm Griffith's experimental findings, indicating that the critical stress cr. is proportional to the square root of the crack length. 29 Chapter 2 Thermodynamics of Continuous Media Example 2 (Elliptical crack embedded in an infinite solid) a/b Figure 2.3 Normalized critical stress versus the ratio of the semi-axes of an elliptical crack [153] Let us now consider the case of an infinite plate endowed with an elliptical crack of semi axes a and b and subjected to a uniform stress cr perpendicular and sa parallel to the crack plane. The change in elastic strain energy U' and the increase in surface energy T due to the presence of the crack are given by [95], u t 2x(\-v)abW = 3pE(k) Y = 2nahy (2.69) Adopting the simplifying assumption of Kassir and Sih [95] stating that the ellipse of semi axes aandb grows into another ellipse with the same foci, application of the Griffith criterion becomes possible and leads to 30 Chapter 2 Thermodynamics of Continuous Media — {U -r)=0 e (2.70) A s a result, the following expression for the critical stress is derived b(l-v) 3/M (Hk' ) -<J. — T. 2(\ + k )E(k)-k 2 2 u u E\k) » K(k) (2.71) A graphical representation of E q . (2.71) is given by Figure 4.3. It is seen that the critical a —>1 stress increases as the ratio — b shaped crack is decreases. Its maximum value for the penny- °c =— a(l-v) 31 (2-72) Chapter 3 Homogenization of Heterogeneous Media HOMOGENIZATION OF HETEROGENOUS MEDIA INTRODUCTION The behavior of fiber reinforced concrete exhibits several characteristics that could be interpreted in terms of the presence of microcracks and fibers. The aim of this chapter is to describe few techniques on how to design a mathematical model of the behavior of concrete that takes into account those characteristic properties of the material. A natural way of describing the response of F R C seems to lie in the idea of considering the macroscopic behavior of matter as a consequence of its microscopic behavior (i.e. at the scale of heterogeneities). If we admit the validity of such an assumption, it seems to be very natural to approach the behavior of F R C by studying friction at crack lips and at the fiber-matrix interface, propagation of cracks, fiber pull out, etc. To design a model that would behave like F R C we can adopt a macroscopic approach [41,101,103] or a microscopic one [2,29,76,112]. A typical result of the former is to postulate the existence of an a priori parameter along with its influence on the macroscopic mechanical properties, and calibrate its evolution from experimental results. In the latter case, the method is less direct, but provides a direct link between the damage parameter and the microstructure. Although homogenization methods are well established for the elastic case, these are in their infancy when local dissipative phenomena such as crack propagation and fiber pull-out need to be included. The aim of homogenization techniques is to replace a highly heterogeneous material by an equivalent homogeneous one whose behavior approximates the macroscopic behavior of the heterogeneous composite. T w o essential reasons explain this approach: • Some characteristic aspects of the mechanical behavior of highly materials can be explained only at the level of the microstructure. 32 heterogeneous Modeling of such Chapter 3 Homogenization of Heterogeneous Media behaviors takes into account the mechanisms taking place at the level of the heterogeneity. • From the numerical point of view, the cost of structural analysis when these materials are used would be extremely high if all the heterogeneities account have to be taken into explicitly. These two reasons define the scale at which we should put ourselves in order to study the behavior of the equivalent homogeneous material, the so-called representative volume element (rve). The rve should be large enough with respect to the microstructure in order to have global information on the behavior of the equivalent homogeneous material , and small enough with respect to the structure so that the equations that govern its behavior could be written in a local form. The different homogenization methods differ in the sense of the equivalence which allows the transition from the behavior of the heterogeneous material to the behavior of the equivalent homogeneous material. MICROSCOPIC AND MACROSCOPIC SCALES In principle, the equations that describe the behavior of each phase are known and may be written at the microscopic level, where we focus our attention on what happens at a point within a considered constituent present in the domain. Sometimes it is even possible to know the conditions that prevail on the surface that bounds a given phase. However, due to the fact that the geometry of the surface that bounds the phase is not observable and/or is too complex to be described, the equations cannot usually be solved at this level. Since the general description and response of a heterogeneous material to a given stimuli at a microscopic level is impractical and may even be impossible, another level of description, is needed, the so-called macroscopic level. A t this level, measurable, continuous and differentiable quantities may be determined and boundary value problems can be stated and solved. Accordingly, the main objective of this chapter is to introduce 33 Chapter 3 Homogenization of Heterogeneous Media the continuum approach that leads to the macroscopic level of describing the constitutive laws of heterogeneous media. Each phase (gas, liquid or solid) is regarded as a continuum, overlooking its molecular structure. A multi-phase medium is called homogeneous with respect to a macroscopic parameter characterizing the geometrical configuration of any phase within the material, if that parameter has the same value at all points of the domain. If not, the domain is called heterogeneous with respect to that parameter. For example, a porous medium domain, D, is homogeneous with respect to porosity, 0, if inZ) V<9 = 0 (3.1) As a consequence of the microscopic configuration of any phase within the rve, in the neighborhood of a point in a heterogeneous medium domain, certain macroscopic properties at that point may vary with direction. A continuous medium is said to be if that property varies with direction at that point. anisotropic with respect to a property A medium is said to be isotropic at a point with respect to a property, if that property does not vary with direction at that point. The advantages of an overall continuum description of a heterogeneous medium are a- It circumvents the need to specify the exact configuration of the interphase boundaries, b- It describes processes occurring quantities, thus enabling mathematical in heterogeneous media in terms of the solution of problems by employing differentiable the methods of analysis, c- The macroscopic quantities are measurable, field problems of practical and can therefore be useful in solving interest. These advantages are however, at the expense of the loss of detailed information concerning the microscopic configuration of interphase boundaries and the actual variation of quantities within each phase. However, as we shall see later, the macroscopic effects of these factors are still accounted for through some averaging technique. 34 Chapter 3 Homogenization of Heterogeneous Media REPRESENTATIVE V O L U M E E L E M E N T Suppose we wish to describe the motion of the fluid inside the void space of a porous medium. For example, we wish to determine the velocity of the fluids at points within the fluid that occupies the void space. At first, we see no difficulty. To achieve our goal, we have to solve the momentum balance equation, which for a Newtonian Fluid (and water is considered as such fluid) takes the form of the well known Navier-Stokes equation. This is a partial differential equation, written in terms of the state variable "fluid velocity". To solve this equation for the fluid domain, i.e., the void space domain, we need also the boundary conditions. The information on boundary conditions contains two parts: (a) information on the geometry of the boundary surface, and (b) conditions that have to be satisfied on the boundary. In the case considered here, the boundary surface is the solidfluid interface. Suppose we know the condition on the boundary surface, e.g., a condition of "no flow through the boundary". Is there any practical way to describe the detailed geometry of this boundary surface? The answer is, obviously, NO. This procedure for obtaining a description of the flow is not possible. However, do we really need to know what happens at every point within the matrix phase, or within any given inclusion? Furthermore, suppose we could predict the behavior at such points, is there any way to measure this behavior, say, in order to verify our prediction? Can we really measure the velocity or the pressure at a point within a fluid occupying the void space in an FRC material? 35 Chapter 3 Homogenization of Heterogeneous Media To overcome the above difficulties which stem from the heterogeneity of the domain (which is a consequence of the presence of solid and void space subdomains), a smoothing operation is usually performed. By the averaging over an rve, we can transform the microscopic description into a macroscopic one, the latter being a continuum description. We must still select the appropriate size and shape for the rve for a heterogeneous material. To determine this size, let us recall what we shall use the rve for. Once rve-averaged values of some property, say of stress or strain, have been introduced into the macroscopic model, we should construct or use instruments that measure these averaged values, i.e., that take averages of the corresponding microscopic values over the selected size of the averaging volume. If our model really describes the considered phenomenon, the predicted and measured values must be the same, within the range of error introduced by the modeling and measuring processes. For each selected averaging volume, the averaged values of the state variables will, in general, be different (see Figure 3.1). Hence, every averaged value must be accompanied by a label that specifies the volume over which this average is taken. But, How could averaged values be compared? Instead of employing an rve of an arbitrary size, the size of an rve should be selected such that the average value of any geometrical characteristic of the microstructure of the region occupied by an inclusion at any point in a heterogeneous medium, is a singlevalued function (or nearly so within an acceptable error) of the location of that point only. The averaged value should also be independent of small perturbations in the size of the rve, which means that it should remain more or less constant over a range of volumes that correspond to the range of variation in instrument sizes. Actually, this process of averaging quantities over a certain volume element should not be new to us, and is very well illustrated by the definition of the density of a continuum. Because of the molecular structure of matter, and its discontinuous structure at that level, 36 Chapter 3 Homogenization of Heterogeneous Media if "point" means "a mathematical point", the answer to the question of density is not unique. The point may "hit" a molecule, or be in vacuum. Under such conditions, there is no meaning to the term "density". T o overcome this difficulty, let us enclose the considered point within a small spherical volume V centered at the latter, and then try to determine the ratio MTV within that sphere, where M is the mass of matter enclosed in the sphere. If the volume, V is very small, the result w i l l be unstable and strongly fluctuating. A s we increase the volume, and the number of molecules increases, the fluctuations in the ratio of M / V w i l l gradually attenuate, reaching some plateau, with ever decreasing fluctuations. This process is shown on Figure 3.1. In this Figure, we note that the rve, may be selected in the range between V = U point, is the ratio MN for V = m and V = . The density, p, at the . Density Figure 3.1. Variation of density with the size of the rve 37 Chapter 3 Homogenization of Heterogeneous Media By using this procedure of averaging to determine p at every mathematical point within a domain, the domain and its behavior have been transformed from the molecular level to the microscopic one. The latter is now a continuum. The density is defined at every microscopic point. The real molecular world is transformed into the fictitious continuum world. Denoting the characteristic dimension of an rve by / and the length characterizing the microscopic structure of inclusions by d, a necessary condition for obtaining non-random estimates of the geometrical characteristics of the heterogeneity at any point within a polyphase medium domain emerges, / >> d. (3.2) In other words, we require that the size of the rve, say the diameter of a sphere, be much larger, say 100 times, the scale of heterogeneity at the microscopic level, indicated by the dimension of an inclusion. This is analogous to the requirement, in passing from the molecular level to the microscopic one, that the size of the rve be sufficiently large, so as to contain a very large number of molecules. Here we require a very large number of grains, pores and fibers, so that if the rve contains a few inclusions more, or a few inclusions less, the average will remain practically unchanged. Thus, the size of the rve must be much larger than the scale of the heterogeneity at the microscopic level, resulting from the presence of solids and voids. We have thus determined the lower bound to the size of the rve. However, as can be seen on the previous Figure, if the domain is heterogeneous at the macroscopic level, an upper bound has also to be determined, as we wish to select an rve that will reflect the properties of the porous medium domain at the considered point. Therefore, another condition, that sets an upper limit to the size of the rve, is required. The selection of the size of the rve is also constrained by the requirement that l « L . (3.3) where L is a characteristic length of the heterogeneous medium domain, such as the size of the structure being analyzed. 38 * Chapter 3 Homogenization of Heterogeneous Media DESCRIPTION OF FEW HOMOGENIZATION TECHNIQUES MACROSCALE CONTINUUM MICROSCALE Figure.3.2. Macroscale vs microscale representation of an rve For each type of medium, knowledge of the most important microscopic mechanisms is necessary for the modeling of materials with microdefects, because a number of characteristic properties of such materials cannot be explained at any scale other than the scale of the microstructure. But at the same time, from an engineering point of view, precise knowledge of the continuum fields such as displacement, stress, dissipation, discontinuities, ...etc., which in the case of a finite element analysis for instance, would require a refined mesh around each type of discontinuity, and for the whole solid, is not of significant interest. Therefore, the scale at which we are interested is located at an intermediate level between the scale of the local mechanisms and the one of global mechanisms of the considered structure. A s we saw earlier, the macroscopic scale is defined as the scale of the material point of the solid and is usually represented by the 39 Chapter 3 Homogenization of Heterogeneous Media rve. This volume is denoted by Q , X is the macroscopic space variable, and x is the microscopic space variable. The microscopic scale on the other hand is the local scale within the rve Q . It is at this level that the displacement, stress, and strain fields are defined. The problem of homogenization is related to the question of whether it is possible to model the global behavior of a heterogeneous material, without explicitly using what happens at the microscopic level, and i f that is possible how could it be done?. Hence we should first of all define the macroscopic strains E(X) and stresses E(X) as functions of the local fields u(x), e(x), tr(x); and then determine the macroscopic constitutive law relating S(X) to E(X). v.\V> >.\V> £2 dn S-\\\\\\\\\\S\\\\\V \\\\\\\\\\\V\V\.NSV \\S\S\S\\\\V\.\'1\\VS\\\N\\\\\V \\\\\\\\\\\\\\\V. .. .wwwwv Figure 3.3. Example of an rve containing a cavity-type inclusion If the rve contains void type defects such as cavities and microcracks we denote by Q the geometric volume occupied by matter, T the volume of cavities and by dD. the boundary of Q . Therefore a = a - r dh = dQvjdr ; df2ndr = 0 (3.4) The strain tensor in small perturbations is defined as a function of the d i s p l a c e m e n t « by: u £ =-(d u +d u ) i 40 j j i (3.5) Chapter 3 Homogenization of Heterogeneous Media The macroscopic strain and stress tensors E and 27 which are intuitively related to averages of the local fields e, u and o~ permit the computation of the macroscopic strain energy of the solid S. Since our objective is to hide the microstructure of Q when describing the behavior of S, it is natural to impose the following condition on the fields E and Z : (3.6) (I) (II) where v is the outer normal to dCl and n is the outer normal to dT. The preceding equation consists of 2 parts; the average strain energy in the element Q , (I), and the work of the forces applied on the internal boundary of Q , (II). In the case where Q does not contain any cavity-type defect (voids, cracks ...etc.), the above mentioned strain energy reduces to the average of the strain energy in the rve. O n the other hand when such defects exist, the strain energy has an extra term due to the work of the forces applied on the boundary of such defects (contact between the lips of a crack, pressure of a fluid on the surface of cavities, etc.). Taking this work in the strain energy expression shows that during the process of homogenization these forces w i l l be considered as internal forces in the equivalent material. To establish homogenized constitutive laws, different approaches have been used since Poisson [138]. A m o n g the most popular of these techniques we have, the variational methods, theory of effective moduli, periodic homogenization, and the self consistent method. In the elastic case, Kroner [98] has proposed a statistical approach that gives a clear distinction between the periodic homogenization method and the self consistent method by clarifying the relative assumptions from which they are derived; that is perfect order for periodic homogenization which utilizes a periodic microstructure, and complete disorder for the self consistent method. 41 Chapter 3 Homogenization of Heterogeneous Media In general, it is not possible to determine the fields 2 and E satisfying equation (3.6) without an extra assumption. The assumption is either related to the microstructure of the material the solid is made of, or to the form of the local fields <y,e oru. The various homogenization techniques are different from each other because of the different assumptions they adopt and the definition of the macroscopic fields that follow. VARIATIONAL METHODS Initially introduced by Hashin and Shtrikman [71], H i l l [75], this method was reconsidered by Bergman and Kantor [31] and in a pure mathematical form by Tartar [163], Francfort and Murat [60], Golden and Papanicolaou [67], and others. In the case of linear problems, these methods yield bounds on the coefficients of the equivalent homogeneous material. The utilization of this method is possible without any prior knowledge of the geometry of the heterogeneities, but the obtained results get better as a better knowledge of the flaws geometry is taken into account. This approach is based on the use of energy theorems along with very suitable choices of the stress and strain fields which have to be statically and kinematically admissible. METHOD OF EFFECTIVE MODULI The key assumption in this method (approach) is related to the distribution on the boundary dQ. of either of the following fields: . the stress vector (stress approach) or, . the displacement field (displacement approach). To fix ideas, we consider the case of a homogeneous, linearly elastic material endowed with cavities and cracks. A s long as the cavities do not propagate and the cracks do not close back, the behavior of the altered material remains elastic. The proposed approach is based on the determination of the compliance or rigidity tensors of the equivalent material. The equivalence between the two materials is then 42 Chapter3 Homogenization of Heterogeneous Media \ 0_ \\\ s = > defined by equating their elastic potentials. The stress approach allows the determination of the effective compliance tensor S, whereas the displacement approach allows the determination of the rigidity tensor C. Z = C:E E = S:I (3.7) Stress Approach Hypotheses - The stress vector is constant on any face ofd - Stresses on opposite sides are of opposite sign Equilibrium of the rve Q shows the existence of a symmetric tensor S(X), satisfying the following equation on dQ: Mxedfi; E(X)\s CT(JC) cr(x)v = Z(X).v (3.8) the macroscopic stress tensor associated with the microscopic stress tensor . It could be shown that Z(X) Z{X) is the average (or) of <J(X) over Q, given by: = (<r) (3.9) =\kio a(x)dv The macroscopic strain E(X) associated with the field u(x) is then deduced from equation (3.6) which represents the strain energy. And the final result could be put under the form: where (u ® v ) = — {u v s i 43 j + u.v ) Vs (3.10) Chapter 3 Homogenization of Heterogeneous Media After imposing the macrostress 27 , the macrostrain E is expressed as a function of 27. The equations of the problem are as follow: div(cr) = 0 in Q a.v = 27.v on d£2 an on e(u) =0 -s:a (3.11) dr s : local compliance tensor Due to its linearity the problem can be reduced to the resolution of 6 elementary problems. The final result is typically a compliance tensor that looks like the following expression: S =T-x\-S <rl <Ti dv i 2 *" PI M mma t (3.12) l where cT*' is the stress localization tensor. For other stress distributions on the boundary, n different compliance tensors are obtained. Displacement approach Basic hypothesis There exists a symmetric tensor E(X) such that the displacement field satisfies the following equation on the faces ofdQ. Vxedn, u {x) = E (X).x i ij The relationship between the macroscopic strain E(X) is the same as before and the macroscopic stress E(X) stress a(x) j (3.13) and the local displacement u(x) is related to the microscopic by: (3.14) Z..= a .X v ds ik 44 r k Chapter 3 Homogenization of Heterogeneous Media After imposing a uniform macroscopic strain E by prescribing a linearly varying displacement field on 8Q, the rigidity tensor is deduced by expressing S as a function of E. The equations of the problem are summarized by: in Q div(cr) = 0 u,. (x) = E u on dO .Xj (3.15) cr.n =0 on dr cr = c : e(u) c : local rigidity tensor The procedures for determining C are basically similar to those used to determine 5, and the final results are analogous, like we can see on the following expression: 'pi • c „ „ / ™ (v m )s „ (v' ) J n dv (3.16) A fundamental limitation of this method is the fact that for the same geometry of pores, the tensors 5 and C are not inverse of each other. SELF-CONSISTENT M E T H O D The preceding approximation is no longer realistic when the inclusion's concentration (porosity in this case) becomes important in which case the interaction between the defects is important. In order to include this interaction in the compliance tensor S for an arbitrary concentration of the defects the self-consistent method can be used. This technique takes into account in some average sense the interactions among all the defects present in the rve. Just like in the case of the method of effective moduli, we distinguish a stress approach and a displacement approach, and equations 3.8-3.10 for the stress approach and 3.133.14 for the displacement approach respectively are still valid. This method was originally developed for polycrystalline metals. The microstructure of these materials consists of grains having different shapes, crystalline orientation and 45 Chapter 3 Homogenization of Heterogeneous Media mechanical characteristics. Hence it is impossible to isolate each and every grain and consider one of the phases forming the grain as a matrix containing all the other phases. To get back to the original matrix-inclusion problem, we adopt as a matrix the equivalent homogenized medium we are looking for. Ignoring local effects, we can imagine that each grain "sees" the whole material to which it belongs as a homogeneous material (the equivalent material we are trying to determine). The other assumptions on which the method is based are basically the same as for the method of effective moduli. Furthermore we assume that the inclusions can be grouped into n "phases", where each phase i is characterized by its compliance tensors,.. Let | Ti| = u l r " ! be the total volume of inclusions in the phase i, and be the average of the q u a n t i t y / i n the phase i: l<,i<.n (/> = j4Z \fds (3-17) If one of the phases contains cavities or cracks the foregoing equations are still valid with the following definitions of averages: (3.18) Let n denote the index of the phase consisting of voids (if it exists). Under the assumption that each phase is homogeneous, (e) is expressed as a function of (a) 46 by the following Chapter 3 Homogenization of Heterogeneous Media relation (e) =s,.(o") (for n^i). This way the macroscopic strain E could be expressed by the following relations: \E E *,+ £ ( ' , • " * / ) A'(S) + A"(S) :2- (3.19) E S s +^ -s )A'(S) 1 i + I A (S) n Consequently, the determination of S reduces to the determination of A . In the classical 1 self-consistent method, these tensors are computed for inclusions with simple geometries such as ellipsoids under the assumption of an infinite medium. This way we have to solve as many matrix-inclusion problems as the number of inclusions in the rve. The main difficulty of this method resides in the fact that the compliance tensor of the matrix used in the elementary calculations is the one we are trying to determine. This leads to the following implicit equation: (3.20) Contrary to the method of effective moduli, the self-consistent method overestimates the interactions between the heterogeneities. This could lead to unexpected results. To overcome this difficulty and take into account local effects in a realistic way, other techniques have been introduced to the self-consistent technique. One of these techniques is the so-called multiple scale homogenization [170]. The basic idea of this technique is that an inclusion of a given size d sees in a homogenized form only the medium made of the matrix and the inclusions of sizes smaller than d, the larger inclusions being ignored. The self-consistent model has been extended by H i l l [77] to nonlinear behaviors. This extension has been applied to elastic-plastic and viscoplastic polycristals, and the results were found to be in agreement with the Voight-type upper bound. A more restrictive upper bound has been established recently by Ponte Castaneda [139]. 47 Chapter 3 Homogenization of Heterogeneous Media PERIODIC HOMOGENIZATION This method is based on the assumption that the microdefects or inclusions are distributed in a periodic fashion throughout the solid under consideration. Hence, the microscopic stress and strain fields are periodic. It has been shown that in this case the microscopic strain field can be written as u (x) = E (X)x +v (x) i ij j i ( 3 21) E : a symmetric second order tensor ij Vj: a periodic function In the case where Q , does not contain any cavity-type defect or crack, the tensor E(X) is interpreted as the average strain of the ce\\E(X) = {e(u)). Hence, E(X) is the macroscopic strain we are looking for. In the case of porous or cracked media, the relation between E(X) and e(u) is given by: (3.22) The average stress tensor obtained by the eventual prolongation of the stress field a by 0 over the cavities, is given by the following expression. (3.23) The tensor (a) satisfies the following equation relative to the strain energy: ( H i l l ' s lemma) 48 (3.24) Chapter 3 Homogenization of Heterogeneous Media Consequently the macroscopic variables are defined in a consistent way with Equation (3.6) by: (3.25) The assumption of periodicity provides a natural way for defining the macroscopic mechanical properties and allows one to reduce the problem of studying the whole solid to the study Of the unit cell. Alike the case of effective moduli there is a displacement and a stress approach: Displacement Approach Starting from equation (3.21): u (x) = E (X)x +v (x) i ij j i we impose Ey which represents the linear part of the displacement, the unknown being the periodic displacement field, v ., solution of the following problem: ( div(a) =0 in Q an-0 on dr a = c:e(Ex tr.v v + v) (3.26) antiperiodic periodic The field r . n is said to be antiperiodic if it takes opposite values on opposite sides of the cell. This condition is due to the continuity of tr.v through dQ. In the case of elastic inclusions or cavities, the existence of the rigidity tensor has been proven when the size of the period tends to zero. In the presence of cracks, however, the convergence problem is still open (Sanchez-Palencia [145]). Nevertheless, we can 49 Chapter 3 Homogenization of Heterogeneous Media compute the equivalent rigidity tensor using averages of the elementary solutions on the unit cell analogous to the ones used in the method of effective moduli. Stress Approach The average stress is imposed on the cell. W e obtain the macroscopic compliance tensor using a procedure similar to the one used i n the displacement approach. It can be shown that the two tensors C and S obtained through these two approaches are inverse of each other. The major difference between the method of effective moduli and the periodic homogenization technique is related (due) to the difference i n the boundary conditions imposed to the unit cell. IV- COMPARATIVE EXAMPLE: Case of an elastic porous medium The differences among the three methods presented earlier are illustrated here through the determination of effective moduli i n the case of a porous medium. Since this chapter is restricted to techniques providing closed form solutions to the problem, the assumption of an infinite medium is practically unavoidable. F o r more realistic geometries, numerical methods are usually invoked (see chapter 8). W e consider the case of a linearly elastic, isotropic two dimensional porous medium where pores are represented by circular holes (see Figure 3.4). 1- Method of effective moduli: A s mentioned earlier, use of the infinite medium approximation requires the replacement of the elementary solutions u i n kl ijkl ~ ijkl S 2\n\ *r 50 ' 1 1 ' (3.27) Figure 3.4. Schematic representation of the porous medium and choice of an rve by the corresponding solutions in the infinite medium. The symmetry of the problem being exactly that of the boundary dr, the equivalent material is necessarily isotropic. The compliance tensor S needs only two coefficients to be defined; namely, the effective Young Modulus E and effective Poisson ratio v. These ones are related to 5 as follows Sun -' (3.28) e °//22 — 51 Chapter 3 Homogenization of Heterogeneous Media 2 X A 11 27 Figure 3.5. Study of the rve associated to the porous medium subjected to Z" loading The isotropy of S allows the determination of its other non zero components in the following way \S 2222 — SJJJJ 1+v 1^1212 ~ ^1111 $1122 (3.29) ~ Thus, the computation of 5 requires the knowledge of only the elementary solution u" This latter one is provided on dr by Muskhelishvili [126]: u" \u'J = R =—(3Cos 0-Sin 0) 2R 2 2 (3.30) Sin20 52 Chapter 3 Homogenization of Heterogeneous Media and hence, the expressions for S 1U1 and S U22 are given by (3.31) Introduction of the porosity of the material, defined as the ratio of the cavity area to the surface of the rve (3.32) leads to l + 3d v = (3.33) v +d(l-3v ) 0 0 It is interesting to notice that, in the case where the Poisson coefficient of the matrix is equal to 1/3, the effective Poisson coefficient is also equal to 1/3 regardless of the porosity of the material. This result, also mentioned by Vavakin [170], for a three dimensional medium, but for v = v = 0.2. 0 It is equally possible to characterize the equivalent medium by its coefficient of compressibility K and its shear modulus ju which are defined as functions of E and v by K = 2(1-v) (3.34) E 2(1+ v) 53 Chapter 3 Homogenization of Heterogeneous Media Hence, the expressions for the considered porous medium 1 K = K n (l + 3d) \M = Moil* l-v 3d) l-v n (3.35) n 2- Self-Consistent Method S is computed using a similar equation as the one used in the previous section with the only difference residing in the fact that the two dimensional medium characterized by (E ,v ) 0 0 is replaced here by the unknown medium we are trying to determine (E, v). The implicit equations in the case of the considered porous medium are given by 3d_ J__J_ E~ l E + n V _ V (3.36) d_ ~E~~E~ n E leading to E = E (l-3d) 0 v =v + 0 d(l-3v ) (3.37) 0 According to this method, - The effective Y o u n g modulus vanishes for a porosity d -1/3, which corresponds to pore radius equal to 0.303 for an rve of a unit width. - For v =l/3, 0 the effective coefficient v does not vary with d and remains constant and equal to its original value v. 0 54 Chapter 3 Homogenization of Heterogeneous Media 3- Periodic Homogenization W e make use in this section of a numerical analysis from the literature [51a]. A s can be seen on Figure 3.6, a regular hexagone has been chosen as a unit cell or rve in order to ensure the isotropy of the effective medium. Figure 3.6. Unit Cell for the isotropic medium Figure 3.7, presents the different predictions of the ratio k/k 0 for the coefficient of compressibility of the porous material to the one of the matrix. A significant difference is noticed in the predictions as soon as the porosity becomes greater than 0.2 which suggests that these predictions, made for the particular case of circular pores, are very sensitive the fundamental assumptions of each method. 55 Chapter 3 Homogenization of Heterogeneous Media Figure 3.7. Comparative chart for the predictions of three homogenization techniques for a porous medium 56 Chapter 4 Continuum Damage Mechanics CONTINUUM DAMAGE MECHANICS INTRODUCTION Many engineering materials, such as concrete, wood and composites exhibit a nonlinear behavior, called damage, and which could be described in terms of microcrack initiation, growth and coalescence leading to the creation of macrocracks. This latter one could be viewed as a sudden localization of damage, the propagation of which would ultimately lead to the ruin of the structure. Even though the idea of reinforcing relatively brittle building materials with fibers has been known and practiced since the ancient times [30], modeling the mechanical behavior of such materials is by no means a trivial task. This is because the nonlinear behavior of these materials depends on the type, size, distribution and orientation of microdefects, fibers and other inclusions within the material. n Figure 4.1. A damaged element To model this behavior in a continuum sense, two different approaches are usually available. The phenomenological approach [101, 102, 103], where relations between global fields are based on the relationship between macroscopic quantities obtained from experiments, and the micromechanics-based approach [33,72, 76, 67] where the mechanical behavior of the material is considered to be a consequence of the behavior of its constituents and their interactions. 57 Chapter 4 Continuum Damage Mechanics In the phenomenological continuum damage mechanics PCDM, a set of continuous field variables that represent the damage of the material,'called damage parameters or internal variables, is introduced without consideration of microscopic phenomena. The mechanical properties of the material at a certain state of damage are prescribed by the damage parameters. The evolution law that determines the value of the damage parameters for a given macroscopic stress or strain is usually deduced from the general framework of the thermodynamic of irreversible processes, and its parameters identified by simple fitting of experimental data. Micromechanics-based continuum damage mechanics MBCDM, offers a natural way of studying the mechanical behavior of heterogeneous materials by considering the global behavior as a consequence of the local one at the scale of the heterogeneities. The MBCDM models thus make an explicit use of local mechanisms and need good knowledge of them to get good predictions. However, this is not always apparent since simplifying assumptions are usually made in order to make the problem manageable or just because the local mechanisms are not well understood. PCDM models, on the other hand, identify their parameters directly from the experiments where the local mechanisms are implicitly accounted for through their final effect on the experimental data. Due to their low tensile strength and strain capacities, fibers are added to cementitious matrices to improve both strength and "ductility" of the material. However, this improvement makes the overall response highly nonlinear and adequate constitutive modeling needs to be developed for such materials. Usefulness and potential of PCDM approach in modeling the complete response of CFRC beams is demonstrated in this chapter through comparison of numerical predictions and experimental results. Ultimate failure of a solid is the result of a succession of complex phenomena, schematically represented in Figure 4.2, that could be divided into two major stages creation of a macrocrack stable or unstable growth of the crack in the solid 58 Figure 4.2. Schematic Representation of the failure process of materials A macroscopic crack is a cut or discontinuity in matter, that is sufficiently large compared to the microscopic heterogeneities but still sufficiently small so that it can be contained in a representative volume element (rve) in the sense of continuum mechanics. Studying the propagation of cracks larger than the one introduced in this definition requires the consideration of modification in the stress and strain fields due to the important perturbations introduced by these geometric defects. This is the realm of fracture mechanics (see chapters 5 and 6). Predicting and analyzing the creation of macrocracks is based on the so-called damage theories. Damage represents the progressive deterioration of the material's cohesion under the action of loads leading to the fracture of the rve. This deterioration process is a complex phenomenon in the case of FRC materials and a full understanding of the mechanisms underlying it has not been reached yet. As can be seen in Figure 4.3, the pre-peak tensile behavior of FRC is characterized by the process of microcrack propagation in the matrix prior to the formation of a continuous crack system across the critical section at the peak load. The process by which load is 59 Chapter 4 Continuum Damage Mechanics Multiple Cracking . HighV f \ Strain Softening Strain Crack Opening Strain Figure 4.3. Typical stress-strain curves for lowfibervolume and highfibervolume FRC transferred from the matrix to the fibers, and the bridging effect of the fibers across the matrix cracks which occurs at a later stage of loading are believed to be the main mechanisms controlling the effectiveness of fibers in enhancing the mechanical properties of the brittle cementitious matrix. The pull-out action of fibers would be mobilized when a crack tends to widen at a critical section; this starts at peak load and tends to dominate the post-peak behavior. A significant amount of energy is consumed in the process of debonding and fiber pull-out. It is generally accepted that for FRCs with Vj<2%, the major contribution of the fibers is after the matrix strain localization, which occurs around the peak of the tensile stress strain curve. However, new processing techniques have helped in the manufacture of thin-sheet products with fiber volume fractions (Vy) as high as 15%. In this type of materials fibers not only improve the ductility of the material, but lead to a significant increase in the strength of the composite. This increase in strength termed pseudostrain hardening is associated with the appearance of multiple cracking in the specimen which requires a higher energy input to open the microcracks. Figure 4.3 provides a typical stress strain response for a high Vf with discontinuous fibers. The response consists of 3 regions: I - elastic region (up to matrix cracking) II - multiple cracking region (up to max. post cracking pt) III- Failure region (crack opening localization) 60 Chapter 4 Continuum Damage Mechanics The complexity of the mechanisms underlying the characteristic behavior of FRC composites suggests that modeling the damage phenomenon, in order to predict the entire response of structural elements can only be a schematic and a crude approximation of the physical reality. In this section a brief review of the theory of damage mechanics before presenting the selected model is given. The theory of damage mechanics goes back to the end of the fifties, with Kachanov's work regarding the creep of metals [90]. But it took the seventies and the early eighties for this theory to receive a sufficiently general base to allow its application to other types of material [102]. At the beginning of the eighties, it was established that damage mechanics could model the strain softening of concrete in a very satisfactory way [114]. The main tool to achieve this goal is the thermodynamics of irreversible processes which provides the definition of a variable representing the matter deterioration while considering this latter one as a continuum. This justifies the use of the classical continuum mechanics equations. Moreover, the hypothesis of a dissipation potential provides the framework leading to damage evolution laws as a function of the loads. Hence, by considering the material as a system described by a set of variables and a thermodynamic potential, the constitutive laws are systematically deduced along with conditions on the kinematics of damage. But an adequate choice of the potential and the damage variable (scalar, tensor, etc.) remains to be made. The notion of hidden or internal variables, introduced by the thermodynamics of irreversible processes, allowed the solution of a certain number of fundamental problems in elasticity, elastoplasticity and elastoviscoplasticity of materials. These hidden internal variables representing the internal state of matter are formally introduced within the general framework of thermodynamics. Unfortunately, this approach cannot take into account the quantitative aspect of the physical phenomena, and often leads to a number of internal variables that are not compatible with their experimental identification possibilities or procedures. To overcome this inconvenience, a more pragmatic approach 61 Chapter 4 Continuum Damage Mechanics is usually adopted. In this alternative approach, one proceeds from the concrete needs towards some form of abstraction according to the following scheme, Describe in a precise manner the phenomena governing the response of the material for the particular type and/or range of loading we are interested in. This is the case because the same material may demonstrate more than one constitutive law to model its behavior (static, dynamic, elastic, plastic etc.) - To capture the effects of these phenomena, one needs to define representative variables and conceive the relationships among them. Usually, this is achieved through a choice of a damage variable and its effects on the mechanical properties of the material. The logical framework for this part is provided by the pre-established thermodynamics concepts. The necessary and sufficient conditions for a complete description of the macroscopic thermo-mechanical process evolutions are given by u: displacement field e : elastic strain tensor field e e' : inelastic strain tensor field n cr: stress tensor field T: absolute temperature f : field of body forces v e\ specific internal energy q: heat flux s: specific entropy X: set of internal variables These quantities are determined for each material and/or specific case of loading through the conservation and balance equations of thermodynamics of continuous media together with the constitutive equations expressing in terms of other quantities, the following functions: 62 Chapter 4 Continuum Damage Mechanics CT S" X h s <f> = (f>(e ,T,%) Helmholtz free energy e The principles of thermodynamics allow the elimination of a certain number of variables in the constitutive relations. If the scope of the analysis is limited to quasi-static cases, the thermomechanical couling is neglected, and the mechanical behavior of the material can be described without the constitutive laws of the three thermodynamic variables (heat flux, entropy and Helmholtz free energy).This leads to the following system of tensorial evolution equations a = F (e ,T, ) e e X e =F { \T, ) in in X= F ,F e in and F are h £ X F (e\T, ) h X the elastic, inelastic and hidden (internal variable) functionals respectively, to be specified for each material. The first equation describes the elastic behavior of materials and is identified with the generalized Hooke's law, whereas the second equation is a plasticity or viscoplasticity constitutive relation. Existence of a dissipation pseudo-potential associated to the normality rule provides a systematic scheme for determining the functional F . The last equation is a formal representation of !n the evolution laws of internal variables. The choice of internal variables is determined by the phenomena that need to be modeled. A hidden variable is associated with every phenomenon that cannot be represented by a measurable variable, such as the internal degradation of matter. The internal variable should not contradict the principles of thermodynamics, and have a sense that is compatible with the acquired knowledge regarding the mechanisms 63 Chapter 4 Continuum Damage Mechanics governing the phenomena. Keeping only the phenomena o f damage and elasticity, the two obvious variables are the strain tensor and a damage variable (to be defined). The damage and strain variables represent physically these phenomena only i n some average sense. This is because a material point is usually identified with a characteristic volume element the material is made of, the so-called representative volume element (rve). F o r instance the characteristic dimension in the rve for metals is of the order o f 0.1-1mm, whereas for concrete this length varies from 1 to 10 cm. Hence trying to give a sense to the strain or damage below that scale would be meaningless. Hypothesis of the effective stress The damage variable is classically defined from the hypothesis of effective stress to which it is easy to give a physical sense in the unidimensional case. Consider a tensile test where the specimen is under a stress cr. If that specimen possesses cavity-type microdefects such as a microcrack and/or microvoid, then only a fraction S of the apparent surface S is resisting the load. The stress cr is the stress acting on the effective section S : 5 = — \a(M) a ds (2) where (M) is the distribution of stress concentration coefficients that exist at crack a vicinities. Since the crack mechanisms and internal crack profiles and densities are not sufficiently understood for evaluation of (M) and S , the following definition of damage a has been introduced. The damage variable is the operator A which when applied to the usual stress cr gives the effective stress & ,cr = zl(cr). The hypothesis o f the effective stress implies that the constitutive law of any damaged material is the same as the constitutive law of the plain material where the usual stress is replaced by the effective one in the equation. When a given load creates a damage characterized by the same loss of resistance in all directions, damage is called isotropic. In this particular case, the 64 Chapter 4 Continuum Damage Mechanics damage variable, which is usually a tensor, is reduced to a scalar and therefore the effective stress is related to the nominal stress by, a = a \-D This way, damage in the sense of deterioration which leads to fracture could be represented by a single scalar variable D such that D = 0 corresponds to the undamaged plain material, and D = 1 for a fully fractured volume element. The parameter D, introduced here as a hidden thermodynamic variable, has the physical sense of a section or a volume weakened due to the development of microcracks. D represents the phase of creation and evolution of the process of deterioration. Notice that the variable D is not directly related to the geometry of microcracks or microvoids; but is related to these quantities only through the equivalence between the behavior of the plain material (written in terms of the stressor) and the behavior of the damaged material (written in terms of the effective stress a). In the case of an elastic behavior, the 3D constitutive law of the material takes the following form a = H :s for the plain material e 0 cr = H(D): s for the damaged material e with H = H(D = 0), H being the elastic rigidity tensor and e is the elastic strain tensor. e 0 In the general case of an anisotropic damage, a scalar representation is no longer possible, because it involves an anisotropic evolution of the mechanical characteristics and therefore a tensorial representation is necessary. Among the difficulties that are raised when working in this direction, one can cite: the respect of compatibility with the isotropic model, the symmetry of the operator H(D) necessary condition for the existence of an elastic potential and the possibility of non radial loads in which the directions of the 65 Chapter 4 Continuum Damage Mechanics principal stresses are not fixed anymore. In the general case of behavior modeling, this leads to the introduction of a fourth order damage tensor, often difficult to be experimentally identified [42]. D A M A G E M O D E L S F O RCEMENTITIOUS M A T E R I A L S [103,116,117] A t constant temperature, concrete can be described by the elastic strain tensor e , e permanent strains e p tensor and the damage tensor = e = £ +e e p the D. ^\e e"\dt p where e is the tensor of plastic strain rates. Each equilibrium state is p characterized by a scalar thermodynamic potential py/ function of e , e E and D where P p is the mass density of the material. A common choice for the function ^ i s the quadratic form represented by the free energy. Assuming that damage affects only the elastic characteristics of the material, one can write py/ = py/ +pyy f where y =y/ (e ,D) e e y e =y/ (£ ,D) p p p and CT = the elastic stress de e d(p ) tensor p CT = Y = ¥ ds ——: the plastic stress tensor the strain energy release rate 3D The phenomena of damage and creation of permanent deformations are irreversible processes leading to a dissipation of the mechanical energy in the form of surface 66 Chapter 4 Continuum Damage Mechanics creation and heat. Application of the two fundamental principles of thermodynamics impose that this dissipation be positive. This is usually represented by the ClaussiusDuhem inequality: <f> = a: e - p\j/ - py/ e >0 p It is easy to distinguish in this relation between the damage dissipated energy rate if> and d the plasticity dissipated energy rate <fi . p <j> =<r:£ -py/ e e d (f> =o: e - py p p p A sufficient condition to satisfy the Claussius-Duhem inequality is given by j> >0 d A s introduction of plasticity in the models follows the classical rules [57] we w i l l focus our attention on the introduction of damage into the elastic constitutive laws. W e adopt the elastic energy as our potential energy py/ e =^H(D):e e where H(D) py/ e :e e is a fourth order tensor interpreted as the secant rigidity matrix. A t this stage we notice that the choice of py/ e has been reduced to the choice of the dependence of the rigidity matrix on the damage variable. The effect of damage appears in H(D), through the evolution of mechanical properties of the material. Damage evolution is governed by a loading surface of the form f(e,H,K )=0 0 67 Chapter 4 Continuum Damage Mechanics where K is the initial damage threshold. A n y loading path that stays within the surface 0 does not lead to any damage evolution, D=0, whereas any loading path that tends to leave the surface w i l l pull it with it. This latter one becomes larger and damage evolves, D>0. Damage Model with Anisotropy and Permanent Deformations Choosing a thermodynamic potential of the form pyf = ±[H :e'}\e +8>] e D with e is the elastic strain tensor e is the permanent strain tensor e p H D is the rigidity matrix of an orthotropic material leads to a formulation incorporating the effects of induced anisotropy. H D according to the following formula H~': or = H~': [L D : a] where H L 0 is the initial rigidity tensor of the material D is a fourth order symmetric tensor representing the damage variable This model divides the total response of a material into three stages: • Stage 1 A t this stage the material behaves in an elastic isotropic manner s = H„' 68 :tr on damage Chapter 4 Continuum Damage Mechanics • Stage 2 Beyond a first threshold, damage appears by affecting in an anisotropic way the mechanical properties of the material e=H a 1 D • Stage 3 Beyond a second threshold, permanent deformations are induced s = H- :a 1 + e p In the most general case, nine scalar damage variables are necessary to describe the model hi ^12 113 hi h2 _hi hi h3 _ In the axisymmetric case, however, only 4 scalar variables are needed to represent damage h h2 ^23 ^12 h hi If h=h=l h2 = hs ^23 we arrive at the case of isotropic damage with / = hs (1-D) Experimental tests on cylinders and cubes have shown that there is a significant drop in the modulus of elasticity in the direction parallel to the loads along with a substantial increase of the Poisson ratios in the perpendicular directions. This led to the following evolution of characteristics: 69 Chapter 4 Continuum Damage Mechanics = o y v _B[l-(l-v )A](£-K ) 0 + v 0 0 B(s-K ) +1 0 with A and B denoting material parameters to be calibrated from experimental tests. Identification of this model requires compressive and axisymmetric triaxial test data. Unilateral Model This model assumes that the material is elastic damageable and that damage is isotropic. T w o independent scalar damage variables are used. These two variables are assumed to be uncoupled in order to account for the closing effects of cracks when the signs of stresses are changed (unilateral character). The constitutive law of the damaged material then takes the following form cr, + 1-D, H,o-i 1~D C where H,o' :cr =^ [u + t '0 r '0 with 70 v )(<T) -v (Trcr) l] 0 + 0 + Chapter 4 Continuum Damage Mechanics In the Claussius-Duhem inequality we have two damage energy release rates, each of which is related to the corresponding scalar damage in the following way 0 =Y D +Y D d t t c c with dD, y = d(py/ ) e 3D C A necessary condition to satisfy this equation is given by D >0 t D >0 c A damaged threshold is defined for each damage type in terms of damage energy release rate f (e,H,K ) = f (e,H,K ) = t 0 c 0 Y -K,(D,) t Y -K (D ) c c c Evolution laws for each damage variable are given by D =F (Y ) C C C Thus, the behavior of the material is given by a = E (1 -D ) s in tension a = E (1 -D ) s in compression 0 0 t c The assymetric behavior displayed by cementitious materials is automatically accounted for, since the damage thresholds and variables are different in compression and in tension, and the memories of D and D are conserved when a change of sign occurs. t c 71 Chapter 4 Continuum Damage Mechanics Scalar Damage Model [114] Keeping only the two phenomena of damage and elasticity, one has a natural choice of variables consisting of the elastic strain tensor and a damage variable (to be defined). The damage and strain variables however physically represent these phenomena only in some average sense. This is because a material point is usually identified with a characteristic volume element the material is made of, the so-called representative volume element (rve). For instance, the rve for metals is of the order of 0.1-1mm , 3 whereas for concrete rve is of the order of 1-10 c m . Hence trying to give a sense to the 3 strain or damage below that scale would be meaningless. The physical phenomena of nucleation, propagation and coalescence of cracks are represented in the framework of continuum mechanics by a scalar variable called damage. The choice of this model is mostly motivated by finding a compromise between the needed accuracy and the cost of calculations. The constitutive law for an elastic material that damages according to the scalar model is given by <r = (l-D)H :e with 0 0<D< 1 The damage evolution law is piloted by the equivalent strain s defined by with where Sj represents a principal strain. Since concrete exhibits an asymmetric behavior in tension and compression, the damage evolution laws in tension and compression are different, the evolution being slower in compression. D,=D (£,A,,B ,K ) t t 0 D =D (s,A ,B ,K ) c c c c 0 72 Chapter 4 Continuum Damage Mechanics where D is the damage due to tension and D the damage due to compression. K t c 0 is the initial damage threshold. At and Bf are parameters modeling the shape of the tensile curve in its nonlinear part, and finally A c and B c are parameters modeling the shape of the compression curve in its nonlinear part. At, Bf, A , c B c are considered to be material characteristics. The case of damage due to a multiaxial loading is defined by an intermediate value between D c and D t which is a function of the two contributions v i a the following linear relation: D = a,D +oc D t c c where ct : Coefficient accounting for tensile damage t a : Coefficient accounting for compressive damage c - C a s e of pure tension a =0 and - Case of pure compression a =1 and - Combination a +a,=l Determination c c (D = D ) t a,=0 (D = D ) C c of a and t a c a: tensor of principal stresses with <T = <7 +<J_ + tx : tensor of positive principal stresses + <y_ : tensor of negative principal stresses e =H .a_ l c with H" : Tensor of mechanical characteristics of the material 1 Hence and 73 Chapter 4 Continuum Damage Mechanics [1 with =I if 0 e i "Rrv >0 Otherwise The evolution laws of D and D t c are identified from the corresponding uniaxial tests. Direct tensile and compressive tests for two volume fractions (V^ = 2%, and V f and two fiber lengths ( L f =3%) = 3 mm, and L = 10 mm) have been performed on pitch based f carbon reinforced mortars. Figure 4.4 shows the ability of the model to represent in a very accurate way the uniaxial behavior of fiber reinforced cementitious materials. The threshold function is expressed by: f = s-K(D) The damage evolution is piloted by the "loading function"/as follow: Y'CASE if f=0 and f =0 Then Dt = Dt (e, At, B , Ko); t _ x (l-AQKp A, e Exp[B (6-K )] t 0 D = D (e, A c , B , K ) c _ 1 c c (l-Ae)Ko 8 2 nd CASE The random A c Exp[B (8-K )] c /<0 , or if \ \f = 0 and distribution 0 0 . f =0 Then of local strengths due AD = 0 to local defects microvoids, etc.) leads to a volume influence on the behavior. 74 (pores, cracks, Chapter 4 Continuum Damage Mechanics • Experiment - - Model • \ >• • / 1 w • / t I 1 > i t t i i / — Reloadii * A • t 0 / ^ — Jnloading 1 Tens inn it— (Vf=37c) 0.000 0.002 0.004 0.006 0.008 Strain (%) 70 E x p e r i m ent M 60 odel # V 50 / | / 40 / <> * 4s 3 0 to A / f / / Rel oading a 20 0 * nloading 10 Compression (Vf=3%l 0 0.000 0.002 0.004 Strain 0.006 0.008 0.010 (%) Figure 4.4. Model's capabilities to represent the tensile behavior and the compressive behavior of FRC 75 Chapter 4 Continuum Damage Mechanics One way to treat the scale effect is to consider it as a consequence of the assumptions used to model the material and its behavior, which allows one to include this effect in the model by taking into account the particular properties of the material and the mechanisms of its rupture. The main influence of the material heterogeneities on the damage evolution is taken into account v i a a multiscale approach which consists in the simultaneous use of information coming from two different scales [115,116]; the scale of volume element and the scale of the rve. Inside a volume element, stresses and strains are computed in a traditional way. The equivalent strain e is evaluated by considering it as a function of the state of deformation in a certain neighborhood of the considered point which is its rve. M a n y previous studies [55,8a,69a] have shown the efficiency of multiscale approaches in the prediction of fracture behavior of materials and the convergence of numerical results with mesh refinement. Different studies have shown that use of local models in numerical applications presents some problems [115,125]. For instance, in the case of a large gradient of stresses, there is a strong dependence of the structural response upon the mesh. Damage evolution is influenced simultaneously by the size and orientation of the elements. The elements located right under the load get damaged first and usually this damage is confined to a narrow region, generally of the size of one element, such that a progressive refinement of the mesh gives results that do not agree with reality. T o determine the origin of these problems, different studies have been conducted on the mathematical aspect as well as on the physical one. On the mathematical side, some authors [27,132] have shown the existence of a certain critical level of deformation beyond which the solution is not unique anymore. This result suggests that beyond that critical point, the finite element solution is not objective and depends strongly on the chosen mesh. O n the other hand, from the physical point of view, studies on the uniaxial loading of specimens [21,22] have shown that the stability in the case of materials exhibiting strain softening is possible only i f the size of the localization zone is larger than a certain critical value. This led to different suggestions on the way to define that critical zone, where the main contributions are: 76 Chapter 4 Continuum Damage Mechanics - limitation of the size of elements [24] - introduction of strain tensors of higher order gradients [26] - use of a nonlocal approach; either for the constitutive law [23] or for the damage evolution law [116,136]. The last possibility has been chosen for the adaptation of the local model, since the only nonlocale quantity is the damage evolution law, which does not require any special modification nor the introduction of supplementary boundary conditions in a finite element analysis. This led Mazars et al. [136] to develop the following multiscale approach. In order to stay in the continuum mechanics framework, the considered damage evolution is nonlocal. Therefore the evolution of the local state of damage is a function of the materials state in a certain neighborhood. The retained approach consists of evaluating the average of the piloting damage variable s representative av on a certain volume V called volume of damage. In the three dimensional case this volume is represented by a sphere of radius Lp, where LQ is the so-called characteristic length which is supposed to be a material property [20,21]. This approach is termed multiscale because it makes use of information coming from different scales; on the one hand the scale of the representative volume element where damage is defined relative to the material heterogeneities and crack interaction and, on the other hand, the scale of the elementary volume element (in the sense of continuum mechanics) where stresses and strains are evaluated in a traditional way. Therefore, instead of computing the traditional local equivalent strain, we compute the averaged quantity defined by 77 Chapter 4 Continuum Damage Mechanics where g(x-s) is the weight function representing the intensity of interaction between the considered point with coordinates x and its neighboring points having coordinates s contained in the representative volume element surrounding the point x. Results from micromechanical studies [20,32,136] led to the following relations: 2(x-sj^ r g(x — s) - - Exp V g(x-s) L if \x-s\< L D J D if =0 \x-s\>L D APPLICATION OF PCDM TO FRC MATERIALS The one parameter damage model described earlier, has been implemented in a finite element code to predict the uniaxial and biaxial bending behavior of C F R C materials. 1- BEAM BENDING The model has been incorporated into the constitutive law of a 4-noded 2 D elasticity isoparametric element and used to predict the bending behavior of pitch-based carbon fiber reinforced cementitious mortars, see Figure 4.5. The characteristics of the fibers used in this study are given in Table 4.1. Diameter Fiber Type Carbon Fiber Pitch-based Mortar Specific Gravity (pm) 18 Tensile Strength E (MPa) (GPa) 1.7 590 W/C S/C SF/C L = 3 mm 0.35 0.5 0.2 L = 10 mm 0.35 1 0.15 f f Table 4.1. Mix proportions and fiber characteristics 78 Chapter 4 Continuum Damage Mechanics 300 mm Figure 4.5. Experimental setting of the three point bending test Two fiber lengths, L =3 mm and L =10 mm have been considered. For each fiber f f length, two volume fractions, V = 2% and V = 3% have been incorporated into the f f cementitious composite. A summary of the comparison between the model predictions and the experimental results is given in Table 4.2. A s can be seen in Figures 4.6 and 4.7 and on Table 4.2, the model predictions of the bending response of C F R C beams are fairly close to the experimental data. One important shortcoming of this continuum model that needs to be mentioned is related to the failure pattern of the specimens. A s shown in Figure 4.8, the model predicts a diffused damage within the material up to failure of the specimen, whereas experimental observations, show that failure is usually a result of the creation and growth of a macrocrack. Thus, P C D M can give a good representation of the physical reality only up to peak load, range where microcrack nucleation takes place in a continuum with no significant growth is the primary deterioration process. After peak load, the response of the specimen or structural member is mainly controlled by the opening and growth of a macrocrack due to propagation and coalescence of microcracks and the phenomenon of localization is said to occur. A t this stage of loading, use of stress-strain functionals to describe the behavior of a quasi-brittle material is physically inadequate since the 79 Chapter 4 Continuum Damage Mechanics observed deformation is rather due to the opening of a crack than to a deformation in a continuum sense. Consequently, PCDM approaches do not capture the physics of the deterioration process after localization and use of existing nonlinear fracture mechanics theories should be more appropriate for a description of the response of the material at this stage. Nonetheless, CDM can still give a fairly accurate prediction of the overall behavior (load vs displacement) of structures. Finally, developed within the general framework of thermodynamics of irreversible processes, PCDM offers a great potential for extending its application to include the effects of durability issues on the long term behavior of cementitious composites. Such an extension should facilitate the development and implementation of an integrated durability-structural design procedure [144] L f Peak Load = 3 mm (KN) Peak Area under Displacement curve (mm) (J) V, = 2% 362.37 0.515 293.58 V,=3% 404.12 0.401 309.25 V = 2% 366.91 0.482 278.64 V =3% 414.17 0.363 294.55 Peak Area under Displacement curve (mm) (J) 363.60 0.601 227.29 421.78 0.890 579.88 2% 367.03 0.483 220.68 V = 3% 413.90 0.740 545.36 Experiment f Model f Peak Load Lj = 10 mm (KN) V = f 2% Experiment V,=3% V = f Model f f Table 4.2. Comparison between model predictions and experimental results 80 Chapter 4 Continuum Damage Mechanics 400 I l l - I l l J\ f i j 350 300 250 — Lf=.3mm •-2% 1 - i ': 230 / 150 Model -Bperirrent \ \ 100 3D L QO 1.2 Q6 1.8 30 24 Deflation (mm) 450 400 * * \ 350 1 300 9-250 J< •« 3200 150 100 50 7 I 1 1 I Lf=3mm , Vf-3% ^ \ • \ i \\ -Experiment •Nfcdel | [ 0.0 0.6 24 1.2 3.0 Deflection (mm) Figure 4.6. Bending behavior of CFRC materials with L =3 mm f Chapter 4 Continuum Damage Mechanics 400 \ 350 300 200 r 1 - Model -Bperinent \\ If ij 100 r \ V 250 150 I Lf=i Omtr> , vj=2% 50 0.0 0.6 12 3.0 24 It Deflection (mil) 450 400 350 1j —I if HUH I v— 300 ?250 & - Experimental -Numerical | 200 150 100 50 0 0.0 0.5 1.0 1.5 20 25 3.0 3.5 4.0 4.5 Deflection (mm) Figure 4.7. Bending behavior of CFRC materials with Lj =10 mm 82 5.0 Chapter 4 Continuum Damage Mechanics CONCLUDING R E M A R K S • A one scalar damage variable continuum model has been implemented into a 4-noded 2D elasticity isoparametric finite element and used to predict the complete response of 4 specimens of CFRC beams up to failure • In spite of its relative simplicity (one scalar damage variable), the model is still capable of predicting most of the characteristic behavior of CFRC beams • One major shortcoming of the model resides in the incompatibility of a stress strain constitutive law with the physical reality after peak load where the response of the specimen becomes controlled by the opening of a major crack rather than by a deformation of matter in a continuum sense Damage Scale • • 0<D<0.3 0.3<D<0.5 Numerical EJ 0.5 < D < 0.7 El 0.7 S D< 1 Figure 4.8. Typical damage distribution as predicted by the model and experimental failure patterns 83 Chapter 4 Continuum Damage Mechanics • A better way of modeling the response of C F R C structures could be achieved by combining C D M and nonlinear fracture mechanics concepts ( N L F M ) through the introduction of an equivalence between a diffused damage and a localized crack so that continuity of the response is guaranteed. This way, C D M could be used from the unloaded stage up to peak load, stage at which damage is converted into a localized crack, and a suitable N L F M theory can be used from then on up to failure P L A T E BENDING In this section we recall briefly the Reisner-Mindlin theory of plates along with the way of introducing damage into a finite element for thick plates called Q4y. Use in a strict sense of the variational forms that have been developed in the literature for the theory of plates with transverse shear and having a reduced number of nodes and nodal variables presents many difficulties. This led to the formulation of some very useful elements (no shear locking and no parasite modes). But this is possible only at the price of a loss of rigor in the application of the variational forms. Q4;K[18] belongs to this category of elements. First, we choose a reference set of coordinates such that the xy plane coincides with the midsurface, the z axis being transverse. The selected theory is based upon the following assumptions. - the normal stresses are negligible - the normal strains are negligible s » 0 with respect to the other stresses; u « 0 - plane sections remain plane after deformation, without necessarily conserving a right angle. - the induced stresses in the midsurface transverse displacements corresponds to a first order are by flexural of the order approximation 84 strains are neglected. If the of the thickness, this approximation Chapter 4 Continuum Damage Mechanics This latter assumption does not present any particular difficulty. The first order theory we chose to model the plate uses the following 5 independent variables - u and v: membrane displacements (in the xy plane) - w : displacement along the z axis - B : rotation of the normal to the mid-surface in the plane xz x - B : rotation of the normal to the mid-surface in the plane yz y The displacement throughout the plate is given by u(x, y) < v(x,y,z) =. v(x, y) ' + Z-P (x,y)w(x,y,z) w(x, y) 0 u(x, y,z) > y Application of the formula for small perturbations, s.. = -^-(w, . +u ), ji allows one to define the state of strains as a function of the 3 variables w, B and B by the following x y kinematic relations wwith u v, {eh 2s ftx,x ,y .x U +V where {%} represents the curvatures, z {^Jthe flexural strains, {e}the membrane strains, and {/} the shear strains. 85 Chapter 4 Continuum Damage Mechanics Since our element has a new constitutive law, the rest being basically the same as the element described in reference [18], we w i l l focus our presentation on the issues related to the implementation of damage into the original law. The stress-strain relationship could be written in a matrix form as follows where {o" }, {r } are the initial stresses. 0 0 For an isotropic material we have {H} =1-v 2 1 v 0 v 1 0 0 0 l^- 2 1 0 2(1 + v) 0 1 The generalized stresses are related to the generalized strains by {M}=[H ]{e} [H ]{ }+{M } mf + F X 0 {T} = [H ]{y} {T }{z} c where {N}, + {M}, receptively. [H ], m 0 { r } represent the normal force, bending moment and shear force [// ] / , and [ # ] m/ are 3x3 symmetric matrices and represent the elastic characteristics of membrane, flexion and membrane-flexion coupling. Carbon fiber reinforced cementitious composites are considered to be isotropic materials, since the fibers are randomly distributed in the bulk of the matrix. Moreover that material was used in plates that are symmetric with respect to the mid-surface. This allows one to take \H mf ] = 0, since membrane and flexural characteristics are decoupled. 86 Chapter 4 Continuum Damage Mechanics The stiffness matrix [K] is given by [*]=W+fcl where +t [ i ^ ] and [Z? ] are matrices relating the curvatures and the shear strains, respectively, to c the nodal variables. Integration in the xy plane is performed using a 2*2 Gauss scheme. If no damage has occurred yet, [H ] / and [H ] are evaluated as in the elastic case: C 1 V 0 V 1 0 0 0 1-v h E 3 n 12 (1-v ) 2 1 0 12 (7 + v) 0 1 2 Introduction of the scalar damage variable is a particularly easy task, since the only elastic characteristic affected by damage is Young's modulus. E = Consequently \H \ F and [H ] C E (1-D) 0 are given by: 87 Chapter 4 Continuum Damage Mechanics 12(l-v ) 1 v 0 v 1 0 0 o !=!L +1 \z (l-D)dz 2 2 1 {H } = ° 12{l + v) 0 5 0 \(1-D)dz E c 1 The parameter D depends on the considered point and the level of loading since it is directly related to the state of deformation. Since damage distribution is not known, we used a simple numerical integration scheme through the depth of the plate at each Gauss point. The depth of the plate was divided into a number of equal intervals, then damage is computed at the center of each interval and assumed to be constant over that interval Figure 4.9 shows a flow chart that summarizes the different steps of the adopted numerical scheme. If first iteration and first step (D=0), adopt E = E n • Compute the global stiffness matrix K[D = 6] and the left hand load vector {F}, • Solve for {u} ; and get the principal strains e . This allows one to compute the { average equivalent strain e av and the damage variable i f the initial threshold is exceeded. If this is the case compute the new global stiffness matrix using the new value of D • Compute the residual {R}=[K]{AU}-{AF} • Test i f there is convergence with the prescribed condition where e is the prescribed precision, and ||/?|| is the norm of the residual vector {R} p • If convergence Then - add a load increment - go to the next step (j=j+l on the flow chart) 88 Chapter 4 Continuum Damage Mechanics RESULTS AND DISCUSION Results of few predictions of the model for plate bending are compared to the experimental results. Use of P C D M is motivated by the fact that prior to peak load, the deterioration of the material is mainly due to the creation, propagation and coalescence of microcracks. After the peak load there is a creation of a macrocrack due to the coalescence of propagating microcracks which would ultimately lead to the failure of the specimen (or structural component). The creation of such a macrocrack could be viewed as a sudden localization of damage. Once that localization occurs, the behavior of the specimen is usually very well predicted by existing fracture mechanics theories. However, when interested in that part of the behavior where deterioration of the material is still considered distributed over a certain volume of material and governed by the local effects of microcracking, damage mechanics provides a more rational framework for predicting the behavior of specimens and predicting the areas where macrocracks w i l l form. That is why some authors [119] have suggested to use damage mechanics to predict the response of structural components up to the peak load, and the use of fracture mechanics from then on to the failure of the structure or specimen. Accordingly, fracture mechanics and damage mechanics are correlated theories. In the numerical implementation of the non local constitutive law, the characteristic length L D which is considered by some researchers [20,21] to be a material characteristic has to be determined. Arguments that were suggested for the use of nonlocal quantities, is the fact that they allow to account not only for the presence of heterogeneities including cavity type defects, but for their interactions as well. For concrete specimens, it is suggested that L »3d D , where d is the dimension of the largest aggregate in the mix [20]. A direct application of this value is not valid since we are dealing with carbon fiber reinforced cement pastes. So in this study, the process of trial and error was used to get an estimate of L . The best results for our problem were obtained for L D D = 5 mm. This suggests that for this particular case, the nonlocal aspect is mainly due to the interaction between microcracks rather than to the presence of large aggregates. 89 Chapter 4 Continuum Damage Mechanics INPUT DATA MESH, B.C'S ...ETC. GLOBAL MATRIX K(D) = 0 LEFT HAND VECTOR / \ SOLVE SYSTEM j = j+ l i = i+l GLOBAL STIFFNESS MATRIX K{D ) M NO NO Yes END Figure 4.9. Flow chart of the main steps followed by the finite element code 90 Chapter 4 Continuum Damage Mechanics Circular support Figure 4.10. Schematic representation of the test setup Considering the symmetry in geometry, boundary conditions and loading (see Figure 4.10), only a quarter of the plate needs to be modeled when using the local approach. Unfortunately, this is no longer valid when using the multiscale approach, unless we account for the "fictitious" integration points that are not present in the considered portion of the structure but which should exist i f the whole specimen is modeled. Therefore, instead of writing a procedure that would take into account those "fictitious" integration points, we decided to choose the easiest solution of considering the whole structure, despite the increase in the computing time. A s we can see in Figures 4.11 and 4.12 , despite of its relative simplicity, the scalar damage model allowed us to predict the strain hardening exhibited by the circular plates for both volume fractions of carbon fibers. Just as in the case of beam bending in the previous section, the model predicts a diffused damage around the center of the plate under the concentrated load and at lesser extent around the support reactions, whereas an examination of experimental failure patterns reveals that propagation of macrocracks is the major mechanism governing the failure of 91 Chapter 4 Continuum Damage Mechanics the plates. A t low fiber volume fractions, failure is usually controlled by the formation and propagation of a single macrocrack. A t high fiber volume fractions, on the other hand, one would observe a lot of microcracking around the center prior to peak load and the formation of a system of macrocracks with a smaller total width than the width of the single macrocrack in the case of unreinforced matrices or with low fiber volume fractions. This suggests that at high volume fractions, the failure pattern gets closer to the damage distribution as predicted by the model. L f Peak Load Peak Displacement (KN) (mm) 0.83 0.54 V =3% 1.39 0.80 V = 2% 0.92 0.51 V = 3% 1.47 0.74 = 6 mm V = 2% f Experiment f f Model f Table 4.3. Comparison between model predictions and experimental results C O N C L U D I N G REMARKS • Based on the Reisner-Mindlin theory, the four-noded thick plate element called Q4y has been extended to include the effect of a scalar damage variable on its constitutive law • A s one would expect, the discrete failure pattern observed in experimental tests and the diffused damage distribution as predicted by the model are fundamentally different. Nonetheless, as the fiber volume fraction increases, the patterns of the two modes becomes closer to each other • The new element has been able to predict the response of C F R C plates subjected to a concentrated center loading up to peak load where the assumption of a diffused damage within the specimen is physically acceptable for both volume fractions. 92 Chapter 4 Continuum Damage Mechanics Figure 4.12. Comparison of the model predictions to the experimental results for V 93 Chapter 4 Continuum Damage Mechanics Nowadays, a lot of efforts have been spent i n developing modeling techniques that are more and more sophisticated, but unfortunately, in the case of cementitious materials, such models are never taking into account the effect of durability issues [144] on the long term behavior of concrete. Due to its general thermodynamic framework, P C D M seems like a potential candidate for such an extension. Damage Scale • 0<D<0.3 • 0.3 < D < 0.5 • 0.5<D<0.7 0 0.7 < D < 1 a) Numerical Low V, f High V, f b) Experimental Figure 4.13 . Failure pattern as predicted by the model and observed experimental failure 94 patterns Chapter 5 Fracture Mechanics Concepts FRACTURE MECHANICS CONCEPTS I- INTRODUCTION The objective of design in civil engineering is the determination of the geometry and dimensions of structural elements along with materials selection in such a way that the elements perform their expected operating function i n a safe, efficient and economic manner. This is the main reason as to why the results of displacement and stress analyses are coupled with an appropriate failure criterion, which is nothing else but a postulate for predicting the event of failure itself. Traditional failure criteria (stress criteria) have failed to give an adequate explanation for number of structural failures that occur at stress levels considerably lower than the ultimate strength of the material and with much less than expected ductility. The quest for a rational explanation to these phenomena has led to the creation of fracture mechanics. Experiments performed by Griffith in 1921 [65] on glass fibers l e d to the conclusion that the strength of real materials is much smaller, typically by two orders of magnitude, than their theoretical strength. Fracture mechanics studies the load-bearing capacity of structures in the presence of initial defects, where a dominant crack is usually assumed to exist. Therefore, it is based on the assumption that all materials contain crack-like defects, which constitute the nuclei of failure initiation. Flaws can appear i n a structure at least in three ways: 1. they can exist in a material due to its composition such as debonding in composites, second phase particles etc. 2. they can be introduced inadequate 3. into a structure during construction, such as in welds, compaction of concrete etc. they can be created during the service life of an element, like fatigue, environment assisted. In the case of cementitious materials, microcracks are usually present even before loading at regions of high material porosity near the interface between the inclusion 95 Chapter 5 Fracture Mechanics Concepts (aggregate or fiber) and the mortar. They are due to shrinkage of the mortar. Cracks are also present in the mortar matrix. Under an applied load, both types of cracks start to increase and new cracks are formed. The interface cracks extend inside the mortar and connect with the matrix cracks. Aggregates and fibers act as crack arrestors. Fracture mechanics is searching for parameters which characterize the onset of a crack to propagate. Such parameters should be able to relate structural performance to laboratory test results, so that the response of a cracked structure can be predicted from laboratory test results. This is determined as a function of material behavior, crack size, structural geometry and loading conditions. Fracture toughness, determined from laboratory tests and considered to be a property of the material, constitutes such a parameter. It expresses the ability of the material to resist fracture in the presence of cracks. B y equating this parameter to its critical value, a relation is obtained between the applied load, crack and structure geometry which gives the desired information for structural failure predictions. A m o n g the critical quantities proposed in the literature, one may mention the stress, the strain, the stress intensity factor, the crack opening displacement, the J-integral, and the strain energy release rate. Each criterion assumes a quantity that has to be related with the loss of continuity and has a critical value that serves as a measure of the resistance of the material to an unstable propagation of crack. During the process of fracture of solids, new surfaces are created i n the medium in a thermodynamically irreversible manner. Fracture mechanics is a branch of mechanics dedicated to the study of problems involving crack initiation, propagation and arrest. The phenomenon of fracture of a solid is complicated and depends on a large number of factors, including the macroscopic effects, the microscopic phenomena which take place at the locations where the fracture nucleates or grows, and the composition of the material. The study of a fracture process depends on the scale level at which it is considered. A t one extreme is the rupture of cohesive bonds in the solid and the associated microscopic phenomena. For such studies quantum mechanics principles should be invoked. A t the other extreme, the material is considered as a homogeneous continuum and the phenomenon of fracture is studied within the framework of continuum mechanics and classical thermodynamics. A m o n g the fracture processes that 96 Chapter 5 Fracture Mechanics Concepts take place at scale levels between these two extremes we can mention, microcracking of concrete, grain inclusions and voids, formation of subgrain boundary precipitates, and slip bands, movement of dislocations. The purpose of this chapter is to present, a clear, straightforward and unified interpretation of the basic problems of fracture mechanics with particular emphasis given to fracture mechanics criteria and their application in cementitious materials, with or without fiber reinforcement. II- LINEAR ELASTIC FRACTURE MECHANICS In this section, L E F M concepts are reviewed. The limitations and applicability of those concepts to model arbitrary, quasi-static crack evolutions are discussed in the framework of modern L E F M . M o r e specifically, the basic methods for determining the linear elastic stress field in cracked bodies, with particular emphasis on the local behavior around the crack tip. In 2 D , the basic concepts necessary to accomplish the modeling are the stress intensity factors and mixed mode-mode interaction theories. The most important crack interaction theories for mixed mode are reviewed It is well known that structures fail due to a variety of causes: excessive elastic deformation, generalized plastic deformation, buckling and fracture. N o structure is perfect and some initial defects w i l l always exist. The presence of flaws or cracks may cause a structure to fail earlier than expected, even though the material is still in the elastic range. This is mainly due to stress concentrations at some particular points, like the presence of flaws or re-entrant corners for example. The classical formula developed by Inglis [84] for elliptical flaws can be used to illustrate the effect of such called stress raisers. CT = CT {\ + 2alb) 0 where cr is the stress at the tip of the elliptical hole, a ( 0 5 is the far field stress, a is the major axis and b is the minor axis as illustrated in Figure 5.1. A s can be seen, when b tends to zero or a » b the elliptical hole approximates a crack and the stress at the tip becomes singular. Cracks result in high stress elevation in the neighborhood of the crack tip, which should receive particular attention since it is at that point that further 97 J ) Chapter 5 Fracture Mechanics Concepts crack growth takes place. The principal idea is that the presence of irregularities w i l l considerably reduce the strength of the material causing the structure to fail far below the expected load capacity of the structure. Figure 5.1. Stress distribution close to an elliptical hole in an oo plate Since a real material cannot sustain an infinite stress, inelastic deformation and other nonlinear effects in the neighborhood of the crack tip are expected to occur, except for the case of ideally brittle materials [85]. This zone is usually referred to as the process zone (Figure 5.2). There are, However, situations in which the extent of the nonlinear effects and the inelastic deformations are very small compared to other dimensions of the problem like crack length, crack ligament (distance from tip to closest boundary), etc., a small scale yielding condition is said to occur. This means that, despite the nonlinear effects, the elastic field surrounding the process zone governs the stability of the crack. This concept forms the basis of the usefulness of L E F M . Situations where the extent of inelastic deformation is pronounced necessitate the use of nonlinear theories and w i l l be dealt with in the next section. 98 4 Chapter 5 Fracture Mechanics Concepts m fc, , j' (•) (b) Figure 5.2. Fracture process in concrete A s we saw earlier, the first mathematical solution of a stress field in a linear infinite flat plate subjected to uniform tension and weakened by an elliptical opening which could be degenerated into a crack was provided by Inglis in 1913 [84]. However, Irwin [86] was the first to recognize the general applicability of the singular stress field in 1957 and introduced the concept of the stress intensity factor that measures the "strength" of the singular stress field. According to Irwin, there are three independent kinematic movements of the upper and lower crack lips with respect to each other as shown on Figure 5.3 Figure 5.3. The three basic modes of crack growth 99 Chapter 5 Fracture Mechanics Concepts A n y deformation of the crack surfaces can result from a superposition of these basic deformation modes, which are defined as follows CL- Opening mode, I. The crack lips move symmetrically b- Sliding mode, II. The crack surfaces slide relative to each other skew-symmetrically respect to the plane xz and symmetrically C- with respect to the planes xy and xz with with respect to the plane xy Tearing mode, III. The crack lips slide relative to each other skew-symmetrically with respect to both planes xy and xz In L E F M , the effects of the process zone are usually neglected. Thus, in the neighborhood of the crack tip a singular stress develops and, consequently, a safety or stability criterion based on stresses is inappropriate. Based on the first law of thermodynamics, an energy criterion was originally proposed by Griffith [65]. H e proposed that a condition for instability occurs when the elastic strain energy released during crack growth by a small amount is equal to the energy necessary to create the new crack surfaces. Thus, Griffith resolved the paradox arising from the Inglis solution of a sharp crack in an elastic body according to which an infinite stress occurs at the crack tip and, consequently, a body with a crack could sustain no load. For the classic, infinitely large plate containing a crack of length 2a oriented perpendicular to an applied uniform tensile stress field (Figure 5.3), the critical applied stress cr. that would cause crack propagation is given by [65,94] (5.2) where y is the surface energy per unit area and E is a function of the Y o u n g ' s modulus E, the Poisson's ratio v and the problem type: E'-E EII. = in plane strain 1 r- (lV) • / in plane stress Later, the critical potential energy release rate G c ( 5 - 3 ) has been introduced in order to incorporate not only the surface energy, but also the energy due to plastic deformations in the process zone [85, 94]. Equation 5.2 can be restated as: 100 Chapter 5 Fracture Mechanics Concepts <y.. G..E (5.4) 1 TUl The energy based stability criterion is then given by: G =G (5.5) C where the potential energy release rate G is defined more generally as: G= (5.6) ——g E' where g is a correction factor due to the geometry and loading of the structure. For Griffith's problem (Figure 5.1), g is equal to unity. The parameter G i s a measure of the c material resistance to crack propagation, usually referred to as critical fracture energy. I" - Figure 5.4. Material element near a crack tip The typical stress distribution near a crack tip when L E F M is assumed, is given by [94] K f(&)+ Negligible terms (5.7) where r is the radial distance from the crack and 9 is the angle with respect to a line parallel to the crack. The coefficient K is the stress intensity factor. The stress distribution is invariant with respect to the geometry and loading of the specimen. 101 Chapter 5 Fracture Mechanics Concepts Therefore, the stress intensity factor must incorporate these effects. It can be recognized that the stress intensity factor fully characterizes the stress field in the vicinity of the crack tip. For this reason, it can be used as a crack propagation criterion. In this case, fracture toughness is represented as a material constant K C K = K (5.8) C For the infinite plate \K - cr -Jm o K (5-9) C Comparing (5.4) and (5.8), it becomes obvious that both measures of toughness are related (5.10) Stress intensity factors can be derived in closed form for a variety of crack geometries, specimen geometries and loading. For general configurations, however, numerical methods must be used for their computation. / / '\A° / / / / / / / Figure 5.5. Elliptical crack at an angle to a uniform applied stress in an infinite medium U p to this stage, only the problem of a plane crack extending through the thickness of a flat plate was solved using the two dimensional theory of elasticity. However, many embedded cracks or flaws have irregular shapes and are three dimensional in nature, but 102 Chapter 5 Fracture Mechanics Concepts for the purpose of analysis, are usually idealized as planes of discontinuities bounded by smooth curves such as penny-shaped and elliptical embedded cracks. A great amount of effort has been spent on the determination of stress distribution in bodies with three dimensional cracks [95]. The fruit Of such studies was the expression of the local stress field near the crack front in a form analogous to the 2 D case in terms of three stress intensity factors, which are independent of the local coordinates, and depend only on the crack geometry, the form of loading and the location of the point along the crack border. This result is fundamental in analyzing the fracture behavior of cracks, and provides uniform expressions for the local stresses under various geometrical and loading conditions, where only the values of stress intensity factors differ. A great variety of stress solutions for internal and external cracks in 3D, including the effects of material anisotropy and non-homogeneity, is provided by Kassir and Sih [95]. These analytical solutions, however, are mainly concerned with bodies of infinite extent, because of the mathematical difficulties that one is faced with, when considering bodies of finite dimensions in which case, numerical and/or experimental methods are used. II.l- MIXED M O D E INTERACTION THEORIES In the previous sections, where crack growth was analyzed mainly within the framework of energy balance, it was assumed that the crack extends in a self-similar manner, which occurs only when the applied loads are directed normal to the crack plane. Unfortunately, the crack follows a curved path in general, and mixed-mode fracture can occur in the plane of the cracked specimen when a- load and crack are not symmetrically b- in the thickness direction aligned when ductile fracture modes are present such as cup and cone failure c- development of shear lips near the specimen surfaces The mixed-mode crack growth problem cannot be ignored as being academic. It is real and must be taken into account in any thorough structural design. Three mixed-mode interaction theories are described for determining crack stability and predicting direction of propagation in cases where crack growth direction is not known a priori. These 103 Chapter 5 Fracture Mechanics Concepts theories are based on the stress intensity factors and have shown good agreement with experimental results for some materials [58,81,154]. II.1.1- Maximum circumferential stress theory This theory [58] assumes that the crack propagates in the direction that maximizes the function cr^, The circumferential stress around the crack tip. The analytical solution of this problem is given by K, Sin 6 + K u (3 Cos 6 -1) = 0 (5.11) where 9 is the angle of propagation measured from the tangent to the crack path, and K,, K are values of stress intensity factors. The crack is assumed to be unstable when n K, C o s ( 0 / 2 ) - - K„ 2 where K Ic Sm(o)Cos{012)>K Ic (5.12) 1 is the critical stress intensity factor for static loading and plane stress conditions. Other versions of the stress criterion based on the maximum tangential principal stress and the maximum tangential strain have been proposed in the literature [57]. The stress criteria used in mixed-mode crack growth are inadequate for the following reasons J. The location of fracture may not always be governed by only one of the six independent stress components. The combination may play a role, 2. A stress quantity cannot be used to describe the fracture resistance of a material, 3. The contradiction that occurs for the case of a moving crack, where the normal stress parallel to the crack is greater than the stress perpendicular which case the maximum experimental stress criterion is in direct to the crack in disagreement with observations. II.1.2- Maximum potential energy release rate theory This theory states that the crack propagates in the direction that maximizes the potential energy released as the crack grows. G{0) is determined from G(G) = - l i m Aa->0 7t{a + Aa)-7t(a) Aa 104 (5.13) Chapter 5 Fracture Mechanics Concepts where 7r(a) and n{a + Aa) represent the potential energy of the body prior to and after crack extension. The critical energy release rate at crack growth is G(0 ) = max[G{e)] = l i m - l i m c n{a + Aa)-n(a) (5.14) Aa Aa->0| The above considerations lead to the following expression G(0) = ^[Kf (0)^1(0)] (5.15) where G is the shear modulus, A: is a function of the Poisson ratio v, Ik = 3 - 4 v in plane stress 3-v . , in plane strain \K = (5.16) 1 +v G(0) is the energy release rate per unit length of crack front, 0 is the angle of propagation, and K, (0) and K (0) are the stress intensity factors for the tip of the u propagation branch, in the limit as the branch vanishes, given by \i* e K,(0) = 2 f 3 Kj COS0 + -K,, 3 + Cos 0, 2 > Sin0 (5.17) \<>' " ( 2 i \K„ Cos0--K, 3 + Cos 0J 2 A Sin0 where K,, K„ are values of stress intensity factors for the tip of the main branch in the absence of the propagation branch. For the analytical computation of the crack extension angle 0, a stress analysis of the branched crack problem has first to be performed. A closed form solution for this 105 Chapter 5 Fracture Mechanics Concepts problem, possessing two different stress singularities, is very difficult to obtain even for the most simple geometry. Even i f an analytical closed form solution were to exist, the proofs for the existence and uniqueness of the strain energy release rate would be a very difficult task. The maximization of the function G(0) is performed numerically using a line search algorithm. A s there is not necessarily only one local maximum of G(o), it is important that the initial trial value be close enough to the solution to ensure convergence. The maximum circumferential stress theory is computationally simple and has proven to give adequate trial values for a sufficiently wide range of stress intensity factor ratios K„/K,. Propagation of the crack is likely to occur when the value of G(c9 ) corresponding to the direction of maximum energy release rate 0^ max G(0^) = ^Kf reaches (5.18) c II.1.3- Maximum strain energy density theory Introduced by S i h [155] in an attempt to circumvent many of the shortcomings that could not be overcome by other criteria, mainly related to the combined effect of specimen size, complex geometries and loadings, this theory states that crack extension takes place in the direction along which the strain energy density possesses a minimum value. The strain energy as a function of the propagation angle is computed by the expression S(0) = a Kf n + 2a l2 K, K +a Kl u (5.19) 12 where a u =—— [(\ + 16nG 7 COS0)(K-COS0)] (5.20) ^2=7^v[2Ow0-(*-l)] \6nG '22 = — r [ ( l - COS0)(K +1)+(l + Cos0){3 Cos0 l6nG 106 -1)] Chapter 5 Fracture Mechanics Concepts The same procedure described for the maximazation of G(#)is used for the minimization of S(d) , using results from the maximum circumferential stress theory to evaluate a first trial sufficiently close to ensure convergence. Propagation is likely to occur when the value of 5(f5 ) corresponding to the direction of minimum strain min energy density 0^ reaches S(^ h^K? n III- (5.21) c NONLINEAR FRACTURE MECHANICS When a cracked concrete structure is subjected to loading, the applied load results in an energy release rate G at the tip of the effective quasi-brittle crack that could be split up q nto two parts: the energy rate consumed during material fracturing in creating two surfaces, G Ic i- and the energy rate to overcome the cohesive pressure o~(vv) in separating the surfaces, G . C T G =G +G q Ic (5.22) <7 According to its definition, the value of G can be computed as follows a G = — t f cr(w) dx dw Aa * * (5.23) a where cr(w) is the normal cohesive pressure and w c is the crack separation corresponding to cr(vv) = 0 . In the case where the crack opening profile (shape) does not depend significantly on the crack length, the above equation becomes G„ = %'a(w)dw (5.24) Eq. (5.22) is a general energy, balance condition indicating that for quasi-brittle fracturing, the energy release rate due to the applied load G is balanced by two fracture q energy dissipation mechanisms. The Griffith-Irwin energy dissipation mechanism, 107 Chapter 5 Fracture Mechanics Concepts represented by the fracture energy release rate G lc and the Dugdale-Barenblatt energy dissipation mechanism represented by the material traction term G . a {Kobayashi et al.[96] , Jenq and Shah [87], Cook et al. [47], C o x and Marshall [49]} have proposed fracture models including the two energy dissipation mechanisms. Another alternative to this approach, is to use only a single fraction energy dissipation mechanism, either the Griffith-Irwin mechanism by assuming G - 0 , leading to the so-called fictitious lc crack models, or the Dugdale-Barenblatt mechanism by assuming CT(W)=0, leading to effective-elastic crack models or equivalent-elastic crack models. III.l- FICTITIOUS C R A C K MODELING The fictitious crack approach assumes that the energy needed to create new surfaces is negligible compared to that required to open them and that the energy dissipation during crack propagation can be completely characterized by the cohesive stress-COD relationship cr(vv). Hence (5.25) The fictitious crack is assumed to initiate and propagate when the principle tensile stress reaches the tensile strength of the material. Hillerborg et al. [78] proposed a fictitious crack model for fracture of concrete. A s can be seen on Figure 5.6, it is assumed that strain localization appears only after the maximum load is reached, whereas, the postpeak fracture behavior or softening can be characterized by a stress-COD curve. The material fracture toughness G , representing the energy absorbed per unit area of crack F and considered to be a material fracture parameter is given by the entire area under the softening stress-COD curve, cr(w). Or analytically (5.26) 108 Chapter 5 Fracture Mechanics w Concepts Slrass, a Crack opening, w (c) Figure 5.6. Principle of the fictitious crack model, after Hillerborg [78] In the fictitious crack model, the softening stress-COD curve, which is considered to be a material property, independent of structural size and geometry, is completely determined when the tensile strength / , the fracture toughness G , and the profile of ( F the cr(w) curve are known. Hillerborg [78] has combined / , and G F to obtain l , ch the so-called characteristic length, which is a purely material property proportional to the length of the fracture process zone. In concrete, the length of the fracture process zone at complete separation of the initial crack tip is of the order of 0.3 l ch to 0.5 l ch according to this model. It is important to keep in mind that the crack propagation in this model is governed by the principal tensile stress exceeding the tensile strength of the material, instead of an energy criterion. The fictitious crack model may be combined with a finite element analysis to predict the fracture behavior of concrete structures. But an appropriate selection of fracture parameters is required in order to obtain good predictions. This approach may physically make sense in the case of metals where the crack tip fracture process zone has the same fracture mechanisms as the crack wake process zone (same yielding process) on which case the original Dugdale's model was based. But in the case of concrete, many different fracture mechanisms may exist in the crack tip and the crack wake process zones, and this model cannot be more than an approximation of the fracturing process. One major problem of this model is the need to determine, experimentally, a size-independent value of the tensile strength. 109 Chapter 5 Fracture Mechanics Concepts III.2- T W O - P A R A M E T E R F R A C T U R E M O D E L It is interesting to note that this model is a member of Equivalent-Elastic crack approach family. In this alternative method, the fracture process zone of concrete is modeled through the replacement of the actual cohesive crack by an equivalent traction-free elastic crack and the single Griffith energy dissipation mechanism is used, by assuming cr(vv) = 0 . This effective crack is governed by i- a LEFM-based criterion, and ii- explicit prescription of an equivalence between the actual and the corresponding effective-crack A s a result, the energy release rate for a medium containing an effective-elastic crack which governs the propagation of such a crack due to change of applied load under opening mode conditions, is given by G =G q (5.27) IC where G i s a function of applied loads, structural size and geometry as well as the ? effective-elastic crack length. In the case of a stable crack propagation, the crack length w i l l increase as the applied load increases. Hence, an additional equation is required to determine the crack length before equation (5.27) is used. Unfortunately, the effectiveelastic crack length cannot be used directly as an independent fracture criterion because it has been shown empirically that it is dependent on the structure's size and geometry. Thus, another quantity should be introduced as a fracture criterion. This is the case of most of the effective-elastic crack models that have been proposed in the literature, where virtually all of them use two fracture parameters to define the inelastic fracture process and to govern crack propagation. It is along these lines that Jenq and Shah [88] proposed a two-parameter model based on the elastic fracture response of structures. According to this model, crack growth w i l l take place i f the following two conditions are satisfied K, = K f, ' CTOD = CTOD (5.28) , c c where K, and CTOD the stress intensity factor and the crack tip opening displacement, 110 Chapter 5 Fracture Mechanics 4 Concepts respectively, whereas, K lc and CTOD are their critical counterparts. CTOD c and K,, are determined according to the following procedure. CMOD = CMOD*+ CMOD' CTOO = CTODV CTOD" CMOD (a) CMOD' (b) Figure 5.7. Test for the determination A s we can see on Figure 5.7, the CMOD component CMOD e c of K s !c = CMOD' c The critical stress intensity factor K 5 lc CTOD c at peak load has been split up into an elastic and an inelastic component CMOD'", CMOD and so that + CMOD' " c (5.29) together with the critical effective-elastic crack length a are computed from the following system of two equations and two unknowns c after measurement of CMOD' value and the maximum stress cr. V CMOD! = o j (5.30) 4cra v u J The value of the crack tip opening displacement CTOD" is then computed CTOD" = CMOD" Note that g i g 3 b 'ac J (5.31) (i = 1,2,3), are geometric functions that can be found in an L E F M handbook [160]. Just like in the fictitious crack model by Hillerborg, Jenq and Shah Chapter 5 Fracture Mechanics <£ 7 | ^ Concepts [88] introduced a material length, Q, which is proportional to the size of the fracture process zone E CTOD ( ^ t (5.32) The material length Q can be used as a brittleness index for the material. The higher the value of <2, the more ductile is the material. In concrete, a one parameter representation of the fracture process, such as K lc or G , is misleading because one would observe F that the fracture toughness becomes proportional to the compressive strength or the strain rate, whereas, experimental evidence shows that concrete becomes more brittle as its compressive stiffness increases. According to the two-parameter model [89], the critical crack extension decreases with increasing compressive strength. Even though, the stress intensity factor evaluated at the effective crack tip, has been shown to be a good candidate as a material parameter through an extensive study initiated by R I L E M [93], the validity of the CTOD c as a material parameter has not yet been confirmed. This is mainly because of the inherent difficulty related to the measurement of such a small parameter. Since the material parameters Kj c and CTOD c are defined for the critical fracture of a structure, the effective-crack approach can only predict the maximum load and the corresponding displacement (or crack length) of a structure. In order to predict the entire stress-strain response of a structure, corresponding R-curves are usually introduced. III.3- R-CURVE BASED MODELS According to the energy principles discussed earlier regarding the Griffith balance approach of crack growth, the critical load is determined from „ 8W dU e „ which resulted from the conservation energy in the entire body. Crack growth is considered unstable when the system energy at equilibrium is maximum and unstable when it is minimum. A sufficient condition for crack stability is 112 Chapter 5 Fracture Mechanics Concepts < 0: unstable a (n+r) fracture 2 SA 2 > 0: stable fracture = 0: neutral (5.34) equilibrium where n is the potential energy of the system, defined by Tl = U -W . This stability condition can be rewritten in terms of the crack driving force G and the crack e resistance R. dJG-R) 8A > 0: unstable fracture < 0: stable fracture = 0: neutral (5.35) equilibrium For the particular case of an ideally brittle material R = 2y = cste and the R - term vanishes from E q . (5.35). Hence the stability criterion for crack growth may be expressed in terms of the stress intensity factor as follows > 0: unstable fracture < 0: stable fracture = 0: neutral (5.36) equilibrium In this case, the fracture resistance is described by the critical stress intensity factor K. Under such conditions, fracture of the material is sudden and there is virtually no crack growth before final instability. For quasi-brittle materials such as fiber reinforced cementitious composites, on the other hand, the process zone is no longer negligible and the final instability is preceded by some slow crack growth. In such circumstances, it was observed experimentally that the fracture resistance increases with increasing crack growth (Figure 5.8). The crack growth resistance curve, or R-curve method, is a one parameter method for the study of fracture in situations where small, slow, stable crack growth -usually accompanied by inelastic deformation- is observed prior to global instability. 113 Chapter 5 Fracture Mechanics Concepts R, G„ I Failure criteria: G = R q dGg dR da da^ = Quasi-brittle materials Gq-curve R-curve Perfectly brittle materials ..../. ao a , a" Figure 5.8. Typical R-curve representation The R - curve describes the changing resistance to fracture with increasing crack size. The slope of the R - curve, known as the tearing modulus, expresses the increase of fracture toughness as the crack length increases prior to instability. The theoretical basis for the R - curve, can be provided by the energy balance equation, which applies during crack growth. For situations in which the energy dissipated to plastic (inelastic) deformation U m is not negligible, the energy balance equation takes the form G =R (5.37) with dW dU e G= dA dA (5.38) dr du in n R =— + dA dA R which represents the rate of energy dissipation during crack growth is composed of two parts 1. the energy consumed in the creation of new material surfaces, 1. the energy dissipated in plastic (inelastic) deformation. During stable crack growth, E q . (5.37) and inequality (5.35) should be satisfied. Both parts of Eq. (5.37), are represented graphically on Figure 5.8 in G - a space. The points of intersection of the G and R-curve refer to stable crack growth, since E q . (5.37) and Eq. (5.35) are satisfied. Stable crack growth continues up to point P, at which the G ( « , c r . ) - c u r v e corresponding to the value cr. of the applied stress is tangent to the R- 114 Chapter 5 Fracture Mechanics Concepts curve. Beyond point P, instability occurs. Thus, the critical stress cr and the critical c crack length a c are determined from the knowledge of point P. The 7?-curve concept was first introduced by Irwin [96] using the energy criterion. Due to the suitability of this approach to the study of quasi-brittle materials, many attempts have been made recently to extend its use to describe the fracture of cementitious materials and ceramics [78,88], both theoretically (construction of the curve) and experimentally (measurement of the crack resistance as a function of the crack growth). The major difficulty of Rcurve behavior is that R- curves depend on both specimen geometry and material parameters. Consequently, the assumption that the R-curve is a material parameter for a given thickness, temperature and strain is open to question. For a given structural geometry and a material fracture property, R-curves can be interpreted as an envelope of a series of Gq-curves for certain type of structures in the following three ways 1. an R-curve may be defined as an envelope of Gq-curves for a series of structures with the same size but different initial crack lengths, or 2. an R-curve is an envelope of Gq-curves for a series of structures of different sizes but the same initial crack length, or 3. an R-curve may be defined as an envelope of Gq-curves for a series of geometrically similar structures with increasing sizes. Based on the energy dissipation mechanisms in the vicinity of the crack tip, R-curves for quasi-brittle materials can be categorized into fictitious crack approach, where a cohesive pressure is introduced on the crack surfaces, and effective-elastic crack approach, where a traction free equivalent-elastic crack is assumed, and R-curves are derived based on L E F M with at least two fracture parameters. It is interesting to note that the R-curve for a given structure is not unique due to the above mentioned different definitions. Bazant et al. [25] Proposed an R-curve based on their size effect model. They considered the maximum loads of a series of geometrically similar structures with different size. The R-curve was established by using failure criteria of the structures with different size. The R-curve obtained is interpreted as an envelope of G - curves q for a series of geometrically similar structures with increasing sizes. Ouyang et al. [133] 115 Chapter 5 Fracture Mechanics Concepts W have developed an R-curve using the effective-elastic crack approach. A governing differential equation for R-curves was derived based on the failure criteria of quasi- brittle materials. The R-curve was then obtained by solving this differential equation. This R-curve is interpreted as an envelope of the energy release rates for a series of structures with the same specimen geometry and the initial crack length but different specimen sizes. Using the fictitious crack approach, Foote et al. [61], proposed a fracture resistance curve in terms of the stress intensity factor Kr. In this model the stress intensity factor at the crack tip has been split into two parts K =K, q where K lc c K + a is the stress intensity factor to create new cracked surfaces, and K stress intensity factor needed to overcome the cohesive pressure cr(w). a P P P x Figure 5.9. Superposition Since K q of stress intensity factors at crack tip should be balanced by K , the following expression is obtained R 116 (5.39) the Chapter 5 Fracture Mechanics Concepts + (5.40) R ~ ic K = T'0 qj(a,x)a(w)dx K K <r K with a (5.41) where <p(a,x) is the appropriate Green's function that represents the stress intensity factor at the crack tip due to a unit force acting at a distance x from the crack tip. Foot et al. [61] suggested the following stress-COD curve f ... \ m (5.42) where f t and w c are the tensile strength and the critical opening displacement. The parameter m is a softening index, m = 0 corresponds to ideally brittle materials and m = oo corresponds to the case of ideally plastic materials. In order to be able to perform the integration in equation (5.36), Castiliano's theorem is used to estimate the crack profile w(x). (5.43) A K - curve is obtained by substituting Eqs. (5.41), (5.42) and (5.43) into E q . (5.40). R Since the crack opening displacement w(x) obtained from E q . (5.43) is a nonlinear function of x, a numerical iterative scheme is usually needed to obtain VV(JC). This K - curve is geometry and size dependent, since the Green's function in E q . (5.41) R and the crack opening displacement calculated from E q . (5.43) are dependent on the size and geometry of the specimen. This R-curve has been further simplified by M a i et al. [120] by introducing a Green's function for the semi-infinite plate and a certain crack opening profile. However, the simplified R-curve (Insensitive). 117 becomes geometry independent Chapter 6 Numerical Modeling of Fracture Growth NUMERICAL MODELING OF FRACTURE GROWTH BACKGROUND AND HISTORY The simplest technique for numerical simulation of a discrete crack involves symmetrical cracks [44]. Schematics of two different configurations that can be modeled with quarter symmetry are shown in Figure 6.1. The configurations shown in the figure are relatively easy to mesh, and the propagation can be represented by merely changing the boundary conditions. The obvious shortcoming of this approach is that crack propagation along a plane of symmetry is the exception, not the rule. Nonetheless, these configurations have been the benchmark problems for many innovations in computational fracture mechanics. They have been used as a vehicle for investigations in such areas as convergence studies [166], crack-tip elements [18,73,167], and techniques for evaluating stress intensity factors [19a]. Figure. 6.1. Example of crack problems that can be represented using quarter symmetry 118 Chapter 6 Numerical Modeling of Fracture Growth A more general methodology for evaluating crack propagation was clearly necessary. In the I960's, two independent approaches or school of thought, nodal release and smeared cracks, evolved. With both approaches, little or no mesh modification was necessary. The first application of crack propagation simulation for non-trivial problems was in the cracking of plain concrete and reinforced concrete using the nodal release method [128]. The nodal release approach is essentially what its name implies. A t each step, one node is released (divided into two nodes) and the crack extends along one element length. A n example of this is shown in Figure 6.2. W i t h this technique, cracks that do not follow symmetry planes can be modeled, and the crack trajectory is not necessarily constrained to remain a straight line. There are problems with this approach, however. Crack paths may become highly idealized, and considerable effort may be necessary to generate a mesh with element edges along the expected crack trajectory, which is usually not known a priori. *********** *********** ****^****** ***£ Sip***** ***jg 5***5 S*i 1* \ / &7 Steel Reinforcement Cracks Figure 6.2. A discrete crack simulation using the technique of nodal release in reinforced concrete beam (after Ngo and Scordelis [128]) At the time, most researchers considered discrete crack modeling to be too difficult for practical simulations (due to necessary mesh modifications).. In 1968, Rashid introduced a fundamentally different approach, the smeared crack method [140]. Here, the idea is quite simple. Rather than representing the crack explicitly, it is modeled implicitly 119 Chapter 6 Numerical Modeling of Fracture Growth > through modifications to the material constitutive relations for elements along the "crack" path. In short, when it is determined that the tensile stress in an element has reached a critical value, the constitutive relationship for that element is changed so that the element has little or no stiffness in the direction of tension. The presence of the crack is effectively "smeared" throughout the element. This effectively intermixes the representation of the crack with the physics governing the fracture process. These are uncoupled in the discrete crack approach SMEARED MODELING OF CRACK PROPAGATION Based on the concept of replacing the crack by a continuous medium with altered physical properties, this approach quickly became the predominate method for performing crack propagation simulations in some areas of mechanics. Not only did it eliminate the necessity to remesh at each simulation step, but with smeared cracks, it was not necessary to treat a crack tip in a special manner, (use of singular elements at the crack tip) since no tip actually exists. A l s o , for many materials (such as concrete), a band of small cracks appealed to one's intuition of how damage accumulates in front of a crack tip and eventually forms new fracture surfaces [22] (however, the physical interpretation of a smeared crack in materials with a relatively finer grain structure, such as metals and glasses, is less obvious). A n illustration of the concept underlying this method is given in Figure 6.3 using the crack band model developed by Bazant and O h [24] In this model, the crack band width, w c must be several times larger than the largest aggregate size in the concrete mix. w includes the actual crack or group of cracks. c Within this band, average properties are used based on orthotropic behavior which depends on the crack's orientation. This band width, assumed to be a material property, must be limited to ensure stability of the solution. The element stiffness matrix is obtained from 120 Chapter 6 Numerical Modeling of Fracture Growth Cracking within the element can be represented through either one of the two ways: Anisotropic model, in which we modify D, Rashid [140] Discontinuous shape functions, in which we modify B, Droz [54] t t (c) (d) Figure 6.3. Random microstructure, scatter of microstress, and crack band or sharp crack models 121 Chapter 6 Numerical Modeling of Fracture Growth In the anisotropic model, assuming we have a crack aligned with the y axis, we maintain B fixed and alter D such that K 0 0 0 0 E 0 0 0 a E 2(1 + v) As can be seen, the normal stresses can only be transmitted across the y axis, and the shear coefficient has been reduced by a factor a which is a function of aggregate interlock. The issue of proper determination of a has been the subject of numerous controversies [40]. Although Rashid's model has been widely used in numerous finite element programs, it has been recently shown by Droz [54] that certain undesirable coupling remains within the stiffness matrix. This coupling causes perturbation of the shear stresses of the element and is unable to correctly model structural response under mixed mode loading cases. It has been numerically shown that this perturbation increases with mesh density. In the new alternative, the matrix D is kept constant whereas B is modified. The idea in this approach is to break any connection between two adjacent nodes across the crack. Among the advantages of the smeared crack approach one can mention [82]: No remeshing is required, making it very convenient from the computational point of view The crack generally is not straight but tortuous. Such tortuosity is readily approximated using the smeared crack approach Distributed damage and cracking have been observed in reinforced concrete structures In the case where parallel cracks in concrete are densely and uniformly distributed, such cracks are well represented by the smeared crack model with a minimum crack band width equal to the actual spacing of the parallel 122 cracks. Chapter 6 Numerical Modeling of Fracture Growth Despite its popularity, there are three significant difficulties inherent in using smeared cracks. The first is the lack of an unambiguous description of a crack's geometry and a crack openning profile. Without explicit geometry it is not possible to model crack propagation along the interfaces between regions of differing materials [see chapter 5]. Also, for some types of fracture problems (hydrofracture for oil well stimulation, for example) cracks are driven by fluid pressure within the crack. The local fluid pressure is a function of the local crack opening. Without an expilicit crack opening profile, these types of analyses are not possible. The second difficulty with smeared cracks is that, if used naively, the results are highly sensitive to the finite element mesh used [22, 23]. What is worse, not only are results not objective, but in the limit as the elements in the crack tip region become small, the solution diverges in a manner such that the load necessary for crack propagation approaches zero. This lack of convergence is counter to one's intuitive understanding of how discrete methods should approach a non-trivial solution in the limit, as the mesh is refined. It has, however, been shown that the divergent behavior can be prevented if care is taken when formulating the constitutive relationships for the cracked elements, and when certain other (sometimes stringent) restrictions are observed when meshing [23]. The third difficulty with smeared cracks has been reported quite recently [142]. This has to do with a fundamental kinematic incompatibility inadvertently introduced by the smeared-crack softening when the direction of propagation is not parallel to the element edges. In this case, significant stresses are observed in some of the elements along the crack face (where stresses should be very small). Rots [142], calls this phenomenon "stress-locking". The overall displacements in the model tend to be accurate, but stresses and strain-energies computed within the finite elements will not be accurate in the region of the crack. This is an artifact of linking crack representation with the physics controlling the fracture process. Stress locking can be illustrated with the following Figures [142]. Figure 6.4 shows a test configuration analyzed by Rots. The physical experiment was performed by Kobayashi et 123 Chapter 6 Numerical Modeling of Fracture Growth al. [97]. Two detailed views of the upper right portion of the plate are shown. The left portion of Figure 6.4 shows the deformed mesh, while the right portion illustrates vectors of principal stress. Regions of "locked-in" stresses are clearly evident along the crack path. This material should be essentially unloaded. These three problems may be more or less significant, depending on one's purpose in performing a simulation and the nature and size of the object being simulated. There is one class of problems for which smeared cracks are the obvious choice. These are problems where a large number of cracks are present in a single structure. In this case, it would be difficult to model all the individual cracks explicitly. Figure 6.4. Test configuration analyzed by Rots together with a representation of the deformed mesh and principal stresses 124 Chapter 6 Numerical Modeling of Fracture Growth DISCRETE MODELING OF CRACK PROPAGATION Initial attempts at discrete crack propagation simulation were made over twenty years ago. Since then, techniques have improved significantly, but such simulations are still performed routinely by only a small number of practicing engineers. This is likely to change in the coming years, as propagation simulations become an integral part of the total life cycle (design, analysis, and maintenance) of many products. The factors expected to drive the increase in popularity of these methods include the following: • The computer processing and graphics performance necessary to perform these simulations is increasingly becoming available on the desks of most practicing engineers, due to the continued exponential increase in computer processing power. • New popular materials, especially composites and ceramics, are more susceptible to failure by brittle fracture than conventional plastic and metallic components. The increased use, and subsequent failure, of these materials will continue to increase interest in discrete crack propagation simulations, which are particularly well suited these types of materials. • Software for crack propagation simulation continues to improve. This is primarily due to new techniques ad algorithms in the areas of computational mechanics and solid modeling. Because crack propagation simulations are currently not widely used, some readers may be unfamiliar with the topic. Part one is included to provide the necessary background information. The specific focus of the chapter is three-fold. First, it serves as an introduction for readers unfamiliar with the topic. Second, it contains a review of the applicable literature. Third, it introduces simulation strategies for crack propagation. 125 Chapter 6 Numerical Modeling of Fracture Growth The purpose of this section is to introduce the topic of discrete modeling of crack propagation, and to present a brief history. There are two defining characteristics of discrete crack modeling: Cracks are represented explicitly and extended discretely. That is, an explicit representation of the geometry of the crack (or cracks) is always maintained. A s a crack grows in an object, the created surfaces are explicitly incorporated into the model. In addition, when cracks are also extended, it is done in a discrete manner: crack increments have finite length, and represent a piecewise, rather than continuous, model of crack evolution. This approach is necessary when performing a propagation simulation with finite simulation-resources and is consistent with other types of discrete mechanics. Discrete modeling of crack propagation is an approximate technique similar to other approximate techniques used for computational mechanics (eg. Finite elements, boundary elements, finite differences, boundary collocation). A s with most approximate techniques, there is a strong notion of convergence towards a "correct" solution as the numerical model is refined and the computational cost is increased. W i t h discrete crack modeling, there are two principal types of approximation in addition to the usual assumptions about idealized material models, geometry, and boundary conditions. First, discrete crack modeling relies on an underlying approximate technique for performing stress analysis (usually the finite or boundary element method). Second, the continuous evolution of a fracture surface is idealized by discrete extensions. A discrete crack simulation w i l l converge to an asymptotic solution only i f both types of approximations are reduced. That is, both more (or better) elements must be used for the stress analysis, and shorter crack increments (i.e., more simulation steps) must be taken. Discrete modeling of crack propagation is an iterative process with four major tasks in each simulation step. A flow chart of the process is shown in Figure 6.5. The input to the process is the initial geometry of an object to be idealized. The first task in the computational procedure is to perform a stress analysis. The purpose of the stress analysis step is to compute the response of the object (stresses, strains, and displacements) to the imposed boundary conditions (loads and restrains). 126 Chapter 6 Numerical Modeling of Fracture Growth Figure 6.5. Flow chart of the iterative scheme 127 Chapter 6 Numerical Modeling of Fracture Growth The second task in the simulation loop is the fracture analysis. This is the evaluation of the parameters that govern the fracture process. The most commonly used parameters are stress intensity factors, energy release rates, J values, crack tip opening displacements ( C T O D ' s ) , and crack mouth opening displacements (CMOD's). It is with these parameters that predictions about the nature of the crack's growth is made. The third task in the simulation process is to perform an assessment to determine whether a condition of failure or arrest has occurred. This is the test that allows an exit from the simulation loop. The notion of crack arrest is only meaningful in the context of unstable crack propagation. It is defined as the point at which the crack-tip velocity drops to zero. The definition of failure, however, is more subtle. The actual definition of failure w i l l vary according to the type of simulation being performed and the ultimate objective of the simulation process. Typical failure criteria are the onset of dynamic crack propagation, a crack tip reaching a free surface, or a critical value of stress intensity factor or crack opening being obtained. In most cases, once failure or arrest has been attained, the simulation process is complete. If failure or arrest is not attained, the fourth task in the simulation process is the finite extension of the crack or cracks. This implies three sub-tasks. First, the fracture parameters, evaluated in step two, are used to predict the new crack direction and length. Second, the geometry of the object is explicitly modified to reflect the new crack configuration. Third, i f necessary, the object is remeshed so that a new stress analysis can be performed. The simulation loop is now complete, and begins again with a new stress analysis. There are three general classes of crack propagation: sub-critical crack growth, quasistatic crack growth, and unstable crack growth. The distinction among the three is the ratio of the generalized forces driving and resisting crack propagation. 128 Chapter 6 Numerical Modeling of Fracture Growth In the case of sub-critical crack growth (for example, fatigue or stress corrosion cracking), the driving forces are less than the toughness of the material (or resistance to cracking). However, cracks w i l l still grow due to the degradation of the material at the crack tip caused by corrosion or repeated loading and unloading. In general, sub-critical crack growth takes place over relatively long periods of time. Crack-growth rates are typically measured in terms of days, months, or even years. Quasi-static crack propagation takes place when three conditions are met. First, the force driving the crack must exactly equal the material's resistance to cracking. Second, the rate of change of the driving forces with respect to the crack length must be negative. Third, the rate of crack growth must be slow enough that inertial forces can be neglected. This situation is attained for some configuration when the loads or displacements applied to the object are monotonically increased at a relatively slow rate. These events are typically measured in the terms of seconds, but may occur over much longer periods of time (as in cracking driven by swelling or long term temperature changes). The final crack propagation is unstable, or dynamic, crack propagation. In this case, the crack driving force is greater than the toughness of the material. The excess of energy, which cannot be absorbed by the fracture process, is converted to other forms such as kinetic and thermal energy. These types of propagation take place very rapidly (typically measured in milliseconds), and inertial effects must be considered. From the point of view of a numerical simulation, there are many similarities among the three classes of propagation. In particular, one set of data structures and algorithms can be used to represent and modify the geometry, meshes, and boundary conditions for all three types of propagation. These data structures and their corresponding algorithms w i l l be called the representation of an object. What differs among the classes is the strategy that is used to dedicate the nature of the modifications to the geometry, meshes, and boundary conditions, during crack propagation. 129 Chapter 6 Numerical Modeling REMESHING of Fracture DURING Growth CRACK PROPAGATION In the early 1980's, the first attempt to model crack propagation without regards to an existing finite element mesh was developed by Saouma and Ingraffea [146, 147]. Here, a very pragmatic approach was taken, splitting elements in the crack path into a number of smaller elements. This approach left most of the mesh intact, making only local modifications in the region. The remeshing process then became one of identifying the particular topology present, and using the appropriate template to remesh the problem. The remeshing process is illustrated for three different crack configurations in Figure 6.6. Using this procedure, arbitrary curvilinear crack propagation was modeled and the remeshing process was automated. The major disadvantage of this method, however, is related to that the resulting meshes are often of poor quality. When poor aspect ratio elements are generated at the crack tip where stress gradients are very high, the accuracy of the near crack-tip solution becomes questionable. Figure 6.6. Three different crack configurations using the splitting element approach 130 Chapter 6 Numerical Modeling of Fracture Growth A more sophisticated approach has been introduced by Swenson [158,159]. Here, the strategy is: to delete a group of elements in a certain neighborhood of the crack tip the crack is extended into this area well shaped elements are introduced around the crack tip create a "transition" mesh that will fill the area between the crack-tip elements and the original mesh With the development of the delete-and-fill approach to remeshing, and robust implementations of the associated algorithms, it is finally possible to analyze arbitrary curvilinear crack propagation for planar and axisymmetric problems when L E F M is applicable [83]. As we saw in chapter 5, the introduction of stress intensity factors for the first time gave engineers an effective tool for investigating the fracture of structural components. Unfortunately, prior to development of numerical methods, this capability was severely limited by the relative scarcity of stress-intensity factor solutions which need to be determined separately for every different geometry and loading configuration, and were only available for a small number, of often, idealized configurations. Handbooks of the known stress intensity factor solutions were compiled. Sometimes, however, the engineer finds himself faced with a problem with no known solution and is forced to use a known solution that has the closest geometry and loading to a known solution (or to use a good dose of judgement to make appropriate modifications). INTERFACE ELEMENTS As seen in chapter 5, the assumption of singular stresses at the crack tip is mathematically correct only within the framework of L E F M , but physically unrealistic. In cementitious materials, a fracture process zone develops ahead of the crack tip. The most popular model simulating this behavior is Hillerborg's fictitious crack model ( F C M ) . 131 Chapter 6 Numerical Modeling of Fracture Growth Due to their ability to model the behavior of geometrical discontinuities in multi-phase materials, interface elements play a very important role in computational mechanics. T w o finite element formulations of interface elements can be made: the nodal (lumped) interface formulation and the continuous (integrated) interface formulation [147]. In order to model the behavior of the stress transfer through a discontinuity, such as cohesive cracks in F R C materials or substrate-repair material in structural repairs, use of interface elements can be invoked because of their ability to represent and model both geometrical discontinuities and material nonlinearities [10]. Even though the nodal interface formulation (lumped) is equivalent to the continuous interface formulation (integrated), they lead to different element stiffness matrices. Figure 6.7. Interface idealization and notations 132 Chapter 6 Numerical Modeling of Fracture Growth Lumped Interface Elements This is a generalization of the nodal approach where rudimentary nodal interfaces are introduced. In the 2D case, this latter method, consists in the introduction of two-noded spring like elements with two degrees of freedom per node. r - - ~ - Figure 6.8. 2D nodal interface elements The element nodal displacement vector is given by v = {v> ,v/,v } 2 (6.1) 2 n t where superscripts indicate element nodes and subscripts n and t represent the normal and tangent to the interface respectively (Figure 6.8). The relative displacement vector is related to the nodal displacements through a matrix B \Au Au = • I Au, n = Bv with B = -1 +1 0 0 0 0 -1 +1 (6.2) The traction-relative displacement constitutive relation is given by a matrix D defined as follows t = = DAu = d n 0 133 0 \Au d, I Au, n (6.3) Chapter 6 Numerical Modeling of Fracture Growth The linear elastic stiffness matrix K is obtained through the principle of minimum total potential energy by setting STI = 0, where p n=U+V = - f Au tdA T vf (6.4) T — is the total potential energy, U is the internal energy, V the potential of external work, and / i s the external force vector. This leads to Kv-f with K=\B DB (6.5) T JA In the case of a 2 D problem, the above equation yields ~d 0 0 ~d d 0 0 0 0 d, -d, 0 0 ~d, d, n K =A n n (6.7) A represents a contribution factor that depends on the type of interpolation and the dimension of the surrounding elements [150]. A" 4 3 4 o -6 O1 1 4 2 -o- 6 -Q -6 3 Quadratic Linear Figure 6.9. Linear and quadratic 2D interface elements A generalization of the nodal interface is provided by the lumped interface formulation where continuous elements are introduced instead o f interface points. In this case, continuous two dimensional elements are used in planar stress or strain or axisymmetric conditions. A n illustration of an r-noded line interface where each node has two degrees 134 Chapter 6 Numerical Modeling of Fracture Growth of freedom is given by Figure 6.9 for the case r=4 and r=6. The element nodal displacement vector is given by v = {v^v„ ,...,v„ ; v / , v , , . . . , v ; } 2 r (6.8) 2 The traction-relative displacements relation is evaluated at the node sets and the matrix B is similar to the one defined in the case of nodal interface elements. Consequently, in the case of a non distorted six-noded line element, B is given by B = -1 0 0 0 0 0 +1 0 0 0 0 0 0 0 0 0 0 -1 0 0 +1 0 0 0 (6.9) A summation over the contribution of each node set is performed i n order to obtain the interface stiffness matrix K and the equivalent nodal loads / in the following manner 6 K. 0 (6.10) 1=1 with d, 0 0 ~d 0 0 0 4d, 0 0 -4d, 0 0 0 d 0 0 -d, -d, 0 0 d, 0 0 0 -4d, 0 0 4d 0 0 0 -d 0 0 d, t ; i { (i — n,t) where A, is the surface nodal contribution factor and nns is the number o f node-sets. Continuous Integrated Interface Elements This formulation is more sophisticated than the lumped formulation. Instead of using relative displacements at isolated node sets, like in the previous case, the continuous integrated interface formulation makes use of interpolated relative displacements at integration points and therefore, interpolation functions are introduced like in the case of 135 Chapter 6 Numerical Modeling of Fracture Growth continuum elements. Since the numerical integration automatically takes into account the correct contributions, the introduction of surface nodal contributions is no longer needed. A n r-noded element nodal displacement vector could be defined as before v = {v' ,v ,...,v ; 2 n r n n v/,v, ,...,v } 2 r (6.11) ( A continuous displacement field vector is defined for the lower and upper edges of the interface under the following form u = {u" ,u' ,u?,u } (6.12) l n n t / and u denote the lower and upper edges of the interface, respectively. A matrix H containing the interpolation polynomial functions defines the relationship between the displacement field vector u and the element nodal vector v as follows U =H v where n = {N ,N ,...,N }. 1 1 r/2 with H = n 0 0 0 0 n 0 0 0 0 n 0 0 0 0 n (6.13) r/2 denotes the number of pair of nodes (node sets). The continuous displacement field u is related to the relative displacement Au through a matrix L defined by Au-Lu where L= ~-l 0 + 1 0 (6.14) 0~ 0 - 1 + 1 The relative displacements are related to the nodal displacements through the following relation Au=Bv with B-LH. (6.15) A n example of a matrix B is given below i n the case of non distorted interface line element 136 Chapter 6 Numerical Modeling B = of Fracture Growth -n n 0 0 0 0 -n n -Nj -N 2 0 0 where n = {N ,N ,N ). l 2 3 (6.16) -N 3 0 Nj N N 0 0 0 0 0 0 0 0 0 -Nj -N -N NJ N N 2 3 2 3 2 3 For line interfaces, the element stiffness matrix K and the equivalent nodal f o r c e s / c a n be computed as follows K=b 1 \ B DB —— T *S (6.17) where b is the out of plane thickness of the specimen in plane strain or plane stress applications. The traction-relative displacement relations need to be evaluated at the location of integration points. In the case of a distorted interface, a transformation to the local reference is necessary in both formulations. The continuous numerically integrated interface formulation gives exactly the same results as the lumped formulation i f the integration points are chosen so that they coincide with the node-set points. However, if another integration scheme is used, the results would be different in the two formulations. In this work, emphasis is on 2 D (line) elements and the continuous (or numerically integrated) interface formulation has been adopted. Linear and nonlinear constitutive laws have been used. The constitutive models retained for the interface elements are the Normal traction-Displacement models. DYNAMIC RELAXATION SOLVER Dealing with nonlinear interface models requires a reliable algorithm to determine accurate solutions. The dynamic relaxation solver has been selected for this purpose because of its simplicity, ability to cope with highly nonlinear material and geomaterial behavior [169]. Intended for the static solution of structural mechanic problems, the 137 Chapter 6 Numerical Modeling of Fracture Growth method is based on the key idea that the static solution is actually the steady state part of a transient response to an instantaneous load step. Finite Element Formulation Discretization of the balance equations by the finite element method leads to a vector equation (system of scalar equations), which i n the case of a linear behavior relate the internal f o r c e s / t o the external forces F through a linear stiffness matrix K f(v) = Kv=F (6.18) where v is the vector o f nodal displacements. For nonlinear problems, on the other hand, the relationship between the internal forces and the displacements is of an incremental nature Af(v) = K (v)Av T where K T (6.19) is the tangent stiffness matrix. Another alternative for treating nonlinear problems, consists in the introduction of a pseudo-force S(v) to account for the nonlinear behavior f(v) = Kv + S(v) (6.20) A s we saw earlier with the finite element formulation of interface elements, nonlinearities are concentrated along the interface elements which contain a small number of degrees of freedom and that tractions derived from the constitutive models can be transformed into equivalent nodal internal forces. The resulting vector form is very suitable for a nonlinear dynamic analysis by an explicit method such as the one presented here. The motion equation that governs the dynamic response of a structure is given by Mv + Cv + f(v) = F(t) 138 (6.21) Chapter 6 Numerical Modeling of Fracture Growth where M is the mass matrix, C the damping matrix, F the external excitation, and t is time. The transient solution is obtained through a direct integration approach using the central difference approximation formulae for the time derivatives •j-rn (-v'-'+v') h h V where / denotes the i' (6.22) time increment, and h a constant h time increment. A n approximation of the velocity at the i' time increment is given by h v' =-(v ~ i 1/2 +v i+1/2 ) (6.23) A formulae for computing the next velocity and displacement vectors is obtained by simple introduction of equations (6.22) and (6.23) into (6.21). 1 ' --C _~h~ 2 1 ' M + -C 2 . h M [F- -f M — C + (6.24) — 2. |v =v ' + hv ,+7 where F' =F(t') and / ' = f(v'). i+112 , Since our physical problem is only a static one, the mass M and damping C matrices are fictitious and therefore are chosen to be diagonal with C = c M. This leads to ^ { £ z £ J ! l ^ = + 2 [2 + ch] v ul A M - [ f ' - / ; ] [2 + ch] = '+hv i+1/2 v 139 (6.25) Chapter 6 Numerical Modeling where M 1 of Fracture Growth denotes the inverse of M. It should be noticed that this assumption preserves the form of the central difference integration formulae. In order to use the algorithm, v~ l! 2 must be defined. Introduction of equation (6.21) together with the following initial conditions \v°*0 v°=0 (6.26) leads to v ' =hM 1 2 (6.27) A summary of this explicit central difference integration algorithm is v 1/2 ifi = 0 \v> =hM =v°+hv 112 (6.28) •MI2 _ =j V [ - \.:i-i/2 2 2 ch i v'-" +2hM~' [2 + ch] J£—Cr- [2 + ch] ifi*0 v M = +hv " i i+ 2 v Stability Considerations A s mentioned earlier, the D R solver generates a fictitious dynamic problem in order to compute the static solution of the physical problem which corresponds to the steady state part of the dynamic solution. Thus, the dynamic properties of the problem (the damping c, time increment h, and the mass matrix M), can be selected so that fast convergence to static equilibrium is achieved. A t every step, the condition f(v) = F is checked to determine whether static equilibrium has been attained, where only F and f(v) are related to the physical problem. 140 Chapter 6 Numerical Modeling of Fracture Growth The specification of appropriate values for the mass matrix elements m, u damping coefficient c, and time increment h determines whether there is convergence to the steady state of the dynamic solution and how fast it is attained. The explicit scheme presented above is conditionally stable. Thus, the time increment must be chosen so that stability conditions are satisfied. The time increment should represent a fraction of the lowest period. The time increment h is set to unity. Use of Gershgorin's theorem, provides a lower bound of the mass values [169] that guarantees numerical stability n^l Y\ ij\ m h2 K (- > 6 ^ where K tj 29 j denote the elements of the stiffness matrix K, and h is the time step. The mass matrix and time increment are chosen having in mind convergence as the main goal, whereas the damping coefficient is selected so that an almost critical damping behavior provides fast convergence to the steady state part of the dynamic solution. The damping factor can be adjusted at every iteration step from Rayleigh's quotient | (v f f l c,=2.\ ' J , (6.30) where superscript T denotes the transpose COHESIVE CRACK ANALYSIS A s seen before, L E F M is limited to the small scale yielding condition; the process zone surrounding the crack tip is small with respect to the zone of K dominance and any other geometric dimension. For some materials and specimen configurations, this condition is impossible. Therefore, nonlinear aspects of the fracturing process need to be examined. The goal of this section is to provide an introduction to such aspects. A more detailed description of the topics discussed here can be found in Kanninen [94]. 141 Chapter 6 Numerical Modeling of Fracture Growth Strain Softening Zone Strain Crack Opening ! Figure 6.10 Schematic representation of stress-strain associated with fracture behavior of strain behavior softening materials A perfectly brittle material follows a straight stress-strain path until fracture is achieved. However, this type of behavior only occurs for a few materials associated with some special geometries. In general, materials behave as shown schematically in Figure 6.10. Because of the high stress concentrations ahead of the crack, the material around the crack tips is highly deformed. Three different zones can be identified around the tip (Figure 4.1). Zone I is the linear elastic zone or the K-dominant zone. The stress at a distance from the tip is low enough to allow the structure to behave elastically. Zone II 142 Chapter 6 Numerical Modeling characterizes of Fracture Growth a region where inelastic effects take place with stresses increasing nonlinearly with increasing strains. Finally, zone HI represents the fracture zone where unloading occurs with strain localization or, in other words, the stress decreases as the strain increases. The process zone, defined as a region where nonlinear effects take place, comprises zones II and IH. The appropriate approach to study the fracturing process depends on the importance of one zone relative to the others. Thus, i f the inelastic and fracture zones are small compared to k-dominant region and the specimen dimensions, including the crack length, L E F M is used to model crack propagation and stability. This situation is usually referred to as brittle fracture behavior. However, i f the inelastic zone is large relative to other dimensions, elastoplastic fracture mechanics approaches are necessary. In yielding materials, like most metals, crack initiation is usually followed by stable crack growth, high plastic deformation occurs with blunting of the crack tip. This type of phenomenon is called ductile fracture. Due to blunting effects the order of strain singularity at the crack tip may change considerably. The exact order of singularity depends on the hardening parameters, reaching sometimes an order of 1/r. In other types of materials like concrete and other ceramics, zone II is characterized by extensive micro-cracking. Actually for these materials, the influence of the fracture zone (zone III) is much more pronounced. For a matter of distinction, this type of behavior is denoted here by cohesive fracture. In this circumstance, models capable of representing the effect of the fracture zone (zone U l ) on the fracture process must be used. The fictitious or cohesive crack is an example of such a model. In the case of concrete and cement-based composites, a cohesive crack model ( C C M ) is usually employed to take into account the softening strain localization phenomenon associated with this dissipation mechanism. A model, proposed by Hillerborg et al [77] considers the energy dissipation in the wake of the crack path, while neglecting the other two toughening mechanisms. C C M is the state of the art tool to analyze the fracture behavior of quasi-brittle materials, including concrete, cement-based composites, rocks, ceramics, and ceramic composites. 143 Chapter 6 Numerical Modeling 0k of Fracture Growth The applicability o f C C M to general mixed-mode conditions has been the focus of recent attention. Throughout the 1980's, Ingraffea and coworkers have employed a strategy of singularity near-cancellation to handle mixed-mode crack propagation, with optional load or crack length control. Crack stability was controlled by iterating on the load or crack length until residual values of stress intensity factors were obtained at the tip of the fictitious crack. Direction of crack propagation was computed using these same residual values. This criterion sometimes failed to predict accurately the critical load to produce crack extension because the singular elements at the crack tip perturbed the stress field in the cohesive zone. Moreover, use o f these elements was inconsistent with the assumption of singularity cancellation inherent in the C C M . Recently, an extension of the fictitious crack model to mixed mode propagation has been proposed by Bocca et al [36]. They employed a crack length control scheme with a stress based crack stability criterion. However, only linear softening behavior is possible in this approach. The Cohesive crack problem The principle of virtual work is used as the integral statement of equilibrium in the formulation of the cohesive crack problem: (6.31) Where S is the geometric volume occupied by the body, dS defines the external surface of the body, a the stress tensor, Se is the incremental virtual strain tensor, b is the body force vector acting per unit volume, Su is the incremental virtual displacement vector, and p is the vector of applied tractions. Expression (6.31) is the weak form of the equilibrium equations and is valid for any stress-strain constitutive law. The cohesive crack model assumes that the process zone can be represented by closing tractions pz (Figure 6.11), acting on both crack faces (pz + 144 = - pz~). A t the tip of the Chapter 6 Numerical Modeling of Fracture Growth process zone, the normal component of pc equals the ultimate tensile strength of the material. Noting that S ac =S a+ =S One may rewrite (6.31) as (6.32) From (6.32) it is clear that the problem can be treated by superimposing the effects of the body forces, the applied external loads and the cohesive loads when the body volume behaves elastically. The cohesive problem can be stated as: a- for a given load factor A, find the process zone size so that the ultimate tensile strength is achieved at the tip of the process zone; or alternatively, b- for a given process zone size, find the load factor X that would allow the ultimate tensile strength to be achieved at the tip of the process zone. Approach (b) is used in this thesis. Due to the arbitrary nonlinear nature of the softening along 8S (Figure 6.11), a numerical technique is usually employed to solve the cohesive C crack problem. A finite element or boundary element approach is most frequently used to solve the equations of equilibrium. u Figure 6.11. Cohesive crack problem 145 4 Chapter 6 Numerical Modeling of Fracture Growth The Influence method The fictitious crack model has to be treated numerically in general. Both the finite element method and the boundary element method can be applied to solve the equilibrium equations. However, due to the nature of the softening model that involves tractions and relative displacements, the boundary element method seems to be a more natural way to approach the problem. Petersson [134], however, used the finite element method to define his procedure called here the influence method. In fact, the influence method is independent of the numerical technique used to solve the equations of equilibrium. Crack displacements are considered through superposition. The contribution of each node along the crack path is taken individually as well as any contribution of the applied loads. The complete crack path is considered initially (Figure 6.12). This is a strong limitation since the crack trajectory needs to be known a priori. General boundary conditions are assumed with proportional loading. f2, u2 Figure 6.12 Discretized cohesive crack problem 146 Chapter 6 Numerical Modeling of Fracture Growth The opening of each node set along the crack is computed for applied opposite tractions or forces (depending on the numerical technique involved) at each node set on the crack path. The opening of each node set for the applied loads is also computed: w = n :p + h F where w(n) is the vector of crack openings, Q{n:n) (6.33) is the influence matrix of crack tractions or forces, p(n) is the crack traction or force vector associated with the cohesive model, h(n) is the influence vector of the applied loads, and F is a load factor (scalar). A n y response of the body can be obtained by superposition of the influence quantities. For example, the displacement d is obtained as d =d.p + d F F where d is the influence vector of the cohesive tractions or forces, and d (6.34) F is the displacement due to the applied loads. Proportional loading is controlled by the factor F. The main advantage of this method is that the equilibrium problem is decoupled from the crack problem. One can model propagation and compute load factor associated to each crack configuration simply by moving the fictitious crack tip along the defined crack path. Computation of the load factor must satisfy a stability criterion for the fictitious crack tip. Assuming that an initial notch is introduced from node 1 to k (Figure 6.12), node k defines then the true crack tip. The fictitious crack tip is then allowed to go from node k+1 to node n. If one moves the fictitious tip to node i while the true tip is still defined by node k, the problem can be solved i f the cohesive pressure acting on the crack surface is known. The objective is to find the load factor F, the traction or force and the opening at each crack node associated with the specific fictitious crack position. The boundary element method is used here to illustrate the solution of the problem. The same ideas could be 147 Chapter 6 Numerical Modeling of Fracture Growth formulated in an equivalent finite element approach i n a similar manner. If the boundary conditions are introduced, equation (6.33) could be modified as follows Pi.k = Pi.k = sK Pk i.j-i + p, = f t Pi+i... a n r e a + n d Wj k is known /..,w) defined by the cohesive and w,. = 0 unknowns and pressure ( 6 ' 3 5 ) w _, „ =0 tJ where g is the function relating the closing pressure on the crack lips to the crack opening. W i t h the introduction of crack boundary conditions, the problem unknowns are reduced to n the openings of the initial notch (l,...,k), the openings at the cohesive zone (k+l,...,i-l), the tractions of the elastic zone, the load factor F. A nonlinear system of equations is usually generated due to the nonlinear nature of the function g. A Newton-Raphson or a modified Newton-Raphson procedure could be used to solve such a nonlinear system. The traction vector p, the opening vector w and the load factor F are then computed for the specified fictitious tip position. A n arbitrary response is retrieved by superposition of p and F. The displacement d for example, is computed through expression (6.34). A s the fictitious crack tip advances, the specimen's response can be recorded. This is a very useful procedure, but unfortunately cannot handle arbitrary crack propagation. The crack path has to be known a priori. It is interesting to notice that the loading may be controlled by imposing surface tractions or forces as well as prescribed displacements. Defined Crack Path Strategy Interface elements provide an alternative solution for the fictitious crack problem. Since a stress-opening displacement constitutive model is associated to interface elements, they hence can be incorporated along a known crack path and use the dynamic relaxation 148 Chapter 6 Numerical Modeling of Fracture Growth method to solve the cohesive crack problem. In this procedure, the fracture process zone is no longer controlled by the user, it is closely related to the constitutive behavior of the interface elements. The location of the fictitious crack tip may be estimated from the opening profile. n o r m a l s t r e s s > Crack Opening Displacement (COD) Figure 6.13. Example of a softening model of interface elements force excessive force displacement control displacement Figure 6.14. Typical force displacement responses 149 Chapter 6 Numerical Modeling of Fracture Growth The compressive stiffness of the interface model (Figure 6.13) provides continuity for the constitutive model. This property is crucial for the Dynamic Relaxation solver. Furthermore, overlapping of the crack faces may be prevented by using a very large value for k . The elastic behavior ahead of the fictitious crack tip can be handled for small c values of crack opening. The load capacity of the specimen or structural member should not be exceeded i f a load control procedure is used to deform the specimen. In this case, the solution is computed under force control and only the increasing load path is recovered because equilibrium may be impossible (see Figure 6.14). Therefore, it is preferable to use the alternative solution consisting in loading the specimen through a displacement control procedure. This procedure is known to handle fairly well the force softening with increasing displacements (Snap-through). Moreover, this procedure can allow detection of a snap back condition. A sudden jump in the force for a small displacement indicates the presence of Snap-back. Arbitrary Cohesive Crack propagation A more general alternative in modeling the growth of arbitrary cohesive cracks is given by the A C C method [35]. This approach combines fundamental capabilities for mesh and geometry representation, stress analysis and nonlinear fracture mechanics into a single procedure where interactive computer graphics and sophisticated data structures to visualize and control geometry evolution are used. It is well accepted that the load-displacement response of quasi-brittle structures is size dependent. Strain softening, as well as snap-back instabilities, may occur for the same material depending on the structure sizes. Whereas softening can be handled with a displacement control procedure, snap-back cases are more difficult to follow. Snap-back branches can be numerically captured [36] i f the loading process is controlled by an increasing function of the fracture process zone (FPZ) length. This procedure whish assures uniqueness of the solution, can handle both snap-through and snap-back instabilities. Traditional 2 D elasticity finite elements are used in conjunction with zero thickness interface elements in the fictitious crack, the behavior of which is usually 150 Chapter 6 Numerical Modeling of Fracture Growth specified as a nonlinear stress-displacement law. The dynamic relaxation ( D R ) solver has been shown to give good results in the case of highly nonlinear problems such as softening and snap-back [169]. The D R although robust, may be slow in some situations such as the case where finite elements with large aspect ratios are used in the mesh, but on the other hand this scheme has a great potential for implementation in parallel computations. In the A C C approach, the aim is to determine a loading condition that satisfies both global equilibrium and a crack stability criterion for a given fictitious crack length. A criterion based on the fictitious crack tip opening profile is used by A C C scheme. N o singular elements are used at the fictitious crack tip, nonlinear analyses are performed until a condition of zero slope of the C O D at the tip is reached. If linear strain isoparametric elements are used, the condition of zero-slope can be written dw (4w -w ) B A ~ds~i=0 I = 0, or (6.36) 4 where w is the opening of the fictitious crack, and A and B are nodes of the element adjacent to the crack tip. A n y one of the stability criteria that have been widely used in the literature, such as the tensile strength [63c], the fictitious stress singularity cancellation [147, 158] could be used in this approach. Monotonic and proportional loading is assumed in the computation of the load factor that satisfies the stability criterion for the current fictitious crack length. A n overestimate of the loading factor is initially assumed, so that the D R solver can be used to generate an equilibrated solution for that load factor. A n important issue in the load factor computation is the initial configuration used in the dynamic relaxation scheme. The fictitious crack tip is always considered completely open at the beginning of the iterations. The closing pressure is introduced within the dynamic iterative analysis. The 151 ^HiChapter 6 Numerical Modeling of Fracture Growth last positive gradient configuration is stored and used as the initial configuration for the next candidate load factor. This procedure allows the F P Z length control scheme to produce a unique solution. Figure 6.15. Zero-slope condition for the opening profile at the fictitious crack tip When the global equilibrium and stability criterion are satisfied for a given crack or fracture process zone length, the fictitious crack tip is propagated to a new position. Thanks to the introduction of automatic mesh generation, application of computational mechanics in the case of moving boundary problems has become much easier than before. This the case of using the finite element method for a discrete analysis of the crack growth. A m o n g the reasons for using the discrete crack model with interface elements and singular crack tip elements, one can mention [35] Interface elements provide a consistent and natural way of describing the fictitious crack model Interface elements are economical since the crack is modeled with the minimum of degrees of freedom Mesh regeneration accomplished associated with crack propagation, automatically 152 and direction change, is Chapter 6 Numerical Modeling The formation of Fracture Growth of a discrete crack is viewed as a problem of changing geometry rather than, as in the smeared crack approach, a change in mechanical properties. The most serious disadvantage of the discrete crack modeling is the computational cost. The method is computing intensive since large numbers of degrees of freedom are involved. This problem is more crucial when many cracks are present, like in the case of shear wall panels and concrete shells. In such cases of distributed damage, where thousands of cracks may be involved, the smeared crack approach is more practical. In the next section we w i l l see an applications of interface modeling in conjunction with F R C materials. In this application, an assessment of C F R C composites as potential candidates for concrete repairs is performed and both linear and nonlinear constitutive laws have been used. APPLICATION TO FRC MATERIALS Concrete structures under severe conditions often start showing symptoms of surface deterioration long before the end of their design lifetimes. This could become very detrimental to their main reinforcement and to the integrity of the whole structure i f not repaired fast and properly. Therefore, proper and timely rehabilitation is a crucial issue for preventing further deterioration of such structures. The inadequacy of the present repair materials in terms of durability and mechanical behavior renders such a task difficult and a new generation of high performance repair materials need to be developed. This section presents the evaluation of C F R C (carbon fiber reinforced cementitious composites) as potential candidates for such repairs. Both interfacial and bulk properties of the material are presented. Addition of fibers leads not only an enhancement in the bulk properties of the matrix by the exhibition of pseudo-strain hardening and development of a high toughness for high volume fractions of fibers, but also to a significant improvement of both the bond strength and bond fracture energy requirements, properties which are fundamental to any good repair material. A series of finite element analyses have been performed 153 to quantify the major interfacial Chapter 6 Numerical Modeling of Fracture Growth mechanisms underlying this kind of interface behavior and to predict the debonding behavior. Overview of durability-related repairs in concrete The need for repairs to concrete structures has increased dramatically over the last few decades. A disturbing fact is that in many cases the structures have less than 25 years of existence especially in the case of structures operating under severe conditions of loading and aggressive environments. Once the diagnosis of the damaged structure has been made, a choice of the repair material must be made. This choice is based usually on which property or combination of properties one wants to restore or improve: structural strength, durability, appearance, fitness for use, etc. A n important feature in the long term durability of a reinforced concrete structure is the protection of steel reinforcement. Figure. 6.16. Spoiled concrete beam. 154 Chapter 6 Numerical Modeling of Fracture Growth A very high percentage of the deterioration which has occurred arises from a failure to provide this protection. This paper is intended to deal mainly with non-structural repairs, or in other words, repairs which are intended to restore long-term durability, but which w i l l not increase to any significant degree the load bearing capacity of the structure. While no structure is absolutely maintenance free, large scale deterioration should not occur. This is why any repair material is usually required to satisfy the following properties: /. Impermeability to aggressive liquids and gases in order to prevent or slow down the corrosion of the reinforcing steel if present, 2. Adequate durability in severe climatic conditions, 3. Ability to bond properly with old concrete and restore structural integrity It is usually well accepted that executing a thick repair is easier than a durable thin repair tends to debond easily. The inadequacy of the present repair materials in terms of durability and mechanical behavior strengthen the need for a new generation of high performance repair materials which should mitigate the effects and slow down the rate of deterioration. C F R C materials (Carbon Fiber Reinforced Composites) are evaluated here as potential candidates for such repairs. But in order to obtain high performance repair materials, a clear understanding and quantification of the major mechanisms governing the repair material-substrate is needed to guide the optimization procedure. T o this end, a combination of numerical modeling to experimental measurements should allow one to quantify phenomena or mechanisms that would be hard or impossible to access experimentally. Furthermore, numerical simulations are usually cheaper and faster than the corresponding experiments providing more cost-effective optimization procedures. However, it is imperative to realize that numerical analysis is not a substitute for experiments since only experimental results allow one to conclude. The present work can be divided in two parts. The first part deals with the identification of an overall tensile constitutive law of the interface between a concrete substrate and a C F R C repair material. This kind of constitutive relations could be useful when assessing the suitability of a C F R C as a repair material for a given concrete structure. In the second part, a fracture- 155 Chapter 6 Numerical Modeling of Fracture Growth based finite element analysis has been carried out to quantify the influence of interfacial cracking on the overall behavior of the repair composite with and without fibers in the system. Mechanical properties of CFRC materials The use of fibers to strengthen materials which are much weaker in tension than in compression goes back to the ancient times. For the purpose of this work, F R C is defined as a material made with hydraulic cement, aggregates of various sizes, incorporating discrete, discontinuous fibers. Traditionally macrofibers are not added to improve the strength, though modest increases in strength may occur [30]. Rather, their role is to control the cracking of F R C , and to alter the behavior of the material once the matrix has cracked, by bridging these cracks and thus providing some post-cracking ductility, or "toughness". Over the last few years, active research has been directed toward the development of cement-based composites with high volume fraction of microfibers (carbon, steel P V A etc.). Under an applied tensile load, these composites display a pseudo-strain hardening capability which was not possible in conventional fiber reinforced concrete, and develop strengths significantly higher than the parent matrix [13]. Moreover, such composites possess significantly improved tensile strengths, tensile strain capacities and toughness. In the context of repairs, it is usually desirable that the mechanical properties of the repair material resemble as closely as possible those of the structure being repaired. The elastic moduli and compressive strengths of F R C materials are not very different from the substrate, which makes them desirable for thin repairs [8, 150]. The effectiveness of fibers in enhancing the mechanical properties of the brittle cementitious matrix is controlled by the process by which load is transferred from the matrix to the fiber (before matrix cracking) and the bridging effect across the matrix cracks at a more advanced stage of loading. A n understanding of the mechanisms 156 Chapter 6 Numerical Modeling of Fracture Growth responsible for stress transfer may permit the prediction of the stress strain curve of the composite and its mode of fracture (ductile vs brittle) and may also serve as a basis for developing composites of improved performance through modification of the fibermatrix interaction (changes in the fiber shape, treatment of the fiber surface, etc.). In brittle matrix composites, the stress transfer between fiber and the matrix should be considered for both the pre-cracking and the post cracking stages. Before any cracking has taken place, elastic stress transfer is the dominant mechanism. The longitudinal displacements of the fiber and the matrix at the interface are geometrically compatible. The stress developed at the interface is a shear stress which is required to distribute the load between the fibers and the matrix (since they differ in their elastic moduli) so that the strains of these two components at the interface remain the same. A t advanced stages of deformation, debonding across the interface usually takes place, and the process controlling stress transfer becomes one of frictional slip. Relative displacements between the fiber and the matrix take place. The frictional stress developed is a shear stress, usually assumed to be uniformly distributed along the fiber-matrix interface. This process is of great importance which controls the efficiency of fibers. Properties such as the ultimate strength and strain of the composite are determined by this mode of stress transfer. The transition from elastic stress transfer prior to debonding to frictional stress transfer after debonding is a gradual process. Interface Modeling The above mentioned mechanisms of fiber reinforcement, which are known to be true in the bulk of F R C materials, cannot take place at the interface between an F R C repair material and an old concrete substrate, since the latter one is in hardened state and no fiber bridging can occur through the interface. A n d yet, previous studies [169] have shown that addition of fibers still significantly improves the tensile behavior of the concrete-CFRC repair composite. Even though the tensile strength and the tensile strain capacity did not improve much with an increase in the fiber volume fraction, the 157 Chapter 6 Numerical Modeling of Fracture Growth toughness of the repair composite, or more specifically its ability to absorb energy, improved in a significant way. In order to model the behavior of the interface, use of interface elements has been invoked because of their ability to represent and model both geometrical discontinuities and material nonlinearities [7]. Even though the nodal interface formulation (lumped) is equivalent to the continuous interface formulation (integrated), they lead to different element stiffness matrices. In this work, emphasis is on 2 D (line) elements and a continuous (or numerically integrated) interface formulation has been adopted. Linear and nonlinear constitutive laws have been used. The constitutive models retained for the interface elements are the Normal traction-Displacement models. A n examination of the failure pattern of the specimens [7] suggests that the interface between the substrate and the repair material is the weakest link in the specimen. Addition of fibers to the system has been shown to improve significantly the post peak tensile behavior of the repair composite. One reason for such an improvement is the reduced initial cracking due to shrinkage where the presence of microfibers has been shown to play a major role in the bulk of the FRC material. In order to quantify the influence of initiation, propagation and coalescence of microcracks on the overall behavior of the repair composite, a series of finite element analysis are being performed. To account for the presence, weakness and degradation of the interface, interface elements are introduce. The interface assumed to be of zero thickness, its presence is modeled through the constitutive law of interface elements in terms of normal stressdisplacement and/or shear stress-displacement behavior. 158 Chapter 6 Numerical Modeling of Fracture Growth Dealing with nonlinear interface models requires a reliable algorithm to determine accurate solutions. T o this end, the Dynamic Relaxation solver has been adopted due to its simplicity and the potential to cope with the nonlinear behavior of materials. Smeared A p p r o a c h (Identification of interfacial constitutive laws) Assumption: All mechanisms leading to the observed behavior of the repair composite can be implicitly lumped into the constitutive law of the interface. A t the microscopic level, these mechanisms could be related to the nucleation, propagation and coalescence of microcracks. A n analysis of Figures 6.17 and 6.18 shows that the experimental results and the numerical predictions are in very good agreement. This leads to the conclusion that the adopted interfacial constitutive laws are fairly accurate and that the adopted solver is robust enough to predict with a very high accuracy the softening behavior of CFRC-Concrete Repair composites. The constitutive laws used for this simulation are: Vf=0% (No fibers) Vf=l%. (Linear elastic then exponential decay) 8<8, 8>S, (Linear elastic then linear softening) if mO if mO where a Vf=2% '• stress transferred through the interface & '• interface opening kj: constant of proportionality (/ = 0,1,2) a : material characteristic 159 Chapter 6 Numerical Modeling of Fracture Growth 8^ : interface displacement associated with the peak load , i = cT; : a - intercept of the softening branch a; : slope of the softening branch The different values of coefficients retained in this investigation are as follow: • Smooth interface: (no surface preparation) k =123, 8 0 =0.006, cr =1450MPa, mo o k, =68.5, 5 =0.012, ml k =113.65, S 2 • m2 a = -43 cr =1610MPa ] =0.01, a = 2230 MPa 2 Rough interface: (with surface preparation) k = 95, 8 mo = 0.004, CT = 2600 MPa, a = -25 k,=85, ml =0.012, a = 2000 MPa o 8 0 k =76.5, 8 } 2 m2 =0.016, a 2 =2400MPa Fracture-based a p p r o a c h Even though the smeared approach allows one to make very accurate predictions of the overall tensile behavior of the CFRC-Concrete Repair composite, it does not give any explanation as to why such a behavior occurs at the interface. What are the major mechanisms underlying such a behavior? H o w does the presence of fibers affect it? In an attempt to shed some more light on these questions, another finite element analysis has been carried out to evaluate the effect of interfacial cracking between the substrate and the repair material on the overall behavior of the repair composite. Since fibers are not expected to physically cross the interface, one explanation is the reduced shrinkage cracking in F R C materials as compared to the plain counterpart [7]. 160 Chapter 6 Numerical Modeling of Fracture Growth Fig. 6.17. Direct Tensile behavior of a CFRC Reinforced Repair composite (No Surface Preparation) Fig. 6.18. Direct Tensile behavior of a CFRC Reinforced Repair composite (With Surface Preparation) 161 Chapter 6 Numerical Modeling of Fracture Growth Assumptions: • The actual and complicated system of interfacial microcracks can be replaced by an "equivalent" crack along the interface • Interfacial behavior: (along all the interface but excluding the crack) - Linear elastic (No Strain softening): <r = k S ( (i = 1,2,3) B y equivalent crack, we mean a crack that would have the same influence on the overall behavior of the specimen. A combined numerical-experimental procedure is introduced to quantify the effect of cracking on the observed behavior. This procedure consists in determining the length of the equivalent crack that would correspond to every point on the load-displacement experimental curve, which is provided to the finite element code as input. For a given crack length, the algorithm iterates on the displacements until the corresponding load is attained while the uncracked portion of the interface remains linear elastic. Once equilibrium is attained for the specified point, the crack is extended and another point is determined. The smaller the crack increments, the more accurate are the crack evolution curves. The properties of the interface away from the crack are determined according to the linear response of the material prior to strain softening. The constants of proportionality k t retained in this analysis to characterize the interface behavior are the same as the ones used in the previous analysis. A n analysis of Figures 6.20 and 6.21 suggests that interfacial cracking is the major mechanism controlling the response of the repaired specimen and shows how addition of fibers and surface preparation positively affect crack formation and evolution at the interface. Concluding Remarks • A n attempt is made to model the behavior of a fiber reinforced cement-based repair material through a smeared and a fracture based model. 162 Chapter 6 Numerical Modeling • of Fracture Growth The smeared model, where all mechanisms leading to the observed behavior are lumped into the constitutive behavior of the interface, appears to predict the descending branch of the tensile load displacement curve quite well. • The numerical-experimental fracture-based approach gives a good quantitative insight into the mechanism of interface fracturing in the repaired specimen, thus providing a powerful tool for the analysis of experimental measurements. • The presence of fibers reduces the total width of the cracking and the rate of crack propagation Fig. 6.19. An example of the finite element mesh used in the fracture-based approach 163 Chapter 6 Numerical Modeling Fig. of Fracture 6.20. Numerical Growth Prediction of Crack Evolution (With Surface Fig. 6.21. Numerical Prediction in the Tensile test preparation) of Crack Evolution (No surface preparation) 164 in the Tensile Test Chapter 7 Micro-mechanical Modeling of FRC Materials MICRO-MECHANICAL MODELING OF FRC MATERIALS INTRODUCTION In this section we shall present briefly some experimental observations that reveal the most important micromechanical mechanisms underlying the characteristic behavior of F R C composites. For a detailed account, the reader is referred to A . Bentur & S.Mindess [30] or N . Banthia and S. Mindess [6]. In this work the term F R C refers to any kind of fiber reinforced cementitious composite. A natural way of describing the response of F R C seems to lie in the idea of considering the macroscopic behavior of matter as a consequence of its microscopic behavior (i.e. at the scale of heterogeneities). If we admit the validity of such an assumption, it seems to be very natural to approach the behavior of such materials by studying friction at crack lips and at the fibber-matrix interface, propagation of cracks, fiber pull out... Plain concrete is a material that requires reinforcement because of its law tensile strength and strain capacities. Historically, continuous reinforcing bars have been used to provide that reinforcement by placing them at the appropriate locations in the structure to withstand the imposed tensile and shear stresses. O n the other hand fibers are discontinuous and usually randomly distributed throughout the cementitious matrix. Therefore they tend to be more closely spaced, and consequently better at controlling cracking. However they are not as efficient in withstanding tensile stresses, even though modest increases in strength may occur. The fibers are not added to improve the strength, rather, their role is to affect the behavior of the material once the matrix has cracked and control the cracking of F R C . This is possible due to the bridging of the fibers across the cracks which helps improving the energy absorption capacity of the material. The structure of fiber reinforced cementitious materials has a great influence on the properties of the composite. Hence, the analysis of these composites, and prediction of 165 Chapter 7 Micro-mechanical Modeling of FRC Materials their performance in various loading conditions, require the characterization of their internal structure. The most important ingredients that must be considered include: the shape and distribution of the fibers, the structure of the bulk cementitious matrix, and the structure of the fiber matrix interface. The fiber reinforcing array can be classified as I D , 2 D or 3 D depending on their dispersion in the cementitious matrices. A n important parameter in controlling the performance of the composite is the spacing between the fibers. However the bulk cementitious matrix is not significantly different from that in other cementitious materials. The vicinity of the reinforcing inclusion in cementitious composites is characterized by a transition zone in which the microstructure of the paste-matrix is considerably different from that of the bulk paste away from the interface. Several effects which should be taken into account with respect to the fiber-matrix bond and the debonding process across the interface are exerted by these fiber-matrix interface characteristics. The matrix in the vicinity of the fiber is much more porous than the bulk paste matrix. The weak link between the fiber and the matrix is not necessarily at the actual fiber-matrix interface, it can be at the porous layer as well. M ultiple Cracking Zone # ^ L High V «•» * \ «»» * I f / / ' \ V : \ Low V f .»' f \ Strain Softening Zone *\ ^ x f III I •< w Strain ) • Crack Opening Strain Figure 7.1. Typical stress-strain curves for low fiber volume and high fiber volume FRC 166 Chapter 7 Micro-mechanical Modeling of FRC Materials The pre-peak tensile behavior of F R C is characterized by the process of microcrack propagation in the matrix prior to the formation of a continuous crack system across the critical section at the peak load. The process by which load is transferred from the matrix to the fibers, and the bridging effect of the fibers across the matrix cracks which occurs at a later stage of loading are believed to be the main mechanisms controlling the effectiveness of fibers in enhancing the mechanical properties of the brittle cementitious matrix. The effects of cracking can be quite different in the precracking case and in the post cracking case. Therefore they should be considered separately. Before any cracking has taken place, the dominant stress transfer mechanism is mainly elastic, and a geometrical compatibility exists between the longitudinal displacement of the fiber and the matrix. The stress developed at the interface is a shear stress which is required to distribute the external load between the fibers and the matrix so that the strains of these two components remain the same at the interface. The pull-out action of fibers would be mobilized when a crack tends to widen at a critical section; this starts at peak load and tends to dominate the post-peak behavior. A n important amount of energy is consumed in the process of debonding and fiber pull-out. The pre-peak tensile behavior of S F R C (Steel Fiber Reinforced Concrete) deviates from linearity when microcrack propagation has already occurred. The elastic shear transfer is the major mechanism to be considered for predicting the limit of proportionality and the first crack of the composite. A t advanced stages of loading, the process controlling stress transfer becomes one of frictional slip since debonding across the interface usually takes place. The frictional stress developed is a shear stress, and relative displacements between the fiber and the matrix take place. The post-peak tensile behavior in S F R C is marked by the opening of one crack at the critical section which transfers practically all tensile stresses to the fibers bridging the crack. From then on, it is the debonding and fiber pull-out that largely provide the post-peak tensile resistance. The matrix residual tensile strength has a relatively small effect in the post-peak regime. This process is of greatest importance in 167 Chapter 7 Micro-mechanical Modeling of FRC Materials the post cracking case. Properties such as the ultimate strength and strain of the composite are controlled by this mode of transfer The transition from elastic stress transfer prior to debonding to frictional stress transfer after debonding is a gradual process during which both processes are present. A comprehensive approach to stress-transfer problems requires simultaneous treatment of the elastic shear transfer, frictional slip, debonding and normal stresses and strains. It has been shown that the tensile behavior after peak strength of quasi-brittle materials such as concrete and rock and certain brittle matrix composites reinforced with short fibers or whiskers can be represented by a tension softening curve that describes the decreasing traction as a function of crack opening displacement and provides an estimate for the material's toughness through the area under the curve. Development of the fracture process zone and hence the R-curve behavior have been shown to be controlled by this softening curve. Different fiber types induce different failure mechanisms, leading to different shapes of the tension softening curve. In composites with good fiber dispersion, the failure mechanism is controlled by fiber pull-out or rupture depending on fiber geometry and interfacial bond strength. Such behavior is exhibited by low fiber volume fractions (<2%) of polypropylene fibers, and steel fibers. T w o major motivations for the development of constitutive models relating the microstructural parameters to the mechanical behavior of fiber composites: guide the optimization of material behavior by tailoring the types and forms of the constituent components - predict the mechanical response of end products made of such materials While pre-peak stress-strain behavior and associated mechanical properties of fiber composites together with the post-peak tension-softening behavior have been extensively studied [6,30,108], the case of pseudo-strain hardening in F R C materials has not received enough attention. 168 Chapter 7 Micro-mechanical Modeling of FRC Materials It is generally accepted that for F R C s with Vf<2%, the major contribution of the fibers is after the matrix strain localization, which occurs around the peak of the tensile stress strain curve. However new processing techniques have helped in the manufacture of thinsheet products with fiber volume fractions (Vf) as high as 15%. In this type of materials fibers not only improve the ductility of the material, but lead to an important increase in the strength of the composite. This increase in strength termed strain hardening is associated with the appearance of multiple cracking in the specimen which requires a higher energy input to open the microcracks: See Figure 7.1 for a typical stress strain response of high V f with discontinuous fibers. The response consists of 3 regions: I - elastic region (up to matrix cracking) JJ - multiple cracking region (up to max. post cracking pt) III- Failure region (crack opening localization) MODEL DERIVATION The previous section shows the difficulty associated with the modeling of F R C composites. In an attempt to propose a model for these materials we w i l l consider the use a multi-scale approach in two steps. A micro-mechanics based model which allows the determination of a stress-crack opening displacement constitutive law for the composite after first cracking of the matrix is constructed. The model assumes that deformation of the material after matrix first cracking is mainly controlled by the microcrack nucleation and propagation even in the pre-peak zone. In a second step the F R C constitutive law is implemented in an interface element and used to compute the response of F R C structures. Such a procedure captures the essence of F R C materials, on the microscopic as well as on the macroscopic scales. The multi-scale approach leads one to consider three different scales: Microscopic Scale: fiber, matrix and fiber-matrix interface - Macroscopic Scale: Scale of structural components or material specimens 169 Chapter 7 Micro-mechanical Modeling of FRC Materials - Mesoscopic Scale: Intermediate scale between the microscopic and macroscopic scales B y nature of the problem, there exists an infinite number of mesoscopic scales. It is at this level the a choice of a representative volume element of the material is made. Transition from the microscopic scale to the macroscopic one is achieved through integrating the behavior of the microscopic scale. Since fibers are randomly distributed ib orientation as well as in location, we introduce a probabilistic density distribution to account for the random distribution of fibers in location and make use of a correction factor to account for the random orientation of fibers throughout the matrix. Experimental tests suggest that two fundamental mechanisms are behind the observed nonlinear behavior of F R C composites: Cracking of the cementitious matrix and debonding of fiber matrix followed by - Fiber slippage. MICROMECHANICAL MODELING STRAIN HARDENING ZONE MODEL ASSUMPTIONS a- Fibers have a 3 D random distribution in location and orientation b- Fibers are straight with cylindrical geometry c- Fibers behavior: linear elastic d- Fiber ruptures when its ultimate strength is reached e- Matrix crack is planar f- Fiber-Matrix bond: elastic prior to debonding frictional shear in the debonded area g- Average embedment length is — 4 170 interfaces, Chapter 7 Micro-mechanical Modeling of FRC Materials Figure 7.2. Single fiber bridging a matrix crack The last assumption implies that the embedment lengths are assumed to be uniformly L distributed between 0 and - j - . A t this stage, the model only predicts the tensile stressf C O D curve, no assumptions regarding matrix bulk properties are made. Some of the above assumptions hold true only for certain fiber types, while others, at best, approximate the real behavior. For instance, the model should not be used for fibers which undergo extensive yielding prior to rupture. For jce[fl,/], let us study the following free body diagram er + do f 171 { Chapter 7 Micro-mechanical Modeling of FRC Materials Equilibrium of forces together with application of boundary conditions at x = I and x = 0 cr, = — a f at x = l c v allows the determination of the debonding length ~ fE r 1 l = 1 v H v E 2 x f and the elastic strain in the fiber is given by sAx) = -^—(x-l)z ^ + ° rE ™ EV m 1 f f f The crack opening, or the relative displacement of the fiber with respect to the matrix is given by 5= \{s (x)-e {x))dx f m and can be written under the form 5 =Ev f cr with f l<—!4 and consequently the expression of the stress as a function of crack opening for the strain hardening regime is given by 2v)EE r^S f cr = V-Vf)rE m 172 Chapter 7 Micro-mechanical Modeling of FRC Materials STRAIN SOFTENING ZONE a) Case of a single fiber in the loading direction A- Assumptions a-fthat have been adopted in the previous model are assumed to be valid B- Fiber is assumed to be completely debonded at least on the shortest embedment length C - The frictional shear bond strength T(S) between the fiber and the matrix is allowed to vary with the slippage distance. Figure 7.3. Schematic representation of a fiber pull out Assumption C is made i n order to account for phenomena such as fiber surface abrasion and accumulation of wear debries observed with synthetic fibers and/or breakdown of the cement at the fiber-matrix interface due to the stiffness and hardness of the steel fiber [107]. The frictional bond is however, assumed to be constant over the entire embedded length of the fiber. A s can be seen in Figure 7.4 : 173 Chapter 7 Micro-mechanical Modeling of FRC Materials r(x) = T[S] r(x) = 0 for x<l f for x>l -S f -8 (inside matrix) (outside matrix) BEFORE LOADING: DURING PULL-OUT : P Figure 7.4. Pull out analysis of a single fiber bridging a crack The x coordinate system in Figure 7.4 is with respect to the undeformed fiber geometry (prior to slippage) and is attached to the fiber end. Thus, in the pull out process the origin, x=0, w i l l be moving with the fiber so that it can remain on its end. The axial force P(x) at a given fiber section is therefore given by X P(x) = P + JT(S)TT d ds 0 f 0 174 Chapter 7 Micro-mechanical Modeling of FRC Materials where P is a constant representing the fiber end anchorage effect when fiber end slips. 0 In the case, where P is neglected, the pull-out force on a fiber with embedment l 0 f is given by i,-e P(l 8)= \x(8)nd dt r f x=0 = r(8) n d (l - 8) f f It is interesting to note that the graphical representation of this equation would be a straight line only i f the shear bond strength r(8) is assumed to be a constant during the entire pull-out process. Most previous models in the literature were based on such an assumption which usually leads to discrepancies with experimental observations. A model in the same spirit of this one [107] has been proposed where the bond strength at a given cross section of the fiber is assumed to depend on the displacement of its center, including the part due to fiber stretching. In this model an a priori law for the dependence of r on the displacement of its application point s is usually assumed and its parameters are determined such that the theoretical predictions are as close as possible to their experimental counterparts. In the model proposed here, which assumes a constant bond strength for a given slippage distance 8, the theoretical prediction can be chosen as close as desired, provided one chooses enough identification points, where the model and the experiment are forced to have the same results. The accuracy o f this identification leads to the accuracy of the shear bond strength variation with fiber slippage. b- Case of a cracked cementitious composite Based on the above described fiber pull-out representation, a model is developed here for approximating the case of stress transfer across a crack i n a cementitious composite containing short and randomly distributed fibers in both location and orientation. T o this end the effects of fiber embedment length and fiber orientation have to be accounted for. 175 Chapter 7 Micro-mechanical Modeling of FRC Materials EFFECT OF EMBEDMENT LENGTH AND INCLINING ANGLE To take into account the effect of fiber embedment length we assume that these fibers follow a random uniform distribution given by p(l ) f =— f with l e[0,L /2] f f L This distribution assumes that the fiber embedment lengths are uniformely distributed between 0 and L /2. f Thus, the number of fibers dN located at z to z+dz bridging across the crack may be calculated from dN = N, p(l ) dl f f for l e [(9, L f f j2\ where JV, is the total number of fibers in the matrix of volume A L c f which contains fibers bridging the matrix crack plane of Area A . Consequently, c total fiber volume AL v = volume per fiber L/A f N= t c A„ =—v A f f f f f where A f f denotes the fiber cross-sectional area and v the fiber volume fraction. The / load transmitted across the crack through all fibers having embedment lengths l f the slippage distance is equal to S SF(l ,S) f = = P{l ,5)dN f 2Ay f —^- P(l ,S)dl L f f ff A L Thus the total load transmitted by all fibers ridging the crack is given by h ——s F{8)= \SF{l ,5) f / o / = 176 when Chapter 7 Micro-mechanical Modeling of FRC Materials Consequently, the total load transmitted by all fibers across the crack can be evaluated as a function of 8 -i-S F 2v A, L, A 2 JV(s) n cr (8) = ^^-T(8)n k L d c f [L 2 f d f (l f - 8) dl f - 8 L , 8 + 128 ] 2 f f It is interesting to note that, even in the case where the frictional shear bond strength r is a constant (does not depend on 8), the transfer of stress across a crack in the composite is not a linear function of the crack opening. In a random distribution of fibers, only a fraction of the volume of fibers bridging the crack are actually parallel to the load direction, and many fibers contributing to enhancing the tensile strain and stress capacities bridge the crack at an angle with respect to the crack flanks. Figure 7.5. A fiber bridging a matrix crack at an angle 177 Chapter 7 Micro-mechanical Modeling of FRC Materials EFFECT OF FIBER INCLINATION: The probability density function, p(0), corresponding to a 3 D random orientation of fibers is given by p(0) = Sind This density function is derived by assuming that a fiber or its extension has an equal likelihood of crossing any point on a hemisphere centered at the point where the the fiber intercepts the crack. The stress transmitted through a fiber with angle 9 with respect to the loading direction is given by CT(8,9) = CT{8)COS 6 2 Thus, the stress transmitted through fibers with orientation 0 to 0 + d0 is given by dcr (8,0) = cr (8)Cos 0 p(0)d0 2 c c and the total stress transferred across the crack can be obtained as cr (8,0) = | c r (8) Cos 0 p(0) d0 2 c c 0 «/2 = cr (8) \Sin0Cos 0d0 2 c 0 er (8) = c — T(8)xd 3A, L d f f [L f -SL 2 f J J f f J 8 + \28 ] 8 + \28 ] 2 J SUMMARY 1- Linear Elastic Zone E = E v +E m m v f f m m 2- J Strain Hardening Zone 2vl CT = 3- J EE, T / / max 8 Strain Softening Zone <J (8) = C 16v —r(8)7td 3A, L, f f 178 d [L , -SL 2 f f 2 Chapter 7 Micro-mechanical Modeling of FRC Materials EXPERIMENTAL PROGRAM This program has been established for the double purpose of model identification and validation in the case of carbon fiber reinforced cementitious composites. Only the direct tensile test data are needed for the identification of the above described micromechanical model, whereas, both compression and direct tension tests are used to identify the parameters appearing in the continuum damage mechanics model. The bending tests presented here, together with other results relative to a double cantiliver beam [62] and plates under central loading [14], have been used for model validation by comparing the model's predictions to the experimental results. Two lengths (3mm and 10mm) of pitch-based carbon micro-fibers have been used in this investigation. Characteristics of these fibers along with the m i x proportions are given in Table 7.1. Fiber Type Diameter (fjm) Pitch-based 18 Carbon Fiber Mortar L f L f Specific Gravity Tensile Strength (MPa) 1.7 590 E (GPa) 30 W/C S/C SF/C = 3 mm 0.35 0.5 0.2 =10 mm 0.35 1 0.15 Table 7.1. Mix proportions and fiber characteristics Only Portland cement of type 10 has been used along with aggregates having a maximum size of 2.38 mm. These aggregates have been used as delivered with no particular pretreatment. Carbon fibers have a uniform cylindrical shape. Various fiber volume fractions have been incorporated into the mixes. For fiber length L =3mm, 0%, 3%, 5% and 7% fiber volume fractions have been used, whereas for L =10mm, 0%, 2%, 3% and 5% f f 179 Chapter 7 Micro-mechanical Modeling of FRC Materials have been used. Since we are dealing with high fiber volume fractions (V. > 1%) of very fine fibers and hence very high surface area, silica fume has been used as a dispersing agent and considerable amounts of superplasticizer (more than 3% by weight of cement in the case of V > 7%) have been incorporated into the mix. Use of the omni-mixer has f greatly facilitated a uniform mixing and dispersion, especially in the case of the short fiber length (L f =3mm). For volume fractions as high as 7%, the mix was very uniform with very good fiber dispersion suggesting that in the case of L =3 mm, even higher f volume fractions could be achieved, whereas in the case of L = 10 mm, carbon fibers f started to ball and show a non uniform dispersion at a volume fraction of 5%. To avoid microcracking due to plastic shrinkage, molds are covered with a plastic material and kept in the moisture room. After 24 hours, the specimen were demolded and put into water for one year. So the results presented in this chapter are relative to a one year old C F R C specimens. Direct Tensile Test 2.5 I [cm] •CT 19.5 •I Figure 7.6. Shape and dimensions of the specimen In order to force cracking to take place in the central area of the specimen and away from the grips, the reduced section specimen shown on Figure 7.6. has been adopted. This way failure takes place in an area of uniform stress where the L V D T ' s have been placed to measure the displacement of the specimen (see Figure 7.7). The acquisition system is ISO Chapter 7 Micro-mechanical Modeling of FRC Materials composed of a micro-computer connected to a 2 K N load cell and two L V D T ' s recording the displacement and load values. The readings of the two 5mm L V D T ' s are later averaged to avoid possible disturbances due to the rotation of the specimen. The tests were performed at a loading rate of 0.1 mm/mn to simulate quasi-static loading conditions under displacement control. The displacement control of the loading is necessary in order to obtain the strain softening response of the specimen. The direct tensile test is one of the most difficult tests to carry and hence requires a lot of care during test preparation. A m o n g the most important problems related to this test, one can mention Alignment of the specimen is difficult Correct positioning of the L V D T ' s Slippage at the grips A s can be seen on Figure 7.8, there is a significant increase in both the strength and strain capacities as the fiber volume fractions become higher. A summary of the main results is presented in Table 7.2. 181 Chapter 7 Micro-mechanical Modeling of FRC Materials Chapter 7 Micro-mechanical Modeling of FRC Materials Direct Tensile^Test Lj=Jm » V >'/,-. * \ V. 7M * * •Vf=0% % if \ •Vf=3% 1 •Vf=5% I i * Vf=7% — f 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 S f r a i n (%) Direct Tensile Test I.f-lOm m v. /s' «, is •s •» ••- V "V \ Vf=4% "*"•« " ^ f1 % * Vf=5% V 1 0.1 0.2 0.3 Strain Figure 7.8. Stress-Strain For L =3 mm f 0.4 (%) Curves in Uniaxial and 183 Vf=0% -Vf=3% \ \ 0.0 • *— L , = 10 mm Tension 0.5 Chapter 7 Micro-mechanical Modeling of FRC Materials L f - Peak Stress (MPa) Peak Strain (%) 3.48 4.50 5.61 6.99 0.031 0.082 0.11 0.14 Peak Stress (MPa) Peak Strain (%) Vf=4% 3.34 4.10 5.04 0.035 0.11 0.14 Vf=5% 5.78 0.16 3 mm Vf=0% Tensile Test Vj=3% Vf=5% Vf=7% L f = 10 mm Vj=0% Tensile Test Vf=2% Table 7.2. Effect of fiber volume fraction on the tensile behavior BENDING TEST One of the simplest ways of studying fracture in CFRC materials is in flexure. The specimen sizes along with a schematic representation of the test setup are summarized on Figure 7.9. 300 mm Figure 7.9. Schematic representation 184 of the test setup Chapter 7 Micro-mechanical Modeling of FRC Materials The test has been carried out in the Instron machine where an L V D T has been instrumented under the beam at mid-span. The L V D T has a maximum displacement of 5 mm. The data acquisition system is composed of a computer connected to the L V D T and a 150 K N load cell, hence allowing the recording of both displacement as measured by the L V D T and load as measured by the load cell. A loading rate of 0.1 mm/mn has been used to simulate quasi-static loading conditions. The load displacement curves for the composites previously characterized in uniaxial tension are given on Figure 7.10. Four specimens have been used to estimate the bending response, as described in Figure 7.9, for each fiber length and fiber volume fraction of the F R C beams. A summary of the results is presented in Table 7.3. L f Bending Test L f Peak Load (KN) Peak Displacement (mm) Vf=0% 350.072 0.331 Vf=3% 416.274 0.905 Vf=5% 464.451 1.239 Vj=7% 520.342 1.274 Peak Load (KN) Peak Displacement (mm) Vf=0% 350.318 0.342 Vf=2% 361.845 0.540 Vj=4% 461.276 0.867 Vj=5% 427.438 1.223 = 3 mm -10 mm Bending Test Table 7.3. Effect offibervolume fraction on the bending behavior 185 Chapter 7 Micro-mechanical Modeling of FRC Materials Figure 7.10. Load deflection Curves for 3 point Bending Test For L =3 mm and L , = 10 mm f 186 Chapter 7 Micro-mechanical Modeling of FRC Materials COMPRESSIVE TEST Load cell LVDT holder LVDT Steel plate Figure 7.11. Schematic representation of the compressive test set up 50 mm x 100 mm cylindrical specimens have been tested under quasi-static compressive loading conditions. This test has been used for the identification of the parameters appearing in the scalar variable continuum damage mechanics model presented in chapter 4. The acquisition system is composed of a micro-computer connected to a 300 K N load cell and two LVDT's recording the displacement and load values. The readings of the two 5mm LVDT's are later averaged to avoid any possible local disturbances on the deformation of the specimen. The tests were performed at a loading rate of 0.1 mm/mn to simulate quasi-static loading conditions under displacement control. The displacement control of the loading has been adopted in order to capture the strain softening response of the specimen. 187 Chapter 7 Micro-mechanical Modeling of FRC Materials Results of this test are summarized in Table 7.4 and Figures 7.12 and 7.13. It is interesting to notice that unlike the case of the tensile response, the compressive strength of the C F R C materials investigated in this study did not increase in any significant way as the fiber volume fraction increases. Moreover, a decrease in the compressive strength has been noticed beyond a certain fiber volume fraction threshold. This threshold has been identified as being V = 7% f for L f = 3mm and V f = 5% for L f =10mm. The pseudo- ductility exhibited by the carbon fiber reinforced cement composites investigated in this study has not been observed to improve in a significant fashion as higher fiber volume fractions are incorporated into the mix proportions. L f Compressive Test L f Peak Stress (MPa) Peak Displacement (%) Vf=0% 41.276 0.66 Vf=3% 39.342 0.759 Vj=5% 40.362 1.171 Vf=7% 39.341 0.837 = 3 mm = 10 m m Compressive Test Peak Stress (MPa) Peak Displacement (%) Vj=0% 36.704 0.448 Vf=2% 32.317 0.801 Vf=4% 32.319 0.803 Vf=5% 28.079 0.924 Table 7.4. Effect of fiber volume fraction 188 Chapter 7 Micro-mechanical Modeling of FRC Materials Figure 7.12. Stress-Strain Curves for L = 3 mm f Chapter 7 Micro-mechanical Modeling of FRC Materials 40 ^v1 v"^ Gitpvxsiai 1 lf=10mt 35 1/ > v 30 25 '20 CO 15 I 1 iiiiiiiiliSlI *w V- - - - -I \ Ij^pil • \ \ lllpllls * • 1 10 S S a \ ~ 0 S ^^^^^^^^^^^ | \ » » il t -\l A? \* J 0.0 1.0 20 3.0 4.0 Figure 7.13. Stress-Strain Curves for L = 10 mm f 190 5.0 60 Chapter 8 Numerical Implementation and validation NUMERICAL IMPLEMENTATION AND VALIDATION NLFM-BASED PREDICTIONS FOR FRC STRUCTURES The micromechanical models developed in the previous chapter are used here to predict the complete behavior of C F R C made structural components. T o illustrate the combined use of a numerical analysis scheme for discrete cohesive crack modeling and the micromechanical constitutive laws for C F R C materials, the responses of a three point bending beam and contoured double cantiliver beam have been analyzed. The analysis, was able to predict the entire response up to ultimate failure of the specimens accounting for the effect of the fiber volume fraction as well as the type of fiber. The simulations were performed for two fiber types (pitch-based carbon microfibers and steel microfibers) and three volume fractions of fibers (vf=l%, 2%, and 3%). The crack path being known a priori in both cases, interface elements simulating cohesive cracks, have been introduced along the expected crack path. This has the advantage - of avoiding the explicit use of a stability criterion during the analysis - avoids remeshing as the crack propagates gives accurate predictions with a relatively low number of degrees of freedom. The main disadvantage of this approach is that, it can be used only when we have a previous knowledge of the crack path, which is the case in the two examples presented here. Since the simulation was related to the entire response, and therefore the performance of the F R C specimens up to failure, a displacement control scheme has been used to pilot the numerical test. This is the case, because a load controlled scheme would not be able to detect the strain softening exhibited by F R C composites in their post-peak behavior. T o this end, the dynamic relaxation solver introduced in chapter 6, has been used again to solve the boundary value problem at hand. 191 Chapter 8 Numerical Implementation and validation BEAM BENDING A 3 point bending test has been carried out with high volume fractions of pitch-based carbon microfibers incorporated in a mortar matrix. T w o fiber lengths were used (3 m m and 10 mm). 300 mm Figure 8.1. Geometry and boundary conditions for the three point bending test A linear elastic stress analysis is performed in a first stage. This allows one to determine the most critical zones and the direction of the principal tensile stress. Once this information is available a series of nonlinear interface elements are inserted along the expected crack path. This information could be available from theoretical considerations or simple examination of failure patterns of tested specimens. Indeed i f there is strong experimental evidence that the crack path follows a certain curve or direction it is perfectly natural to force the crack to follow such a path without paying too much attention to the various criteria for determining the direction of propagation. Unfortunately, in many practical situations such an information is not available and propagation criteria still have all their usefulness. A s mentioned earlier, the presence of the nonlinear interface elements simulating the effect of fibers bridging the crack w i l l automatically open according to the cr(w) equivalent to a crack opening and propagation. 192 stress-COD constitutive law which is Chapter 8 Numerical Implementation and L f validation Peak Load (KN) = 3 mm Experimental Model f = 3% 416.274 0.905 V f = 5% 464.454 1.239 V f = 7% 520.343 1.274 V f = 3% 445.476 0.875 V =5% 503.944 1.227 V 549.884 1.282 f =7% Peak Load (KN) = 10 mm 0.540 V 461.276 0.867 V =5% 427.438 1.223 V =2% 375.283 0.507 V = 4% 473.927 0.622 V = 5% 424.036 0.922 f = 4% f Model f f Table 8.1. Comparison Displacement (mm) 361.845 f j Peak V =2% f Experimental Displacement (mm) V f Lf Peak between model predictions 193 and experimental results Chapter 8 Numerical Implementation and validation 600 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 Deflection (mm) Figure 8.3. Experimental results versus numerical predictions for Lf=3mm 194 0 4 0.0 : M 0.5 1 1.0 1 1.5 1 1 2.0 2.5 1^- 3.0 1 3.5 1 = 4.0 Deflection (mm) Figure 8.4. Experimental results versus numerical predictions for Lf=10mm 195 1 4.5 1 5.0 4 Chapter 8 Numerical Implementation COUNTOURED and DOUBLE validation CANTILIVER BEAM The 300 mm shown in Figure 8.5 (Banthia and Genois 1995), has been selected in order to assess the effectiveness of two types of fibers on the overall performance of the F R C composite. Pitch-based carbon microfibers as well as steel microfibers have been incorporated into the composite mixes Genois (1995). The test consisted in propagating a crack in a pre-notched specimen while recording the load required to open the crack and the corresponding displacement of the loading points. A wedge is driven between two pairs of wheels to split the specimen into two halves. The experiment is performed in displacement control. Figure 8.5. Specimen geometry and a schematic of the loading procedure (after Genois [62]) 196 Chapter 8 Numerical Implementation and validation Before a cohesive crack analysis is performed, a linear elastic stress analysis is performed to obtain the location and direction of the maximum tensile stress. This is usually a good indication about the area where a crack is likely to nucleate and/or propagate. Due to the symmetry of our specimen and strong evidence from experimental tests, the path of the major propagating crack is along the symmetry line of the specimen. Thus interface elements have been introduced along the crack path to simulate a cohesive crack propagation. The effect of the fiber type and volume fraction is accounted for through its specific CT(VV) constitutive law as presented in chapter 7. Figure 8.6. Numerical prediction of the crack growth at an intermediate stage 197 Chapter 8 Numerical Implementation and validation Steel fiber - Cement Paste Experimental f = 1% 1186.59 87.728 V f = 2% 1485.51 131.593 = 3% 1884.06 167.06 V, = 1% 1313.41 79.985 V = 2% 1585.14 130.026 = 3% 1992.75 142.559 f f V f Table 8.2. Comparison between model predictions (Steel fiber-Cement V results Peak Displacement (um) f = 2% 1616.07 97.488 V f = 3% 2125.01 155.023 V f = 2% 1651.79 84.703 V, = 3% 2312.50 151.826 Model Table 8.3. Comparison and experimental paste) Peak Load (N) Steel fiber - Mortar Experimental Peak Displacement (um) V V Model Peak Load (N) between model predictions and experimental (Steel fiber - Mortar) 198 results Chapter 8 Numerical Implementation and validation '••mm* Exper.Vf=3% DCB: Steel-Paste Experiment V S Numerical • - -- - -Model.Vf=3% Exper.Vf=2% - -*- - Model.Vf=2% Exper.Vf=J% Model.Vf=l% V \ \ \ 1 \ •iiiliil 0 100 200 300 400 500 CMOD (um) Figure 8.7. Numerical Prediction VS Experimental results (Steel fiber-Cement paste) 199 600 700 800 Chapter 8 Numerical Implementation and validation J Exp.Vf=3% Sim.Vf=3% / I 1 f • 1 r• : :// t'i Exp.Vf=2% i1CB: Steele•tartar •' Sim.Vf=2% v. • »* i • W ;.'/ \ IKIlBillllI •X \ \ * \ /"/ 111ll|l|lHiH!ll!iti!i!iiW!l!lIiliIlilHt ^^JIM^^^^h^^^^^^H^flfl l l l l i i l l l l l l f illlpliMplMiili^MMi lliilllllpiMSiPiii \ \ iiiiiiiiiiiiiiiiiiiiiiiiiti • - V v \ ii! \ Hi NjiiyhisllllS;!;;;;:;:;;:.:... | !J ! ;•;! ii; ^ I i M.i H! NI LL^Lf:S^[MHin^ jjiJllllljMllllBiBM ISpiiilll^Bil • - . - \ :-i;' : \ 1 s ••iiMl \ Siiii\ii iiiiiiSi iiiiiiiiiiiiiiiiiii:! ' \ • p i i a i p i i i ^ i i H \sm$smsmsmm \ V 0 100 200 300 CMOD 400 500 (um) Figure 8.8. Numerical Prediction VS Experimental results (Steel fiber-Mortar) 200 600 700 I Chapter 8 Numerical Implementation and validation Carbonfiber- Cement Paste Experimental Peak Displacement (N) (urn) V f =1% 1337 163.061 V f = 2% 1684.98 310.29 V = 3% 2106.23 351.45 = 1% 1419.41 128.232 V =2% 1813.19 240.633 V,=3% 2274.9 301.544 f V f Model Peak Load f Table 8.4. Comparison between model predictions and experimental results (Carbonfiber-Cementpaste) Carbon fiber - Mortar (urn) V =l% 1675.72 157.859 Vf = 2% 2010.87 197.722 V = 3% 2318.72 301.367 =1% 1784.40 132.346 V =2% 2137.68 172.21 V 2841.88 234.396 f V f Model Peak Displacement 1 (N) f Experimental Peak Load f f = 3% | Table 8.5. Comparison between model predictions and experimental results (Carbon - Mortar) 201 Chapter 8 Numerical Implementation and validation 0 100 200 300 400 CMOD 500 (urn) Figure 8.9. Numerical Prediction VS Experimental results (Carbon fiber-Cement paste) 202 600 700 800 Chapter 8 Numerical Implementation and validation Vf=]% li('H (\uh<»i-\1i>nai Experimental results Vf=2% Vf=3% - — - SimVf=2% ----- SimVf=3% - -+ - SimVf=0% \ \\ 0 4 100 200 300 400 500 CMOD (urn) Figure 8.10. Numerical Prediction VS Experimental results (Carbonfiber-Mortar) 203 600 700 Chapter 8 Numerical Implementation and validation A COMBINED CDM-NLFM APPROACH From our previous chapters it is tempting to consider continuum damage mechanics and fracture mechanics as two correlated theories. The equivalence between the two approaches can be established i f one considers the situation where damage is equal to one at a material point, or in a small region defining the size of an initial flaw in the theories of fracture. Physically, this usually is related to the occurrence of strain softening, most of the time associated with the transition from a diffused damage of microcracks, in the case of cementitious materials, to localization into one macrocrack. A s we saw in chapter 4, the local version of damage mechanics has been shown to be inadequate to describe this advanced stage of loading. A m o n g the major problems that have been faced with the numerical implementation of such models is the fact that failure occurs without dissipation of energy. Indeed, damage models predict that the dissipated energy density is finite at each material point, but i f damage localizes into a region of zero volume, the total amount of energy dissipated to a crack viewed as a line along which damage is equal to one at each material point, vanishes. Most fracture theories use energy considerations or quantities as crack propagation criteria. Consequently, there is a major inconsistency between fracture mechanics which assumes that this quantity is finite and greater than zero, and classical damage theories which predict that this same quantity is zero. On the other hand, progressive cracking, is usually introduced in fracture mechanics as cohesive zone along the crack lips which results in the so-called fictitious crack models. This is actually a way of lumping the volume distribution of microcracks into a surface [57], and the equivalence between damage and fracture is obtained through energy considerations: the energy dissipated during fracture is equal to the work of the pressure introduced on the crack surfaces. Since the fracture process zone is collapsed into a line, the distribution of damage in this zone cannot be obtained and this approach is inadequate if one is interested in predicting the size and shape of a diffused damage. A bridge between damage mechanics and fracture mechanics theories can be established in a consistent way through nonlocal damage models or gradient dependent theories. 204 Chapter 8 Numerical Implementation and validation Choice between the gradient-type theories making use of higher gradient of relevant variables and nonlocal theories where some variables are made dependent on what is happening in an entire neighborhood surrounding a given material point, is most of the time a matter of preference of the model builder. Use of higher order theories usually leads to the introduction of "nice" (although of higher order) differential equations, whereas use of nonlocal theories yields integrodifferential equations. Physically, the two approaches are still expected to yield the same finale solution. Nonlocal damage theories have solved the problem of spurious mesh dependency due to damage localization and have the capability of capturing the diffused deterioration of materials and predicting the complete response of a specimen or structural member up to failure. Even though C D M - b a s e d and fracture theories are usually applied to the same problems and expected to yield similar predictions, it is physically more appealing or more realistic to use damage mechanics when dealing with pre-peak response where no major initial flaw is supposed to exist until peak load where usually the material has undergone some deterioration and flaws become apparent. Once the damaged areas have been identified and quantified, equivalent crack-like flaws can be introduced into the structure and the concepts of fracture mechanics can be invoked to resume the response of the specimen or structural component up to ultimate failure. This is physically more sound, since after peak load the response of a given structure is usually controlled by the opening and propagation of a major crack rather than a diffused damage. This idea of combining the two theories is illustrated in this section through the prediction of the complete response of a fiber reinforced double cantiliver beam. Similar results have been performed in the previous section performing a pure numerical N L F M approach. A theoretical justification for the connection between the two approaches is provided here in the framework of thermodynamics for the simple case of linear elasticity. BRIDGE BETWEEN A unified DAMAGE framework AND FRACTURE for presenting thermodynamics of continuous media damage [102], and fracture is provided by since this theory deals with the energetic considerations, which allow a consistent link between local damage variables and global 205 Chapter 8 Numerical Implementation and validation fracture variables. These considerations are usually based on a particular choice of the free (reversible) energy stored in the material during the deformation process. The state equations are deduced from the free energy defined as <F = U-TS where U is the internal energy, T temperature and S the entropy. For a given state of damage D, the strain energy of an elementary volume element and its associated specific free energy are given by ~^Cijkl ij kt U = £ £ i// = u-Ts where C j]kl is the local stiffness matrix at a given state of damage and s is a component tj of the strain tensor. If a structural component or a finite continuum is subjected to a load Q, the strain energy, associated with the damaged and/or cracked material, stored in the entire body is given by U=-Kq 2 2 where K is the global stiffness and q the displacement corresponding to the applied load Q. In the simple case of linear elasticity and isotropic damage, the relationship between the current stiffness matrix C° u and the initial stiffness of the undamaged material Cf jkl is given by C° =C° (l-D) t kl Assuming a constant and uniform temperature the state laws for the damaged and cracked bodies are respectively given by: The damaged material ^=f^- = dD C° {l-D)e kl 2 206 1 kl 1 Chapter 8 Numerical Implementation and validation The cracked structure dW dq „ G dW = I dK = —q — dA 2 dA 2 Wher e A is the actual area of the crack and -G is the fracture energy release rate. A s we saw in chapter two, the first and second principle are respected i f the Clausius-Duhem inequality is satisfied. For the two cases in which we are interested, this inequality yields YD>0 -C IIU GA>0 s a s D>0 hl dA J 2< L Noticing that K is a quadratic form and that K decreases when the fractured surface A increases, these equations lead to D>0 A>0 Clearly showing that irreversibilities correspond to increase in either diffused (damage mechanics) or localized damage (macrocracking in fracture mechanics). When a cracked concrete structure is subjected to loading, the applied load results in an energy release rate G at the tip of the effective quasi-brittle crack that could be split up into two parts: q i- the energy rate consumed during material fracturing in creating two surfaces, G ii- the energy rate to overcome the cohesive pressure [c G =G, q Ic +G CT(VV) in separating the surfaces, G . a (5.22) a According to its definition, the value of G can be computed as follows a dxdw 207 and (5.23) Chapter 8 Numerical Implementation and where normal cohesive pressure and CT(W) is the validation w c is the crack separation corresponding to cr(w) = 0 . In the case where the crack opening profile (shape) does not depend significantly on the crack length, the above equation becomes G = %'a(w)dw (5.24) a E q . (5.22) is a general energy balance condition indicating that for quasi-brittle fracturing, the energy release rate due to the applied load G is balanced by two fracture q energy dissipation mechanisms. A bridge from damage mechanics to fracture theories can be established by defining the way of transforming a given damage zone into an equivalent crack. Assuming that the equivalent infinitesimal crack extension dA of a crack to a given damage increment dD(x) at point x is given by dA = ^YdD(x)dx where V is the volume occupied by the damaged zone. For the total evolution of damage, from D(x) = D =0 o to Dj(x) >0 at point x, the equivalent crack is ^' A e ( X ) Y dD(x)dx G = i A e denotes the crack, the growth of which would have consumed the same energy as the energy consumed during damage propagation in the structure. T o illustrate, this combined C D M - N L F M approach, the case of the contoured double cantiliver beam has been reconsidered in this section. Figure 8.11.bis, shows that the numerical predictions agree fairly well with the experimental results. In conclusion it is interesting to keep in mind, that C D M approaches do not assume the presence of any pre-existing flaws and can help identifying areas where cracking is likely to occur setting the stage for a fracture analysis after peak load which is usually a more appropriate approach once microcracks have localized into a given major crack. 208 Chapter 8 Numerical Implementation and validation A. = True Crack Interface Elements Figure 8.11. Schematic representation of the equivalence 209 CDM-NLFM Chapter 8 Numerical Implementation and validation (b) Figure 8.11.bis. Numerical Prediction VS Experimental results for Carbon fiber- Mortar and steel fiber-Mortar mixes 210 Chapter 8 Numerical Implementation and validation ON THE OVERALL BEHAVIOR OF FRC MATERIALS A s we saw in chapter 2, the macroscopic behavior of a damaged material may be completely known once the specific Helmholtz free energy y/ and the damage evolution laws D are known. Indeed, given the specific free energy y/ = y/(e ,T,D ) ij ij one can compute the associated thermodynamic forces and the stiffness tensor dy/ de , s s= - dy/ ~8T dy/ 5D, m =P c a y de ds i} a The rate of change of the internal state of the solid is governed by the damage evolution equations as independent equations of evolution for every internal damage variable (6). D. = D (e ,T,D ) Another alternative for writing the kinetic equations describing damage evolution may be as derivatives of a suitably chosen potential function [141], making use of the generalized normality postulate. The Clausius-Duhem inequality which ensures that the first and second principle of thermodynamics are satisfied leads to Y.D.. >0 IJ IJ U p to this point we have been dealing with phenomenological continuum damage models which do not reflect explicitly any of the dissipative or energy transfer taking place at the mesoscale since they do not account for the specific arrangement and organization of 211 Chapter 8 Numerical Implementation and validation damage entities. However, they are still able to predict the macroscopic stiffness in a rather accurate way. In a micromechanics-based formulation, on the other hand, damage variables are chosen in such a way that the most important aspects of damage morphology are explicitly incorporated into the model. A m o n g the many damage descriptors that have been proposed in the literature one can mention the second rank tensor proposed by Kachanov [91] for the characterization of deterioration due to cracking 7 M ^=TE i>»« X fl (8 () B) 1Q a=l where A is the area of the 2 D R V E , a index denoting the number of the crack under consideration (0 <a <M ), M being the total number of cracks, n ( a ) t h e unit vector normal to crack number a, a is the half length of crack a. Similar damage variables have been developed to account for different mechanisms taking place in a whole variety of single phase materials and composites [1,158,161, 162] A major shortcoming in this approach, resides in the fact that macroscopic damage variables obtained through a spatial average of individual damage entities, assumes that damage is uniformly distributed throughout the R V E . This reprents a serious limitation in formulating damage evolution laws that include s u b - R V E length scale interaction effects. The nonlinear response of engineering materials depends on the type, size, distribution and orientation of microdefects in the materials. W i t h an increasing load, microdefects in a material evolve, and the effects of interaction among the microdefects become dominant, governing the macroscopic failure. Horii's theoretical work, based on micromechanics, has revealed that the dominant mechanism of localization is the effects of interaction among microdefects. Practically, the damage growth kinetic equations are usually postulated based on empirical evidence or computational suitability. A number of simplifying assumptions to make the problem tractable are often made. These have usually the inconvenient of restricting the validity of the analysis to the dilute damage case under a single specific type of loading. In an attempt to estimate some of the effects of the simplifying assumptions on the determination of effective moduli and damage 212 Chapter 8 Numerical Implementation and validation evolution in classical closed form homogenization techniques a numerical analysis is performed in this section. Figure 8.12. Numerical visualization of the main difference between dilute and interacting cracks in an RVE The effect of initial crack patterning on the RVE-effective stiffness and thermodynamic force associated with the chosen damage variable has been investigated through parametric studies for a number of periodic crack distributions of parallel cracks in a twodimensional case. Both brittle and fiber reinforced cementitious materials considered. Overall moduli and thermodynamic forces associated with damage have were calculated at each damage increment within a numerical simulation of evolving cracks. The aim of this section is an attempt to illustrate the effect of the relative size of cracks length with respect to the RVE size on the effective stiffness and thermodynamic forces. The inability of a low-order damage measures to capture the difference crack distributions that yield markedly different global 213 behaviors. between Chapter 8 Numerical Implementation and validation A.schematic representation of a typical R V E used in the numerical study is shown on Figure 8.13. where a distribution of traction free or cohesive cracks parallel to the x, direction in a R V E consisting of a linearly elastic, homogeneous, isotropic matrix is presented. The traction-free cracks were used to simulate the case of unreinforced matrices, whereas the cohesive cracks were introduced to account for the presence of fibers in the cementitious composite. Periodic boundary conditions were applied to simulate a repeating mesostructure. It is interesting to note that periodicity of the R V E does not necessarily mean periodicity of inclusions distribution within the R V E . One can easily immagine the genesis of a macroscopic material by joining together a large number of R V E ' s containing randomly distributed defects or inclusions. A displacement u 2 applied in the x 2 was direction to the upper R V E boundary. Crack whose propagation criterion is satisfied are extended in a self-similar manner. The propagation criterion is K, = K IC (Tj = f t for traction-free cracks for the cohesive case, with CT being the tensile principal stress and f ; t the tensile strength of the material The effective moduli and damage thermodynamic force together with local driving forces, were computer at each stable damage configuration. Damage distribution is characterized by the crack density tensor in equation (10) which gives a typical low order damage representation [92]. For the particular configuration considered in this study, this tensorial damage tensor has only one non-zero component given by 7 D = D 2 2 - - E«, N j=l and the damage thermodynamic force, also has only one non-zero component, oD 22 214 Chapter 8 Numerical Implementation and validation A n expression for Young's modulus E in the x 2 direction can be derived using Hooke's 2 law for orthotropic media and associated symmetry properties [102] to yield e 22 where E and v 0 E 2 -a 22 E +v 2 0 EG 0 (15) =0 22 are the initial Young's modulus and Poisson's ratio. e 0 22 and a 22 are the average normal strain and stress, acting on the R V E boundary in the x , direction. Figure 8.13. A typical RVE prior to strain localization Continuous damage evolution in a damaged R V E may be numerically simulated through a sequence of N steps of damage increments associated with N+l stable damage states, D, U) where i=0, threshold, W , 0) ,N. Each stable damage state, D , U) has a strain energy density necessary for damage evolution to take place under a given load. If the strain density is chosen to be the thermodynamic potential, then equation (14) can be approximated using the three-point formula [79] y(0 _ D (>+/) 215 •D (i-I) (16) Chapter 8 Numerical Implementation and validation where i denotes the i' stable damage state, Y the damage thermodynamic force. T w o - point formulas analogous to Equation (16) can be used to evaluate Y U) and 7 [79] ( J V ) Equations (15) and (16) may therefore be used to compute the averaged modulus and the damage thermodynamic force. Uniformly distributed crack patterns are used to illustrate the effect of the relative distribution of crack lengths on the effective moduli and damage thermodynamic force. Figure 8.14 shows three uniform crack distributions over the R V E consisting of one, four and sixteen cracks, respectively. Each distribution has the same initial normalized horizontal and vertical spacing between neighboring cracks, whereas the initial crack length of the second and third are about Vi and A of that of the first configuration. l Figure 8.15 reveals that, whereas the evolution of the normalized stiffness E /E 2 0 is practically the same for all configurations over the entire range of damage variation, the damage thermodynamic force is different for each case. A s can be seen on Figure 8.17, the normalized thermodynamic force, Y/Y REF , necessary for damage evolution decreases as the characteristic crack size of the distribution increases, where Y is the damage thermodynamic force and Y REF is the thermodynamic force associated with the initial damage (crack density) shown in pattern HI chosen as a reference. Despite the difference in the value of the normalized damage thermodynamic force, the three configurations show a similar tendency in that, YJY REF decreases rapidly with increasing damage, then assumes an asymptotic value which is crack size dependent. This might be caused by the increasing interaction between parallel cracks as they evolve within the R V E . A physical interpretation of the aforementioned results could be achieved i f the damage thermodynamic force is thought of as the critical energy threshold per increment of damage that must be exceed for deterioration to occur. Recall that in chapter 4, the thermodynamic force associated with the continuum scalar damage variable has been shown to play the equivalent role of the critical strain energy release rate in fracture mechanics. Having this analogy in mind, one could say that configurations with larger cracks reach a critical value and evolve at lower load levels and lower strain energy 216 Chapter 8 Numerical Implementation and validation release rates than configurations with smaller cracks for the same damage state (crack density), despite the fact that they have the same macroscopic stiffness. U p to this stage, we have been dealing mainly with the effect of relative crack size on the R V E response. Another important issue influencing the overall behavior is the spatial distribution of flaws. Crack configurations that tend to prefer a horizontal or vertical distribution may display a sensibly different overall behavior. Four crack distributions with a non uniform distribution are used to investigate this effect (see Figure 8.18). Unlike the case of uniformly distributed cracks, Figure 8.18 suggests that the stiffness deterioration strongly depends on the initial configuration of cracks. Distributions with a preference towards horizontal patterns (/ and III) exhibit a relatively fast decrease in stiffness as damage evolves compared to those leaning towards a vertical pattern (// and IV). Compatible results with the ones obtained in this section were observed by Lacy et al. [100] and by Deng and Nemat-Nasser [52]. h/w AJW RVE Pattern I Pattern I Pattern II Pattern III 1 1 1 0.3 0.14 0.14 Pattern II Figure 14. Influence of uniform crack distributions on the overall response of the RYE 217 Pattern III Chapter 8 Numerical Implementation and validation Figure 8.15. Effect of uniform crack distribution on the effective moduli 218 Chapter 8 Numerical Implementation and validation Figure 8.16. Effect of fiber volume fraction on the effective moduli (uniform crack distribution) 219 Chapter 8 Numerical Implementation and validation Figure 8.17. Effect of uniform crack distribution Thermodynamic on the force associated with damage 220 Chapter 8 Numerical Implementation and validation h/w /w a RVE Pattern I Pattern II Pattern III Pattern IV 2 0.5 8 1/8 0.1 0.1 0.05 0.08 II Figure 8.18. Crack patterns used to assess the influence of crack distributions on the overall response of the RVE 221 Chapter 8 Numerical Implementation and validation 120 r 1 1 —O— Pattern —•— Pattern II - - Y/Yref vs Damage (LEFM) r— 1.00 iV i r P a t t e r n l - * - Pattern — III 0.80 § 0-60 0.40 0.20 0.00 4 0 1 5 I 1 1 1 1 1 10 15 20 25 30 35 Damage (%) Figure 8.19. Effect of non uniform crack distribution on the effective moduli 222 1 40 Chapter 9 Conclusions and Recommendations CONCLUSIONS AND RECOMMENDATIONS SUMMARY High fiber volume fraction cementitious composites are a new construction materials which exhibit high performance when subjected to a judicious mechanical tailoring. They can be optimized for either a high ductility or a high strength. The effect of fiber addition is no longer limited to the softening stage, but there is a significant pseudo-strain hardening phase prior to strain localization. A t this stage of loading, it becomes easier to open a new crack than to propagate an existing one, and the phenomenon of multiple cracking occurs leading to the observed pseudo-strain hardening. Extensive experimental and analytical studies have been conducted on understanding the behavior of fiber reinforced cementitious composites at the materials level. Even though, much work is still needed in order to understand the fundamental behavior of the material, there is a noticeable lack of analytical models that close the link between the material behavior and the structural response of components making use of those materials. Although, it is possible to experimentally investigate the performance of fiber reinforced structural members, analytical modeling remains the most attractive approach due to the difficulty of carrying out test at the structural level. Based on the understanding of F R C ' s fracture behavior and in the light of the above background, an analytical model has been proposed to give some basic understanding of the material's behavior. The model is later injected into a finite element code to provide a link between the material behavior and the structural performance. Numerous experimental studies, show that under a tensile load, cracking is initiated on faces normal to the maximum principle stress when its magnitude reaches the first crack strength. The crack opens then with increasing bridging stress, due to the high volume fraction of fibers bridging the crack, which leads to a pseudo-strain hardening. Once a certain crack 223 Chapter 9 Conclusions and Recommendations opening displacement is reached, localization occurs and the crack response becomes one of strain softening. The strain hardening zone where there is multiple cracking is accounted for through the introduction of an "equivalent" cohesive crack which would exhibit the same behavior as a representative volume element of the microcracked material at this stage of loading. In other words, the tensile behavior of an rve containing several microcracks and short fibers randomly distributed within a cementitious matrix is assumed to be equivalent to the behavior of a cohesive crack embedded into the same matrix. The problem is thus reduced to the determination of the cohesive pressure in terms of the crack opening displacement. The post-peak behavior is characterized by localization of strains into a major crack, due to the propagation and coalescence of microcracks that have been developed in the previous phase (multiple cracking). This phase exhibits a strain softening, and the response of the material is mainly controlled by the opening of the crack that w i l l eventually lead to the ultimate failure of the specimen. Following the traditional approach for this stage of loading, a micromechanical model is proposed for the determination of a cohesive pressure. Using the above mentioned micromechanical models, a nonlinear boundary value problem is defined and solved numerically using a discrete crack analysis. The cohesive cracks are modeled trough the introduction of nonlinear interface elements along the crack path. The dynamic relaxation solver is used to solve the problem. This procedure was used to analyze the effect fiber content on the response of a three point bending CFRC beam the response of a countoured double cantiliver beam the performance of CFRC composites as repair 224 materials Chapter 9 Conclusions and Recommendations Moreover, the effect of crack patterns and fiber volume fraction on the effective moduli of FRC materials the damage variable and its associated thermodynamic force have been analyzed CONCLUSIONS The major conclusions obtained from this investigation are as follows. A continuum damage mechanics approach has been considered for the analysis of the three point bending tests. A s far as the overall response is concerned, one could say that both approaches yield fairly accurate results. However, one point needs to be mentioned, that is, the use of a stress-strain functional to characterize the post peak behavior is a questionable formulation. This is because after the peak load, the diffused strain due to microcracking localizes into a major macrocrack, and any subsequent deformation of the structure or the specimen is more related to the opening of the crack rather than to a bulk deformation of the material. - Micromechanical-based constitutive models have been developed for both the strain hardening and the strain softening zones for cementitious composites containing randomly oriented and distributed microfibers. These models have been implemented in a finite element code and a discrete cohesive fracture analysis has been performed. The obtained results are shown to be in a very good agreement with experimental data for all fiber volume fractions and the fiber types investigated in this study. - A more physically sound procedure would be to use a combined CDM-NLFM approach, where damage mechanics would be used up to peak load and N L F M for the prediction of the subsequent response. Such a procedure has been proposed through the introduction of an energy equivalence between a diffused deterioration like that predicted by a continuum damage model and a localized cut through the material like it would be predicted by a fracture mechanics analysis. A s expected, the results obtained here were in very good agreement with the experimental data. 225 Chapter 9 Conclusions and Recommendations - Although either N L F M or C D M concepts can be used to predict the entire structural response of F R C materials, the choice of which one to use depends on the expected results from the analysis. If one, is only interested in the global response of a structural component, either one is good, and i f the knowledge of the damaged zones has to be known at any stage of loading then, the combined approach provides the best alternative. If on the other hand one is interested in optimizing the reinforcing act of the fibers within the composite, then the micromechanics-based N L F M approach is the best alternative since it provides some quantification of the fiber efficiency through its effect the stress-COD constitutive law governing the opening of the cracks. A numerical homogenization procedure has been used to compute the evolution of the effective moduli and damage together with its associated thermodynamic force in the case of F R C materials. The generality of the finite element analysis allowed the investigation of parameters which are usually ignored or oversimplified in closed form solutions. Crack interaction and presence of fibers are two examples of such parameters. Such an investigation allows one to quantify the a priori assumptions used in the closed form solutions and help in the proper selection of potential damage variables for the construction of micromechanics-based continuum damage models. In summary one could say that damage evolution is far more sensitive to crack patterns than the effective moduli due to the nature of the two quantities. RECOMMENDATIONS FOR FUTURE RESEARCH - Refinement of the micromechanical constitutive laws would certainly yield a better understanding and quantification of the fiber reinforcing mechanism which could eventually lead to an accurate optimization of the fibers at the industrial level. A s we saw earlier, despite the use of relatively simple mechanisms to develop such models, the overall response is in very good agreement with the experimental data. 226 Chapter 9 Conclusions and Recommendations - In the N L F M predictions of the F R C beam structures it has been assumed that the opening and propagation of a cohesive crack occurring in the lower part of the beams is the main mechanism controlling the degradation of the structure while any deterioration due to the compressive loading has been neglected. This assumption, clearly more suitable for the case of the double cantiliver beam than for the beam bending case, has led to fairly accurate results in both cases. A little improvement, however, could be expected in the case of beam bending i f one incorporates the damage in the compressive part of the beam by considering the bulk F R C material to be elastic damageable through the use of some kind of P C D M instead of a simple linear elastic behavior as it has been assumed in this study. - A n interesting extension would be to determine the complete stiffness tensor for F R C materials. This actually should not be very difficult, since experimental results in the literature have shown that the shear behavior of F R C materials can be reduced to their tensile behavior since the cracking is controlled by the maximum principal tensile stress just like in the pure tensile case. Thus it should be only a matter of bases rotations and some mechanical judgement to get the full response. This could be useful for somebody who wants to use a traditional finite element code for instance to compute the behavior of an F R C structural component. - A consistent underprediction of the post peak response has been observed in the case of beam bending as well as the case of the contoured double cantiliver beam for all fiber types and fiber volume fractions investigated in this study. It is thought that such a consistent discrepancy is due to the introduction of errors during the identification of the model parameters characterizing the material at hand. Indeed, despite the use of a very small displacement loading rate (0.1 mm/mn) to record the direct tensile response of cementitious composites, the true post-peak behavior is very difficult to capture due to a lack of control over the propagation of cracking within the material leading to an underprediction of the material's response. This tensile response is subsequently used to identify the model parameters for both the P C D M and the N L F M approaches. Better predictions should be achieved with the use of more 227 Chapter 9 Conclusions and Recommendations accurate tests such as closed loop testing procedures, where crack opening within the material is used as a feed-back to control the loading rate to avoid any sudden fracture of specimens. - Over the last few years, durability of concrete has been identified as one of the major parameters controlling the long term behavior of concrete structures. Unfortunately, at the present stage of knowledge, durability design is practically independent from the structural design. Moreover, there is no rational design method for durability and only few rule of thumb prescriptions are being used and a rational integrated structural-durability design method needs to be developed. Such a method should lead to a performance-based design procedure. Durability of concrete structures in harsh environments has been shown to be controlled by the phenomena of mass transport within concrete, such as acids, sulfates etc. A n integrated durability-structural design procedure could be established with the development of a long term constitutive law that couples reactive flow in porous media with continuum damage mechanics. W i t h such a law, the combined effect of deterioration due to loading and that due to the environment could be estimated, thus allowing to track the evolution of the mechanical properties of the material with time. In the P C D M model proposed in this thesis for the modeling of F R C materials, it has been shown that the state of strain within a certain neighborhood of any material point can be used as the quantity piloting the evolution of damage due to loading. A similar quantity needs to be identified to characterize the evolution of damage due to the environmental attacks. It is believed that permeability of the material could be considered as a potential parameter for the characterization of damage due to environmental attacks. The combination of C D M with reactive flow in porous media should allow monitoring the evolution of mechanical properties of concrete structures with time. 228 Bibliography [1] D . H . A l l e n , C E . Harris, and S.E. Groves A Thermomechanical Constitutive Theory for Elastic Distributed Damage. I: Theoretical Development Int. J . Solids and Struct., V o l . 23, 1987, pp. 1310-1318 Composites with [2] S. Andrieux, Y . Bamberger, J.J. Marigo A model of micro-cracked material for concretes and rocks Journal of Theoretical and Applied Mechanics, V o l . 5, N 3, 1986 [3] J . H . Argyris, L . E . V a z and K . J . W i l l a m Numerical Techniques for the Finite Element Analysis of Elastic and Inelastic Material Nonlinearties in: Physical Non-Linearities in Structural Analysis, Symposium Senlis, France M a y 27-30, 1980, Published by Springer Verlag, Berlin Heidelberg N e w Y o r k 1981, J . Hult and J. Lemaitre (eds.), pp 6-16 [4] M . Arminjon, T. Chambard et S. Turgeman Homogeneisation d'un mortier avec une repartition aleatoire de fibres C. R . Acad. S c i . Paris, t. 316, Serie II, p. 1505-1510, 1993 [5] M . Arminjon, T. Chambard et S. Turgeman Variational M i c r o - M a c r o Transition, with Application to Reinforced Mortars Int. J. Solids Structures V o l . 31, N° 5, pp. 683-704, 1994, Elsevier Science L t d [6] N . Banthia, S. Mindess Fiber Reinforced Concrete: Modern Developments Second University-Industry Workshop on Fiber Reinforced Concrete, Toronto 1995 [7]. N . Banthia, and C . Y a n , " M i c r o - F i b e r Reinforced Cement-based Composite for Repairing," Journal of Concrete International, 1997. (in press) [8]. N . Banthia, C . Y a n , and S. Mindess, "Restrained Shrinkage Cracking i n Fiber Reinforced Concrete: A N o v e l Test Technique," Journal of Cement and Concrete Research, V o l . 26, N o . 1, 1996. pp. 9-14. [9]. N . Banthia, C . Y a n , and W . Y . Lee, "Restrained Shrinkage Cracking i n Fiber Reinforced Concrete with Polyolefin Fiber," Proceedings, Fifth International Conference on Structural Failure, Durability and Retrofit, Singapore, N o v . 1997. (in press) 229 [10]. N . Banthia, and C . Y a n , M i c r o - F i b e r Reinforced Cement-based Composite for Repairing," Proceedings, Rehabilitation and Repairing of Concrete Structures, American Concrete Institute, 1997. (in press) sv [11]. N . Banthia, N . Y a n , C . Y a n , C . Chan, and A . Bentur, ^ B o n d - s l i p Characteristics of Steel Micro-Fibers Reinforced Cement Composite Symposium on Bonding and Interfaces in Cementitious Materials, Materials Research Society, Boston, U . S . A . , V o l . 370, 1995, pp. 539-543. [12] N . Banthia, e t a l . Micro-fiber reinforced cement composites. I. Uniaxial tensile response. Cdn. J. O f C i v i l Engg, V o l 21, N o 6, 1994, pp. 999-1011. [13] N . Banthia A n d S. Dubeau Carbon and Steel Microfiber-Reinforced Cement-Based Composites for Thin Repairs. A S C E , J. O f Mat. In C i v i l Engg, V o l 6, N o . 1, 1994, pp. 88-99 [14] N . Banthia, J. Sheng and Y . Ohama Carbon Fiber Reinforced Cement Plates Under Transverse Loading [15] G . B a o a n d Y . Song Crack Bridging Models for Fiber Composites with Slip-Dependent Interfaces J. Mech. Phys. Solids V o l . 41, N o 9, pp 1425-1444,1993 [16] G.I. Barenblatt Micromechanics of Fracture Theoretical and Applied Mechanics 1992, S.R. Bodner , J. Singer, A . Solan & Z . Hashin (Editors) Elsevier Science Publishers B . V . 1993 I U T A M [17] R . S . Barsoum On the use of Isoparametric Finite Elements in Linear Fracture Mechanics Int. J. N u m . Meth. Engng., V o l . 10, 1976, pp. 25-38 [18] J . L Batoz and G . Dhatt "Modelisation des structures par elements finis" Volume 2: poutres et plaques, Les presses de l'universite Laval, Quebec, 1990 [19]Z.P. Bazant,J. Mazars France-U.S. Workshop on Strain Localization and Size Effect due to Cracking and Damage Journal of Engineering Mechanics, V o l . 116, n° 6, June, 1990 230 [20] Z . P Bazant and G . Pijaudier-Cabot "Measurement of characteristic length of non local continuum" J. Eng. Mech., A S C E , V o l . 115, n° 755-767 (1989) [21] Z.P. Bazant W h y Continuum Damage is Nonlocal: Micromechanics Argument Journal of Engineering Mechanics, vol. 117, n° 5, M a y , 1991 [22] Z.P. Bazant Instability, ductility and size effect in strain softening concrete J. E n g . M e c h . D i v . ( A S C E ) V o l . 102, 1976, pp. 331-344 [23] Z.P. Bazant "Numerical simulation of progressive fracture in concrete structures" :Recent developments, Computer A i d e d Analysis and Design of Concrete Structures, (Ed.F.Damjanic), Pineridge Press, Swansea, Part 1, pp 1-17 (1984) [24] Z . P Bazant & B . H . O h "Crack band theory for fracture of concrete" Materials and Structures, R i l e m 16 (1983), 155-177 [25] Z.P. Bazant and M . T . Kazemi Determination of Fracture Energy, Process Zone Length and Brittleness Number from Size Effect, with Application to Rock and Concrete International Journal of Fracture, V o l . 44, 1990, pp. 111-131 [26] T. Belytschko and D . Larsy "Localization limiters and numerical strategies for strain softening materials" Damage and Cracking, Strain Localization and Size effect, (Eds J. Mazars & Z.P. Bazant) Elsevier, Amsterdam, pp. 349-362 (1989) [27] A . Benallal, R. Billardon, G . Geymonat "Some mathematical aspects of the damage softening rate problem" Damage and Cracking, Strain Localization and Size Effect (Eds J. Mazars and Z.P. Bazant), Elsevier, Amsterdam, pp 247-258 (1989) [28] A . Benallal, C . C o m i , and J. Lemaitre Critical Damage States at Crack Initiation A M D - V o l . 1 4 2 / M D - V o l . 34, Damage Mechanics and Localization A S M E 1992 [29] A . Bensoussan, J.L. Lions and G . Papanicolaou Asymptotic Analysis for Periodic Structures Studies in Mathematics and its Applications, North Holland, 1978 231 [30] A . Bentur, S. Mindess Fibre Reinforced Cementitious Composites Elsevier Applied Science, 1990 [31] D . J . Bergman and Y . Kantor Improved Rigorous Bounds on the Effective Elastic M o d u l i of a Composite Material J. M e c h . Phys. Solids, vol 32, no 1, 1984, pp 41-62 [32] Y . Berthaud " A n ultrasonic testing method" : A n aid for material characterization, Int. Conf. Fract. Concrete Rock, Houston (1987) [33] M . Berveiller, J. Zaoui Methodes self-consistantes en mecanique des solides heterogenes C R . 15 C o l l . du groupe Francais de Rheologie, Paris 1980 [34] J. Betten Damage tensors in continuum mechanics Journal de Mecanique theorique et applique, V o l . 2, n° 1, 1983, pp. 13-32 [35] T . N . Bittencourt Computer simulation of linear and nonlinear crack propagation in cementitious materials Ph.D. Thesis, Cornell University, M a y 1993 [36] P. Bocca, A . Carpinteri, and S. Valente Effects in the M i x e d M o d e Crack Propagation: Softening and Snap-Back Analys Engng Fract M e c h , V o l . 26, 1987, pp. 185-201 [37] P. Bompard, D . W e i , T. Guennouni and D.Frangois Mechanical and fracture behavior of porous materials Engineering Fracture Mechanics, V o l . 28, N ° 5/6, pp. 627-642, 1987 [38] A . F . Bower and M . Ortiz A Three-Dimensional Analysis of Crack Trapping and Bridging by Tough Particles J. M e c h . Phys. Solids V o l . 39, N o 6, pp 815-858,1991 [39] B . Budiansky On the Elastic M o d u l i of Some Heterogeneous Materials J. M e c h . Phys. Solids, vol 3, 1965 [40] L . Cedolin and S. D e i Poli Finite Element Studies of Shear Critical r/c beams Engng Mechanics, A S C E , 103 (3), 1977, pp. 395-410 232 [41] J . L . Chaboche Phenomenological Aspects of Continuum Damage Mechanics Theoretical and Applied Mechanics. P. Germain, M . Piau and D . Caillerie (eds.) Elsevier Science Publishers B . V . (North-Holland), I U T A M 1989 [42] J . L . Chaboche, J. Lemaitre, D . Marquis, and S. Savalle Discussion on Problems of Models Identification in: Physical Non-Linearities in Structural Analysis, Symposium Senlis, France M a y 27-30, 1980, Published by Springer Verlag, Berlin Heidelberg N e w Y o r k 1981, J. Hult and J. Lemaitre (eds.), pp 37-51 [43] J.L. Chaboche, N . E l Mayas et P. Paulmer Modelisation fhermodynamique des phenomenes de viscoplasticite, restauration et vieillissement C . R . Acad. S c i . Paris, t. 320, Serie II, p. 9-16, 1995 [44] S.K. Chan, I.S. Tuba, and W . K . W i l s o n O n the Finite Element Method i n Linear Fracture Mechanics Engng. Fract. Mech., V o l . 2, (1970), pp. 1-17 [45] Y . C Chiang, A . S . D . Wang and T . W . Chou On Matrix Cracking in Fiber Reinforced Ceramics J. M e c h . Phys. Solids V o l . 41, N o 7, pp 1137-1154,1993 [46] R . M . Christensen Mechanics of Composite Materials John W i l e y & Sons, 1979 [47] R . F . Cook, C J . Fairbanks, B . R . L a w n , and Y . W . M a i Crack Resistance by Interfacial Bridging: Its role in Determining Strength Characterization Journal for Materials Research, V o l . 2, N o . 3, 1987, pp. 345-356 [48] J.P. Cordebois et F . Sidoroff Endommagement Anisotrope en Elasticite et Plasticite Journal de Mecanique theorique et appliquee, V o l . 3, n° Special, 1982, pp. 45-60 [49] B . N . C o x and D . B . Marshall Concepts in the Fracture and Fatigue of Bridged Cracks A C T A Metallurgica et Materialia, 1994 [50] W . A . Curtin Fiber Pull-Out and Strain Localization in Ceramic Matrix Composites J. M e c h . Phys. Solids V o l . 41, N o 1, pp 35-53, 1993 233 [51] P. De Buhan et S. Maghous Une methode numerique pour la determination du critere de resistance macroscopique de materiaux heterogenes a structure periodique C. R. Acad. Sci. Paris, t. 313, Serie II, p. 983-988, 1991 [52] H . Deng and S. Nemat-Nasser Microcrack Arrays in Isotropic Solids Mech. Mater., V o l . 13, 1992, pp.15-36 [53] F. Devries, H . Dumontet, G . Duvaut and F . Lene Homogenization and Damage for Composite structures International Journal for Numerical Methods in Engineering, V o l . 27, pp. 285-298 (1989) [54] P. Droz Numerical Modeling of the Nonlinear Behavior of Massive Unreinforced Concrete Structures Doctoral Thesis N o . 682, Swiss Federal Institute of Technology, Lausanne, 1987 (in Freeh) [55] D . C . Drucker Yielding, F l o w and Failure in: Inelastic Behaviour of Composite Materials, ed. Carl T. Herakovich, 1975, pp 1-15 [56] G . Duvaut Homogeneisation et materiaux composites in: Trends and application of mathematics to mechanics, eds. P . G . Ciarlet et M . Roseau [57] L . Elfgren Fracture Mechanics of Concrete Structures, from Theory to Applications Chapman and H a l l , London, 1989 [58] F. Erdogen, and G . C . Sih On the Crack Extension in Plates under Plain Loading and Transverse Shear J. Basic Engineering, V o l . 8 5 , 4 , pp. 519-525 (1963) [59] J . D . Eshelby The Determination of an Elastic Field of an Ellipsoidal inclusion and related problems. Proc. R o y . Soc. London, 1957 234 [60] G . Francfort and E . Murat Homogenization and Optimal Bounds in Linear Elasticity A r c h . Rat. M e c h . A n d Analysis [61] R . M . Foote , Y . W . M a i , and B . Cotterell Crack Growth Resistance Curves in Strain Softening Materials Journal of Mechanics and Physics of Solids, V o l . 34, N o . 6, 1986, pp. 593-607 [62] I. Genois Fracture Resistance of Micro-Fiber Reinforced Cement Composites M A . S c . Thesis, University of British Columbia, 1995 [63] P. Germain Sur Quelques Concepts Fondamentaux de la Mecanique in: Large Deformations of Solids: Physical Basis and Mathematical Modelling Esevier Applied Science London and New York, edited by. J. Gittus, J. Zarka, S. Nemat-Nasser pp 3-25 [64] P. Gilormini Insuffisance de l'extension classique du modele autocoherent au comportement non lineaire C. R. Acad. S c i . Paris, t. 320 Serie II, p. 115-122, 1995 [65] A . A . Griffith The phenomena of rupture and flow in solids Philosophical Transactions of the Royal Society of London A 2 2 1 , 163-198 (1921) [66] A . A . Griffith The theory of rupture Proceedings of First International Congress of Applied Mechanics Delft, pp. 55-63 (1924) [67] K . Golden and G . Papanicolaou Bounds for Effective Parameters of Heterogeneous M e d i a by Analytic Continuation C o m m . Math. Physics, vol. 90, 1983, p 473 [68] B . Halphen et N . Q . Son Sur les materiaux standards generalises Journal de Mecanique, V o l . 14, n° 1, 1975, pp.39-63 [69] Z . Hashin The Elastic M o d u l i of Heterogeneous Materials Journal of Applied Mechanics, M a r c h 1962, pp. 143-150 235 [70] Z . Hashin Extremum Principals for Elastic Heterogeneous M e d i a with Imperfect Interfaces and their Application to Bounding of Effective M o d u l i J. M e c h . Phys. Solids V o l . 40, N o 4, pp 767-781,1992 [71] Z . Hashin, B . W . R o s e n The Elastic M o d u l i of Fiber Reinforced Materials J. Applied Mechanics Trans of A S M E , vol 31, (1964) [72] Z . Hashin, S. Shtrikman A variational Approach of the Theory of Elastic Behavior of Multiphase Materials. J. M e c h . Phys. Solids, vol 32, n 1, 1984 [73] R . D . Henshell and K . G . Shaw Crack tip elements are unnecessary. Int J. N u m . Meth. Engng., V o l . 9, 1975, pp. 495-509 [74] E . Herve and A . Zaoui Modeling the effective behavior of non-linear matrix-inclusion composites Eur. J. M e c h . , A/Solids, 9, n« 6, 505-515, 1990 [75] R. H i l l "Theory of Mechanical Properties of Fiber-Strengthened Materials: I Elastic Behavior" J. M e c h . Phys. Solids, V o l . 12, p 199 (1964) [76] R. H i l l A Self Consistent Mechanics of Composite Matrials J. M e c h . Phys. Solids, 13, 1965 [77] R. H i l l Continuum Micro-Mechanics of Elastoplastic Polycristals J. M e c h . Phys. Solids, V o l 13, 1965 [78] A . Hillerborg, M . Modeer, and P . E . Peterson Analysis of Crack Formation and Crack Growth in Concrete by Means of Fracture Mechanics and Finite Elements Cement and Concrete Research, V o l . 6, N o . 6, 1976, pp. 773-782 [79] J.D. Hoffman Numerical Methods for Engineers and Scientists New York, M c G r a w - H i l l , 1992 236 [80] J . H . Huang, R . Furuhashi and T. M u r a Frictional Sliding Inclusions J. Mech. Phys. Solids V o l . 41, N o 2, pp 247-265,1993 [81] M . A . Hussain, S. L . Pu and, J. H . Underwood Strain Energy Release Rate for a Crack under Combined M o d e I and M o d e II Fracture Analysis, A S T M , STP560, pp. 2-28 (1974) [82] A . R. Ingraffea, V . Saouma Numerical Modeling of discrete crack propagation in reinforced and plain concrete [83] A . R. Ingraffea, T . N . Bittencourt, and J . L . A . O Sousa Automatic Fracture Propagation for 2 D Finite Element Models Proc. M E C O M 90, R i o de Janeiro, 1990 [84] C . E . Inglis Stresses in a plate due to the presence of cracks and sharp corners Transactions of the Institution of Naval Architects, v o l . 55, 1913, pp. 219-241 [85] G . R . Irwin Plastic zone near a crack and fracture toughness In Proceedings of the Seventh Sagamore Ordnance Materials Conference Syracuse University, 1960, pp. IV-63 [86] G . R . Irwin Analysis of stresses and strains near the end of a crack traversing a plate J. A p p l . M e n . , 24, (1957), pp. 361-354 [87] Y . S . Jenq and S.P. Shah Fracture Toughness Criterion for Concrete Engng Faract. Mech., V o l . 21, 1985, pp. 1055-1069 [88] Y . S . Jenq and S.P. Shah A T w o Parameter Fracture M o d e l For Concrete Journal of Engineering Mechanics, V o l . I l l , N o . 4, 1985, pp. 1227-1241 [89] R. John and S.P. Shah Fracture Mechanics Analysis of H i g h Strength Concrete Journal of Materials in C i v i l Engineering, A S C E , V o l . 1, N o . 4, 1989, pp. 185-198 [90] L . M . Kachanov "Time of rupture process under creep conditions". Tzv, A k a d . Nauk. S S R Otd. Teck. Nauk., n° 8 1958, pp. 26-31 237 [91] L . M . Kachanov A Continuum M o d e l of M e d i u m with Cracks J. Engng Mech., V o l . 106, 1980, pp. 1039-1051 [92] L . M . Kachanov Effective Elastic Properties of Cracked Solids: Critical Review of Some Basic Concepts A p p l . M e c h . Rev., V o l . 45, 1992, pp. 304-335 [93] B . L . Karihaloo and P. Nallathambi Notched Beam Test: M o d e I Fracture Toughness In RJJJEM Report 5, Fracture Mechanics Test Methods for Concrete, ed. S.P. Shah, and A . Carpinteri, Chapman & H a l l , London 1991, pp. 2831-2836 [94] M . P. Kanninen and C . H . Popelar Advanced Fracture Mechanics Oxford University Press, N e w Y o r k , 1985 [95] M . K . Kassir and G . C . Sih Grifith's theory of brittle fracture in three dimensions International Journal of Engineering Sciences 5, 899-918 (1967) [96] A . S . Kobayashi, N . M . Hawkins, and R . C . Bradt Facture Process zone in concrete and Ceramics: A matter of scaling In Toughening Mechanisms in Quasi-Brittle Materials, ed. S.P. Shah, Kluer Academic, Dordrecht, (1991), pp. 35-46 [97] A . S . Kobayashi et al. Fracture Process Zone of Concrete In Application of Fracture Mechanics to Cementitious Composites, ed. S.P. Shah, Martinus Nijhoff Publ., Dordrecht, 1985, pp. 25-50 [98] E . Kroner Linear Properties of Random Media: The systematic Theory C R . 15 C o l l . D u Group Francais de Rheologie, Paris 1980 [99] E . Kroner The Statistical Basis of Polycristal Plasticity in: Large Deformations of Solids: Physical Basis and Mathematical Modeling Elsevier Applied Science London and N e w York, edited by J. Gittus, J. Zarka, S. Nemat-Nasser pp. 27-40 [100] T . E . Lacy, D . L . M c D o w e l l , P . A . W i l l i c e , and R. Talreja On Representation of Damage Evolution in Continuum damage Mechanics Int. J. Damage Mech., V o l . 6, January 1997, pp. 62-94 238 [101] F . A . Leckie Tensorial Nature of Damage Measuring Internal Variables in: Physical Non-Linearities in Structural Analysis, Symposium Senlis, France M a y 27-30, 1980, Published by Springer Verlag, Berlin Heidelberg N e w Y o r k 1981, J. Hult and J. Lemaitre (eds.), pp 140-155 [102] J. Lemaitre, J . L . Chaboche Mechanics of Solid Materials Cambridge University Press, 1990 [103] J. Lemaitre et J. Mazars Application de la theorie de l'endommagement au comportement non lineaire et a la rupture du beton de structure Annales de l'institut technique du batiment et des travaux publics, n° 401, Janvier 1982, Theorie et Methodes de Calcul 24, pp. 114-138 [104] F. Lene and P. Paumelle Micromechanisms of Damage in Woven Composite Composite Material Technology, P D - V o l . 45, A S M E 1992, pp. 97-105 [105] C . K . Y . Leung and V . C . L i Effect of Fiber Inclination on Crack Bridging Stress in Brittle Fiber Reinforced Brittle Matrix Composites J. Mech. Phys. Solids V o l . 40, N o 6, pp 1333-1362, 1992 [106] A . J . Levy The Debonding of Elastic Inclusions and Inhomogeneities J. M e c h . Phys. Solids V o l . 39, N o 4, pp 477-505, 1991 [107] V . C . L i Performance Driven Design of Fiber Reinforced Cementitious Composites in: Fiber Reinforced Cement and Concrete. Edited, by R . N . Swamy [108] V . C . L i , Y . Wang, S. Backer A micrmechanical model of tension softening and bridging toughning of short random fiber reinforced brittle composites J. M e c h . Phys. Solids V o l . 39, N 5, 1991 [109] J . L Lions Perturbations singulieres dans les problemes aux limites et en controle optimale Lecture notes in Mathematics 323, Springer Verlag (1973) [110] A . Litewka, A . Sawczuk, J. Stanilawska Simulation of oriented continuous damage evolution Journal de Mecanique theorique et appliquee, V o l . 3, n° 5, 1984, pp. 675-688 239 [ I l l ] J.J. Marigo Formulation d'une loi d'endommagement d'un materiau elastique C . R. Acad. S c i . Paris, t. 292, Serie II, p. 1309-1312, 1981 [112] J.J. Marigo U n modele d'endommagement de materiau elastique fragile par croissances de microvides via l'homogeneisation H I 4426.07 E D F - D E R 1983 [113] S.J. Matysiak and C z . Wozniak On the Microlocal Modeling of the Thermoelastic Periodic Composites Journal of Technical Physics, J. Tech. Phys., 29, 1, 85-97, 1988 [114] J. Mazars Application de la mecanique de l'endommagement au comportement non lineaire et a la rupture du beton de structures These de doctorat d'etat, Universite Paris V I , 1984. [115] J. Mazars A Description of M i c r o - and Macroscale Damage of Concrete Structures Engineering Fracture Mechanics, V o l . 25, n° 5/6, pp. 729-737 , 1986 [116] J. Mazars, G . Pijaudier-Cabot and C . Saouridis Size effect and continuous damage in cementitious materials International Journal of Fracture 51: 159-173, 1991 [117] J. Mazars and G . Pijaudier-Cabot Continuum Damage theory - Application to Concrete Journal of Engineering Mechanics, vol. 115, n° 2, February, 1989 [118] J. Mazars Damage Models For Concrete and their Usefulness for Seismic Loadings in: Experimental and Numerical Methods in Earthquake Engineering, J. Donea and P . M . Jones (eds.), 199-221. [119] J. Mazars and G . Pijaudier-Cabot From damage to fracture mechanics and conversely - A combined approach in: Proceedings of Engineering Mechanics, v l 1995, A S C E , N e w Y o r k , N . Y . , U S A 1995, pp 231-234 [120] Y . Y . M a i , B . Barakat, B . Cotterell, and M . V . Swain R-curve Behavior in a Macro-Defect-Free Cement Paste Philosophical Magazine A , V o l . 62, N o . 3, 1990, pp. 347-361 240 [121] G . A . Maugin The Thermomechanics of Plasticity and Fracture Cambridge Texts in Applied Mathematics Cambridge University Press. [122] S. Mindess The Fracture of Fiber Reinforced and Polymer impregnated Concretes: A Review In Fracture Mechanics of Concrete, ed. F . H . Wittmann. Elsevier Science Publishers, Amsterdam, 1983, pp. 481-501 [123] S. Mindess Fiber Reinforced Concrete: Challenges and Prospects In Fiber Reinforced Concrete: Modern Developments, Second UniversityIndustry Workshop on Fiber Reinforced Concrete, Toronto 1995, pp. 1-11 [124] W . Muschik Thermodynamical Constitutive Laws - Outlines in: Constitutive L a w and Microstructure: Proceedings of the Seminar, Wissenschaftskolleg - Institute for Advanced Study, Berlin West Germany, February 23-24, 1987 Axelrad and Mushik (eds.), pp. 3-25 [125] S. Murakami, M . K a w a i and H . Rong Finite Element Analysis of creep crack growth by a local approach Int. J. M e c h . S c i . , 30 1988, pp. 491-502 [126] N.I. Muskhelishvili Some basic problems of the mathematical theory of elasticity Nordhorff, 1953 [127] S. Nemat-Nasser, M . H o r i Micromechanics: Overall Properties of Heterogeneous Materials North-Holland, 1993 [128] D . N g o and A . C . Scordelis Finite Element Analysis of Reinforced Concrete Beams J American Concrete Inst, V o l 64, 1967, pp. 152-163 [129] L . N i and S. Nemat-Nasser Interface Cracks in Anisotropic Dissimilar Materials: A n Analytic Solution J. M e c h . Phys. Solids V o l . 39, N o 1, pp 113-144, 1991 [130] S. Nirmalendran, H . Horii Analytical modeling of microcracking and bridging in the fracture of quasi-brittle materials. J. M e c h . Phys. Solids V o l 40, N 4, 1992 241 [ 131 ] Y . Okui, H . H o r i i and N . A k i y a m a A Continuum Theory for Solids Containing Microdefects Int. J. Engng S c i . V o l . 31, N ° 5, pp. 735-749, 1993 [132] M . Oritz " A n analytical study of the localized failure modes of concrete" M e c h . Mat., 6, pp. 159-174 (1987) [133] C . Ouyang, B . Mobasher, and S.P. Shah A n R-curve Approach for Fracture of Quasi-brittle Materials Engng. Fract. Mech., V o l . 37, N o . 4, 1990, pp. 901-916 [134] P . E . Petersson Crack growth and Development of fracture zones in concrete and similar materials Technical Report T V B M 1006, L u n d Institute of Technology, L u n d Sweden, December 1981 [135] G . Pijaudier-Cabot, Z . P . Bazant Localization Limiting Properties of Nonlocal Damage Models A M D - V o l . 1 4 2 / M D - V o l . 34, Damage Mechanics and Localization A S M E 1992 [136] G . Pijaudier-Cabot & Z . P . Bazant "Non local damage theory". J. of Eng. Mech., A S C E , V o l . 113, n°10, pp. 287-293, 1988 [137] G . Pijaudier-Cabot, J. Mazars, J. Pulikowski Steel-Concrete B o n d Analysis with Nonlocal Continuous Damage Journal of Structural Engineering, V o l . 117, n° 3, March, 1991 [138] S. D . Poisson Second mem. Sur la theorie du magnetisme M e m . De l ' A c a d . De France 5, 1822 [139] P. Ponte Castaheda The Effective Mechanical Properties of Nonlinear Isotropic Composites J. M e c h . Phys. Solids V o l . 39, N o 1, pp 45-71, 1991 [140] Y . R . Rashid Ultimate Strength Analysis of Prestressed Concrete Pressure Vessels Nuclear Engng & Design, V o l . 7,1968, pp. 334-344 [141] J.R. Rice Inelastic Constitutive Relations for Solids: A n Internal-Variable Theory and its Application to Metal Plasticity J. Mech. Phys. Solids, V o l . 19, 1971, pp. 433-455 242 [142] J . G . Rots Computational Modeling of Concrete Fracture Ph.D. Thesis, The technical University of Delft, September 1988 [143] A . A . Rubinstein and K . X u Micromechanical M o d e l of Crack Growth in Fiber-Reinforced Ceramics J. M e c h . Phys. Solids V o l . 40, N o 1, pp 105,125, 1992 [144] K . Sakai Integrated Design and Environmental Issues in Concrete Technology Proceedings of the International Workshop "Rational Design of Concrete Structures under Severe Conditions" Hakodate, Japan, 1995 [145] E . Sanchez-Palencia Comportement local et macroscopique d'un type de milieux physiques heterogenes Int. J. Engng Sci., 1974, V o l . 12, pp. 331-351, Pergamon Press. [146] V . E . Saouma Interactive Finite Element Analysis of Reinforced Concrete: A Fracture Mechanics Approach Ph.D. Thesis, Cornell University, 1981 [147] V . E . Saouma and A . R . Ingraffea Fracture Mechanics Analysis of Discrete Cracking Proc. I A B S E Colloquium on Advanced Mechanics of Reinforced Concrete, Delft, June 1981, pp. 393-416 [148] C . Saouridis and J. Mazars Prediction of the failure and size effect in concrete v i a bi-scale damage approach Engineering Computations, V o l . 9, 1992 [149] C . Saouridis Identification et numerisation objectives des comportements adoucissants: Une approche multi-echelle de l'endommagement du beton These de Doctorat, Universite Paris V I , 1988 [150] J.C.J Schellekens Interface elements in finite element analysis TU-Delft report N o 25-2-90-5-17, T N O institute for building materials and structures, Delft, the Netherlands 1990 243 [151] E . Shrader Mistakes, Misconceptions, and Controversial Issues Concerning Concrete and Concrete Repairs. Concrete International, Sept. 1992, pp. 52-56. [152] F . Sidoroff Description of Anisotropic Damage. Application to Elasticity in: Physical Non-Linearities in Structural Analysis, Symposium Senlis, France M a y 27-30, 1980, Published by Springer Verlag, Berlin Heidelberg N e w Y o r k 1981, J. Hult and J. Lemaitre (eds.), pp 237-244 [153] G . C . S i h and H . Liebowitz Mathematical Theories of Brittle Fracture In Fracture-An Advanced Treatise, V o l . II, Mathematical Fundamentals (ed. H . Liebowitz), Pergamon Press, pp. 67-190 (1968) [154] G . C . S i h Strain Energy Density Factor Applied to M i x e d M o d e Crack Problems Int. J. Fracture Mech., V o l . 10, 3, (1974), pp. 305-321 [155] G . C . Sih and B . MacDonald Fracture Mechanics Applied to Engineering ProblemsStrain Energy Density Criterion Engng Fracture M e c h , V o l . 6, (1974), pp. 361-386 [156] C . Stolzet A . Zaoui Analyse morphologique et approches variationnelles du comportement d'un milieu elastique heterogene C . R . Acad. Sci. Paris, t. 312, Serie II, p. 143-150, 1991 [157] P. Suquet On Bounds for the Overall potential of Power L a w Materials Containing Voids with Arbitrary Shape Mechanics Research Communications V o l . 19(1), pp. 51-58, 1992 [158] D . V . Swenson Modeling M i x e d - M o d e Dynamic Crack Propagation Using Finite Elements Ph.D. Thesis, Cornell University, Dec. 1985 [159] D . V . Swenson and A . R . Ingraffea A Finite Element M o d e l of Dynamic Crack Propagation with an Application to Intersecting Cracks Proc. 4 Int Conf on Numerical Methods in Fracture Mechanics, ed. Luxmoore et al., 1974, pp. 493-507 th 244 [160] H . Tada, P . C . Paris, and G . R . Irwin The Stress Analysis of Cracks Handbook 2 ed., Paris Productions, St. Louis, M O , 1985 n d [161] R. Talreja Internal Variable Damage Mechanics of Composite Materials In Yielding Damage and Failure of Anisotropic Solids, ed. J.P. Bohler, London: Mechanical Engineering Publications, 1990, pp. 509-533 [162] R. Talreja Continuum Modeling of Damage in Ceramic Matrix Composites M e c h . Mater., V o l . 12, 1991, pp. 165-180 [163] L . Tartar Estimations fines des coefficients homogeneises in Ennio de Giorgi Colloquium, 1985, Pittman Press, P. Kree E d . [164] L . Tartar Etude des oscillations dans les equations aux derivees partielles non lineaires in: Trends and application of mathematics to mechanics, eds. P . G . Ciarlet et M . Roseau [165] L . Tartar Memory effects and Homogenization Archive for Rational Mechanics and Analysis, v o l . 111, n° 2, 1990, pp. 121-133 [166] P. Tong, and T . H . H . Plan On the Convergence of the Finite Element Method for Problems with Singularities Int. J Solids & Struct., V o l . 9, 1973, pp. 313-321 [167] D . M . Tracy Finite Elements for Determination of Crack T i p Elastic Stress Intensity Factors Engng Fracture Mech., V o l . 3, 1971, pp.255-256 [168] V . Tvergaard M o d e l Studies of Fiber Breakage and Debonding in a Metal Reinforced by Short Fibers J. M e c h . Phys. Solids V o l . 41, N o 8, pp 1309-1326,1993 [169] P . G . Underwood Dynamic relaxation In computational methods for transient analysis, (ed. T. Belytshko and T.J.R. Hughes), elsevier science publishers, 1983, pp. 245-265 245 [170] A.S. Vavakin, R.L. Salganik Effective Elastic Characteristics of body with Isolated Cracks, Cavities and Rigid Non-Homogeneities. Izu and USSR Mekhanik Tvergoda Tela, vol 13, no 2, 1978 [171] Y. Wang, V.C. Li, and S. Backer Analysis of Synthetic fiber Pull-out from a cement matrix In, Bonding in Cementitious Composites, eds. S. Mindess ans S.P. Shah, Material Research Society, Pittsburg 1988, pp. 159-65 [172] Y. Weitsman and H. Zhu Multi-Fracture of Ceramic Composites J. Mech. Phys. Solids Vol. 41, No 2, pp 351-388,1993 [173] Y. Weistman Damage Coupled with Heat Conduction in Uni-Axially Reinforced Composites Eds. V. Stokes and D. Krajcinovic, AMD-Vol. 85, New York: ASME, 1987, pp. 161-174 [174] J.R. Willis On Methods for Bounding the Overall Properties of Nonlinear Composites J. Mech. Phys. Solids Vol. 39, No 1, pp 73-86,1991 [175] J.R. Willis Theoretical Determination of the Overall properties of Composite Materials 1UTAM (International Union of Theoretical and Applied Mechanics) 1980 [176] J.R. Willis Bounds and Self-Consistent Estimates for the Overall Properties of Anisotropic Composites J. Mech. Phys. Solids, vol 25, 1977 [177] Yingqing Lawrence Cui Interaction of Fiber and Transformation Toughening J. Mech. Phys. Solids Vol. 40, No 8, pp 1837-1850,1992 246
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Constitutive modeling of fiber reinforced cement composites
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Constitutive modeling of fiber reinforced cement composites Boulfiza, Mohamed 1998
pdf
Page Metadata
Item Metadata
Title | Constitutive modeling of fiber reinforced cement composites |
Creator |
Boulfiza, Mohamed |
Date Issued | 1998 |
Description | The role of fibers in the enhancement of the inherently low tensile stress and strain capacities of fiber reinforced cementitious composites (FRC) has been addressed through both the phenomenological, using concepts of continuum damage mechanics, and micromechanical approaches leading to the development of a closing pressure that could be used in a cohesive crack analysis. The observed enhancements in the matrix behavior is assumed to be related to the ability of the material to transfer stress across cracks. In the micromechanics approach, this is modeled by the introduction of a nonlinear closing pressure at the crack lips. Due to the different nature of cracking in the pre-peak and post peak regimes, two different micro-mechanical models of the cohesive pressure have been proposed, one for the strain hardening stage and another for the strain softening regime. This cohesive pressure is subsequently incorporated into a finite element code so that a nonlinear fracture analysis can be carried out. On top of the fact that a direct fracture analysis has been performed to predict the response of some FRC structural elements, a numerical procedure for the homogenization of FRC materials has been proposed. In this latter approach, a link is established between the cracking taking place at the meso-scale and its mechanical characteristics as represented by the Young's modulus. A parametric study has been carried out to investigate the effect of crack patterning and fiber volume fractions on the overall Young's modulus and the thermodynamic force associated with the tensorial damage variable. After showing the usefulness and power of phenomenological continuum damage mechanics (PCDM) in the prediction of FRC materials' response to a stimuli (loading), a combined PCDM – NLFM¹ approach is proposed to model (predict, forecast) the complete response of the composite up to failure. Based on experimental observations, this approach assumes that damage mechanics which predicts a diffused damage is more appropriate in the pre-peak regime whereas, NLFM is more suitable in the post-peak stage where the opening and propagation of a major crack will control the response of the material and not a deformation in a continuum sense as opposed to the pre-cracking zone. Tensile and compressive tests have been carried out for the sole purpose of calibrating the constitutive models proposed and/or developed in this thesis for FRC materials. The suitability of the models in predicting the response of different structural members has been performed by comparing the models' forecasts with experimental results carried out by the author, as well as experimental results from the literature. The different models proposed in this thesis have the possibility to account for the presence of fibers in the matrix, and give fairly good results for both high fiber volume fractions (v[sub f] > 2%) and low fiber volume fractions (v[sub f] < 2%). Use of interface elements in a finite element code has been shown to be a powerful tool in analyzing the behavior of concrete substrate- FRC repair materials by the introduction of a zero thickness layer of interface elements to account for the interface properties which usually control the effectiveness of the repair material. |
Extent | 17778605 bytes |
Genre |
Thesis/Dissertation |
Type |
Text |
File Format | application/pdf |
Language | eng |
Date Available | 2009-06-02 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
DOI | 10.14288/1.0050167 |
URI | http://hdl.handle.net/2429/8580 |
Degree |
Doctor of Philosophy - PhD |
Program |
Civil Engineering |
Affiliation |
Applied Science, Faculty of Civil Engineering, Department of |
Degree Grantor | University of British Columbia |
Graduation Date | 1998-05 |
Campus |
UBCV |
Scholarly Level | Graduate |
Aggregated Source Repository | DSpace |
Download
- Media
- 831-ubc_1998-271110.pdf [ 16.95MB ]
- Metadata
- JSON: 831-1.0050167.json
- JSON-LD: 831-1.0050167-ld.json
- RDF/XML (Pretty): 831-1.0050167-rdf.xml
- RDF/JSON: 831-1.0050167-rdf.json
- Turtle: 831-1.0050167-turtle.txt
- N-Triples: 831-1.0050167-rdf-ntriples.txt
- Original Record: 831-1.0050167-source.json
- Full Text
- 831-1.0050167-fulltext.txt
- Citation
- 831-1.0050167.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0050167/manifest