1*1 National Library of Canada Acquisitions and Bibliographic Services Branch 395 Wellington Street Ottawa, Ontario K1A0N4 Bibliotheque nationale du Canada Direction des acquisitions et des services bibiiographiques 395, rue Wellington Ottawa (Ontario) K1A0N4 Your We Voire rdlironco Our Ha Noire leliirence NOTICE The quality of this microform is heavily dependent upon the quality of the original thesis submitted for microfilming. Every effort has been made to ensure the highest quality of reproduction possible. If pages are missing, contact the university which granted the degree. Some pages may have indistinct print especially if the original pages were typed with a poor typewriter ribbon or if the university sent us an inferior photocopy. AVIS La qualite de cette microforme depend grandement de la qualite de la these soumise au microfilmage. Nous avons tout fait pour assurer une qualite superieure de reproduction. S'il manque des pages, veuillez communiquer avec I'universite qui a confere le grade. La qualite d'impression de certaines pages peut laisser a desirer, surtout si les pages originates ont ete dactylographiees a I'aide d'un ruban use ou si I'universite nous a fait parvenir une photocopie de qualite inferieure. La reproduction, meme partielle, de cette microforme est soumise a la Loi canadienne sur le droit d'auteur, SRC 1970, c. C-30, et ses amendements subsequents. Reproduction in full or in part of this microform is governed by the Canadian Copyright Act, R.S.C. 1970, c. C-30, and subsequent amendments. Canada TWO-DIMENSIONAL BOUNDARY ELEMENT ANALYSIS OF SEISMIC CRACKING IN CONCRETE GRAVITY DAMS Vinod Batta A Thesis in The Department of Civil Engineering Presented in Partial Fulfillment of the Requirements for the Degree of Doctor of Philosophy at Concordia University Montreal, Quebec, Canada February 1995 © Vinod Batta, 1995 1*1 National Library of Canada Acquisitions and Bibliographic Services Branch 395 Wellington Street Ottawa, Ontario K1A0N4 Bibliotheque rationale du Canada Direction des acquisitions et des services bibliographiques 395, fue Wellington Ottawa (Ontario) K1A0N4 your tih Voti* uttmnct Our tile Notre relertoca THE AUTHOR HAS GRANTED AN IRREVOCABLE NON-EXCLUSIVE LICENCE ALLOWING THE NATIONAL LIBRARY OF CANADA TO REPRODUCE, LOAN, DISTRIBUTE OR SELL COPIES OF HIS/HER THESIS BY ANY MEANS AND IN ANY FORM OR FORMAT, MAKING THIS THESIS AVAILABLE TO INTERESTED PERSONS. L'AUTEUR A ACCORDE UNE LICENCE IRREVOCABLE ET NON EXCLUSIVE PERMETTANT A LA BIBLIOTHEQUE NATIONALE DU CANADA DE REPRODUIRE, PRETER, DISTRIBUER OU VENDRE DES COPIES DE SA THESE DE QUELQUE MANIERE ET SOUS QUELQUE FORME QUE CE SOIT POUR METTRE DES EXEMPLAIRES DE CETTE THESE A LA DISPOSITION DES PERSONNE INTERESSEES. THE AUTHOR RETAINS OWNERSHIP OF THE COPYRIGHT IN HIS/HER THESIS. NEITHER THE THESIS NOR SUBSTANTIAL EXTRACTS FROM IT MAY BE PRINTED OR OTHERWISE REPRODUCED WITHOUT HIS/HER PERMISSION. L'AUTEUR CONSERVE LA PROPRIETE DU DROIT D'AUTEUR QUI PROTEGE SA THESE. NI LA THESE NI DES EXTRAITS SUBSTANTIELS DE CELLE-CINE DOIVENT ETREIMPRIMES OU AUTREMENT REPRODUITS SANS SON AUTORISATION. ISBN 0-612-01273-5 Canada &ATTA Nome. VIWOO DfsMribrion A&tlmds International is arranged by broad, general subject categories. Please select the one subject which most nearly describes the content of your dissertation. Enter the corresponding four-digit code in the spaces provided. CIVIL £N5tlOEERH>i6 SUUECTTERM r° l5hM U'M'I SUBJECT COCK Sufejext Categories THE HUMANITIES AND SOCIAL SCIINCIS CMjMMNMOmONS AND THl ARTS ArcWt_taf» ArtHiJlory Cinema,,, Dance Fin. Art. Information Science Joumoliim Ubrory Science Mail Communications Muiic Speech Communication Theater 0729 0377 0900 0378 0357 0723 0391 0399 .0708 0413 0459 0465 Psychology Reading;.. Religious ScieTices Secondary Social Sciences .„ Sociology of Special:.. Teocnei' Training Technology Tests ana Measurements , Vocational 0525 0535 0527 0714 0533 0534 0340 0529 0530 0710 ,...0288 0747 PHILOSOPHY, RELIGION AND THEOLOGY Philosophy Religion General Biblical Studies Clergy History of Philosophy of Theology American Studies Anthropology Archaeology. Cultural Physical Business Administration General Accounting Banking... Management Marketing Canadian Studies Economics General Agricultural Commerce-Business Finance History Labor Theory Folklore Geography Gerontology History General 0422 0318 0321 0319 0320 0322 0469 0323 ,.0324 0326 0327 0310 0272 0770 0454 0338 0385 0501 0503 0505 0508 0509 0510 0511 0358 0366 .....0351 0578 SOCIAL SCIENCES EDUCATION General 0515 Administration ,.,....0514 Adult and Continuing., ,.. 0516 Agricultural 0517 Aft , , 0273 Bilingual and Multicultural 0282 Business , 0688 Community Coilegr ..,.„.,. 0275 Curriculum and Instruction ,0727 Early Childhood 0518 Elementary 0524 Finance , ,..0277 Guidance and Counseling 0519 Health, 0680 Higher 0745 History °r 0520 Home Economics ...0278 Industrial 0521 language and literature....,, 0279 Mathematics 0280 Music... 0522 Philosophy of 0998 Physical... 0523 LANGUAGE, LITERATURE AND LINGUISTICS Language General 0679 Ancient 0289 linguistics 0290 Modern 0291 Literature General 0401 Classical 0294 Comparative 0295 Medieval 0297 Modern 0298 African 0316 American ..0591 Asian 0305 Canadian (English) 0352 Canadian (French) ...0355 English 0593 Germanic 0311 Latin American 0312 Middle Eastern 0315 Romance 0313 Slavic and East European 0314 Ancient 0579 Medieval 0581 Modern 0582 Black 0328 African 0331 Asia, Australia and Oceania 0332 Canadian 0334 European.. .......0335 Latin American 0336 Middle Eastern 0333 United States 0337 History of Science 0585 Law 0395 Political Science General 0615 International Law and Relations. 0616 Public Administration 0617 Recreation. 0814 Social Work 0452 Sociology General 0626 Criminology and Penology ...0627 Demography 0938 Ethnic and Racial Studies 0631 Individual and Family Studies 0628 Industrial and Labor Relations 0629 Public and Social Welfare.... 0630 Social Structure and Development 0700 Theory and Methods 0344 Transportation 0709 Urban and Regional Planning ....0999 Women's Studies 0453 TNI SCIENCES AND ENGINEERING Agriculture General Agronomy Animal Culture and Nutrition Animal Pathology Food Scionce ana Technology Forestry and Wildlife Plant Culture Plant Pathology „ Plant Physiology Range Management , WoSd Technology Biology Anatomy Biostatistics Botany Cell... Ecology Entomology... , Genetics... Limnology Microbiology Molecular:'. Neurojdonce .., Oceanography Physiology Radiation. Veterinary Science Zoology.. Biophysics General Medical Biogeochemistry Geochemistry IKHOGKAl SCIENCES 0473 .0285 0475 0476 0359 0478 0479 0480 0817 0777 0746 0287 0308 0309 0379 0329 , 0353 0369 0793 0410 0307 0317 ,0416 ,. 0433 0821 .,0778 0472 0786 0760 , .0425 0996 Geodesy Geology. Geophysics Hydrology Mineralogy Paleobotany Paleoecology Paleontology Paleozoology Palynology Physical Geography Physical Oceanography 0370 0372 0373 ,.0388 0411 0345 0426 0418 0985 0427 0368 • 0415 Speech Pathology Toxicology Home Economics ., PHYSICAL SCIENCES Pure Sciences Chemistry General Agricultural Analytical Biochemistry Inorganic Nuclear Organic .' Pharmaceutical Physical Polymer Radiation Mathematics Physics Acoustics Astronomy and Astrophysics Atmospheric Science Atomic Electronics and Electricity Elementary Particles and High Energy Fluid and Plasma Molecular Nuclear Optics Radiation Solid State Statistics Applied Sciences Applied Mechanics Computer Science 0460 0383 0386 Engineering General Aerospace Agricultural Automotive Biomedical HEALTH AND ENVIRONMENTAL SCIENCES Environmental Sciences .„ 0768 Health Sciences General 0566 Audiology 0300 Chemotherapy 0992 Dentistry 0567 Education 0350 Hospital Management 0769 Human Development 0758 Immunology 0982 Medicine and Surgery 0564 Mental Health .....0347 Nursing 0569 Nutrition 0570 Obstetrics and Gynecology ..0380 Occupational Health and Therapy 0354 Ophthalmology 0381 Pathology 0571 Pharmacology 0419 Pharmacy 0572 Physical Therapy .....0382 Public Health.!. 0573 Radiology 0574 Recreation 0575 General 0306 0485 0749 0486 0487 ...0488 0738 0490 0491 0494 0495 0754 0405 General 0605 0986 0606 0608 0748 0607 0798 0759 0609 0610 0752 0756 0611 0463 0346 0984 ICiviL..,.„Z.,,...; • 0543jl-^ Electronic! and Electrical 0i44 Heat and Thermodynamics ...0348 Hydraulic 0545 Industrial 0546 Marine 0547 Materials Science ....0794 Mechanical 0548 Metallurgy 0743 Mining 0551 Nuclear 0552 Packaging 0549 Petroleum 0765 Sanitary and Municipal 0554 System Science 0790 Geotechnology 0428 Operations Research 0796 Plastics Technology 0795 Textile Technology 0994 General Behavioral Clinical Developmental Experimental Industrial Personality Physiological Psychobiology Psychometrics Social Chemical ™__ 0542 0537 0538 0539 0540 0541 PSYCHOLOGY EARTH SCIENCES , 0621 0384 0622 0620 0623 ..0624 0625 0989 0349 0632 0451 ® CONCORDIA UNIVERSITY SCHOOL OF GRADUATE STUDIES This is to certify that the thesis prepared By: Entitled: VINOD BATTA Two-Dimensional Boundary Element Analysis of Seismic Cracking in Concrete Gravity Dams and submitted in partial fulfillment of the requirements for the degree of DOCTOR OF PHILOSOPHY (Civil Engineering) complies with the regulations of the University and meets the accepted standards with respect to originality and quality. Signed by the final examining committee: Chair Dr. O. Schwelb ^AW^^* !______ Dr. J.L. HumaT External Examiner Dr. R. Hail kiml External-to-Program Examiner Prof. C. Marsh AT- /T/~ l^LJ^ Dr. M.M. El-BajUy ^J Examiner i WudUi Dr. V. Gocevski' Dr. O.A. Pekau Examiner Thesis Supervisor ABSTRACT TWO-DIMENSIONAL BOUNDARY ELEMENT ANALYSIS OF SEISMIC CRACKING i:N CONCRETE GRAVITY DAMS Vinod Batta, Ph. D. Concordia University, Montreal A two-dimensional formulation based on the boundary element method is presented to investigate earthquake induced fracture in concrete gravity dams. The cracking is represented by the discrete approach wherein each crack surface is modelled in a separate dam domain employing the concept of multi-domain discretization. The principles of linear elastic fracture mechanics are used to define the stress field in front of the crack tip and the stress singularity occurring at the tip is captured utilizing traction singular quarter-point boundary elements placed on each side. Stability of cracks under dynamic loads is monitored employing a propagation criterion based on the maximum tensile strain theory. The dynamic analysis is carried out by direct integration of the equations of motion. The ability of the formulation to predict correctly the stress singularity at the crack tip is demonstrated. Also verified is the ability of the formulation to simulate and monitor i changing crack profile under dynamic loading through comparison with available experimental data. Seismic cracking of the prototype Koyna dam is investigated in detail employing both single and multiple fracture models. For the single fracture model, a parametric investigation is conducted to determine the influence of various analytical and material parameters on the fracture response of a gravity dam subjected to an earthquake. Refinement of the numerical modelling of the dam by means of internal collocation is studied and the adopted dual-reciprocity approach is shown to be adequate for seismic fracture analysis of dams. In order to obtain an enhanced database, the seismic cracking response of a monolith of the Pine Flat dam, which represents a more standard gravity-Hi type dam, is also studied. It is demonstrated that the final pattern of cracking, as well as the fracture process itself, is largely unaffected by the geometrical properties of the; cross-section. Furthermore, an attempt is made to examine the likelihood of hydrodyntimic pressure build-up in cracks on the water retaining face of a dam. To this end, the corresponding opening/closing data during various phases of the dam response is analyzed. It is found that the dynamic water pressure is more likely to influence response during the post-rupture phase, rather than affect the behavior prior to rupture. Finally, noting that the fracture investigations undertaken herein assume fixed conditions for the dam base, a time-domain boundary element procedure is proposed to allow future dam/foundation coupling in the dynamic crack propagation analysis of concrete dams. IV ACKNOWLEDGEMENTS The author acknowledges with gratitude the guidance, support and encouragement provided by his thesis supervisor, Professor O. A. Pekau, during the entire course of this research. His critical reviews, together with his excellent sense of engineering judgement, have proved invaluable to the quality of this work. The initial introduction to this problem by Professor Zhang Chuhan of Tsinghua University, Beijing, and subsequent fruitful discussions during his visits to Concordia University are also gratefully acknowledged. The financial support for this work v/as provided by Natural Sciences and Engineering Research Council of Canada under Grant No. A8258. The excellent computing facilities provided at the Civil Engineering Department at Concordia University, largely because of the efforts of Professor Pekau, have been very useful in conducting this research. Thanks are also due to SNC-Lavalin Inc. for allowing the author to take time off from work to complete this research. Among his friends, the author wishes to mention Denis, the two Richards and Tony for their friendship and support. The deepest gratitude goes to the family of the author and his Bhua for their inspiration, support and encouragement throughout his studies. Finally, the utmost patience and understanding demonstrated by author's wife Sumeera deserve greatest appreciation. Her help in the preparation of this thesis is also acknowledged thankfully. Of course, these acknowledgements will not be complete without an affectionate mention of author's six month old son Varun, who has been a continuous source of joy during the final stretch of this work. v TO MY PARENTS CHIRANJILAL JI AND KAUSHALYA DEVI TABLE OF CONTENTS LIST OF FIGURES LIST OF TABLES LIST OF NOTATIONS CHAPTER I INTRODUCTION 1.1 OBJECTIVE AND SCOPE 1.2 REVIEW OF EXISTING LITERATURE 1.2.1 Preliminary remarks 1.2.2 Past investigations 1.3 ORGANIZATION OF THE THESIS CHAPTER II IMPLEMENTATION OF BOUNDARY ELEMENT METHOD FOR CRACK PROPAGATION ANALYSIS 2.1 INTRODUCTION 2.2 BOUNDARY ELEMENT FORMULATION 2.3 BOUNDARY DISCRETIZTION AND NUMERICAL SOLUTION 2.3.1 Crack tip boundary elements 2.3.2 Discretized boundary equations 2.3.3 Simulation of crack closure 2.3.4 Solution procedure 2.4 CRACK PROPAGATION 2.4.1 Computation of stress intensity factors 2.4.2 Criterion for crack extension 2.5 REDISCRETIZATION 2.6 COMPUTER IMPLEMENTATION 2.7 APPLICATIONS 9 9 9 12 12 14 15 16 20 20 21 23 27 28 1 1 2 2 3 6 xi xv xvi vii 2.7.1 Free vibration of cantilever beam 2.7.2 Dynamic analysis of centrally cracked plate 2.7.3 Crack propagation in model cantilever structure 2.8 CONCLUDING REMARKS CHAPTER III SEISMIC ANALYSIS OF SINGLE CRACK PROPAGATION IN CONCRETE GRAVITY DAMS 3.1 INTRODUCTION 29 29 30 32 46 46 3.2 DYNAMIC BOUNDARY ELEMENT FORMULATION WITH INTERNAL COLLOCATION 47 3.3 MODELLING OF STATIC LOADS 3.4 PARAMETRIC STUDY OF SEISMIC CRACKING IN KOYNA DAM .. 3.4.1 Boundary element model 3.4.2 Material characteristics and loading 3.4.3 Effect of internal collocation points 3.4.4 Effect of dynamic fracture toughness 3.4.5 Effect of crack length increment 3.4.6 Effect of static loads 3.5 CONCLUDING REMARKS CHAPTER IV ANALYSIS OF MULTIPLE SEISMIC CRACKING IN CONCRETE GRAVITY DAMS 4.1 INTRODUCTION 4.2 IMPLEMENTATION OF BEM FOR MULTIPLE CRACK PROPAGATION 4.2.1 Modelling and monitoring of multiple cracks 4.2.2 Choice of analytical factors influencing fracture analysis 4.2.3 Applications 4.3 MULTIPLE CRACKING OF FONGMAN DAM 76 76 77 77 78 78 79 49 50 51 52 53 54 57 58 59 viu 4.4 SINGLE AND MULTIPLE CRACKING OF KOYNA DAM 4.4.1 Potential cracking models 4.4.2 Assessment of static stability of cracks 4.4.3 Comparison of fracture processes for single and multiple cracking models 4.4.4 Behavior with higher concrete fracture toughness 4.4.5 Approximate post-cracking behavior based on multiple cracking model 4.4.6 Summary of observations concerning hydrodynamic pressure build-up in crack C2 4.5 CONCLUDING REMARKS CHAPTER V ELASTODYNAMIC BOUNDARY INTEGRAL FORMULATION FOR FUTURE DAM-FOUNDATION INTERACTION ANALYSIS 5.1 INTRODUCTION 5.2 MODELLING OF THE FOUNDATION MEDIUM 5.2.1 Literature review 5.2.2 Present approach 5.3 ELASTODYNAMIC BOUNDARY INTEGRAL FORMULATION 5.3 NUMERICAL TREATMENT 5.4 APPLICATIONS 5.4.1 Response of a half-plane to external dynamic stresses 5.4.2 Surface rigid strip foundation 5.4.3 Embedded rigid strip foundation 5.5 NUMERICAL RESULTS (i) Half-space under discontinuous prescribed stress distribution (ii) Half-space under continuous prescribed stress distribution (iii) Surface rigid strip footing under harmonic vibrations (iv) Embedded rigid strip footing under harmonic vibrations 81 82 84 84 88 89 91 92 109 109 Ill Ill 113 113 117 121 121 121 123 123 123 124 125 126 IX 5.6 FUTURE IMPLEMENTATION FOR DAM-FOUNDATION INTERACTION 5.6.1 Time-domain sub-structure analysis of dam-foundation system 5.7 CONCLUDING REMARKS CHAPTER VI CONCLUSIONS AND RECOMMENDATIONS 6.1 CONCLUSIONS 6.2 RECOMMENDATIONS FOR FUTURE WORK APPENDIX A STANDARD FUNCTIONS AND PROCEDURES USED IN THE BOUNDARY INTEGRAL FORMULATION Al FREQUENCY INDEPENDENT FUNDAMENTAL SOLUTION A2 DISPLACEMENT FIELD AND TRACTIONS CORRESPONDING TO FUNCTIONS fj = C - R(^-x) A3 BODY FORCE DUE TO SELF-WEIGHT A4 UPDATING OF STIFFNESS MATRIX IN CRACK PROPAGATION ANALYSIS APPENDDC B SEISMIC MULTIPLE CRACKING ANALYSIS OF PINE FLAT DAM Bl INTRODUCTION Bl.l Past investigations of cracking in the Pine Flat dam B2 NUMERICAL RESULTS B2.1 Linear response analysis and locations of potential cracking B2.2 Analysis of fracture response B2.3 Behavior with higher concrete fracture toughness B2.4 Approximate post-cracking behavior B2.5 Summary of observations concerning hydrodynamic pressure build-up in crack C2 B3 CONCLUDING REMARKS , 127 128 131 147 147 149 158 158 158 159 160 161 161 161 163 164 165 167 167 168 169 X LIST OF FIGURES Figure 2. \\ Figure 2.2 Figure 2.3 Figure 2.4 Figure 2.5 Boundary elements: (a) isoparametric quadratic element; (b) traction singular quarter-point element 35 Typical nodal pair on crack flanks and the normal direction 36 Crack tip coordinate system for: (a) displacement field; (b) calculation of stress intensity factor 37 Rediscretization for movement of crack tip from B to B' (a) One quarter of a plate with central crack; (b) time history of normalized S.I.F 39 40 38 Figure 2.6 Figure 2.7 Figure 2.8 Model structure with initial crack (dimensions in mm) Boundary element discretization of model structure and dynamic characteristics 41 Time histories of response for model structure: (a) support excitation; (b) computed crest displacement; (c) comparison for crest acceleration 42 Time histories of stress intensity factors: (a) comparison of computed and experimental Ky, (b) computed K 43 Comparison for stage-wise development of crack profile Computed crack length as a function of time Tallest non-overflow monolith of Koyna dam 44 45 61 Figure 2.9 Figure 2.10 Figure 2.11 Figure 3.1 Figure 3.2 Figure 3.3 Figure 3.4 Figure 3.5 Time histories of ground acceleration for Koyna earthquake of Dec, 11, 1967: (a) stream direction; (b) vertical 62 Observed pattern of cracking in Koyna dam 63 Koyna dam: (a) analyzed section with the pre-assigned single crack; (b) BE discretization with internal points (scheme 3/3/3) 64 Crack profiles with and without internal points for Kjd = 9.0 MPa.m1/2 65 xi Figure 3.6 Figure 3.7 Figure 3.8 Figure 3.9 Figure 3.10 Figure 3.11 Figure 3.12 Time histories of crest response with and without internal points for KId = 9.0 MPa.m1/2: (a) displacement; (b) acceleration 66 Envelopes of principal tensile stress (MPa) with and without internal points for KI(J = 9.0 MPa.ni1'2 Crack profiles for different fracture toughness Kjd 67 68 Time histories of horizontal projection of crack length for different fracture toughness Kjd 69 Time histories of crest response for different fracture toughness K]d: (a) displacement; (b) acceleration 70 Instantaneous deformation of cracked dam for different fracture toughness Kjd Time histories of principal tensile stresses at selected locations on upstream face of dam: (a) heel; (b) elevation 66.m 72 73 71 Figure 3.13 Figure 3.14 Figure 3.15 Effect of maximum crack extension on crack profile Effect of loading condition on crack profile for Kjd = 2.0 MPa.m . (a) dam self-weight and reservoir ignored; (b) all loads considered. 74 Time histories of combined mode stress intensity factor K for different loading conditions: (a) dam self-weight and reservoir ignored; (b) all loads considered 75 Fongman dam model with initial cracks: (a) boundary element discretization; (b) final cracking profiles 95 Figure 4.1 Figure 4.2 Figure 4.3 Envelope of principal tensile stresses for seismic analysis of Koyna dam without cracks 96 Boundary element discretizations of Koyna dam for different fracture models: (a) single downstream crack CI; (b) single upstream crack C2; (c) multiple cracking 97 Figure 4.4 Envelopes of principal tensile stresses for different fracture models: (a) single downstream crack CI; (b) single upstream crackC2; (c) multiple cracking 98 Fracture process for single downstream crack CI: (a) final cracking profile; (b) time history of crack length; (c) time history of crest displacement; (d) time history of crack mouth opening 99 Figure 4.5 xii Figure 4.6 Fracture process for single upstream crack C2: t'a) final cracking profile; (b) time history of crack length; (c) time history of crest displacement; (d) time history of crack mouth opening . 101 Fracture process for multiple cracking model: (a) final cracking profile; (b) time histories of crack length; (c) time history of crest displacement; (d) time histories of crack mouth openings 103 Final cracking profiles for multiple fracture model and higher magnitudes of dynamic fracture toughness: (a) Kjd = 5.5 MPa.m ; (b) KId = 9.0 MPa.m1/2 105 Envelope of principal tensile stresses for multiple cracking model and KId = 9.0 MPa.m1/2 106 Figure 4.7 Figure 4.8 Figure 4.9 Figure 4.10 Figure 4.11 Time histories of crack mouth openings multiple cracking model: (a) KId = 5.5 MPa.m1/2; (b) KId = 9.0 MPa.m1/2 107 Post-rupture behaviour for multiple cracking model: (a) Kjd = 2.0 MPa.m\"2; (b) KId = 5.5 MPa.m\"2 108 132 133 134 135 Figure 5.1 Figure 5.2 Figure 5.3 Figure 5.4 Figure 5.5 Strip footing on the surface of half-space Strip stress distribution on half-space: (a) applied loading; (b) time histories of vertical displacements Half-plane subjected to discontinuous stress: (a) loading; (b) time-history of displacement at point A Harmonic vibration of rigid strip foundation: (a) boundary discretization; (b) effect of mesh truncation on vertical stiffness Harmonic vibration of surface rigid strip footing: (a) vertical stiffness for relaxed contact; (b) vertical compliance for welded contact; comparison of vertical stiffness for welded and relaxed contact. 136 Harmonic vibration of surface rigid strip footing: (a) horizontal stiffness for relaxed contact; (b) horizontal compliance for welded contact; (c) comparison of horizontal stiffness for welded and relaxed contact Figure 5.6 138 Figure 5.7 Harmonic vibration of surface rigid strip footing: (a) rocking stiffness for relaxed contact; (b) rocking compliance for welded contact; (c) comparison of rocking stiffness for welded and relaxed contact 140 xiii Figure 5.8 Figure 5.9 Figure 5.10 Figure Bl Figure B2 Figure B3 Figure B4 Figure B5 Figure B6 Coupled horizontal-rocking compliance for surface rigid strip footing ,... 142 Harmonic vibration of embedded rigid strip footing: (a) strip footing and BE discretization; (b) time vs frequency domain analysis.... 143 Effect of embedment on the harmonic response of rigid strip footing: (a) rocking stiffens; (b) horizontal stiffness; (c) vertical stiffness 144 Pine Flat dam: (a) tallest non-overflow monolith; (b) envelope of principal tensile stresses for seismic analysis without cracks 171 Boundary element discretization of Pine Flat dam with multiple cracking 172 Fracture process for K^ =2.0 MPa.m1/2: (a) final cracking profile; (b) time-histories of crack lengths; (c) time history of crest displacement; (d) time histories of crack mouth openings 173 Final cracking profiles for higher magnitudes of dynamic fracture toughness: (a) KId = 5.5 MPa.m1/2; (b) KId = 9.0 MPa.m1/2 175 Envelope of principal tensile stresses for Kjd = 9.0 MPa.m .... 176 Post-rupture behaviour for KId = 2.0 MPa.m1/2 177 XIV LIST OF TABLES Table 2.1 Table 4.1 Table Bl Comparison of frequencies for cantilever beam Crack mouth opening data for C2 Crack mouth opening data for C2 34 94 170 xv LIST OF NOTATIONS a0 frequency factor (COD/C2) half-width of strip footing or loading velocity of dilatation wave Rayleigh damping matrix velocity of shear wave spatial functions used to transform inertial domain integral stiffness matrix spring stiffness mass matrix number of nodes on dam boundaiy excluding dam-rock interface unit outward vector normal to T number of nodes on dam-rock interface component of traction traction component of static fundamental solution distance between x and ^ tan8c/2 time component of traction traction component of dynamic fundamental solution component of displacement component of velocity component of acceleration displacement component of dynamic fundamental solution displacement component of static fundamental solution relative displacement b cl [en] C2 f> M ks [m] n (n) nb Pi * Pki r s t l(n)i *(n)ik «i \"« \"i \"ik * ur xvi ..0 u X cimposed support acceleration source point dynamic compliance of the half-space constant coefficients in the boundary integral equation inverse of matrix [F] boundary element system matrix boundary element system matrix combined mode stress intensity factor dynamic stiffness of the half-space opuning mode stress intensity factor sliding mode stress intensity factor fracture toughness of material number of boundary elements in the discretization generalized mass matrix total number of time steps shape functions components of dynamic force under the footing nodal force vector width of discretization on each side of footing; also load vector receiver element source element Rayleigh damping coefficient time-dependent functions used to transform inertial domain integral Rayleigh damping coefficient Kronecker's delta circumferential strain circumferential strain factor ij ki C[E] [G] [H] K Kij Ki Kn KId M [M] N Ni Pi {Q} R Rc Se a of P 5ik ee ee xvii components of traction corresponding to t7.. material constant related to v constant used to obtain spring stiffness ks shear modulus of material Poisson's ratio angle measured from crack tip crack extension angle mass density stress tensor forcing frequency displacement field corresponding to tj. field point; also generalized coordinate length of each boundary element time step size boundary of the region domain of the region ; xviu CHAPTER I INTRODUCTION 1.1 OBJECTIVE AND SCOPE Accurate prediction of the performance of concrete gravity dams under strong ground motions poses significant analytical challenge. The coupled dam-foundation-reservoir system requires a proper modelling of each sub-system including the inherent non-linearities which may arise from cracking of the dam concrete, non-linear behavior of the foundation material and/or cavitation at the dam-reservoir interface. The present study investigates primarily the first of these problems, namely the cracking behavior of concrete gravity dams under external excitation. The problem is of considerable importance because a large number of existing gravity dams are 50 years or more in age and were designed using the uniform seismic coefficient method which is no longer considered adequate (Chopra 1987). Moreover, for many of these dam sites the seismic codes have been scaled-up, which further necessitates a re-evaluation of the safety of these structures. Since the past experience - cracking of Koyna, Hsigfenkiang and Sefidrud dams (Ahmadi and Khoshrang 1992) - has shown that the dam concrete may crack significantly under an earthquake with strong enough ground shaking, such safety evaluations entail proper crack propagation analysis. The question to be addressed is to what extent a crack of known size and location will grow, since this generally will determine the operational safety of a dam under the dynamic loadings resulting from a future earthquake. The principal objective of the current research is to develop an efficient and accurate numerical procedure for analyzing seismic cracking in concrete gravity dams and to perform a detailed investigation of their fracture behavior under earthquakes. 1 Employing the analytical procedure developed in this study, a parametric evaluation involving various analytical and physical parameters is conducted to examine their influence on the fracture response of dams. Since the interaction of reservoir water within a crack on the water retaining face of a dam is of considerable importance, an attempt is made to study the crack opening/closing behavior and relate it to the possible effect of water pressure on fracture process itself. In addition, independent of the above investigation, a formulation is presented also for the two-dimensional soil-structure interaction problem typically encountered at the base of a gravity dam. In this regard, it has been demonstrated (Chopra 1987) that the response of a dam subjected to earthquake induced vibrations is affected significantly by the interaction with the foundation rock. However, since most of the formulations developed for modelling such interaction (Vaish and Chopra 1974; Dasgupta and Chopra 1979; Abdalla 1984) work in the frequency domain, they cannot be incorporated in the analysis of dams undergoing crack propagation. The non-linearity inherent in such analyses necessitates that the stiffness or compliance functions of the foundation medium be derived in the time domain. A procedure is therefore developed in the time domain so that it can be integrated in the seismic crack propagation analysis of dams. The boundary element method, which offers many advantges over the finite element method in the modelling of infinite domains, is adopted for this formulation also. 1.2 REVIEW OF EXISTING LITERATURE 1.2.1 Preliminary remarks In carrying out a dynamic crack propagation analysis three distinct aspects of the problem arise: namely, the choice of an appropriate numerical tool, proper modelling of the crack itself and, finally, the selection of a suitable criterion for crack extension. In order to formulate an accurate and yet efficient methodology to solve this complex problem, a proper blending of these features is required. 2 Current numerical techniques employed in fracture analysis comprise the finite element and the boundary element methods. In either method, a proper representation of the strain discontinuity imposed by the presence of cracks in the domain is imperative. The available alternatives consist of explicitly representing the crack as a physical separation of the nodes or else modelling the formation of a crack by altering the stress-strain relationship in the cracked zone in such a way that no tensile stress develops perpendicular to the crack. These alternative techniques are, respectively, known as the discrete and the smeared crack approaches. While the finite element method is capable of easily integrating either of these approaches in its formulation, the boundary element technique poses certain difficulties when used in conjunction with the smeared crack approach and thus exhibits limitations when applied to structures which have a cluster of small cracks situated in specific zones. On the other hand, in analyzing the propagation of a discrete crack, the finite element method requires complete remeshing of the domain around the original and the propagated cracks. This redefinition of the mesh enlarges the system and destroys its bandwidth. Alternative use of boundary element discretization requires only the addition of elements to represent the extended crack profile. 1.2.2 Past investigations The available studies of seismic cracking in concrete gravity dams have employed mostly the finite element method together with either smeared or discrete modelling of the crack. The smeared crack model combined with tensile strength propagation criterion was employed in early studies of earthquake induced cracking (Pal 1976; Mlakar 1987; Graves and Derucher 1987). However, as demonstrated by Bazant and Cedolin (1979), the metliod lacks objectivity because the force required to propagate the crack after initiation depends upon the size of the finite elements in front of the crack. To overcome this spurious mesh dependence it was suggested that the cracking model must incorporate the principles of fracture mechanics so that the fracture energy can be conserved during the 3 cracking process. This procedure was subsequently employed by Vargas-Loli and Fenves (1989) to study the effect of cracking on the earthquake response of gravity dams. In order to circumvent the limit on the size of the finite elements, the size reduced strength criterion proposed by Bazant and Cedolin (1979) was used. Although objective results were obtained, the analysis could only predict diffused crack patterns, which is the expected limitation of smeared crack procedures. Improvements in the form of localization of the crack band (Bazant and Oh 1983) and the use of a coaxial rotating crack model (Rots 1991) have been attempted and the latter was recently adopted by Bhattacharjee and Leger (1993) to study cracking of concrete dams under static and dynamic loads. The formulation employs nonlinear fracture mechanics and relies on constitutive modelling for crack initiation as well as for the propagation behaviour. The procedure, however, limits the maximum size of the finite elements which can be used in the mesh and thereby introduces a computational disadvantge when analyzing large dams. The above considerations have prompted the use of discrete crack models. Discrete models provide a better physical representation of cracks and, in the case of hydraulic structures such as dams, also permit modelling of the water pressure within the crack. Employed in a finite element procedure the model requires a remapping of the mesh at each stage of crack extension. Employing a semi-adaptive mesh updating procedure and a strength-based criterion for the cracking process Skrikerud and Bachmann (1986) were the first to apply discrete crack model in analyzing cracking in dams. However, results for the Koyna dam showed that the predicted cracking was strongly mesh dependent. The procedure was subsequently refined by Feltrin et al. (1992) who modified the discrete crack model to incorporate strain softening of the concrete but found this refinement to be of minor importance. Discrete crack propagation based on linear elastic fracture mechanics was proposed by Ayari and Saouma (1990). A contact-impact model for crack closure was developed and seismic cracking analysis of Koyna dam was presented. Unfortunately, 4 results were presented for only an unrealisticly high value of the concrete fracture toughness. Chapuis et al. (1985), on the other hand, used a material model derived with particular regard to the effect of strain rate and seismic loading history for analyzing cracking of the Pine Flat dam. However, the hybrid smeared-discrete crack model used in following the crack trajectory requires two sets of analyses. More recently, Zhou and Lin (1992) employed an idealized sliding structural model to study the seismic cracking of the Fongman dam in China. The analytical results were also corroborated by model testing on a shake table. In recent years discrete crack propagation analysis procedures have also been developed using the boundary element method, although for the dynamic case the literature is scarce. The general advantages of adopting a boundary approach, particularly the reduction of the dimensions of the problem, are now well established. Furthermore, the two-dimensional problem posed by the present fracture analysis of a gravity dam monolith requires only the use of line elements. Modelling and propagation of discrete cracks is thus considerably facilitated as the rediscretization process involves only realigning the line elements along the new crack profile. Moreover, since the method requires discretization only of the boundary of the structure and utilizes the exact solution of the governing differential equation in the interior, it offers a more accurate solution for fracture problems than the finite element technique. However, unlike the latter, the boundary element method yiulds a set of equations which is neither symmetric nor sparse, although it is dimensionally much smaller than the corresponding finite element system. Application of the boundary element method to the analysis of cracked structures poses an analytical difficulty resulting from the identical coordinates of nodes along the two crack surfaces, which lead to a singular system of equations. Various strategies have been proposed in the literature (Cruse 1978; Snyder and Cruse 1975; Blandford et al. 1981) to overcome this limitation. Of these, the multidomain discretization proposed by Blandford et al. (1981) and employed by Ingreffea et al. (1989) provides the most accurate 5 mathematical representation of the crack but was limited therein to fracture analysis under static loads. For dynamic fracture mechanics, a brief summary of previous applications of the boundary element method is presented by Manolis and Beskos (1988), which reveals that most of these attempts (Fan and Hahn 1985; Mettu and Nicholson 1988; Dominguez and Chirino 1986) are restricted to the analysis of stationary cracks under dynamic loads and thus do not include crack propagation. To the author's knowledge, the only available boundary-element-based dynamic crack propagation analysis formulation is by Feng (1994). The dynamic analysis procedure, however, employed the mode superposition technique for obtaining the time-step response and was therefore not very efficient for the present non-linear analysis, as the eigenvalue problem was required to be solved after each modification to the crack geometry. The formulation also relied on the assumption of a hypothetical cracked structure wherein overlap of the crack surfaces was permitted during crack closure in order to maintain unchanged the elastic stiffness of the structure. Change in stiffness accompanying crack closure was simulated by load pulses applied at the contact nodes. An alternative crack closure scheme was also suggested employing an iterative force method, which further reduced the efficiency. In addition, the formulation was limited to the treatment of a single crack propagating through the dam. 1.3 ORGANIZATION OF THE THESIS In the present study, the multidomain boundary element method has been adopted for structural idealization. The discrete crack approach is employed for fracture modelling and the crack propagation procedure is generalized to account for the presence of multiple cracks. Principles of linear elastic fracture mechanics are used to conduct a time-history analysis of crack propagation. Although normally nonlinear fracture mechanics, accounting for the fracture process zone in front of the crack tip, is recommended for concrete structures, for large structures such as dams linear elastic fracture mechanics is 6 considered appropriate (Bruhwiler et al. 1991;. It is argued that, in the case of concrete dams, the size of fracture process zone is small compared to the size of the structure. Brittle fracture thus dominates the cracking since the fracture energy contribution to the total energy is not significant (Feltrin et al. 1992; Bhattacharjee and Leger 1993). Implementation of the boundary element procedure for dynamic crack propagation analysis is presented in Chapter II. Techniques for modelling the crack and monitoring its stability are also described in this Chapter, together with the procedure used for remeshing after each crack extension. The accuracy of the dynamic boundary element formulation in capturing the crack tip stress singularity in a dynamic environment is verified using a classical solution. The capability of the developed procedure to conduct a time step-by-step crack propagation analysis and to accurately predict the crack profile is demonstrated by comparing the * analytical results with those available for a laboratory specimen of a gypsum cantilever beam tested on a shake table. Application of the developed procedure to study earthquake induced cracking in a concrete gravity dam is presented in Chapter III. Also discussed is a refinement to the transient boundary element formulation in the form of internal collocation within a domain in order to permit a better representation of the inertia of the system. Necessity of internal collocation for fracture analysis of dams is evaluated and the influence of various analytical and physical parameters on the cracking response of dams is investigated. Results are presented for the Koyna dam with a single crack propagating from the downstream face. In Chapter IV the above procedure is extended to a more general case where the occurrence and simultaneous propagation of multiple cracks can be simulated. A detailed examination of the cracking process for the Koyna dam employing both single and multiple fracture models is conducted. Particular attention is focussed on the patterns of rupture predicted by these fracture models, as well as on the associated crack opening behaviors. 7 Chapter V presents a time-domain boundary element formulation for two-dimensional elastodynamic foundation problems. The developed procedure is verified by computing the response of the half-space, representing the soil medium, to external dynamic stresses and comparing the results with available analytical solutions. The applicability and usefulness of the formulation for future dynamic interaction problems is established through successful analyses of surface as well as embedded rigid strip foundations. Future implementation of the procedure to incorporate foundation interaction effects in the crack propagation analysis of concrete gravity dams is also discussed in this Chapter. Finally, the significant conclusions of the present research are summarized in Chapter VI. Recommendations for future work in the area are also given in this Chapter. CHAPTER H IMPLEMENTATION OF BOUNDARY ELEMENT METHOD FOR CRACK PROPAGATION ANALYSIS 2.1 INTRODUCTION The present Chapter formulates the problem of dynamic crack propagation using boundary elements. Multidomain boundary discretization is adopted to model the structure and a procedure is developed to study the extension of an existing crack under external excitation. Analysis is conducted by means of a time marching scheme and the extent, direction and time of propagation of the crack is determined employing the principles of linear elastic fracture mechanics theory. In addition to the direct integration of the equations of motion, the formulation proposed herein avoids also the hypothetical cracked structure of Feng (1994). To simulate closure of the crack, dimensionless springs are employed to avoid overlap, thus providing closer analytical modelling of the physical situation. Comparison with the available experimental results for a gypsum cantilever beam tested on a shaking table under quasi-harmonic excitation confirms the good performance of the dynamic crack propagation analysis proposed herein. 2.2 BOUNDARY ELEMENT FORMULATION The boundary element formulation, also known as boundary integral equation analysis, redefines a problem in terms of quantities assigned to the boundary of the structure only and thus reduces the integral operations by one step. The governing equations are replaced by a boundary integral operator using Betti's reciprocal work theorem. In the multidomain boundary element method, the structure is divided into subdomains by the crack profile in such a way that the two crack surfaces are represented 9 as traction free boundaries belonging to two separate domains. The integral equations are then formulated for each domain and assembled. Beyond the crack tip these domains are analytically separated by a line extending from the crack tip. Continuity and equilibrium conditions are satisfied along this arbitrary line. For static problems the boundary integral expressions are widely available; however, their use in a dynamic analysis is limited because of the additional requirement of formulating the mass matrix as a function of the boundary nodes only. Earlier attempts to use the boundary solutions in elastodynamics, typically that of Cruse and Rizzo (1968), work in transformed domain and require back-transformation to obtain the time-dependent solution. Alternatively, transient boundary solutions have been proposed by Nardini and Brebbia (1982) and Ahmad and Banerjee (1985) formulated using the frequency-independent functions and particular solution approach. The procedure of Nardini and Brebbia (1982) is ^Hrpted for the present formulation and is briefly summarized below. The boundary integral expression for a 2-D isotropic, homogeneous and linear elastic body £2 with boundary T can be written as -Cki ft) U{(£) + juu(t r x)p,(x) dT- jp^il,x) r -pjui(x,t)u*ki(^,x)dQ a. ut (x)dT = 0 (2.1) where p is the mass density; % and x are the source and field points, respectively; Cki denote constants uniquely defined by the position of point E, with respect to the boundary; U; and p; represent displacements and tractions on the boundary; ui denote the components of acceleration; u*ki and p*ki comprise the frequency independent fundamental solution defined in Appendix A. In equation (2.1) standard tensor notation, with a repeated subscript implying summation, is adopted. A comma denotes spatial differentiation, while 10 temporal differentiation is denoted by an overhead dot. In order to transform the inertial integral, an approximation for acceleration «; within the domain is employed. Expressing the time-dependent displacements uj(x,t) as the sum of m spatial functions f*(x) multiplied by the unknown time-dependent functions al (t) yields ii,. (x, t) = a\\ (t)fj (x); (j = 1 to m) The acceleration ui (x, t) can then be written as */,.(*, 0 = a>t(t)fj(x) (2.3) (2.2) Defining a displacement field \\|fj. such that the corresponding stress tensor is given by iun,m = V' (2-4) where 8]j is the Kronecker delta, and utilizing equation (2.2), permits expressing the inertial domain integral of equation (2.1) as p Jii,. (x, t) u*ki {\\, x) dQ. = po{(0 {- Cki (x) \\p[. (x) + juki (|, x) r^. (x) dT a r -lPlia,x)%{x)dT} (25) where r^. are the tractions corresponding to the stress field x\\im. Combining equations (2.5) and (2.1) yields the following expression: Cki (\\) u, ($) + \\p*u a, x) M. (x) dT - \p":{"h":20.981,"w":94.941,"x":405.163,"y":949.489,"z":77},"ps":{"_opacity":1},"s":{"letter-spacing":"0.291ki {%, x)Pi (x) dT r r -pdi(t) {Cki($)^u(x)+\\pli&,x)viu(x)dr-\p":{"h":20.982,"w":65.405,"x":452.859,"y":1004.929,"z":93},"ps":{"_opacity":1},"t":"worditt,*)%(*)dT} r r = 0 (2.6) 11 This is the governing integral equation which, in the multidomain formulation, is applied to the boundary of each subdomain. It should be noted that the representation of the internal displacement field by equation (2.2) introduces an approximation because the functions f(x) result in displacements \\|f|. and the corresponding stresses xl. which satisfy only static equilibrium given by equation (2.4). On the other hand, the procedure of Ahmad and Banerjee (1985) eliminates this approximation by employing a domain displacement field that satisfies more closely the differential equation of a dynamic problem. This is achieved by considering also the inertial effect in defining the displacement field. The need for this, or similar refinement such as the introduction of interior point collocation, is discussed in the next Chapter. 2.3 BOUNDARY DISCRETIZTION AND NUMERICAL SOLUTION Numerical solution of integral equation (2.6) is obtained by discretizing boundary r into a finite number of elements. Displacements and tractions within these boundary elements may be expressed as interpolation functions using standard isoparametric formulation. Since use of these elements in the vicinity of the crack tip does not account for the true variation of displacements and stresses, special crack tip boundary elements are used. 2.3.1 Crack tip boundary elements Linear elastic fracture mechanics theory suggests that the displacements and tractions, respectively, exhibit Jr and l/Jr variation in the vicinity of the crack tip, where r is the distance from the tip. The required displacement variation is achieved by moving the mid-point node of the quadratic element to the quarter point, as is commonly done in the corresponding finite element formulation. This process is shown in Figure 2.1. Thus, the displacement variation within the crack tip element is expressed by 12 3 u = a0 + axJ% + a2\\ = ]>>,.($)«,. t = i (2.7) where a<), a, and a2 are coefficients; £, is the normalized coordinate; Nj(^) represent the element shape functions; and Uj are the nodal displacements. For the quarter-point element of Figure 2.1(b), shape functions N^(^) are given by N,(%) = 1-371 + 2^ N2(k) =4(J|-0 JV3U) = -V| + 2^ (2-8) For the element tractions, however, the quarter-point boundary elements require additional modification of the interpolation functions. This is because the boundary element method requires the displacements and tractions to be independently represented, unlike the finite element method wher; iU° traction singularity is inherently obtained by differentiating the displacement shppe irnctions. The correct traction shape functions, representing the analytical singularity, are obtained in this case by multiplying the displ.- ..Int shape function by l/«/f. Thus, the traction variation over the crack tip element becomes where t\\ represent the nodal values of traction. The traction singular quarter-point boundary elements defined above are introduced only at the crack tip, whereas the elements adjoining these are of the standard isoparametric type. 13 2.3.2 Discretized boundary equations A discretized matrix form of equation (2.6) can now be written as [M]{ii} + [H]{u}=[G]{p} (2.10) where [H] and [G] are the system matrices of the corresponding static problem; {u}, {u} and {p} are the displacement, acceleration and traction vectors, respectively; and [M] represents the generalized mass matrix which has the form [M]= p([G][r|]-[H][\\|/])E (2.11) with [E] = [F]\"1, where [F] is the matrix of values of functions f\"(x) evaluated at the boundary nodes. Various choices for the class of functions f1 have been proposed in Nardini and Brebbia (1982). Since these functions are used to formulate the inertial matrix only and noting that the mass integral does not involve derivatives, simple interpolation functions provide reasonable accuracy. Herein, functions fJ = C - R(E,j, x) are adopted, where R(^;, x) is the distance between point ^j and x, and C is a suitably chosen constant. The common value for C is unity and this is also employed herein. The corresponding expressions for the displacements \\|^\\ and tractions T^. are given in Appendix A. Equation (2.10) is the governing discretized boundary integral equation for the dynamic system. Since use of nodal forces is preferred over nodal tractions {p} for most engineering applications, a general transformation of this equation is carried out. For distribution matrix [N], with coefficients derived from element interpolation functions, the nodal forces {Q} can be expressed as {Q} = [N]{p} (2.12) Using equation (2.12) and premultiplying equation (2.10) by [N][G]', one obtains 14 [«]{«} + [« {u} = {Q} where [k] = [MEG]\"1!*] [m] = [AT] [GJ-^Mj (2.13) (2.i4) In equation (2.14), [k] and [m] represent the stiffness and mass matrices of the system, respectively. It should be noted that equation (2.13) is valid over the entire boundary of the system. In most practical applications, the forcing functions are in the form of external excitations or support movements acting over only certain portions of the boundary. If subscript '1' is used to identify those portions of the boundary on which no kinematic excitations act, equation (2.13) can be partitioned and rewritten as [mn] {ur} + [*n] {ur} = [m0] {li0} + {Q,} (2.15) in which ii° denotes the vector of imposed support accelerations, [m0] is a suitably transformed inertial matrix; and {ur} represents the unknown relative displacement vector with respect to the support movement. 2.3.3 Simulation of crack closure Discrete cracks require special modelling in a dynamic environment. This is because a cracked structure subjected to external excitation experiences alternate opening and closing of the crack. During the opening phase, the two surfaces separate and become traction free over the open portion of the crack. On the other hand, when the crack is in the process of closing the surfaces come into contact and bear against each other over the closed length of the crack. In modelling the crack by an analytical technique such as the present multidomain boundary element method, the two crack surfaces are represented as 15 traction tree boundaries. Thus, the correct representation of the opening phase is inherent in the analytical model; however, in representing crack closure the model allows overlap of the surfaces. This limitation arises because, once a certain number of nodes on the two crack flanks have come into contact, they are not prevented from moving past one another when the structure undergoes subsequent deformation. The above analytical limitation is overcome by employing a series of non-dimensional and massless springs of high stiffness between adjacent nodes of the two crack flanks. Stiffness properties of these springs are determined such that overlap is prevented. A brief description of the solution procedure used for equation (2.15) and the required modifications to prevent overlap is given below. 2.3.4 Solution procedure Equation (2.15) represents dynamic equilibrium for the boundary element system without damping. With the introduction of viscous damping equation (2.15) becomes [«„] («r) + [c„] K) + [*n] K> = [m0] {u0} + {fij where [c^] is the Rayleigh damping matrix which has the well-known form (2.16) [cn] = a[mn] + pTkn] in which a and P represent the Rayleigh damping coefficients. Although the system matrices [m^] and [k^] are not symmetric, it is nonetheless possible to decouple equation (2.16) into component modes (Feng 1994). Coefficients a and (3 can thus be obtained from standard procedure by relating them to the assumed damping in two selected modes of vibration. Equation (2.16) is solved to obtain the unknown displacement {ur} for a given 16 external excitation {u } and external forces {Qi). The choice of the appropriate integration technique is, however, critical. Direct time step-by-step integration using the Houbolt method is employed in the present study. Not only does the method provide numerical damping for the higher modes it also, as discussed later, avoids the additional assumptions required to carry out the dynamic crack propagation analysis. For convenience, various subscripts in equation (2.16) are dropped and it is re-written as [m] {u} + [c] {u} + [k] {u} = {R} (2.17) for the linear case where no crack propagation is considered. The load vector {/?} = ([m0] {ii0} + {Q}). In the nonlinear analysis involving a propagating crack, the term [k]{u} in equation (2.17) is replaced by the nonlinear restoring force vector {pr}. Thus, equation (2.17) becomes [m] {ii} + [c] {ii} + {pr} = {/?} (2.17a) in which the restoring force vector {pr} is updated after each stage of crack propagation based on the cunent tangent stiffn i matrix of the cracked dam obtained as discussed in Appendix A. At the (n+1)* time step, displacement {un+j} can be written as where {5un} is the incremental displacement at the (n+1)\"1 time step. Using the backward difference expressions for acceleration and velocity and utilizing equation (2.18), equation (2.17a) can be written as (^[»]^[c] + W){5M,} = (A2r++m]+_^[C])K} (2.19) -(A7> [m] +2Al^) K-il <^5 [m] +3Al[C]) {M«-2> + {/?« + i> 17 in which At is the time step size. Equation (2.19) is solved to obtain the incremental displacement {5un}. Thus, 5«„ = [k]~l{Rn + l} (2.20) where the effective stiffness matrix [k] is expressed as [*]=^[»]+SSiW + W and the pseudo load vector {Rn +1} is given by 2 {R„ + l} = {*B+i} + [m](3{Mn}-4{MB_1} + {M„_2})/Ar: + 7 3 1 [c](l{un}-hun_l}+\\{un_2})/At. Since the solution for 8un requires knowledge of un, un_i and un_2, the solution for the first two time steps is obtained by employing the central difference explicit scheme with a fraction of At as the time step. Once the boundary displacement vector {un+1} is obtained from equation (2.18), relative displacements are calculated for all the pairs of nodes situated along the crack. Thus, if 'i' and 'j' denote a pair of nodes belonging to the upper and lower crack flanks (Figure 2.2), overlap is said to occur if uri