Univerza v Ljubljani Fakulteta za gradbeništvo in geodezijo mi 11 i7t| i l il iim I milil 111M M UMimi JAKA DUJC METODA KONČNIH ELEMENTOV ZA RAČUN MEJNE NOSILNOSTI IN LOKALIZIRANE PORUŠITVE PLOSKOVNIH KONSTRUKCIJ DOKTORSKA DISERTACIJA Ljubljana, 2010 To naslovno stran smo dodali pri Digitalni knjižnici Slovenije, začasno, dokler ne prejmemo izvirne različice. BIBLIOGRAFSKO-DOKUMENTACIJSKA STRAN Obseg in oprema: Ključne besede: UDK Avtor: Mentor: Somentor: Naslov: 519.63:624.073(043.3) Jaka Dujc izr. prof. dr. Boštjan Brank prof. dr. Adnan Ibrahimbegovic Metoda konCnih elementov za raCun mejne nosilnosti in lokalizirane porusitve ploskovnih konstrukcij 211 str., 6 pregl., 100 sl., 429 en. konstrukcije, metoda konCnih elementov, armirani beton, jeklo, rezultante napetosti, elastoplasticnost, lokalna porusitev, vkljucena nezveznost Izvlecek V delu obravnavamo račun mejne nosilnosti in mejne duktilnosti konstrukcij z metodo končnih elementov. Ko je konstrukcija na meji svoje nosilnosti, so nekateri njeni deli v neelastičnem stanju, v najbolj kritičnih konstrukcijskih elementih pa se pojavi lokalizirana porusitev materiala, ki je povezana s koncentracijo neelastičnih deformacij na majhnem območju. Pri krhkih materialih predstavlja lokalizirano porusitev nastanek in razvoj večje razpoke, pri duktilnih materialih pa nastanek lokaliziranega nestabilnega plastičnega območja. Zato je račun mejne nosilnosti povezan tako z modeliranjem standardnih neelastičnih materialnih efektov, kot tudi z modeliranjem lokalizacijskih efektov v materialu, ki se jih pogosto označuje kot mehčanje materiala. V tem delu za modeliranje standardnih neelastičnih materialnih efektov izpeljemo in uporabljamo elastoplastične materialne modele, dotaknemo pa se tudi elastoviskoplastičnega materialnega modela za plosče in nelinearnega materialnega modela za armiranobetonske plosče. Skupno vsem izpeljanim in uporabljenim neelastičnim modelom je, da so definirani na nivoju rezultant napetosti. Za modeliranje pojava lokalizacij v materialu obstajajo različni matematični modeli in numerični algoritmi, ki pa so pogosto premalo natančni ali premalo učinkoviti. Zato v tem delu uporabimo novejso metodo, ki temelji na končnih elementih z nezveznimi pomiki znotraj končnega elementa. Dokazemo, da je mogoče v okviru te metode izpeljati tudi kom-pleksnejse končne elemente - kompleksnejse tako v smislu zahtevnejse kinematike samega osnovnega modela, kot v smislu kompleksnejse aproksimacije nezveznih pomikov. Razvijemo več končnih elementov za račun različnih konstrukcijskih elementov. Najprej je predstavljen nelinearni končni element za račun mejne nosilnosti armiranobetonskih plosč, nato pa končni element za neelastično analizo kovinskih plosč, ki temelji na elastoplastičnem in elastoviskoplastičnem materialnem modelu z izotropnim in kinematičnim utrjevanjem. Pri slednjem izpeljemo računski algoritem, ki hkrati zajame obe neeleastični formulaciji. Nadalje izpeljemo končni element za kovinske lupine z elastoplastičnim materialnim modelom z izotropnim in kinematičnim utrjevanjem, ki temelji na rezultantah napetosti in na kriteriju plastičnega tečenja, ki ga določata dve ploskvi. Zadnja dva izpeljana končna elementa sta namenjena modeliranju lokalizirane porusitve materiala v nosilcih in ravninskih telesih. Elastoplastični končni element za Euler-Bernoullijeve nosilce ima vključeno nezveznost v rotacijah, njegovi konstitutivni parametri pa so določeni s posebnim postopkom, ki temelji na predhodni nelinearni analizi izbranega dela nosilca s podrob-nejsim modelom. Zadnji izpeljani končni element je elastoplastični stirivozlisčni končni element za ravninska telesa, ki ima vključene nezveznosti v pomikih. Kinematika nezveznosti je taksna, da omogoča linearne nezvezne skoke v pomikih tako v smeri normale, kot tudi v smeri tangente na črto nezveznosti. Numerični primeri dokazujejo, da je mogoče z izpeljanimi končnimi elementi in s pripadajočimi algoritmi učinkovito in na dokaj robusten način izračunati tako mejno nosilnost konstrukcije, kot simulirati popolno porusitev konstrukcije. Med drugim je predstavljen primer računa nastanka in sirjenje razpoke v krhkem materialu in primer strizne porusitve v duktilnem materialu. Vse računalniske kode za izpeljane končne elemente so bile pripravljene v programskem okolju AceGen, ki omogoča razvoj programske kode v simbolnem zapisu in optimizacijo generirane kode. Numerični primeri so bili izračunani v programskem okolju AceFem. BIBLIOGRAPHIC-DOCUMENTARY INFORMATION Title: UDC Author: Supervisor: Co-supervisor: 519.63:624.073(043.3) Jaka Dujc assoc. prof. Boštjan Brank prof. Adnan Ibrahimbegovic Finite element analysis of limit load and localized failure of structures Notes: Key words: 211 p., 6 tab., 100 fig., 429 eq. structures, finite element method, reinforced concrete, steel, stress resultants, elastoplasticity, localized failure, embedded discontinuity Abstract The dissertation deals with limit load and limit ductility analysis of structures by the finite element method. When structure is at its limit load, several structural components behave inelastically, while in the critical parts of the structure, due to localization of inelastic strains, failure of material appears. Localized effects in brittle materials are related to appearance and formation of a large (macro) crack, while failure in ductile materials is governed by localized shear bands. The study of limit load is thus related to modeling both standard inelastic material effects, as well as modeling of localized failure of material, often reffered to as material softening. Standard inelastic material effects are in this work described with elastoplastic, elas-toviscoplastic and nonlinear elastic material models. All the material models are defined at the level of stress-resultants. Several mathematical approaches and numerical algorithms for modeling localized effects are at hand, but they are often inefficient or inaccurate. Therefor, we use an up-to-date approach, based on a finite element method with embedded discontinuity. We derive new finite element formulations with a quite complex kinematics of the basic elements, as well as rather complex description of discontinuous displacement fields. We derived several finite element formulations for analysis of different structural components. First we present a finite element for limit load analysis of reinforced concrete plates. Stress-resultant elastoplas-tic and elastoviscoplastic plate finite element formulation along with a unified computational procedure that covers both formulations are presented next. Further, a nonlinear shell finite element, based on a two-surface yield function, that includes both isotropic and kinematic material hardening is presented. The last two finite elements derived in this work are intended to model the localized failure in planar beams and 2D solids. The embedded discontinuity in rotations was built into elastoplastic Euler-Bernoulli beam finite element, and a procedure, based on a precomputed analysis of a part of a structure, by using a refined numerical model, is proposed to obtain the beam constitutive model parameters. Finally, we derive an elastoplastic quadrilateral two-dimensional finite element formulation with embedded strong discontinuity, whose kinematics can model linear jumps in both normal and tangential displacements along the discontinuity line. Numerical simulations show, that the derived finite elements, along with the accompanied numerical algorithms, are an efficient and a rather robust tool for limit load and failure analysis of structures. Among other examples, we present a simulation of crack growth in brittle material and a simulation of shear band failure in ductile material. All the computer codes of the finite element formulations presented in this work have been generated through the symbolic programming of the finite element computer code and the expression optimization in AceGen computer program. The performance of these elements has been presented in numerous numerical examples, all performed by the AceFem computer program. INFORMATION BIBLIOGRAPHIQUE-DOCUMENTAIRE Notes: Mots cles: CDU Auteur: Directeur de these: Co-directeur de these: Titre: 519.63:624.073(043.3) Jaka Dujc prof. Boštjan Brank prof. Adnan Ibrahimbegovic Analyse Elements Finis de la charge limite et de la rupture localisee des structures 211 p., 6 tab., 100 fig., 429 eq. structures, methode elements finis, beton arme, acier, efforts generalises, elastoplasticite, rupture localisee, enrichissement discontinu Resume Ce travail a pour objet l'analyse limite des structures par la methode des elements finis. Lorsqu'une structure atteint sa charge limite, certaines de ses composantes sont dans la phase inelastique de leur comportement, alors que dans les parties les plus critiques, du fait de la localisation des deformations inelastiques, se produit la rupture du materiau. Les effets de localisation sont, dans les materiaux fragiles lies a l'apparition et au developpement de macro fissures alors qu'ils sont, dans les materiaux ductiles, gouvernes par les bandes de cisaille-ment localisees. L'etude de la charge limite est ainsi reliee a la modelisation du comportement inelastique standard du materiau mais egalement a la modelisation des effets localises corre-spondant au comportement adoucissant des materiaux. Le comportement inelastique standard du materiau est, dans ce travail, decrit par des modeles elastoplastiques, elastoviscoplastiques ou elastiques non lineaires. Tous les modeles de comportement sont definis en termes d'efforts generalises. Un certain nombre d'approches mathematiques et d'algorithmes numeriques sont disponibles mais sont bien souvent inefficaces et manquent de precision. Ainsi, nous utilisons une approche developpee plus recemment s'appuyant sur une methode d'elements finis enrichis de discontinues. Nous avons developpe de nouvelles formulations d'elements standards prenant en compte des cinematiques et des descriptions des champs de deplacements discontinus complexes. Plusieurs formulations d'elements finis ont ete developpees pour l'analyse de differents com-posants structurels. Nous presentons, dans un premier temps, un element fini dedie a l'analyse limite des plaques en beton arme. La formulation d'un element de plaque elastoplastique et elastoviscoplastique ecrite en efforts generalises associee a une procedure commune d'integration sont presentees ensuite. Un element de coque non lineaire, faisant intervenir une fonction seuil a deux surfaces incluant a la fois un ecrouissage isotrope et un ecrouissage cinematique est ensuite presente. Les deux derniers elements finis developpes dans ce travail sont dedies a la modelisation de la rupture localisee dans les poutres planes et les solides bidimensionnels. L'element de poutre d'Euler-Bernouilli est enrichi par une discontinuite en rotation. Une strategie s'appuyant sur l'analyse prealable, par un modele raffine, d'une partie de la structure est proposee afin d'obtenir les parametres du modele constitutif de la poutre. Enfin, nous presentons la formulation d'un element quadrangulaire a discontinuite forte dont la cinematique permet de prendre en compte des sauts de deplacements lineaires dans les deux directions normale et tangentielle le long de la surface de discontinuite. Des resultats numeriques montrent que les elements developpes ainsi que les algorithmes associes constituent un outil efficace et robuste d'analyse de la charge limite et de la rupture des structures. Parmi les exemples, nous presentons la simulation de la propagation d'une fissure dans un materiau fragile ainsi que le developpement d'une bande de ci-saillement dans un materiau ductile. Les codes numeriques associes aux formulations presentees dans ce travail ont ete genres par l'outil de programmation symbolique et d'optimisation de code AceGen. Les performances des elements sont presentes a travers un grand nombre d'exemples numeriques realises a partir du code AceFem. ZAHVALA Za pomoč pri pripravi doktorske disertacije se zahvaljujem izr. prof. Boštjanu Branku in prof. Adnanu Ibrahimbegovicu. Vajina vrata so mi bila vedno odprta, tudi ko ni slo za reševanje strokovnih tečav. Posebna zahvala gre moji mami Silvestri in sestri Jerneji, ki sta mi vedno nudili podporo. Z besedami ne morem opisati, kako sem vama hvalečen. Za podporo med študijem in številne nepozabne trenutke se iskreno zahvaljujem Ireni, Jozetu, Primožu, Urošu in Vidi. Tekom studija sem se druzil z veliko ljudmi, zaradi katerih sta bila delo na fakulteti in vsakdan prijetna in zanimiva. Seznam ljudi, katerim bi se rad zahvalil, je predolg, da bi ga lahko napisal na to stran. Moram pa izpostaviti Miho, Mateja in Vasjo, ki so mi nudili podporo, ko sem jo najbolj potreboval. Zahvalil bi se tudi DDC svetovanje inzeniring d.o.o, katerega stipendist sem bil v Času dodiplomskega studija. Na koncu bi se rad zahvalil ,se kolegom iz dodiplomskega studija, katere sem pri hitri izdelavi diplomske naloge nenamerno izpustil. Gorazd, Jernej, Marko, Matija, Uros in ze prej omenjena Miha in Vasja, hvala za nepozabna dozivetja in druzenje tekom dodiplomskega studija in v letih po njem. Vesel sem, da smo kljub geografskim razdaljam ohranili stike in da na sedaj ze tradicionalnih srečanjih negujemo prijateljstva. Contents 1 Introduction 1 1.1 Motivation........................................................................1 1.2 Background......................................................................3 1.2.1 Limit load analysis......................................................3 1.2.2 Failure analysis..........................................................5 1.3 Goals of the thesis ..............................................................7 1.4 Outline of the thesis ............................................................8 2 Limit load analysis of reinforced concrete plates 11 2.1 Introduction......................................................................11 2.2 Constitutive model for reinforced concrete plates..............................12 2.2.1 Basic idea................................................................12 2.2.2 Moment curvature relationship ........................................14 2.3 Isotropic and anisotropic reinforcement........................................17 2.3.1 Isotropic reinforcement..................................................17 2.3.2 Anisotropic reinforcement..............................................18 2.4 Numerical examples..............................................................21 2.4.1 Rectangular simply supported plate with anisotropic reinforcement 21 2.4.2 Simply supported square plate with anisotropic reinforcement ... 22 2.4.3 Circular clamped plate with isotropic reinforcement..................23 2.4.4 Plate with two free edges................................................23 2.4.5 Square plate with point supports in the corners......................24 2.5 Concluding remarks and chapter summary....................................26 3 Inelastic analysis of metal plates 29 3.1 Introduction......................................................................29 3.2 Inelastic plate models............................................................29 3.2.1 Plate elastoplasticity....................................................31 3.2.2 Plate elastoviscoplasticity ..............................................34 3.3 Finite element formulation......................................................35 3.3.1 Space discretization......................................................35 3.3.2 Computational issues for plasticity ....................................38 3.3.3 Computational issues for viscoplasticity................................41 3.4 Numerical examples..............................................................41 3.4.1 Limit load analysis of a rectangular plate..............................41 3.4.2 Limit load analysis of a circular plate..................................42 3.4.3 Elastoplastic analysis of a skew plate..................................47 3.4.4 Cyclic analysis of a circular plate......................................48 3.4.5 Elastoviscoplastic analysis of a circular plate..........................50 3.5 Concluding remarks and chapter summary....................................50 4 Inelastic analysis of metal shells 53 4.1 Introduction......................................................................53 4.2 Inelastic geometrically exact shell formulation ................................53 4.2.1 Geometry, kinematics and strains ......................................54 4.2.2 Variational formulation ..................................................57 4.2.3 Stress-resultant constitutive equations for elastoplasticity ............59 4.3 Finite element formulation ......................................................62 4.3.1 Space-domain discretization ............................................62 4.3.2 Computational issues for plasticity ....................................63 4.4 Numerical examples ..............................................................73 4.4.1 Pinched cylinder with isotropic hardening ............................74 4.4.2 Half of a sphere..........................................................76 4.4.3 Limit load analysis of a rectangular plate ..............................76 4.5 Concluding remarks and chapter summary ....................................79 5 Illustration of embedded discontinuity concept for failure analysis 81 5.1 Introduction......................................................................81 5.2 1D finite element with embedded discontinuity................................82 5.3 Computational procedure for failure analysis ..................................90 5.4 Numerical examples ..............................................................93 5.4.1 Tension test of a bar ....................................................94 5.4.2 Tension test of four parallel bars ........................................95 5.5 Concluding remarks and chapter summary ....................................96 6 Failure analysis of metal beams and frames 99 6.1 Introduction................................... 99 6.2 Beam element with embedded discontinuity .................101 6.2.1 Kinematics ...............................101 6.2.2 Equilibrium equations .........................105 6.2.3 Constitutive relations..........................106 6.3 Computation of beam plasticity material parameters ............110 6.4 Computational procedure ...........................114 6.5 Examples ....................................117 6.5.1 Computation of beam plasticity material parameters ........118 6.5.2 Push-over of a symmetric frame .................... 123 6.5.3 Push-over of an asymmetric frame ................... 125 6.5.4 Bending of beam under constant axial force ............. 127 6.5.5 Collapse of a simple frame ....................... 129 6.5.6 Darvall-Mendis frame .......................... 130 6.6 Concluding remarks and chapter summary .................. 135 7 Failure analysis of 2D solids 141 7.1 Introduction ................................... 141 7.2 Family of ED elements for planar problems .................. 141 7.2.1 Kinematics ............................... 142 7.2.2 Equilibrium equations ......................... 147 7.2.3 Constitutive relations .......................... 150 7.3 Computational procedure ........................... 153 7.4 Examples .................................... 159 7.4.1 Tension test ............................... 159 7.4.2 Bending test ............................... 161 7.4.3 Partial tension test ........................... 161 7.4.4 Three point bending test ........................ 162 7.4.5 Four point bending test ........................166 7.4.6 Delamination..............................168 7.4.7 Elasto-plastic tension test.......................170 7.5 Concluding remarks and chapter summary..................173 8 Conclusion 177 Razsirjeni povzetek 181 Bibliography 205 List of Figures 2.1 Positive directions of the transverse displacement and rotations of the plate 12 2.2 Face with tensial stresses (darkened) caused by the positive bending moments 13 2.3 Moment curvature relationship for the principal direction....................14 2.4 Layers of reinforcement and the effective area of the reinforcement in the principal direction ................................................................15 2.5 Equilibrium in the cross-section of the plate..................................16 2.6 Directions of the principal curvatures and the principal moments............18 2.7 Rectangular plate with anisotropic reinforcement under distributed load . . 22 2.8 Simply supported square plate with anisotropic reinforcement ..............23 2.9 Coarse (left) and fine (right) mesh of the circular plate......................24 2.10 Load versus center displacement diagrams for circular plate..................25 2.11 Geometry and loading of the plate with two free edges ......................25 2.12 Bending moment versus displacement at the center curves ..................25 2.13 Force versus displacement diagram for the square plate......................26 3.1 Notation of the used finite element ............................................36 3.2 Load - displacement diagram for simply supported rectangular plate . . . . 42 3.3 Load - displacement diagram for clamped rectangular plate ..................43 3.4 Load - displacement diagram for simply supported rectangular plate; a parameter case ..................................................................43 3.5 Load - displacement diagram for clamped rectangular plate; a parameter case ..............................................................................44 3.6 Meshes used for: (a) rectangular plate - fine, (b) rectangular plate - coarse, (c) circular plate - fine and (d) circular plate - coarse........................45 3.7 Simply supported circular plate - limit load analysis..........................45 3.8 Clamped circular plate - limit load analysis....................................46 3.9 Skew plate - (a) fine mesh, (b) coarse mesh....................................47 3.10 Skew plate - elastoplastic analysis..............................................47 3.11 Spreading of plastic zones ......................................................48 3.12 Clamped circular plate - cyclic load...................... 49 3.13 Loading curve for viscoplastic analyses.................... 50 3.14 Time response for force-prescribed viscoplastic circular plate........ 51 3.15 Time response for displacement-prescribed viscoplastic circular plate .... 51 4.1 A sketch of a two-surface yield function.................... 60 4.2 A sketch of a two-surface yield function and the closest point projection . . 69 4.3 Geometry, loading and boundary conditions for pinched cylinder...... 74 4.4 Load versus displacement curves for pinched cylinder............ 75 4.5 Initial and deformed configuration for pinched cylinder........... 75 4.6 Geometry, loading and boundary conditions of half of a sphere....... 76 4.7 Load versus displacement curves for half of a sphere............. 77 4.8 Initial, deformed fully loaded and deformed unloaded configuration for half of a sphere .................................... 77 4.9 Geometry for rectangular plate ........................ 78 4.10 Load versus displacement curves for the thin plate.............. 78 4.11 Load versus displacement curves of the thick plate.............. 79 5.1 Tension test of an idealized bar........................ 81 5.2 Standard isoparametric 1d solid finite element................ 83 5.3 Isoparametric 1d solid finite element with embedded discontinuity..... 84 5.4 Strain distribution ............................... 85 5.5 Cohesive law at the discontinuity....................... 88 5.6 One point numerical integration scheme.................... 90 5.7 Tension test of one bar............................. 94 5.8 Reaction force versus imposed displacement curve.............. 94 5.9 Tension test of four parallel bars........................ 95 5.10 Reaction versus imposed displacement of structure (left) and structural components (right)............................... 96 6.1 Beam finite element with embedded discontinuity..............101 6.2 Strain-free mode of the element........................103 6.3 Rigid-plastic cohesive law at discontinuity ..................107 6.4 Evaluation of beam material parameters by using results of refined analysis 112 6.5 Uniaxial stress - strain curve..........................119 6.6 Boundary conditions for the shell model analysis...............119 6.7 Failure mode of the representative part of the frame member as computed by the shell model ............................... 120 Doktorska disertacija. Ljubljana, UL, FGG._xvli 6.8 Bending moment versus rotation curves for the end cross-section......120 6.9 Plastic work versus end cross-section rotation curves.............121 6.10 Approximation of the ultimate bending moment of the cross-section .... 122 6.11 Symmetric frame: geometry and loading...................124 6.12 Load versus displacement and dissipated energy versus displacement curves for symmetric frame ............................................................124 6.13 Symmetric frame: deformed shape and locations of softening plastic hinges at utop ~ 60 cm.................................125 6.14 Asymmetric frame: geometry and loading ..................126 6.15 Load versus displacement and dissipated energy versus displacement curves for asymmetric frame..............................126 6.16 Asymmetric frame: Deformed shape and locations of softening plastic hinges 127 6.17 Comparison of results for the bending of the beam under compressive axial force ..............................................................................128 6.18 Simple frame: the beam and the shell model.................130 6.19 Comparison of results for simple frame example...............131 6.20 Deformed shapes of the simple frame.....................131 6.21 Internal forces at the right support of the shell model............132 6.22 Geometry and loading of the portal frame...................133 6.23 Vertical load versus vertical deflection curves.................134 7.1 Quadrilateral finite element with embedded discontinuity..........142 7.2 Different element separation modes......................144 7.3 Neighborhood of a point of interest at the discontinuity re.........151 7.4 Numerical integration scheme.........................153 7.5 Rigid-plastic cohesive law with linear softening in tension..........159 7.6 Tension test on a square block.........................160 7.7 Reaction force versus imposed displacement curves .............160 7.8 Geometry of the bending test of the square block and the imposed displacement versus pseudo-time curves........................161 7.9 Top reaction force versus imposed top displacement curves.........162 7.10 Partial tension test of the square block....................163 7.11 Stress at integration points in and Q- region versus imposed displacement curves ......................................................................163 7.12 Three point bending test of a notched concrete beam............164 7.13 Rigid-plastic cohesive law with exponential softening in tension.......164 7.14 Coarse (top) and fine (bottom) finite element meshes for the three point bending test...................................165 7.15 Reaction force versus imposed displacement curves and scaled (100 times) deformed meshes ................................................................165 7.16 Four point bending test.............................166 7.17 Finite element mesh for the four point bending test.............167 7.18 Load versus crack mouth sliding displacement curves and the corresponding crack paths ...................................167 7.19 Scaled deformed mesh of the "n0 + m0" (left), "n0 + n1" (middle) and "n0" formulation................................168 7.20 Delamination test data.............................169 7.21 Coarse (top) and fine (bottom) finite element meshes for the delamination test ................................................................................169 7.22 Reaction force versus imposed displacement diagram ............170 7.23 Deformed configuration of the coarse (left) and fine (right) mesh......171 7.24 Tension test of a metal strip..........................171 7.25 Rigid-plastic cohesive law with linear softening................171 7.26 Total reaction force versus imposed displacement curves...........172 7.27 Discontinuity paths for several mesh sizes...................172 7.28 Scaled deformed configurations ........................173 List of Tables 3.1 Limit load solutions for circular simply supported plate .......... 44 3.2 Limit load solutions for circular clamped plate ................ 46 5.1 The convegrence results ............................ 95 6.1 Summary of results of the shell model analyses................121 6.2 Comparison between approximations and shell analyses results.......123 6.3 Comparison of the presented formulation with the literature........135 Chapter 1 Introduction The objective of the first chapter is to motivate the reader and to outline the areas, which are aimed to be improved in this work. A brief background on accomplishments in the material nonlinear modeling and modeling of failure is provided. The goals of this thesis are presented and finally an outline of the manuscript is given. 1.1 Motivation From the early days people were faced with a challenge to provide a safe shelter. Once people outgrew the natural shelters and started building man-made structures, the need to acquire the knowledge about materials and structures appeared. For millennia it was only possible to determine the material toughness and the behavior of structures by experimental observations and even nowadays the material properties and estimates regarding structural behavior are determined by experiments. Over the course of human history the understanding of laws of nature increased drastically and the theoretical models started to be used to describe the natural processes. With the uprise of computers in the second half of the twentieth century numerical simulations started playing the leading role in predicting the processes in nature. We are now able to obtain numerical solutions for problems where analytical results can no longer be attained. In particular, we can now simulate the behavior of complex structures, predict their limit load and even describe their behavior until complete failure. Given the short time frame and low cost of getting such computational solutions, when compared to obtaining experimental results, increase its importance ever since. Currently the most widely used tool for analysis of structures is the finite element method. This is a numerical tool to solve differential equations characterizing a certain boundary value problem. The complete three-dimensional analysis of structure, that would have taken into account its exact geometry and all the phenomena related to material geometrically nonlinear behavior (e.g. yielding of steel, cracking of concrete, slip between concrete and steel reinforcement, local failure of material, local and global buckling...) is due to the limited amount of computer power, and due to the limitation of numerical models, describing the above mentioned effects, currently impossible. Therefore a certain degree of mathematical idealization is needed. By considering certain kinematic assumptions we can model the 3D structure as 2D (plates, walls and shells) and 1D (bars, beams) by using plate, 2D solid, shell, bar and bream finite elements. Furthermore we can lower the computational costs by using simplified constitutive models. In this work we derive several finite elements with inelastic material models that are suitable for limit load and failure analysis of different structural components. In order to obtain adequate solutions we derive elasto-plastic material models for analysis of steel beams, plates and shells and a nonlinear elastic model for analysis of reinforced concrete plates. These models are defined at the level of stress resultants, which makes them more robust and computationally cheaper. The novelty of this work is related to representation of material models. For instance, concerning the formulation for nonlinear analysis of reinforced concrete plates, Eurocode 2 recommendations for constitutive behavior of concrete were adopted. When analyzing plastic and viscoplastic stress-resultant plates, a unified computational procedure is presented that covers both formulations. Further, the stress-resultant plasticity for shells was revisited in order to provide a new view on return mapping algorithm for two-surface stress-resultant plasticity. When one wants to determine the behavior of structures when the load produces localized material failure then the elasto-plastic and nonlinear elastic material models no longer suffice. The failure of structure or structure components usually begins due to very localized effects, which the continuum constitutive models are unable to capture adequately. As an example, one can think of a macro crack forming within a concrete wall, where the thickness of the crack is much smaller than the overall dimension of the wall itself. This makes the challenging task for the applied numerical method since the discontinuity in the wall has to be described. In order to simulate such problems, we incorporated singular fields, representing localized zones of failures, into the finite element framework, resulting in a finite element formulation with embedded strong discontinuities. Several contributions to this field are made in this work. The first one provides a new beam finite element with embedded rotational discontinuity, suitable for elastoplastic and failure analysis of planar frames. To determine beam hardening and softening parameters, a sequential multi-scale (shell-beam) procedure is proposed. The next contribution is a new formulation for an elastoplastic 2D solid finite element with embedded discontinuity. A quadrilateral finite element is derived that is suitable for nonlinear collapse analysis of 2D solids and can describe both brittle and ductile failure of the material. 1.2 Background The number of contributions on the limit load analysis, inelastic analysis and the failure analysis in the literature is enormous. In the following a brief overview of the state of the art is given. In Section 1.2.1 we present approaches for the limit load analysis and in Section 1.2.2 we summarize the numerical developments on the modeling of failure. 1.2.1 Limit load analysis Limit load analysis of structures has a significant practical value and there are numerous articles and books written on this topic. The most studied problems are limit load analysis of frames and plates. In the case of simple geometry of a reinforced concrete plate one can determine the limit load by using the plastic lines theory which is explained in numerous books (e.g. [Save and Massonet, 1972], [Nielsen, 1984], [Moy, 1996], [Park and Gamble, 2000], [Radosavljevic and Bajic, 1990]). The basic principle in this method is the assumed plastic (breaking) lines that divide the plate into individual rigid parts that can only rotate around these lines. If we consider that the virtual work done by the fully plastic (breaking) bending moments acting at the plastic lines is equal to the virtual work done by the external load we can derive the limit load of the reinforced concrete plate. Note that this method gives the upper value of the limit load. The limitations of this approach become clear if the geometry of the plate is more complex. Note, that with this approach we only obtain the estimate of the limit load and there are no informations regarding the corresponding limit ductility. For a more complex structural geometry, or when the limit load is not the only relevant quantity (e.g. one wants to know more on the displacement field, load-displacement history, irreversible deformation,...), one can use an inelastic (nonlinear) analysis based on the finite element method. One already available option is the use of commercial programs that allow the nonlinear description of materials, e.g. of the reinforced concrete (e.g. [Hobbit et al., 2007]). When using the commercial software, we divide for example the plate in the thickness direction into layers of concrete and reinforcements, and use the standard elastoplastic material model for steel (see e.g. [Simo and Hughes, 1998], [Kojic and Bathe, 2005], [Lubliner, 1990], [Ibrahimbegovic, 2006], [Ibrahimbegovic, 2009]), to describe the behavior of reinforcement, and a rather sophisticated nonlinear material model to describe the behavior of the concrete. Such a nonlinear analysis requires the use of finite elements for shells and a knowledge of several parameters that describe the behavior of concrete. An alternative to the above approaches was presented in [Ibrahimbegovic et al., 1992], [Ibrahimbegovic and Frey, 1993b] and [Ibrahimbegovic and Frey, 1994]. This method is also based on material nonlinear analysis of reinforced concrete plates, yet it uses a constitutive law that is based on stress resultants instead of stresses, which simplifies the analysis significantly. With this approach we can determine the limit load of the plate in the case of monotonically increasing loading. If the loading is such that includes significant amount of loading and unloading, one should use a different approach to compute the limit load. A similar approach, that is capable of describing the cyclic and even dynamic loading is considered in [Koechlin and Potapov, 2007]. Here, the constitutive relations are also written in terms of stress resultants, in particular a yield function is used, which considers the damage model for concrete and elastoplastic model for steel reinforcement. The standard approach for the inelastic analysis of metal plates and shells is the use of finite elements for shells that require through the thickness integration, see e.g. [Brank et al., 1997]. Such elements are included in most of the commercial finite element based programs (e.g. [Hobbit et al., 2007]). When using the finite elements with through the thickness integration, we divide the shell in the thickness direction into several layers, determine the stress state in each layer by using the standard elastoplastic material model for steel and then by numerical integration along the thickness direction determine the values of stress resultants for a particular cross-section. Such a nonlinear analysis is computational quite expensive, since we need to determine and store the values of internal variables (plastic strain and hardening variables) for all the integration points in the thickness direction. Another way, which is less common, introduces the elastoplastic constitutive equations directly in the 2D stress resultant form; see e.g. [Simo and Kennedy, 1992], [Skallerud et al., 2001], [Crisfield and Peng, 1992] for the shell case. The latter approach is on one hand computationally much faster than the former one, but on the other hand it fails to describe the spreading of the plastification through the plate thickness, see e.g. [Auricchio and Taylor, 1994]. This drawback can be removed (to a certain extent) by a pseudo-time dependent value of the yield parameter associated with the plate bending response. This was first suggested by Crisfield, see e.g. [Crisfield, 1981] and references therein, and was later used e.g. by [Shi and Voyiadjis, 1992] for plates and [Zeng et al., 2001] and [Voyiadjis and Woelke, 2006] for shells. Another way to approximately describe the 3d effect of spreading of plasticity throughout the plate is to use the generalized plasticity model for plates, which is based on two functions (both defined in terms of stress resultants), the yield function and the limit function, see e.g. [Auricchio and Taylor, 1994]. However, if one wants to evaluate only the limit load of the metal plate, the above mentioned modifications are unnecessary, since both the stress resultant formulation and the stress based formulation with the through-the-thickness numerical integration provide the same result. 1.2.2 Failure analysis In the above section we have mentioned some of the approaches that are suitable to determine the behavior of a structure up to the limit load. The question now arises how can one include the effects that are related to material failure. The failure of material is related to the appearance of a discontinuity in the displacement field (e.g. cracks in brittle materials or shear bands in ductile materials) which is in contradiction with the smooth nature of the finite element method. The most trivial approach, that is also incorporated in most commercial finite element programs, is element deletion. The basic idea behind this procedure is simple. Once we detect that the strength of the material has been reached within a single element, e.g. maximum principal stress at the integration point of the element was surpassed, we consider this element to be totally damaged and that it no longer contributes to the resulting system of equations. The discontinuity is in this case represented by the deleted finite element. It is clear that this approach has a drawback of being dependent on the finite element discretization. Instead of deleting the element, once the strength of material is surpassed, one can at the integration points of this element apply a softening stress - strain law in which the stress tends to zero with the growing strain. These procedures are often denoted as softening models. Softening models are problematic since multiple strain states are possible for a certain stress field leading to an ill-posed mathematical problem and similarly as for the element deletion, the finite element solution is pathologically dependent on the mesh size. This drawback can be circumvented if the element size is included in the constitutive relations in terms of the characteristic length, leading to the so called smeared crack approach first proposed in [Rashid, 1968] and is often used to analyze the failure of brittle materials. In this method the continuum model can capture the softening response objectively, since the proper element softening dissipation is achieved by the definition of the softening law being dependent on the element size. The proper energy dissipation distributed over the volume of the element was achieved in [Bazant and Oh, 1983] and [Rots et al., 1985] leading to a band smeared crack approach. Alternative approaches where the characteristic length is also included in the constitutive law are among others: non-local continuum ( [Bazant et al., 1984]), higher-gradient models ( [Coleman and Hodgon, 1985]) and Cosserat continuum ( [de Borst and Sluys, 1991]). In all the above mentioned procedures the actual discontinuity is not present in the model and we only model its influence on the continuum. The first papers that considered the actual discontinuity as a geometric entity in the finite element formulation are [Ngo and Scordelis, 1967] and later [Hillerborg, 1976] introducing the discrete crack approach. Within this formulation the force - displacement relations are enforced along the finite element edges, in particular at the element nodes. Once the nodal force at the node ahead of the discontinuity tip exceeds the strength of the material the node splits into two new nodes and the discontinuity can now propagate to the next node. The position of the discontinuity is thus restricted to the finite element boundary which makes the procedure mesh dependent. A better result is obtained if the crack path is predefined with an appropriate mesh alignment or if the remeshing is used, see e.g. [Ingraffea and Saouma, 1985]. In [de Borst et al., 2004] a comparison of smeared and discrete crack approaches is given. Similar to the discrete crack approach is the cohesive zone model developed in [Needle-man, 1987]. In this approach the boundaries of the finite elements are modeled by the cohesive finite elements that have a certain traction-separation law incorporated. The formulation was further developed in [Tvergaard, 1990] to include both normal and tangential separation and a further improvement in [Xu and Needleman, 1994] allowed to deal with the crack branching. If the cohesive elements are included in the finite element mesh at the beginning of the simulation, the solution is influenced by cohesive elements' stiffness. This can be avoided if the elements are included only after reaching a failure criterion, see e.g. [Camacho and Ortiz, 1996] and [Ortiz and Pandolfi, 1999]. Similarly as for the discrete crack approach the cohesive zone model also suffers to be mesh objective and several techniques to diminish this drawback can be used, see e.g. [Carter et al., 2000], [Bittencourt et al., 1996], [Ortiz and Quigley, 1991], [Marusich and Ortiz, 1995] and [Ingraffea and Saouma, 1985]. In [Ortiz et al., 1987], the first time the discontinuity in the strain field was introduced into the bulk of a finite element to employ the strain localization. The displacement field in this approach remained continuous and the method became known as the weak discontinuity approach. The strong discontinuity approach, in a sense, that the displacement field in the finite element framework is discontinues itself, was for the first time accomplished in [Dvorkin et al., 1990] and [Simo et al., 1993] and it was further developed e.g. in [Simo and Oliver, 1994], [Oliver and Simo, 1994], [Oliver, 1995], [Oliver, 1996a] and [Oliver, 1996b]. The main advantage of the method is in the way it incorporates the discontinuity in the finite element framework. Namely, once the discontinuity propagation criterion was reached in a particular finite element, additional parameters are introduced in this finite element only. These local unknowns, which one can interpret as the jumps in the displacement field, can be condensed out at the element level, thus leaving the global system of equations unchanged. Since the strong discontinuity with all the related parameters is embedded within individual finite elements, this method is often denoted as the embedded discontinuity finite element method (ED-FEM). ED-FEM approach received a lot of attention in the past two decades and numerous authors took part in developing and improving it in various fields of finite elements. The embedded discontinuity approach was successfully included in the beam finite element formulation in the works of [Ehrlich and Armero, 2005], [Armero and Ehrlich, 2006b], [Armero and Ehrlich, 2004], [Wackerfuss, 2008] and [Dujc et al., 2009] among others. The theoretical formulation of the localized failure in the Reisner-Mindlin plate theory was presented in [Armero and Ehrlich, 2006a] and built into triangular and quadrilateral plate finite elements. The theoretical formulation of combined hardening and softening plasticity model was presented in [Ibrahimbegovic and Brancherie, 2003] and applied into 2D solid finite elements. This formulation along with the formulation that includes damage constitutive model for the bulk response and the damage constitutive model for the discontinuity was presented in [Brancherie, 2003]. The modeling of failure in solids was also studied in [Brancherie and Ibrahimbegovic, 2008], [Manzoli and Shing, 2006], [Linder and Armero, 2007], [Dujc et al., 2010] and [Linder, 2007], where also a detailed review on the history of theoretical and numerical modeling of failure is presented. Recently, in [Ibrahimbegovic and Melnyk, 2007] the approach of modeling of localized failure in heterogeneous materials was considered where both strong and weak discontinuity approach was incorporated. For a review on approaches to derivation of embedded discontinuity model we also refer to [Jirasek, 2000] and [Jirasek and Zimmermann, 2001]. Similar to ED-FEM approach is the extended finite element method (X-FEM), whose concepts are based on the partition of unity ( [Melenk and Babuska, 1996] and [Babuska and Melenk, 1997]). Contrary to ED-FEM, the parameters describing the behavior of discontinuity are in X-FEM treated globally. This means that the number of global equilibrium equations is not fixed during the simulation but it grows with the formation of the discontinuity path. A comparison of the ED-FEM and the X-FEM is for example given in [Jirasek, 2000] and [Ibrahimbegovic and Melnyk, 2007]. 1.3 Goals of the thesis The particular goals of the thesis follow the main goal, which is to study and develop numerical tools for analysis of limit load, limit ductility and complete failure of structures and structural components. The particular goals can be summarized as follows: • The first goal is to apply the Eurocode 2 [Eurocode 2, 2004] recommendations to de- scribe the constitutive behavior of reinforced concrete plates presented in [Ibrahim-begovic et al., 1992], [Ibrahimbegovic and Frey, 1993b] and [Ibrahimbegovic and Frey, 1994], and to include the resulting constitutive model into the AceGen [Ko-relc, 2007b] code of the existing Reissner-Mindlin plate finite element presented in [Bohinc et al., 2009]. • The second goal is to derive the small strain stress resultant elastoplastic plate finite element that includes both isotropic and kinematic hardening, further extend the plasticity formulation into the visoplasticy formulation of Perzyna type and include the constitutive relations into the AceGen code [Korelc, 2007b] of the existing Reissner-Mindlin plate finite element presented in [Bohinc et al., 2009]. • The third goal is to derive the inelastic geometrically exact shell formulation and upgrade the existing shell finite element (see [Brank et al., 1995] and [Brank and Ibrahimbegovic, 2001]) programed in AceGen [Korelc, 2007b], to include the return mapping algorithms needed for shell stress-resultant plasticity with isotropic and kinematic hardening. • Fourth goal is to derive a new elastoplastic Euler-Bernoulli beam finite element with embedded discontinuity for push over analysis of steel frame structures. In this context an approach to determine the material parameters needed for an analysis of newly derived Euler-Bernoulli beam had to be derived. • The last goal of the thesis is to derive the family of quadrilateral finite elements for analysis of 2D solids, including both hardening and the localized softening effects, where the opening of the discontinuity is considered as linear in normal and in tangential direction. The element should be suitable for failure analysis of brittle as well as ductile solids. 1.4 Outline of the thesis In addition to Section 1 the thesis consists of seven chapters, where we consider different topics related to the limit load analysis and the failure analysis of structures. The limit load analysis of reinforced concrete plates is considered in Chapter 2. Kinematic quantities with dual strain variables related to plate boundary value problem are presented, followed by a detailed description of nonlinear elastic constitutive relations for reinforced concrete plates. With respect to the pattern of reinforcement we differ between the case of isotropic reinforcement and the case of anisotropic reinforcement. Concep- tual algorithms for both cases are presented, followed by several numerical simulations. Concluding remarks and a short summary end this chapter. In Chapter 3 we derive small strain elastoplastic plate finite element formulation in terms of stress resultants. Nonlinear isotropic and linear kinematic hardening are considered. We further extend the plasticity formulation into the visoplasticy formulation of Perzyna type. Both elastoplastic and elastoviscoplastic stress resultant plate formulations are derived by exploiting the hypotheses of instantaneous elastic response and the principle of maximum plastic dissipation (plasticity) or the penalty-like form of the principle of maximum plastic dissipation (viscoplasticity). The performance of the finite element is shown by several numerical examples. We conclude the chapter with concluding remarks and a short summary. Elastoplastic analysis of thin metal shells is presented in Chapter 4. First we describe the kinematics and the variational formulation of the geometrically exact shell. By the thermodynamics consideration and the use of modified version of the classical Ilyushin-Shapiro yield function we derive the evolution equations for generalized plastic strain, internal variable related to isotropic hardening and the internal variables related to kinematic hardening. Next we present the spatial finite element discretization and computational algorithms related to shell stress resultant multi-surface plasticity. The results of several numerical simulations, accompanied by the results from literature are further presented. The chapter closes with concluding remarks and a short summary. In Chapter 5 we present the finite element method with embedded discontinuity for failure analysis in 1D solids. With the introduction of one additional kinematic parameter into the standard 1D solid isoparametric finite element, we derive the enriched strain field and by considering the principle of virtual work we obtain the mesh related equilibrium equations along with the local element-wise equations that arise due to strain enrichment. The response of the bulk material is linear elastic, while the response of the discontinuity is considered to be rigid plastic, where all the necessary ingredients are obtained by considering the principles of thermodynamics. Next we present the computational procedure that determines the new values of nodal displacements and the new values of local (element) variables. Two numerical simulations are presented. We conclude the chapter with concluding remarks and a short summary. We note that this chapter is included in the thesis only for illustrative purposes. The principles that were presented in Chapter 5 for 1D solids are in Chapter 6 incorporated in elastoplastic Euler-Bernoulli beam finite element. First we consider the kinematics of Euler-Bernoulli beam finite element with discontinuities in the rotation field and in the axial displacement field and derive the operators needed to enrich the generalized strain field. Again we use the principle of virtual work to obtain the global and the local equilibrium equations. Next we considered the elasto-plastic response of the bulk material of the beam and a rigid plastic response of the discontinuity. The remaining ingredients of the elastoplasticity with hardening and the rigid-plasticity with softening are obtained from the consideration of thermodynamics of associative plasticity and the principle of maximum plastic dissipation. Next we discuss computation of the beam hardening and softening parameters by using the shell model. We present the details of the computational procedure and show the performance of our approach by several numerical simulations. Concluding remarks and a short summary close the chapter. The failure analysis of 2D solids is presented in Chapter 7. Within the standard isoparametric quadrilateral finite element we introduce a discontinuity line with four additional parameters representing four modes of separation. We determine the interpolation matrices related to the additional parameters and with the derivation of the enriched displacement field we obtain the strain field and the strain field operators related to the parameters of the discontinuity. By considering that the additional parameters can be viewed as the incompatible degrees of freedom, we obtain the operators related to virtual strain and with them by considering the principle of virtual work we derive the local and the global equilibrium equations. Next we considered the elasto-plastic response of the bulk material and a rigid plastic response of the discontinuity. All the remaining ingredients are obtained from the consideration of thermodynamics of associative plasticity and the principle of maximum plastic dissipation. In Section 7.3 we present the procedure to determine the new values of the nodal displacement, hardening plasticity variables related to integration points of the element and the local softening plasticity variables related to integration points of the discontinuity along with the procedure that determines the discontinuity propagation within the finite element mesh. Several numerical simulations are presented. We conclude this chapter with concluding remarks and a short summary. Concluding remarks of the thesis are given in Chapter 8. Remark 1.1. In all chapters, the notation is specialized to the topic of the particular chapter. Therefore the validity of notation ends with the chapter end. Chapter 2 Limit load analysis of reinforced concrete plates 2.1 Introduction In this chapter we present the approach for a nonlinear analysis of concrete plates by using the finite element method. We use a constitutive law that is based on stress resultants instead of stress, which makes the analysis computationally much cheaper. The idea is taken from [Ibrahimbegovic et al., 1992], [Ibrahimbegovic and Frey, 1993b], [Ibrahimbe-govic and Frey, 1994]; the difference is that we apply the Eurocode 2 [Eurocode 2, 2004] recommendations to describe the constitutive behavior of reinforced concrete. With this approach we can determine the limit load of the plate in the case of mono-tonically increasing loading. If the loading is such that includes significant amount of loading and unloading one should use a different approach to compute the limit load. We consider that the local effects related to the change of stiffness of plate and the redistribution of loading are small and negligible. We also assume that the displacements of the plate are small enough, that we can neglect the membrane forces in the plate. The chapter is organized as follows: in Section 2.2 we present the kinematics of the Reissner-Mindlin plate model, define the stress resultants and describe the details of the nonlinear elastic constitutive relations for reinforced concrete plates. With respect to the pattern of reinforcement we distinguish between the case of isotropic reinforcement and the case of anisotropic reinforcement and present the computational algorithms for both cases, in Section 2.3. In Section 2.4 several numerical simulations are presented. Finally, in Section 2.5 concluding remarks and a short summary are provided. 2.2 Constitutive model for reinforced concrete plates 2.2.1 Basic idea According to Reissner-Mindlin plate theory we model the plate as a surface in the x1 x2 plane that has at each point three degrees of freedom: displacement in the x3 direction and two rotations of the plate normal around axis x1 and axis x2 (see Figure 2.1). We can then express the curvature vector k and the vector of transverse shear deformations Y as r de2 dei dei de2 T r 2 ir K = --TT-1 = [«11, «22, 2k^] , (2.1) dx1 dx2 dx1 dx2 r dw „ dw „ lT . Y = r axr + e2-ax2- e1]T <2'2) With the proposed constitutive law for reinforced concrete we relate the total values of deformations (curvatures and shear strains) with the plate stress resultants m = [m11,m22,m,12]T , q = [q1,q2]T, (2.3) where m is the vector of bending moments and q is the vector of shear forces per unit length. With the chosen coordinate system and the chosen cinematic and constitutive relations we obtain positive values of bending moments, when we have tensial stress on the lower face of the plate, t.i. the face with the normal n in the —x3 direction (Figure 2.2). Figure 2.1: Positive directions of the transverse displacement and rotations of the plate Slika 2.1: Pozitivne smeri precnega pomika in zasukov plosce We treat the bending and the shear part of the constitutive law separately. The shear forces are computed according to standard linear elastic constitutive relations for isotropic material q = C s Y , C 5 Ech s = 6 2(1 + vc) 1 0 0 1 (2.4) Figure 2.2: Face with tensial stresses (darkened) caused by the positive bending moments Slika 2.2: Nategnjena stran plosCe (potemnjena), ko so momenti pozitivni where Ec is the Young's modulus of the concrete, h is the thickness of the plate and vc is the Poisson ratio of the concrete. Therefore, the derived model for reinforced concrete plate is unable to describe the shear failure of the plate. Regarding bending we consider two distinguished states: I the state before cracking of the concrete and II the cracked state, when the cracks due to tensial failure of the concrete are already present in the plate. For the intact concrete (state I) we consider the linear elastic relations between the bending moments and the curvatures. If we neglect the contribution of the reinforcement we get m = C B k , C Ech3 B 12(1 - Vc2) 1 Vc 0 Vc 1 0 0 0 1-V 2 (2.5) The state II starts when one of the moments in one of the prescribed orthogonal directions, determined with the angles 0 and 0 + 2 (see Figure 2.4), reaches the value mcrack = /cthr, where /ct is the tensial strength of the concrete. In other words, we check for the crack initiation in the two orthogonal directions, where one can interpret mcrack as the bending moment that initiates crack in the planar beam with the unit width. The two orthogonal directions can be either in the directions of the principal curvatures,in the directions of the principal bending moments or in the directions of the maximum resistance of the plate. More on the choice of crack directions will be explained in Section 2.3. Even in the subsequent phase of the loading, after the cracking has been initiated (in state II), we check for the progress of cracks along the two orthogonal directions, that are not necessary equal to the directions at the first appearance of the crack. In the cracked state we disregard the influence of Poisson's ratio and assume independent responses for the two orthogonal directions. Further details will be addressed in Section 2.2.2 below. c 2.2.2 Moment curvature relationship In what fallows we present the approach that gives us the relationship between the bending moment m and the corresponding curvature K for the direction that is determined with the angle 0 (Figure 2.4). First we provide the values of m and K for three distinguished states of the reinforced concrete cross-section (Figure 2.3): • appearance of the first crack in the concrete (point A, m = mcrack, K = Kcrack), • start of yielding of the reinforcement (point B, m = my, K = Ky), • and the failure of the concrete in compression (point C, m = m/, K = K/). Once we have points A, B and C determined we simply connect them and obtain the peace-wise linear constitutive relation. The key parameter that determines the shape Figure 2.3: Moment curvature relationship for the principal direction Slika 2.3: Diagram moment ukrivljenost za glavno smer of the diagram in Figure 2.3 is the effective area of the reinforcement in the principal direction determined with the angle 0, which we compute according to a^ = ^^ ai cos2(0 — a), (2.6) where ai is the area of the reinforcement layer i with the orientation angle ai (Figure 2.4). In order to obtain the points A, B and C we adopt the Eurocode 2 recommendations for the constitutive relations for concrete and steel. For the concrete in compression we use the following relation between the stress ac and strain ec Gc = f fck(1 — (1 — %)n) za 0 < £c < £c2 , (2.7) 1 fck za ^c2 < £c < £cu3 a2 Figure 2.4: Layers of reinforcement and the effective area of the reinforcement in the principal direction Slika 2.4: Plasti armature in efektivna kolicina armature v glavni smeri where Sc2 = 0.002, scu3 = 0.0035, n = 2 for concrete fck < 50 and fck is the compressive strength of concrete. The tensile strength and the corresponding strain are according to Eurocode 2 fct = fct(fck) and ect = j- ■ We adopt the standard bilinear relationship between the reinforcement strain ss and the reinforcement stress as ■ Es^s Za S s fy ^ H s s s " j , (2.8) Jy za Ss > e where fy is the tensial strength of steel and Es = 20000kN/cm2 is the Young's modulus of steel. Characteristic points for moment curvature diagrams In order to obtain the characteristic points A, B and C in Figure 2.3 we need the following data: the tensile strength of steel fy, the Young's modulus of steel Es, the effective area of the reinforcement a^[cm2/cm], the compressive strength of concrete fck, the tensile strength of concrete fct or its corresponding strain sct = , the Young's modulus of concrete Ec, the thickness of the plate h and the effective height of the plate d = h — a, where a is the distance from the bottom face of the plate to the reinforcement. In Figure 2.5 Scc denotes the strain at the face in compression, sct is the strain at the face in tension, Ss is the strain at the position of reinforcement, Fc is the resultant force arising from the compressive stress in concrete, Fs is the force in the reinforcement and x is the distance between the neutral axis and the plate mid-plane. Bending moment and the corresponding curvature at the crack initiation are mcracfc fct— j Kcrack 22 ^ t . (2.9) 6 h Ech ecc top face Qt bottom face Figure 2.5: Equilibrium in the cross-section of the plate Slika 2.5: Ravnotežje v prerezu ploSce The strain in the reinforcement that corresponds to the start of yielding is equal to f es = Ey and the force in the reinforcement at that time is Fs = fyBy considering, that the compressive force in the concrete is equal to the force in reinforcement /X ^ acdx3, (2.10) we obtain the position of the neutral axis x that satisfies the equilibrium in (2.10). We can then provide the edge (face) values of the strains in concrete ecc(x,es) and ect(x,es) (Figure 2.5). The values of the bending moment and the curvature at the start of yielding of the reinforcement are then my = j ^ acX3dx3 + (h - a)Fs , Ky = £ct + e =1/2arctan(K12—1K22) + *§ = { 0 ^ T* < 0 • (2.32) 4. Compute the effective area of reinforcement for both directions «1 = a^ = ^2 ai cos2(0 — ai), a2 = 2 = ^ a cos2(0 + 2 — ai). (2.33) 5. Determine the relationships m1(K1) and m2 (k2). 6. Determine the values of principal moments m1 = m1(«1), m2 = m2(«2)- (2.34) 7. Compute the shear strains and shear forces in the X1X2X3 coordinate system Y = Y (u) = [71,72]T, q =[<1,<2]T = CsY • (2.35) 8. Use the principle of virtual work as the basis for the finite element method analysis OT(u, = / (m18n1 + m28n2 + q1 Sj1 + q2č>Y2)dQ — ^^ / pdwdQ. (2.36) 2.4 Numerical examples The presented constitutive model and the corresponding algorithms were compiled with the quadrilateral finite element for plates with quadratic interpolation of transverse displacement w and linear interpolation of rotations 01 and 02, which was presented in [Bohinc and Ibrahimbegovic, 2005] and is described in detail in Section 3.3.1. The computer code was generated by using symbolic manipulation code AceGen and the examples were computed by using finite element program AceFem, both developed by Korelc, see [Korelc, 2002], [Korelc, 2007b] and [Korelc, 2007a]. In the cases with isotropic reinforcement the results obtained with the isotropic model are completely the same as the ones obtained with the anisotropic model with rotating crack. The only difference is in the computation time, which is shorter for the isotropic model, since less numerical operations are needed. It is assumed in the following examples that the positive or negative reinforcement is the same at any point of the plate. 2.4.1 Rectangular simply supported plate with anisotropic reinforcement We consider a rectangular simply supported plate under distributed load. The properties of the plate are: the thickness h = 80 mm, the length l = 300 cm and the width b = 200 cm. The plate reinforcement is orthogonal with the area of the reinforcement in the first direction a1 = 251 mm2/m and in the second direction a2 = 559 mm2/m. The position of reinforcement is at a distance c1 = c2 = 14 mm from the bottom face of the plate. The material properties are: Ec = 24 kN/mm2, fck = 26.5 N/mm2, fct = 2.5 N/mm2, Es = 205 kN/mm2 in fy = 460 N/mm2. The finite element mesh used in our simulations consist of 8 x 8 finite elements. Two numerical simulations were performed; one with the fixed crack approach and the second with the rotating crack approach. In Figure 2.7 we present the total load versus displacement at the center obtained in numerical simulations along with experimental results adopted from [Ibrahimbegovic and Frey, 1994]. We can see from Figure 2.7 that there is a reasonable agreement between the rotating crack approach and the experimental data and that the prediction of the limit load is fairly accurate whereas the fixed crack approach significantly over-estimates the true resistance of the plate. The reason for that is that the directions of principal action do not coincide with the directions of the plate reinforcement. The fixed crack approach would give reasnoble results only in the case, where these two directions would coincide, i.e. in the case of optimally designed plate. Since this approach in general leads to a large over-estimate of the plate resistance it is not suitable for numerical simulations. Displacement at the center [mm] Figure 2.7: Rectangular plate with anisotropic reinforcement under distributed load Slika 2.7: Pravokotna anizotropno armirana plosca pod vplivom ploskovne obtežbe 2.4.2 Simply supported square plate with anisotropic reinforcement In this example we consider a simply supported square plate with anisotropic reinforcement under distributed loading. The plate's thickness is h = 5.1 cm and the side of the plate is a = 183 cm. The plate reinforcement is orthogonal where the area of the reinforcement in the first direction is a1 = 2.81 cm2/m and in the second direction is a2 = 2.35 cm2/m. The position of reinforcement is determined with d1 = 3.9 cm for the first direction and with d2 = 4.39 cm for the second direction, where dj is the distance form the top face of the plate to the position of the reinforcement. The material properties are: Ec = 3242 kN/cm2, fck = 3.5 kN/cm2, fct = 0.379 kN/cm2, vc = 0.18, Es = 20691 kN/cm2 in fy = 37.59 kN/cm2. The numerical simulations were performed with the mesh of 8 x 8 finite elements. In the simulations we assumed that both layers of reinforcement are positioned at the same distance from the top face of the plate: d = di+ d2 = 4.15 cm. In the corner nodes of the plate we did not support the vertical displacement which was in agreement with the experimental boundary conditions, see [Ibrahimbegovic et al., 1992]. Figure 2.8 presents the load versus displacement at the center curves. We can see that there is a good agreement of results in the region of small displacements (up to cca. 20 mm). When the displacements get bigger, the membrane effects significantly increase the resistance of the plate, which is not accounted for in our model. Displacement at the center [mm] Figure 2.8: Simply supported square plate with anisotropic reinforcement Slika 2.8: Kvadratna anizotropno armirana plosca pod vplivom ploskovne obtežbe 2.4.3 Circular clamped plate with isotropic reinforcement We consider a circular plate with isotropic reinforcement. The thickness of the plate is h =1 m and its radius is R =10 m. The material properties are: Ec = 2000 kN/cm2, fck = 3.5 kN/cm2, fct = 0.56 kN/cm2, vc = 0.16, Es = 21000 kN/cm2 in fy = 46 kN/cm2. The plate has an isotropic reinforcement with y = =1% reinforcement in each direction, placed at a = a! = 3 cm from the top and bottom face of the plate. In Figure 2.9 we present the coarse (left) and the fine (right) finite element mesh that we used in the numerical simulations. The plate is built in, i.e. the edge nodes have prescribed all the degrees of freedom (w = 01 = 92 = 0). In Figure 2.10 we present the load versus transverse displacement at the center curves. Along with the results of our simulations we also present the results obtained by the 3-D solid finite elements (see [Ibrahimbegovic et al., 1992]). We can see that there is hardly any difference in results obtained by the course and the fine mesh and despite the simplicity of our approach, we have a good agreement with the results obtained by the more sophisticated approach. 2.4.4 Plate with two free edges In this example we consider a rectangular plate that is simply supported along the shorter edges while the longer edges are free, see Figure 2.11. The geometric properties of the plate are: the length of the shorter edges l1 = 45.72 cm, the length of the longer edges l2 = 76.2 cm, the thickness h = 3.81 cm and the distance from the top face of the plate Figure 2.9: Coarse (left) and fine (right) mesh of the circular plate Slika 2.9: Groba (levo) in fina (desno) mreža koncnih elementov pri krožni plosci to the reinforcement d =3.1 cm. The material properties are: Ec = 2900 kN/cm2, fck = 3.2 kN/cm2, fct = 0.2 kN/cm2, vc = 0.18, Es = 20000 kN/cm2 and fy = 22 kN/cm2. The plate is loaded with the line load P (Figure 2.11). The reinforcement is placed only in the direction parallel to the longer edges where a1 = 2, 74 cm2/m. The numerical simulation is performed with the finite element mesh of 10 x 10 elements. In Figure 2.12 we compare the results of our approach with the results obtained by the commercial program Abaqus [Hobbit et al., 2007], where the behavior of the reinforcement was described with the elastoplastic material model for steel and the behavior of the concrete was described with the damage material model. We can see that there is a reasonable agreement in the results. 2.4.5 Square plate with point supports in the corners We consider a square plate with point supports in the corners. The plate is loaded with the point load at the center of the plate. The plate's properties are: the thickness h = 4.4 cm, the length of the side of the plate l = 91.4 cm, the amount of isotropic reinforcement ^ = 0.85% and the distance from the top face of the plate to the position of reinforcement d =3.3 cm. The properties of steel and concrete are: fck = 3.8 kN/cm2, Ec = 3600 kN/cm2, vc = 0.15, fct = 0.53 kN/cm2, Es = 20000 kN/cm2 and fy = 34.5 kN/cm2. The numerical simulation is performed with a mesh of 20 x 20 elements. Center displacement [cm] Figure 2.10: Load versus center displacement diagrams for circular plate Slika 2.10: Krivulje obtežba - pomik pri krožni plosci P P A" "A a : b =2:6 b a a l 2 Figure 2.11: Geometry and loading of the plate with two free edges Slika 2.11: Geometrija in obtežba plosce z dvema prostima robovoma n e c e h M ni di n e m 1.8 1.6 1.4 1.2 1.0 0.8 0.6 0.4 0.2 0 -ROTATING CRACK --- ABAQUS H-1-1-1-1-1-1-1-1—- 0123456789 Displacement at the center [mm] Figure 2.12: Bending moment versus displacement at the center curves Slika 2.12: Krivulje upogibnega momenta v odvisnosti od precnega pomika na sredini plosce Displacement [cm] Figure 2.13: Force versus displacement diagram for the square plate Slika 2.13: Krivulje sila pomik pri kvadratni plosci In Figure 2.13 we present the transverse displacement of the point that is at a distance 10 from the center point of the plate in the x1 direction. This example also shows that the presented approach gives results that are in reasnoble agreement with the experimentally obtained results (see [Zahlten, 1993]). 2.5 Concluding remarks and chapter summary In this chapter the analysis of reinforced concrete plates is considered where we distinguish between the uncracked state I of the plate and the cracked state II of the plate, activated by the tensile failure of the concrete. In state I we consider the standard linear elastic constitutive relations for isotropic material q = CsY, m = Cb k. In state II we again consider a linear elastic response for the shear while for the bending part we disregard the influence of Poisson's ratio and assume independent responses for the two orthogonal directions defined by 0 and 0 + | of the form m = m (k). The response in those two directions is determined with the effective area of reinforcement in that direction a

b is virtual work of external moments ma and external forces qa acting on the plate boundary, and (0) is virtual quantity that corresponds to (o). We consider displacement w, rotations 9Y, stress resultants ma3, qa and load as functions of position x = [x1,x2]T G Q and pseudo-time t G [0,T], i.e. w = w (x, t), 9j = 9j (x,t), ma/ = ma/ (x,t), qa = qa (x,t), p = p (x, t), ma = ma (x, t), qa = qa (x,t). Equation (3.1) can be written in matrix form as G (m, q; w, 0) = KT mdQ + 0T qdQ — wpdQ — Gext,b = 0, where the following mappings have been defined (3.4) Ka3 Ya ma/ 0 = [91, 92]t , K = [K11 ,K22, 2K12] dw T d92 d91 d91 d9i dx1 dx2 dx1 dx2 T Y = [Y1 j Y2] T m dx1 [m11, m22, m12]T + 92j -d--91 dx2 qa ^ q T (3.5) [q1jq2]T. For further use we also define the following strain and stress resultant vectors £ 0- (315) which is assumed to be non-negative. Note that equation (3.15) can be derived from the second law of thermodynamics, see e.g. [Khan and Huang, 1995], [Simo and Kennedy, 1992]. By assuming that the elastic process is non-dissipative (i.e. the state variables do not change during that process and V = 0) one has a = ^ = Cie. (3.16) die v 7 By further consideration of (3.15) one can define the dual variables, i.e. the hardening variable q and the variables that control kinematic hardening a, as d4 d~(0 d4 2 H 2 H q = -^ = = -- (i) , a = -^ = -3= - 3HkinX, (3.17) where k = [k11, k22, k12, k13, k23]t. By using (3.16) and (3.17) in (3.15) we obtain the reduced material dissipation (i.e. the dissipation of the plastic process) as Vp = aTip + qi + aTk ^ 0. (3.18) The principle of maximum plastic dissipation states that among all the dual variables (a,q, a) that satisfy the yield criteria, one should choose those that maximize plastic dissipation. The problem can be written in the following form: Find minimum of Lp (a,q, a,Y), where Lp (a, q, a, j) = -Vp (a, q, a) + (a, q, a), (3.19) and y ^ 0 plays the role of Lagrange multiplier. From the above minimization problem and (3.11) we obtain explicit forms of evolution equations for the internal variables d CP d 6 = -ip + 7d6 = q=^ ip = Yy2A (a + a), d a d a --v-' V dLP = -£ + # = Q=^ č = 7-fl - (3==1) 7- V(a + a)T A (a + a), (3.2Q) ay V °y P P dLp .06 . -r— = -x + 7^- = Q x = 72A (a + a). d a d a --v-' V Note that k = ip. Remark 3.3. Equation (3.20)2 is generalization of the equivalent plastic work variable Wp = aTep. Namely, by inserting x = ep and (3.20)1 into (3.18), and using (3.11), one gets Dp = 72 - , which implies, see (3.20)2, that £ = ^Dp. The loading/unloading conditions follow from the demands that 7 is non-negative, 0 is non-positive, and the plastic dissipation Dp equals zero for elastic process when 0 < 0 7 ^ 0, 0 ^ 0, 70 = 0. (3.21) In addition to (3.21) we have the condition 0 = 0 if 7 > 0 (the consistency condition). It guarantees the admissibility of the subsequent state in the case of change of state variables. The consistency condition 7>0; 0 = O=(^a- + md, (3.22) pseudo-time derivatives of (3.16) and (3.17), a = C (e — ep), q = S" (£) £, a = — |Hkink and equations (3.20) lead to the following expression for 7 Y = (vT Cv + (£)(32 + 3HkmvtV) vT(3.23) If (3.23) and (3.20) are used in a = C (e — ep), one can write a = Cepe, where f C if 7 = 0 Cep = <^ C__cvvTC if . > 0 (3.24) [ C vtcv+s"(t)fp + 2Hfcinvtv if ' > 0 is elastoplastic tangent modulus of the elastoplastic plate model. By computation of the internal variables (i.e. by integration of (3.20)) and by using (3.16) one recovers the stress resultants a (x, t) = [mT, qT]T appearing in the weak form of the equilibrium equations (3.4). 3.2.2 Plate elastoviscoplasticity A stress resultant viscoplastic constitutive equations for plates of Perzyna type are obtained by a modification of the elastoplasticity model presented in the previous section. The basic difference between the viscoplasticity and plasticity is that in the former model the stress states {a,q, a}, such that 0 (a,q, a) > 0, are permissible, while in the latter are not. The state variables remain the same, except for the viscoplastic strain evp, which replaces ep. The constrained minimization problem (3.19) for plasticity is here replaced by the penalty form of the principle of maximum plastic dissipation (see e.g. [Simo and Hughes, 1998] section 2.7 and [Ibrahimbegovic et al., 1998] for details), which can be stated as: Find minimum of Lvp (a, q, a), where Lvp (a, q, a) = -Dvp (a, q, a) + 1 g (0 (a, q, a)), (3.25) n n G (0, to) is penalty parameter (also called viscosity coefficient or fluidity parameter), Dvp is viscoplastic dissipation of the same form as (3.18) and g is penalized functional. A usual choice for g is ... f 102 if 0 > 0 g (0) = \ 0 if 0< 0 • (3.26) With this choice for g the minimization of (3.25) leads to ^ = -.vp +1 (,) f = 0, ^ = -S +1 (« £ = 0, d-f- = -K + i (« f = 0, (3.27) da n da dq n dq da n da where (0) = g' (0) = dg. Equations (3.27) provide the corresponding evolution equations of the state variables for the viscoplastic plate model. If one defines the viscoplastic multiplier y as y = 1 (0) then the evolution equations for viscoplastic model can be written as those for elastoplastic model, see (3.20). Their integration, which leads to stress resultants a (x,t) = [mT, qT]T, can be performed in very similar manner as for plasticity, as shown below. 3.3 Finite element formulation 3.3.1 Space discretization In the finite element solution of the plate bending problem a domain under consideration Q is discretized by a mesh of finite elements so that Qk = (J^ Qe, where nel is number of elements in the mesh. In this work we use one of the quadrilateral plate elements originally introduced in [Ibrahimbegovic, 1993]. Its geometry is defined by the bilinear mapping | ^ xk (| G [-1,1] x [-1,1]; xk G Qe) with xk (I) |ne =£ N (|) Xa xa [x1a, x2a] T I = [e\e2] T a=1 where xa are coordinates of the finite element node a and (3.28) a1 Na (e\e2) = 1(1+eae ^ (1+eae2) ea -1 1 34 11 ea -1 -111 (3.29) 2 The subscript h is used to denote the discretely approximated quantities. Interpolation of the rotations is based on bilinear polynomials (3.29) (fh) = 0h (I, t) |ne = ^ Na (I) 0a (t) , 0a = [91a, 92a]T ^ 2 ' a=1 (3.30) while interpolation of the transverse displacement is performed in more elaborated way as w (I, t) |ne= £ Na (I) wa (t) + J] Ne (I) JKnTjK (0j (t) — 0k (t)) 8 (3.31) a=1 E=5 The second term in (3.31) is such that the shear distribution along each element edge is constant. In (3.31) Ijk = {(x^ — xw)2 + (x2K — x2j)2)1/2, njK = [cos ajK, sinajK]T (see Figure 3.1) and ne (i) = 1 (1 — e1)2 (1+6e2), e = 5,7 Ne (|) = 1 (1 — C2)2 (1+ 6C1), E = 6, 8 E 5 6 7 8 J 12 3 4 K 2 3 4 1 Interpolation of bending strains follows from (3.5)2 and (3.30) (3.32) X2 Xi X3 ■e Figure 3.1: Notation of the used finite element Slika 3.1: Notacija pri uporabljenem koncnem elementu Kh (I,t) |ne =Y, Ba (I) 0a (t) B n a=1 0N a,xi N N a,x2 (3.33) a,xi N a,x2 0 where notation Nax. = has been used. We further choose a bilinear distribution a,xi dxi over the element for the shear strains ( Yh ) = Yh (I, t) |n = g N (I) Yi (t), Yi = [Y11, Y2i]T , where the nodal shear strains yi are obtained from (3.3), (3.30) and (3.31) as (3.34) 1 Y i = TT- tiJ nIK jK nuWK + iJ niKwj - (jlK niJ + TJ niK) wi + 2 nijnTKOk - 2 niKnJ0j + 2 (nunTK - niKnJ Qi _ I 12 3 4 J 2 3 K 2 3 4 1 (3.35) The notation for strains and transverse displacement can be further simplified as Kh (I,t) | n = Y, B a (I) Ua (t) , Yh (I,t) |n =Y, Ga (I) Ua (t) , a=1 a=1 Wh (I, t) | ne = g M a (I) Ua (t) , Ua = ( ^ ) , (3.36) a=1 where Ba follows from (3.33), Ga from (3.34) and (3.35) and Ma from (3.31). The virtual quantities Kh (I), Yh (I) and Wh (I) are interpolated in the same manner as corresponding quantities in (3.36). One can also introduce more compact notation: ih = [Kh'T, Yh'T] T T T T B a Ua, £ = B a Ua, B a T t B a , Ga When the above interpolations are introduced in the weak form of equilibrium equations (3.4) one gets for an element (e) the following discretized equation (we assume that only load p = p (I, t) is active) G(e) (Wh (I, t) , Oh (I, t) ; Wa, da) =52 Ua^ (Wa, da) ^ (wh (I, t) , Qh (I, t)) , (3.37) a=1 where U ^ da T r(e) ' a A 0, 0 T, ih (I, t) = ih [Oh (I, t) , yh (2 x 2 Gauss integration points are used for the present element) leads to : /n(e) BT (I) a (ih (I, t), ip (I, t), i (I, t), k (I, t)) dQ(e)- , /n(e) MTT (I) p (I, t) dQ(e) (3.38) p = [ph, 0, 0 , ih (I,t) = ih (Oh (I,t) , yh (I,t)), etc. Numerical integration of (3.38) T rr = £ WG[ BaJIc) a(ih (IG,t), ip (IG,t) ,£ (IG,t) ,k (IG,t)) - ) det J (IG) (e) a G=1 MT (Ig) P (Ig, t) G) (3.39) where £G are £ coordinates evaluated at the Gauss point, WG is Gauss point weight, and J is Jacobian matrix of the mapping £ ^ xh. It can be seen that the values of state variables need to be obtained only at the integration points for a particular value of pseudo-time. The component of the element consistent tangent stiffness matrix are K (e) _ d r (e) a ab d Ub 4 £ G=1 'T WgBa (£g) d a (£G,t) d eh (£G,t) dsh 0 and 0n+1 = 0. Backward Euler integration of evolution equations is performed, i.e. en+1 = s« + Yn+1 V n+1, e«+1 = e« + Yn+1/#n+1, Xn+1 = Kn + Yn+1Vn+1 • (3.43) Equations (3.43) and 0n+1 = 0 can be written as the following set of nonlinear equations (with respect to spn+1 = k«+1, en+1 and Yn+1) en+1 - en - Yn+12A (c [s^ - s«^) - | H^e*^ = 0, Cn+1 - Cn - Yn+1(l + = 0, (3.44) V ay J (c (s^ - sn+1) - 2 Hfcmsn+1)T a(C (s^ - en+1) - 2 #^+1) - (1 + 2 = 0, which can be solved iteratively by Newton procedure to get the final values sn+1, en+1 and y«+1. For more effective solution of (3.43) one can write relations (3.16) and (3.17) as, see (3.20) an+1 = C (sn+1 - = aS£ - Yn+1C2A (an+^ + an+1) , Vn + 1 / 2 1 T qn+1 = -S (en+1), en+1 = en + Yn+1-V (an+1 + an+1) A (an+1 + an+1), (3.45) /T.. v ■-V- ^n + 1 2 t • i 2 an+1 = - "Hkinxn+1 = aJ+1 - Yn+17THkin2A (an+1 + an+1) 3 +3 ^-v-' Vn + 1 which leads to (an+1 + an+1) 15 + Yn+1 2CA + 3 HkinA 1 an+l + an+S) , (3.46) v Wn+l(Yn + l) qn+1 = qn+1 (Yn+1)• A closed form expression for the inverse of the matrix in (3.46) can be obtained by using spectral decomposition of c and A; the procedure is very similar to the one at the plane stress situation, see [Ibrahimbegovic and Frey, 1993a], [Simo and Hughes, 1998], [Simo and Kennedy, 1992], [Kleiber and Kowalczyk, 1996], and will not be repeated here, see e.g. [Brank, 1994] for details. Explicit inversion of Wn+1 enables expressing (an+1 + an+1) in terms of single unknown Yn+1. Since Yn+10n+1 = 0, see (3.21), one gets a single nonlinear equation in terms of Yn+1 0n+1((^n+1 + an+1) (Yn+1) , qn+1 (Yn+1)) = 0n+1 (Yn+1) = 0. (3.47) The solution of (3.47) is obtained by Newton iterative procedure j'(k) A (k) = ,(k) Y(k+1) = Y(k) + Ay(k) (3 48) — 0n+1AYn+1 = 9n+1 j Yn+1 = Yn+1 + AYn+1J (3.48) where (k) is iteration counter and 0n+1 = . The final (converged) solution of (3.47) is marked by the bar, i.e. Yn+1. The final values at the end of the pseudo-time increment are also marked with the bar; for example £pn+1 = Kn+1 and £n+1 can be computed from (3.46) and (3.43) by using Yn+1. (k) Remark 3.4. Numerical experiments show that for the present case the function 4>n)1 may be very steep; i.e. very large differences in function value can be obtained for small differences in function argument. However, no difficulties in computation of the numerical examples presented below were observed if the convergence criterion was based on the value of we used < 10-12 for convergence criterion. Remark 3.5. When yield function with a (t), see (3.13), was used in the numerical computations, the equivalent plastic curvature (3.14) of the pseudo-time increment [tn,tn+1] was evaluated with kpn and an. The kpn+1 for the next time step was calculated with numerical integration of equation (3.14) ^n+1 = K + ~Eh- [(AK?1;n+1)2 + (AKp22,n+1)2 + AK^ AK^ + ^7^+1)74] 1 , V 3a y (3.49) where AK11n+1, AKp2 n+1, AK12 n+1 are the first three components of Ae^+1 = Yn+1Vn+1. To complete the elastoplastic implementation issues the consistent tangent matrix d^n+1/den++1 has to be derived for Yn+1 > 0. It is obtained by differentiation of = C-1 cn+1 + epn+1, eqs. (3.43), and consistency condition d0n+1 = 0. After some manipulation one can have i -1 d£(i) den+1, (3.50) dcn+1 — C 1 + 27n+1Hn+1A + -=-—-—-Hn+1vn+1v^+1Hn+1 2 J n+1H kin + 6cn+1 where 4 - _T _ _ 2(1 + 2 <«"+') H-+1 = 15 + -Yn+lHkinA, fn+1 = vn+1Hn+1 vn+1, cn+1 = 2 . • + 3 ^2 — 2Yn+1^" (Cn+1) (3.51) The inverse of H-+ can be obtained in closed form. The matrix in (3.50) can be obtained in closed form as well by using Sherman-Morrison formula, see e.g. [Press et al., 1992]. 3.3.3 Computational issues for viscoplasticity The above discussed viscoplastic plate model allows one to define a unified framework for both stress resultant elastoplasticity and stress resultant viscoplasticity for plates. Namely, the integration procedure for the plate viscoplasticity is essentially the same as for the plate plasticity, except that for 4>n+1 > 0 one looks for Yn+1 _ ^ (0n+1) > 0. Its final value is obtained by iterative solution of nonlinear equation n _ — At Yn+1 + 0n+1 (Yn+1) _ 0 ^ Yn+1. (3.52) The consistent tangent matrix is obtained in the same manner as for plasticity except that one has to replace in its derivation the consistency condition d0n+1 _ 0 by d0n+1 — AtdYn+1 _ 0. The form of the consistent tangent matrix is the same as (3.50) except that Cn+1 is replaced by cn+1 + ^. 3.4 Numerical examples The finite element code for inelastic plate analysis was generated by using symbolic code manipulation program AceGen developed by Korelc [Korelc, 2007b] and implemented into the finite element analysis program AceFEM, see [Korelc, 2007a]. We note that the plate element used in this work is locking-free as shown in [Ibrahimbegovic, 1993] (see results for PQ2 element). 3.4.1 Limit load analysis of a rectangular plate A rectangular plate of elastic-perfectly plastic material under uniformly distributed load is analyzed for two sets of boundary conditions: simple supported (of hard type) and clamped (of hard type). The plate characteristics are: thickness h _ 0.5 cm, length l _ 150 cm, width b _ 100 cm. Material parameters are: Young's modulus E _ 21000 kN/cm2, Poisson's ration v _ 0.3 and yield stress ay _ 40 kN/cm2. Numerical analysis was performed with a coarse mesh of 8 x 8 and with a fine mesh of 60 x 40 elements, see Figure 3.6. We compare our results with those obtained by ABAQUS' [Hobbit et al., 2007] quadrilateral shell element (S4R element) with through-the-thickness stress integration (with 5 integration points) and von Mises yield criterion. Load-displacement Midpoint deflection [cm] Figure 3.2: Load - displacement diagram for simply supported rectangular plate Slika 3.2: Diagram obtežba pomik pri prostoležeci pravokotni plosci curves are presented in Figures 3.2 and 3.3. There is a difference in results of both analyses since the stress resultant formulation does not account for gradual through-the-thickness plastification. However, equal limit load is obtained in both cases. It is interesting to see that the mesh density plays more important role in the accuracy of the limit load computation than the chosen way of definition of elastoplastic constitutive model. Namely, the difference between the coarse and fine mesh in predicting the limit load is around 20 % for the clamped plate, see Figure 3.3. By replacing m0 with am0, we can estimate gradual spreading of plastic zones through the thickness. In Figures 3.4 and 3.5 we present load-displacement curves obtained by using a parameter and coarse mesh. We used a constant value of a across one time increment and therefore the yield criterion is no longer smooth in pseudo-time. To reduce the influence of this effect we used small time increments. In case of simply supported plate (Figure 3.4) the first yield is well predicted, yet the curve in subsequent states is below the ABAQUS' curve. Results for clamped plate are much better since one can hardly distinguish between ABAQUS and stress resultant formulation when using time increment At = 0.0025. 3.4.2 Limit load analysis of a circular plate We analyze a uniformly loaded circular plate of the same elastic-perfectly plastic material as in the previous example. Again we consider simply supported and clamped plates. Midpoint deflection [cm] Figure 3.3: Load - displacement diagram for clamped rectangular plate Slika 3.3: Diagram obtežba pomik pri vpeti pravokotni plosci Midpoint deflection [cm] Figure 3.4: Load - displacement diagram for simply supported rectangular plate; a parameter case Slika 3.4: Diagram obtezba pomik pri prostolezeci pravokotni plosci v primeru uporabe parametra a Midpoint deflection [cm] Figure 3.5: Load - displacement diagram for clamped rectangular plate; a parameter case Slika 3.5: Diagram obtežba pomik pri vpeti pravokotni ploSci v primeru uporabe parametra a The radius of the plate is r = 50 cm and the thickness is h = 0.5 cm. The meshes are shown in Figure 3.6. In Figures 3.7 and 3.8 we plot load-displacement curves. Again we see that the coarse mesh overestimates limit load in the case of clamped plate. Overall correspondence of two formulations is reasonable. In Tables 3.1 and 3.2 we compare our Limit load Yield criterion Reference q = 1.629 ay present present (fine mesh) q = 1.625 ay - [Eurocode 3, 2007] q = 1.500 ay Tresca (analytical solution) [Lubliner, 1990], [Sawcžuk, 1989] q = 1.629 ay Von Mises (analytical solution) [Sawcžuk, 1989] q = 2.000 ay Von Mises (analytical upper bound) [Sawcžuk, 1989] q = 1.500 g ay Von Mises (analytical lower bound) [Sawcžuk, 1989] Table 3.1: Limit load solutions for circular simply supported plate Tabela 3.1: Mejne obtežbe pri prostoležeci krožni krožni plosci results for limit load with analytical solutions found in textbooks on plasticity [Lubliner, 1990], [Sawcžuk, 1989] and Eurocode [Eurocode 3, 2007]. Our result for simply supported plate is in complete agreement with solution based on von Mises yield criterion, see Table 3.1. In the case of clamped plate our result is slightly greater than the von Mises yield based solution, see Table 3.2. Figure 3.6: Meshes used for: (a) rectangular plate - fine, (b) rectangular plate - coarse, (c) circular plate - fine and (d) circular plate - coarse Slika 3.6: Uporabljene mreze: (a) pravokotna plosca - fina, (b) pravokotna plosca -groba, (c) krozna plosca - fina in (d) krozna plosca - groba Midpoint deflection [cm] Figure 3.7: Simply supported circular plate - limit load analysis Slika 3.7: Prostolezeca krozna plosca - analiza mejne obtežbe Midpoint deflection [cm] Figure 3.8: Clamped circular plate - limit load analysis Slika 3.8: Vpeta krožna plosca - analiza mejne obtežbe Limit load Yield criterion Reference q = 3.240 h2 ay q = 3.125 h2 ay q = 2.815 h2 ay q = 3.138 h2 ay present Tresca (analytical solution) Von Mises (analytical solution) present (fine mesh) [Eurocode 3, 2007] [Lubliner, 1990], [Sawcžuk, 1989] [Sawcžuk, 1989] Table 3.2: Limit load solutions for circular clamped plate Tabela 3.2: Mejne obtežbe pri vpeti krožni krožni plosci 3.4.3 Elastoplastic analysis of a skew plate We consider a skew plate of elastic-plastic material with hardening under uniformly distributed load. The plate thickness is h = 0.5 cm, the longer side is a = 150 cm, the shorter one is b =135 cm and the in-between angle is 0 = 45°. All material properties are the same as in the above examples, except for isotropic hardening modulus, which is now Hiso = 0.1E = 2100kN/cm2. The plate is supported along the shorter edges with five equally spaced point supports restraining displacements and allowing both rotations. Mesh is shown in Figure 3.9. Load versus centre displacement diagrams are presented Figure 3.9: Skew plate - (a) fine mesh, (b) coarse mesh Slika 3.9: Romboidna plosca - (a) fina mreža, (b) groba mreža - STRESSRESOWXNT - HUE - STRESSRESOWXNT - COARSE 20 40 60 Midpoint deflection [cm] Figure 3.10: Skew plate - elastoplastic analysis Slika 3.10: Romboidna plosca - elastoplasticna analiza in Figure 3.10. Both curves have similar shapes, yet the curve obtained with the coarse mesh is again above the curve obtained with the fine mesh. We see that the yielding of the plate significantly reduces its stiffness, yet the limit load is never reached because of the isotropic hardening. When using stress resultant plasticity model one can easily track the spreading of plastic žones. In Figure 3.11 we see that the yielding starts in the corners Load = 11.776kN/m2 Load = 13.248kN/m2 Load = 14.72kN/m2 Load = 17.618kN/m2 Load = 18.043140625kN/m2 Figure 3.11: Spreading of plastic žones Slika 3.11: Širjenje plastificiranega obmocja of the shorter diagonal, then it reaches the centre of the plate and spreads in the direction of longer diagonal corners. 3.4.4 Cyclic analysis of a circular plate We consider a clamped circular plate under cyclic loading conditions. The plate is loaded with uniformly distributed load with the amplitude that corresponds to twice of the load at the first yield pmax = 2 (1.5(h)2ay). We examine three hardening cases: (i) isotropic hardening (Hiso = 2100 kN/cm2, = 0 kN/cm2), (ii) kinematic hardening (Hiso = 0 kN/cm2, = 2100 kN/cm2), (iii) combined isotropic and kinematic hardening (Hiso = = 1050 kN/cm2). The remaining material and geometry parameters are the same as those adopted for the limit load analysis, except for plate thickness which is now h = 4 cm. In Figure 3.12 the load-displacement curves for the first two cycles are presented for stress resultant formulation and stress resultant formulation with parameter a. Isotropic hardening enlarges the yield surface which can be seen in Figure 3.12 where the curve appears as a closed loop after the first cycle. After the initial cycle the plate can sustain greater stress resultants and still remain elastic. The shift of the initial curve to the right represents the plastic deformation. Purely kinematic hardening curve is wider and is virtually unchanged from one cycle to another. In this case the siže of the yield surface is unchanged whereas the effect of kinematic hardening changes the position of it. We can look at the combined hardening curve as a combination of the purely isotropic and purely kinematic hardening curves. Isotropic hardening effects prevail and after the first cycle the plate remains in elastic state. 8000 6000 4000 2000 0 o -2000 I—i -4000 -6000 Isotropic hardening 8000 6000 4000 2000 0 0-2000 1—i -4000 -6000 -0.2 0 0.2 0.4 Midpoint deflection [cm] Kinematic hardening - m0 = const. ............ m0 = m0(A) At=0.005 yr yS 8000 6000 4000 2000 0 o-2000 I—i -4000 -6000 -0.6 -0.4 -0.2 0 0.2 0.4 Midpoint deflection [cm] Isotropic and kinematic hardening - m0 = const. ............ mQ = mQ(A) At=0.005 -0.2 0 0.2 0.4 Midpoint deflection [cm] Figure 3.12: Clamped circular plate - cyclic load Slika 3.12: Vpeta krožna plosca - ciklicna obremenitev For further elastoplastic numerical examples we direct the reader to [Dujc and Brank, 2006]. 3.4.5 Elastoviscoplastic analysis of a circular plate In this example we consider a clamped circular plate of elastoviscoplastic material. All material and geometry parameters are the same as in the case of limit load analysis. Three different values of viscosity parameter are chosen, n = 0, n = 1, n = 10, for two sets of loading conditions. In the first set we gradually apply a point load in the center of the plate until it reaches its final value F = 22 kN at time t =1. The second set is displacement driven with a prescribed final value of midpoint deflection w =11 cm. Loading curves Time Figure 3.13: Loading curve for viscoplastic analyses Slika 3.13: ObteZne krivulje za viskoplasticne analize for both loading sets are presented in Figure 3.13. We show the time-deflection curve of the plate under first loading condition in Figure 3.14. The viscosity coefficient n has a significant effect on a nature of inelastic response. The value n = 0 corresponds to plasticity whereas for values n > 0 the inelastic deformations are time dependent. One can note that the strain in elastoviscoplastic material held at constant stress will gradually reach the level of strain in a time independent material. Time response of the plate for the strain driven loading is presented in Figure 3.15. We see a hardening like response in viscoplastic materials (n > 0) but resistance is slowly dropping to the value corresponding to the time independent material. 3.5 Concluding remarks and chapter summary The stress resultant plasticity for plates has been revisited and reformulated in this chapter. The basic ingredients of the constitutive law are the usual additive split of elastic Time Figure 3.14: Time response for force-prescribed viscoplastic circular plate Slika 3.14: Časovni odživ krožne plosce pri predpisani obtežbi Time Figure 3.15: Time response for displacement-prescribed viscoplastic circular plate Slika 3.15: Časovni odživ krožne plosce pri predpisanih pomikih and plastic (viscoplastic) strains £ = £ + £ , the strain energy function including isotropic and kinematic hardening ^ (£e, K) = 2£e,TC£ + S (£) + 1 (2Hkl1)j kTDk, stress resultant yield function 0 (a,q, a) = (a + a)T A (a + a) - - —^J = 0, and the principle of maximum plastic dissipation for the plasticity case Lp (a, q, a, j) = -Dp (a, q, a) + Y0 (a, q, a), which is in the viscoplastic case replaced by Lvp (a, q, a) = -Dvp (a, q, a) + - g (0 (a, q, a)), n where ' 102 if 0 > 0 g (0)^0^ if 0< 0' is the penalized functional. The update of plastic internal variables is carried out by solving only one equation 0n+1((ara+1 + an+1) (Yn+1) , qn+1 (Yn+1)) = 0n+1 (Yn+1) = 0 ^ Yn+1, while in the viscoplastic case this equation is replaced by n — - At Yn+1 + 0n+1 (Yn+1) = 0 ^ Yn+1- By setting n = 0 in the viscoplastic equation we obtain the plastic case, thus both inelastic formulations can be treated within one computational framework. Numerical results of the presented formulation have been compared with the stress formulation (ABAQUS [Hobbit et al., 2007]) as well as with the stress resultant formulation with a parameter that takes into account gradual spreading of through-the-thickness plastification. It has been shown that, regarding the accuracy of the limit load computation, the mesh density plays more important role than the type of elastoplastic formulation. An extension of this work to geometrically nonlinear shells will be addressed in Chapter 4. Chapter 4 Inelastic analysis of metal shells 4.1 Introduction In this chapter we present an approach for an elastoplastic analysis of metal shells by using the finite element method. In our approach we use a constitutive law that is based on stress resultants, thus reducing the number of required numerical integration points in the thickness direction to only one, which makes the analysis computationally much faster with respect to the approach dealing with stress, see e.g. [Brank et al., 1997]. A similar approach was already presented in [Simo and Kennedy, 1992], with the use of a different shell finite element and a different hardening law. In Section 4.2 we first describe the kinematics and the variational formulation of the geometrically exact shell (note that this part is included to ensure the complete description of the derived finite element, while the emphasis of this chapter is on the constitutive relations presented and the return mapping algorithms) followed by the constitutive equations for stress-resultant elastoplasticity for shells. In Section 4.3 we present the spatial finite element discretization and computational algorithms related to shell stress resultant multi-surface plasticity. The results of several numerical simulations are presented in Section 4.4. We close the chapter with concluding remarks and a short summary in Section 4.5. 4.2 Inelastic geometrically exact shell formulation In this section we present inelastic nonlinear shell model that is formulated entirely in stress resultants and can accommodate large displacements and large rotations. 4.2.1 Geometry, kinematics and strains We model a shell as a surface that has an inextensible unit vector (called shell director) attached at its every point. The surface, which will be called midsurface, since it represents the middle-surface of the shell, is embedded in the 3d space. The position vector of a material point in the initial (undeformed, reference) stress-free shell configuration is then defined by x (e1 ,e2 ,e) = ^o (e\e2) + e t (e\e2), (e\e2) g a, e g f . (4.1) Here, e1 and e2 are convective curvilinear coordinates that parametrize the midsurface; A is the domain of that parametrization; T, ||T || = 1, is the shell director that coincides with the normal vector to the midsurface; and e is through-the-thickness convective coordinate defined in the domain F = [—h/2,h/2], where h represents initial shell thickness, here assumed to be constant. In what follows, we always determine the components of the above vectors with respect to the fixed orthonormal basis ei = i = 1, 2, 3, in the 3d space, i.e. X = Xiei, = , T = Tiei. We further define the shell director as T = A0e3, where A0 is a given (initial) rotation tensor, A-1 = A^, det A0 = 1. If one introduces at a point of midsurface an orthonormal basis e = ee as es = T, e ± es, Hej = 1, e = e3 x eb (4.2) the rotation tensor A0 at that point can be represented as A0 = [e1, e2, T]. It is further assumed that the position vector to the material point in the deformed configuration is given by x (e1,e2,e) = (e 1,e2) + u(e1,e2)] + et (e1,e2), Htii = 1. (4.3) vC«1,«2) In (4.3) u is a displacement vector of a midsurface point, and t is new position of shell director at the deformed configuration. We will define t as the following sequence of two rotations t = A0Ae3. The components of newly defined vectors in (4.3) are also determined with respect to the fixed orthonormal basis ei, i.e. x = xiei, ^ = ^iei, u = uiei, and t = tiei. The rotation tensor A is viewed in this work as a function of a constrained rotation vector m, i.e. A = A (m), see e.g. [Brank and Ibrahimbegovic, 2001] and [Brank et al., 1997] for details. Since the rotation around the shell director (i.e. drilling rotation) plays no role in the present theory, the constrained rotation vector has only two components with respect to the basis ei, i.e. m = a = 1, 2. By using the Rodrigues formula for the representation of A (m), one ends up with the following expression for t = t (m), see [Brank and Ibrahimbegovic, 2001] ( sin $ \ t = A0 ( cos$e3 + m x e3j , $ = ||m|| . (4.4) The vectors of the convected basis Gi at the initial configuration are related to the position vector X and to the convected coordinates £a,£ as Ga _ de _ de0+edp G3 _ ^ _ • (5) Similarly, the vectors of the convected basis gi at the deformed configuration are d x d p .81 d x ,A ga _ de _ de+ede' g3 _ de _ *• (46) The corresponding dual base vectors Gi and gi are defined through the relationships Gi • Gj _ 8j and gi • gj _ 8j, where 8j is a Kronecker symbol. Note that G3 _ G3. The identity tensor of the shell reference configuration (or the shell metric tensor) is G _ Gi 0 GJ — Gij Gi 0 GJ, where Gj — G i • G j. The differential volume element is given as dV _ VGdCdC1dC2, where \[G _ G3 • (G1 x G2). The base vectors at the reference midsurface and at the deformed midsurface are obtained by setting £ _ 0 in (4.5) and (4.6), respectively, i.e. Ai _ Gi |^=0, and ai _ gi |^=0. For the reference configuration we have A _ dp0' A3 _ T• (4.7) The corresponding dual base vectors of A% and ai are defined as Ai • Aj _ 5j and ai • aj _ 8j, respectively. Note that A3 _ A3. The identity (or metric) tensor of the shell reference midsurface is A _ Aa 0 Ap _ AapAa 0 Ap _ Aap (Aa)Y (A/3)s eY 0 _ ea 0 e^, where Aa$ _ Aa • Ap and Aap _ Aa • Ap. The differential surface element is given as dA _ \fAdC1 2, where \[A _ ||A1 x A2||. We can further define a tensor, called shifter, which transforms the base vectors of the midsurface to the base vectors of the shell body. The shifter from the shell reference configuration, denoted as Z, shifts Ai and A% (defined at a midsurface point) to Gi and Gi, respectively, i.e. Gi _ ZAi and Gi _ Z-TAi. In what follows, we assume that a shell is suficently thin that Z ~ I. Having defined the base vectors, we can proceed with the expression for the deformation gradient d x dx dX dei dX 1 gi 0 Gi• (4.8) l_d£ iJ In (4.8) we used notation £3 _ £. By knowing F, we can obtain the components of the Green-Lagrange strain tensor with respect to the convected basis A% E _ 2 (fTF - G) _2 [(^ 0 gt)(gj 0 Gj) - Gi • Gj (g1 0 Gj)] (4.9) 1 (gi • gj - G i • G j) Gi 0 Gj _ Eij Gi 0 Gj j^Ejj Ai 0 Aj • zxii By evaluation of the dot products in (4.9) one can get the strains Ej, which are varying quadratically with respect to the £ coordinate Ej = £j + + (£ )2 nij - (4.10) Note, that £33 = k33 = na3 = n3a = n33 = 0. In this work we will truncate the strains Eaf after the linear term, and the transverse shear strains Ea3 = E3a after the constant term, i.e. Eaf ^ £af + £Kaf, (4 n) Ea3 ^ £a3" The Green-Lagrange tensor that we will work with will have the following components in basis Ai E « Eij Ai 0 Aj = £af Aa 0 Af + £Kaf Aa 0 Af + £a3 (Aa 0 T + T 0 Aa). (4.12) Expressions for £af, Kaf and 2£a3 in (4.11) are the classical expressions for the shell membrane, the shell bending and the shell transverse shear strains, respectively. Their explicit forms follow from using gi and Gi in (4.9) £af = 1 (V,a ■ - ^0,a " 0, 6, ^ 0, = 0. (4.42) In addition to (4.42) we have the condition = 0 if j, > 0 (the consistency condition). It guarantees the admissibility of the subsequent state in the case of change of state variables. 4.3 Finite element formulation 4.3.1 Space-domain discretization Let the initial shell midsurface A be discretized by nel nonoverlaping elements with nen nodes such that A & (jn=1 Ae = Ah. Over the element domain Ae the initial shell configuration (the midsurface and the shell director) are interpolated as nen nen (e 1,e2) = £ Na (e\e2) ^ Th (e\e2) = £ n (e ^2) t„, (4.43) a=1 a=1 where (o)a are the corresponding nodal values. In this work we choose nen = 4 and the bi-linear shape functions Na (e1 , e2), defined over the square domain Ae = [—1,1] x [— 1,1]. Note, that Ta is chosen to coincide with the normal vector to a given shell midsurface at that nodal point. However, due to the bi-linear interpolation (4.43) Th is only approximately perpendicular to the base vectors A^ = d^/de". The interpolation of the shell deformed configuration t is performed in a similar fashion as T nen nen ■h (e1, e2) = E Na (e1, e2) «a, ^ = ^h+uh, th (e1, e2) = E Na (e\ e2) ta wo. u a=1 a=1 (4.44) The virtual quantities S^ and St are interpolated in the same manner as ^ and t nen nen (e\e2) = E Na (e1,e2) s^, sth (e \e2) = E Na (e\e2) sta. (4.45) a=1 a=1 Derivations of t, S^ and St with respect to ea coordinates are obtained trivially. To avoid the transverse shear locking, the assumed natural strain (ANS) concept is chosen. The transverse shear strains are evaluated, using (4.13) and the interpolations (4.43) and (4.44), only at element edge mid-points A, B, C and D, where ^0 = 2 (^0J + ) and I = A, B, C, D, J = 1, 2, 3, 4 and K = 2, 3, 4,1. The transverse shear strain field over the 4-node shell element is given by the interpolation suggested by [Bathe and Dvorkin, 1985] 713 = 1 (1 - e2) yA + 2 (1+e2) YC3, 723 = 1 (1 - e1) y23 + 1 (1 + e1) ^ (4.46) The transformation to 713 and 723 is given according to (4.15). 4.3.2 Computational issues for plasticity The solution of the weak form of equilibrium equations (4.26), discretized in space is searched for at discrete pseudo-time points 0 < t1 < ... tn < tn+1... < T. The value of (o) at the pseudo-time instant tn+1 is denoted as (o)n+1. In what follows a pseudo-time increment At = tn+1 — tn will be considered. As a result of space discretization, the evolution equations (4.41) become ordinary differential equations in pseudo-time that are related to each finite element integration point. This enables introduction of operator split method, see e.g. [Ibrahimbegovic, 2006]. This method consists of two sequential solution procedures; one solution procedure (called global) searches for the values of nodal dis-placemens/rotations at pseudo-time instant tn+1 (at frozen values of internal variables), while the other solution procedure searches for the values of internal variables at integration points at tn+1 (while keeping frozen nodal displacements/rotations data). The later will be considered in this section. For numerical integration of evolution equations (4.41) the backward Euler integration scheme will be used. At a pseudo-time increment At = tn+1 — tn, the pseudo-time integration problem at an integration point located at xh (£G) G can be stated as: By knowing the values of the internal variables at tn, i.e. e«, £«, xn, find by integrating (4.41) such values of the internal variables at tn+1, i.e. £«+^£«+1^, Kra+1, which will satisfy the yield criterion. The best guess for the strains at the end of the pseudo-time increment, £«+, is given data. Here (i) is an iteration counter of the global Newton-Raphson solution procedure that searches for nodal displacements/rotations at time tn+1. Prior the integration of evolution equations the following test is performed: (i) Assume that the pseudo-time step from tn to tn+1 remains elastic and evaluate the trial (test) values of strain-like and stress-like internal variables / \ / \ „(i) p,trial _trial /-i ^n+l = C £n+1 £n+1 V trial _ qn+1 = ^/,trial ^ra+l 2 03? = — 2 #*«*«+?, (4.47) (ii) If both yield functions, evaluated with the trial values, ^r,ra+1 = (^«r+ai, 9«r+'i, a«r+a1l) — 0, then, see (4.42), 7M,n+1 = YMAt = 0. The final values at the end ofthe pseudo-time increment (marked with the bar) equal the trial values, i.e. £«+ = £«'+7al, £«+1 = £«+7al and Kn+1 = x«+a_l. The pseudo-time step is indeed elastic. (iii) In the case that one or both yield functions for the trial values are violated, then 7M,n+1 > 0 and 0M,n+1 = 0 for at least one ^ = 1, 2. The backward Euler integration of evolution equations needs to be performed, i.e. ^=1 ^=1 ^=1 (4.48) to get final values of internal variables at the end of the pseudo-time step. Since equations (4.48) contain five unknowns, £«T1,£«+1, xn+1,71,n+1 and Y2,n+1, one needs to solve those equations together with 01,n+1 = 0 and 02,n+1 = 0. Remark 4.1. When dealing with two-surface plasticity we have in general three different sets of equations, where 01,n+1 = 0 and 02,n+1 = 0 is only one of them. The remaining two occur when only one yield surface is active and then instead of solving the above equations we solve 01,n+1 = 0 and Y2,n+1 = 0 or Y1,n+1 = 0 and 02>n+1 = 0. The procedures that produce the right set of equations and their solutions will be further addressed in Section 4.3.2, while for brevity we will here only consider the case when both functions are active. The rest of this section will be related to the solution of eqs. (4.48) constrained by 01,n+1 = 0 and 02,n+1 = 0. There are two ways of handeling this problem, where both of them have their benefits and drawbacks. One option is to reduce the number of unknowns by: (i) expressing stress resultants with eqs. (4.48), and (ii) use those stress resultants in 01,n+1 = 0 and 02,n+1 = 0, while the second option is to solve all the equations simultaneously. We first consider the reduction of the unknowns. Let us consider the case when both yield surfaces are active. By using (4.48), one can express stress resultants at tn+1 as "n+1 .(i) p,trial\ _ trial .. — ^ .. 1 — & .. — C I En+1 Sn+1 i _(i) _p Acp in + 1-£n-A&n+1 n+1 ,= 1 At ip . , n + 1 qn+1 = an+1 ( ih + Y1 Y^,n+1pn+1] , ,=1 2 . i 2 — 3 Hkin Kn+1 = an+1 — 3 Yn+1v ,,n+1- ,= 1 (4.49) Since vn,n+1 = 2A, (an+1 + an+1), one can conclude from (4.49) that ("n+1 + an+1) 24 Is + Yv,n+1y2CA, + 3 HkinA, n=1 1 trial trial "n+1 + an+1 qn+1 = qn+1 ^^n + ^ Yn,n+1Pn+1^ (4.50) (4.51) where the matrix product C An is CA, -1C nP n 5 sign(n) cbp J_Cbp 2\/3n0m0 0 sign(n) c n p 0 2\/3n0m0 0 0 C 12 2 For the linear isotropic hardening case one can from (4.38), (4.41) and (4.51) easily obtain that Pn+1 2(Kh£n + ^y) ^ - 2 En=1 Yn,n+1Kh (4.52) 0 It has been proven that the inverse in (4.50) can be obtained explicitly. Since the shear part of (^n+1 + an+1) is clearly uncoupled from the membrane and bending parts, one can use (4.50) to get the following expression for the shear stress resultants 1 (qfn+1 + a n+1 1 + (Yn+1+Yn+o qr (ks + § h^) trial | q,trial \ qn+1 + an+1 J 9o v" ' 3 The fact that P and C have the same characteristic subspaces is further used, i.e. P = QAP QT , where QT = Q-1 and C = QAC Qt , Q = 722 110 — 1 1 0 00 (4.53) ■ 1 2 0 0 j= 0 3 2 0 , Ac = 0 0 3 ^ 0 1 — - 0 T+- 1+V 00 0 0 E 2(1+-) J By defining n+1. Now, the yield functions can be rewritten by using (4.58) and (4.51) as 0,(71,n+1, 72,n+1) = (^n+1 + an+1)T A, (an+1 + an+1) -11- qn+1 (cn + S,=1 7,,n+1^n+1^ 0, (4.59) V = 1, 2, and further as 0,(71,n+1,72,n+1) trial trial T T trial trial an+1 + an+J W n+1A, W n+1 (an+1 + an+1 ) trial T T trial trial -11- qn+1 cn+e ,=1 7,,n+1^n+1 ) 0, (4.60) V 1, 2. 0 1 2 2 One finally gets only two equations for 71>n+1 and 72>n+1, RP (Y1,n+1, Y2,n+1) = 01 (Y1,n+1, Y2,n+1) 02 (Y1,n+1, Y2,n+1) 0, (4.61) that are highly nonlinear. Once (4.61) are solved for converged solutions Y1,n+1 and Y2,n+1, the strains (4.48) and stress-resultants (4.47) at the end of the time increment can be computed. In the second case we simultaneously solve (4.48) along with 01,n+1 = 0 and 02,n+1 = 0 for five unknowns £Pn+1,i1n+1, Kn+1, Y1,n+1 and Y2,n+1. By considering the first and the last equation in (4.48) we can eliminate Kn+1 since Kn+1 'n+1' In this situation, when all the internal variables are considered as the unknowns, we can obtain the values of stress and stress like variables & = C f6(i) 6p &n+1 = C ^6n+1 - 6n+1 qn+1 = qn+1 (^n+O ' 2 an+1 = — 3 Hkin K n+1, and with them express the yield functions 0^,n+1 = (&n+1 + an+1)T n+1 + an+1) — ( 1 — ) V ) The final system of equations that needs to be solved is then 1, 2. (4.62) (4.63) (4.64) (4.65) RP (6n+1,in +1Y1,n+1,Y2,n+1) -6n+1 + £n + E ^=1 Y^,n+1v^,n+1 — cn+1 + cn + E n=1 YIJ.,n+1Pn+1 01,n+1 02,n+1 (4.66) Note that in this case all the expressions are rather simple. The drawback is that we need to solve a system of 11 equations (8 for plastic strains, 1 for isotropic hardening and 2 yield functions), which means we need to invert a matrix of size 11 x 11. Solution algorithm In the previous section we have presented the equations that one must solve in order to obtain the updated values of internal variables, i.e. (4.61) or (4.66). Both ways of handling the evolution of plastic variables and the yield functions produce completely the same results and for simlicity we will only consider the first case, i.e. (4.61). As already mentioned in the previus section we have chosen the equations (4.61) by assuming that both yield surfaces p = 1, 2 are active 01 (Y1,n+1, Y2,n+1) 02 (Y1,n+1, Y2,n+1) R12 (Y1,n+1,Y2,n+1) 0. (4.67) 0 The two other possibilities occur when only one yield surface is active. Let us first consider the option when only ^ = 1 is active. In this case we replace equations (4.67) with R (Yl,n+l, Y2,n+l) = [ 0l(7l,n+1' 72,n+1M = 0. (4.68) L Y2,«+1 The last option is when only ^ = 2 is active and the equations that need to be solved are then Rp (Yl,«+1, Y2,n+l) 0. (4.69) Yl,n+1 02(Yl,«+1, Y2,n+l) Figure 4.2 depicts three different situations for the trial values of stress. For the trial \<$S-1 = 0 trial, 12 _tnal,1 \&n+1 01riaL = 0 Figure 4.2: A sketch of a two-surface yield function and the closest point projection Slika 4.2: Skica projekcije testnih vrednosti napetosti na funkcijo tecenja doloceno z dvema ploskvama value ^f1 we have 01«+1 > 0 and 02,«+1 < 0. Here we know in advance that the solution of equations (4.68) gives us the closest point projection and the correct update of the internal variables. Similarly we know for ct«+11,2 that we need to solve equations (4.69) since 01«+1 < 0 and 02«+1 > 0. But in the case of CT«+a_,12, when both trial values of yield functions are violated, we can not be sure in advance which set of equations is the right one. In the last case, when we are not sure which equations to use, we have two alternative strategies: • Procedure 1. In this method we look for solutions for all posible sets of equations and the admissibility of the solutions is checked by testing weather the Kuhn-Tucker conditions hold. • Procedure 2. In this method we can change the set of equations during the iterative process by enforcing the admissibility constrain Y^,n+1 > 0 for all active surfaces. A brief descriptions of the algorithms associated with Procedures 1 and 2 are given below. Procedure 1: In this procedure we determine the solutions of three different sets of equations —a,n+1 = [Y 1,n+1,a, Y2,n+1,JT , K (—a,n+1) = 0 for a = 1 2 12, (4.70) where (◦) denotes the converged values of variables. Each solution is obtained with the following iterative procedure. 1. First we define the iteration counter k k = 0, and set the initial values of plastic multipliers to zero y as+ = [0, 0]t . 2. Then we start an iterative loop in which we first determine the current values of Rpa RP,(k) = RP (—(k) \ Ra = Ra I I a,n+1 J , and chech for convergence \\RPa(k)\\ < tol.. 3. If the convergence test is satisfied we are happy with the solution — = (k) Ya,n+1 — a,n+1, and we exit the iteration loop. Otherwise we move on to 4. i i r* —RP d 4. Here we compute the current value of matrix K a = k Rp,(k)=k Rp (— and compute the update for plastic multipliers A—ah=-(k Rp,(k))-1Ra(k). 5. Then we set the iteration counter to k = k + 1, update the values of the plastic multipliers (fc) = (fc-l) . a (fc-l) Y a,«+1 = Y a,«+1 + AY a,«+1, and go to 2. The admissibility of each solution is checked with the Kuhn-Tucker's loading/unloading conditions: Yl,«+1,a > 0 & Y2,«+1,a > 0 & 01 (Ya,«+1) — 0 & 02 (Ya,«+l) — 0, a =1, 2,12. (4.71) N-V-' ? We denote the solution that satisfies all the conditions with a Y«+i = 7s,«+i- (4.72) Note that this procedure is robust and it always provides the solution but it is quite computationally expensive since three independent iteration procedures are performed. Procedure 2: The below algorithm is a variation of a general multi-surface closest point projection iteration procedure, that is presented in [Simo and Kennedy, 1992] and [Simo and Hughes, 1998]. The difference is that in our work we deal with at most two active yield surfaces, whereas the original algorithm allows for an arbitrary number of active yield surfaces. The solution is obtained with the following iteration procedure. 1. First we define the iteration counter k = 0, set a = 12, to determine the starting set of equations and set the initial values of plastic multipliers to zero Y«+i = [0,0]T . 2. Then we start an iterative loop in which we first determine the current values of Rp Rp,(k) = Ra (y«+I) , and chech for convergence ||Rp,(k)|| < tol.. 3. If the convergence test is satisfied we are happy with the solution — = (k) Y n+1 = Y n+1, and we exit the iteration loop. Otherwise we move on to 4. 4. Here we compute the current value of matrix KRP = dRP k Rp ,(k)=K Rp (7n+1 and compute the update for plastic multipliers A7 n+ = —(K Rp,(k))-1 RP,(k). 5. Then we set the iteration counter to k = k + 1, and if a = 1 or a = 2 update the values of the plastic multipliers (k) = (k-1) + a (k-1) Y n+1 = Y n+1 + AY n+1 , and go to 2. If a = 12 go to 6. 6. Here we compute the test values of plastic multipliers (k),test = (k-1) + a (k-1) Tn+1 = Tn+1 + AT n+1 , and check if those values are admissible ?? (k),test -1 n (k),test -1 n Y1,,i+1 > Y2,n)+1 > 0- If both tests are satisfied we update the values of the plastic multipliers (k) = (k-1) + a (k-1) Tn+1 = Tn+1 + AT n+1 , and go to 2. If one of the tests fails, the values of plastic multipliers remain unchanged (k) = (k-1) Tn+1 = Tn+1 , we change the parameter a according to If Y(kn+1 < 0 then a = 2, If Y2k2+1 < 0 then a =1, and go to 2. This procedure is computationally much cheaper and is also rather robust. Theoretically, if the yield surface is convex, then the above algorithm is unconditionally convergent. Practically, it may occur that the above procedure does not produce the desired solution and in that case a step size adjustment is needed, see e.g. [Simo and Hughes, 1998] and references therein. Consistent tanget modulus In order to ensure the quadratic rate of convergence of the global iterative procedure we have to consistently linearize the global system of equations. This requires to compute the implicit dependencies among the state variables and the strain vector in order to obtain the consistent tangent modulus d*n+1/d£n+1. For that purpuse we assume the following functional dependencies * (£„+1) = C (£„+i - <+1 (Yn+1 (£n+l))) • (4.73) The challenging part is to obtain the derivatives of plastic strain with respect to the total strain. By applying the chain rule we have d£n+1 _ d£n+1 n+1 d£n+1 ^Tn+l d£n+1 (4.74) and the only unknown derivative here is . The implicit dependencies are obtained from the consistency condition 0M,n+1 _ 0. We assume that the yield surfaces (the plastic equations) are the functions of the total strain and the plastic multipliers 0M,n+1 _ (£n+1> Yn+1 fcn+0) ^ RP (£n+1, Yn+1 fcn+0) • (4.75) By the chain rule derivation of the above equation we obtain Rp _ dRP d£n+1 + dRP dYn+1 d£n+1 M'n+1 d£n+1 dt dYn+1 d£n+1 dt rRP dR + KRdY^ _ 0, (4.76) d£ n+1 ^ ' =0 and the derivative that we are looking for is then dYn+1 _ -(K Rp )-1 dR , (4.77) d£n+1 d£n+1 where all the derivatives on the right hand side of (4.77) can be obtained by considering functional dependencies. 4.4 Numerical examples In this section we present some numerical examples, computed by the above derived 4-node element. The computer implementation of the element is in complete accordance with the above derivation. The local Cartesian frames are introduced at the element integration points. The computer code was generated by using symbolic code manipulation program AceGen developed by Korelc [Korelc, 2007b], [Korelc, 1997]. The element codes were introduced into the finite element analysis program AceFEM, see [Korelc, 2007a]. 4.4.1 Pinched cylinder with isotropic hardening We consider a short cylinder bounded by two rigid diaphragms at its ends. The cylinder is loaded with two concentrated forces at the middle section. Due to symmetry, only one octant of the cylinder is modeled. The geometry, loading, boundary conditions and Figure 4.3: Geometry, loading and boundary conditions for pinched cylinder Slika 4.3: Geometrija, obtežba in robni pogoji pri cilindru the finite element mesh of the octant are presented in Figure 4.3, where a = 300 cm, bc1 denotes the edge with the rigid diaphragm while bc2, bc3 and bc4 denote the edges with the symmetry boundary conditions. The thickness of the cylinder wall is 3 cm and the material properties are: Young's modulus E = 3000 kN/cm2, Poisson's ration v = 0.3 and yield stress ay = 24.3 kN/cm2. The plastic behavior is characterized by linear isotropic hardening response with hardening modulus Kh = 300 kN/cm2. The load versus displacement curves of our simulation along with the curves obtained in [Simo and Kennedy, 1992] and [Brank et al., 1997] are presented in Figure 4.4. One can see that all the formulations have similar responses at the low levels of loading while the differences grow with the increase of loading. Note that in this example Procedure 2 presented in Section 4.3.2 did not always produce the solution for internal variables. We believe this was caused by the relatively large steps in loading that occur due to local buckling of the cylinder. The robust Procedure 1 had no problems finding the solution. Figure 4.5 depicts the initial and the final deformed configuration of the cylinder. Z -a rt o 1—I 0 50 100 150 200 250 300 Displacement [cm] Figure 4.4: Load versus displacement curves for pinched cylinder Slika 4.4: Diagram obtežba - pomik pri cilindru 25 000 20 000 15 000 10 000 5000 0 Figure 4.5: Initial and deformed configuration for pinched cylinder Slika 4.5: Zacetna in deformirana lega cilindra 4.4.2 Half of a sphere We consider half of a sphere loaded with two inward and two outward forces. Due to symmetry, only one quadrant is modeled. The geometry, loading, boundary conditions Figure 4.6: Geometry, loading and boundary conditions of half of a sphere Slika 4.6: Geometrija, obtežba in robni pogoji pri polovici sfere and the finite element mesh of the quadrant are presented in Figure 4.6, bc1 and bc2 denote the edges with the symmetry boundary conditions and the edge bc3 is free. The radius of the sphere is 10 cm and the thickness is 0.5 cm. The material properties are: Young's modulus E =10 kN/cm2, Poisson's ration v = 0.2, yield stress ay = 0.2 kN/cm2 and the linear hardening modulus Kh = 3 kN/cm2. The load versus displacement curves of our simulation along with the curves obtained in [Simo and Kennedy, 1992] and [Basar and Itskov, 1999] are presented in Figure 4.7. We note that [Basar and Itskov, 1999] used stress-based plasticity for shells. In this example both procedures described in Section 4.3.2 did not exhibit any problems. Figure 4.8 depicts the initial configuration, the configuration under the maximum load and the final unloaded configuration, where irreversible nature of plastic deformation is evident. 4.4.3 Limit load analysis of a rectangular plate A clamped (of hard type) rectangular plate of elastic-perfectly plastic material under uniformly distributed load in the Z direction is considered in this example. We consider Z Figure 4.7: Load versus displacement curves for half of a sphere Slika 4.7: Diagram obtežba - pomik pri polovici sfere Figure 4.8: Initial, deformed fully loaded and deformed unloaded configuration for half of a sphere Slika 4.8: Zacetna, deformirana obremenjena in deformirana neobremenjena konfiguracija polovice sfere two cases, a thin plate with thickness 0.5 cm and a thick plate with thickness 5 cm. The other geometry and material parameters are the same for both cases: length a = 150 cm, width b = 100 cm, Young's modulus E = 21000 kN/cm2, Poisson's ration v = 0.3 and yield stress ay = 40 kN/cm2. The geometry and the finite element mesh of the plate Figure 4.9: Geometry for rectangular plate Slika 4.9: Geometrija pri pravokotni plosci are presented in Figure 4.9. We compare the results of the geometrically nonlinear shell formulation with those obtained by the geometrically linear plate formulation presented in Chapter 3. Load-displacement curves for the thin plate are presented in Figure 4.10. Displacement [cm] Figure 4.10: Load versus displacement curves for the thin plate Slika 4.10: Krivulje obtežba pomik pri tanki plosci There is a big difference in results of both formulations, since the plate is so thin, that the assumptions of the geometrically linear plate formulation are not valid. The influence of geometrical nonlinearity is negligable only when load is small. With the increase of load, the displacements also grow and the axial forces start to play a big role, thus giving a much stiffer response. We can see that limit load of the plate formulation is around 0.097 kN/cm2, while the shell formulation at load level 2 kN/cm2 is still in elastic regime. In Figure 4.11 we plot the load versus displacement curves for the thick plate. We can see Displacement [cm] Figure 4.11: Load versus displacement curves of the thick plate Slika 4.11: Krivulje obtežba pomik pri debeli plosci that in this case we have a much better agreement in results. The limit load computed with the plate formulation is 0.95 kN/cm2, while in the shell formulation axial stiffening occures, thus giving a greater resistance of the plate. 4.5 Concluding remarks and chapter summary An inelastic geometrically exact shell finite element formulated entirely in terms of stress resultants has been presented in this chapter. The basic ingredients of the constitutive law are the usual additive split of elastic and plastic strains e = ee + ep the strain energy function including isotropic and kinematic hardening ^ (ee, e1, *) = 2ee'TCee + S ) + \ktDk, stress resultant yield function defined with two yield surfaces (see Figure 4.1) (ct, q, a) (ct + a)T A (ct + a) - ( 1 V = 1, 2, the plastic dissipation = cttep + qC'1 + aTk > 0, 2 and the principle of maximum plastic dissipation 2 Lp (a,q, a, 71,72) = -Dp (a,q, a) + > 7^ (a,q, a) m=1 which here considers the multi-surface nature of the yield function. In Section 4.3.2 we first present two options of computing the updates of plastic variables. In the first approach we only solve the equations related to yield function RP (Y1,n+1, Y2,n+1) 01 (71,n+1,72,n+1) 02 (71,n+1,72,n+1) 0, and with the newly computed values of plastic multipliers determine the updates of internal variables. In the second approach we determine the values of plastic variables directly, since evolution equations are included in the system of nonlinear equations RP (en+1,^n+171,n+1,72,n+1) -<+1 + ePn + E M=1 1^,n+1v M,n+1 tl tn+1 tn+1 + + E u=1 01,n+1 02,n+1 When dealing with multi-surface plasticity we are in general not a priori sure which of the yield surfaces is active and in the case of two surfaces, we differ between three sets of equations R12 01,n+1 02,n+1 0, R1 01,n+1 Y2,n+1 0, Rp2 Y1,n+1 02,n+1 0, which cover all the possible outcomes, i.e. R1p2 both yield surfaces are active, R1 only the first yield surface is active and R2 only the second yield surface is active. Two procedures for obtaining the proper solution have been presented. In Section 4.4 several numerical examples have been presented, which show a very satisfying performance of the presented approach compared to results from literature. 0 Chapter 5 Illustration of embedded discontinuity concept for failure analysis 5.1 Introduction In this chapter we illustrate the embedded discontinuity concept for failure analysis of solids. The failure of a 1D solid contains all the significant features that are present in the failure of a more complex 2D and 3D solids yet it is simple enough to clearly illustrate the principle of the embedded discontinuity finite element method. F u A= 0 x= L 1 Aui 1 i Au 2 ! ! 2 F 0 " u(x) u 1 " 1 w __Jj 1 ^ L x Figure 5.1: Tension test of an idealized bar Slika 5.1: Natezni test idealizirane palice We consider a tension test of an idealized bar presented in Figure 5.1. The response of the bar is linear elastic up to the point when the ultimate material stress is reached. Then the resistance of the bar reduces linearly with the increasing imposed displacement, see bottom right part of Figure 5.1. The bottom left part of Figure 5.1 depicts one of the key features of the failure of the component, i.e. the resistance of the bar is reduced due to the localized failure that occurs in the small neighborhood of the bar's weakest point (we choose the midpoint of the bar). In the sketch we presented the failure by the necking effect that is common in the tension tests of metal bars. Note that except for the midpoint where irreversible (plastic) strains appeared the rest of the bar is elastic, thus the softening response of the bar is governed only by the behavior of the bar's midpoint. In the top right part of Figure 5.1 we present the distribution of displacements along the length of the bar for the bar that is in softening regime prior to complete failure. One can interpret the jump in displacements u as the localized plastic strain at the position of the discontinuity. The chapter is organized as follows: in Section 5.2 we present the kinematics, equilibrium equations and constitutive relations in a 1D solid finite element with embedded discontinuity. In Section 5.3 we present the computational procedure for failure analysis with the embedded discontinuity finite element method. Two numerical simulations are presented in Section 5.4. Finally in Section 5.5 concluding remarks and a short summary close the chapter. 5.2 1D finite element with embedded discontinuity In this section we will present the formulation that is capable of modeling the failure of the idealized bar presented in the previous section. Kinematics As a starting point we consider a standard isoparametric 1D solid finite element of length Le whose geometry is defined by xh(0 = 5] N*(0xa, xh e ne = [X1,X2], (5.1) a=1 where superscript h is used to denote the discretely approximated quantities, £ e [-1,1] is the element's local coordinate, xa is the global coordinate of the node a, is the element's domain and Na (£) = 1(1+ £a£) , aa -1 1 • (5.2) The element's displacement field is interpolated as uFem (0 = £ Na(£K (5.3) a=1 where da is the displacement of the node a and subscript FEM is used to denote the displacement interpolation according to standard finite element method. Figure 5.2 depicts Le x\,d\ X2,d2 0 Ni N2 d2 "fem Figure 5.2: Standard isoparametric 1d solid finite element Slika 5.2: Standardni izoparametricni koncni element za palico 1 the finite element's geometry along with the graphic interpretation of shape functions and the displacement field. Clearly this finite element is not able to model the jump in displacement. In order to capture the discontinuity in the displacements we consider the same finite element as described above but now enriched with one additional kinematic parameter a, see Figure 5.3. The enriched displacement field can be then written as 2 uh(£) = E N„(£K + MoCe^, (5.4) >a=1 ^ j enrichment uFem X1,d1 d1 Xd ,A X2,d2 1 O !-►X — ?= 1 /2 o =0 ^ A -1/2 d2 M a uh Figure 5.3: Isoparametric 1d solid finite element with embedded discontinuity Slika 5.3: Izoparametricni koncni element za palico z vgrajeno diskontinuiteto where Ma (C) = (C)) - N2®, Hxd (xh(C)) = 1 for xh > xd (5.5) 0 otherwise and xd is the position of the discontinuity within the element (we choose the element's midpoint). The total strain field can be then computed as the space derivative of this displacement interpolation leading to 2 e = = B«(CK + Ga(C)a, (5.6) dx a=1 where Bi(C) = --, B2(C) = -Je, Ga = r*((C) +(C). (5.7) Ga Ga In above expression we have used the following rule when deriving the Heaviside function dH Xd dx = sXd(C) to for xh(C) = xd 0 otherwise (5.8) We further divide the strain field into a regular part e and singular part e e = e + e = ^ £a(C)"« + Ga« + Ga 0 is the Lagrange multiplier. By using (5.28) and (5.24), we get from (5.29) above the following evolution equations: =p = dL .L. d< = — = -a + y~QI = 0=^ a = j, (5.30) =p = dL = — d< = — . . ■ar = ^ = l = j (531) along with the Kuhn-Tucker loading/unloading conditions and the consistency condition j > 0, < < 0, =0, j<<=0. (5.32) 5.3 Computational procedure for failure analysis In this section we present a procedure for solving the set of global (mesh related) equations and the set of local (element related) equations generated by using the finite element with embedded discontinuity presented in Section 5.2. The solution ought to be computed at discrete pseudo-time values 0,t1,t2,..., Tn-1, Tn, Tn+1,... ,T by means of the incremental-iterative scheme. We will consider the solution in a typical pseudo-time incremental step from Tn to Tn+1. We assume that all the variables related to an element e are given at Tn, i.e. given: = [d^d^, atf, t?. (5.33) We will then iterate in the pseudo-time step in order to compute the converged values of the variables at Tn+1, i.e. find: it = ^gr > 0^1, (5.34) The computation of solution (5.34) is split into two phases: (a) In the global (mesh related) phase we compute the current iterative values of nodal displacements at Tn+1, while keeping the other variables fixed global phase: d^? = dgf-^ + AdJ^-^, (5.35) where (i) is the iteration counter. The computation of iterative update Ad^H 1J will be explained further below. (e) =(e) (b) In the local (element) phase we compute the values of the variables an^, fn+1 while keeping dn+'(i) fixed. In the rest of this section we will first describe in detail the phase (b). This will be followed by the description of the phase (a). ŠG= 0, wG= 2 o-X-o Figure 5.6: One point numerical integration scheme Slika 5.6: Enotockovna numericna integracijska shema We first provide the trial value of stress in the integration point of the element (we use a 1 point Gauss integration scheme, see Figure 5.6) * n+i" = *(e«+i ,aie) ,&)), (5.36) where = 0 is the coordinate of the integration point. Since we use the rigid plasticity law at the discontinuity, we are unable to determine the trial value of traction at the discontinuity by using the constitutive law. Instead we use the local equilibrium equation (5.19), i.e. r __i rx2 ttrial = _ G = __*trialdx (-ra+1 = Ga*n+1 dx = re / *n+1 dx Jne L Jx! = L ^ X1) ^fGx*n+1 = *n+1 , (5.37) =L =2 where is the Gauss quadrature point weight. Note that in the 1D setting, when using only one integration point, the traction at the discontinuity is equal to the stress in the integration point. Next we provide the trial value of the failure function at the discontinuity where 0 = c? - (*u - q!?/) < 0, (5.38) q^TTi = min[- ), *u]. (5.39) When computing the traction like softening variable in (5.39) we have to consider that the traction can never be less than zero therefore we set the upper value of this variable to the ultimate tensial strength of the material. If the trial failure criterion (5.38) is satisfied, the values of softening plasticity local variables remain unchanged =trial / \ , . =(e) =(e) 0 < 0 aie|1 = a!e), = Cn . (5.40) In the case of violation of the criterion (5.38) the values of softening plasticity variables are updated by using the evolution equations (5.30) and (5.31) and the backward Euler integration scheme ai+1 = a!e) + 7n+1, £n+1 = Cn + 7n+1, (5.41) where 7n+1 = 7n+1(Tn+1 - rn) is the plastic multiplier. The value of the plastic multiplier 7n+1 is determined from 0 = tn+1 (fn+O - (*u - qn+1 (^+1)) = 0 (fn+J = 0, (5.42) where we considered that the traction and the traction like softening variable are computed as tn+1 = °G+1 = E Bada,n+1 + Gaoq.i (^n+1^ , (5.43) = =(e) Qn+1 = min[-KsCn+1,Ou]- (5.44) For the linear softening case one can determine by considering (5.36) to (5.39) and (5.41) to (5.44) the value of plastic multiplier explicitly as (see e.g. [Ibrahimbegovic, 2009] Section 8.3.4 for more details) Tn+1 =trial 0 EG a - K s Ks =trial K s for Qn+1 < Ou „ p =trial 0 for Qn+1 = O (5.45) whereas for a nonlinear softening case an iteration procedure has to be used. Once the local variables are computed, we turn to the global phase (a) of the iterative loop in order to provide, if so needed, new iterative values of nodal displacements. First, the set of global equilibrium equations (5.16) is checked with newly computed oA ix Figure 5.9: Tension test of four parallel bars Slika 5.9: Natezni preizkus štirih vzporednih palic In the second example we consider a simple structure composed of four parallel bars, see Figure 5.9. We consider that all the bars have the same Young's modulus E = 1 and the same cross-section area A =1 but different material parameters that define the softening response, i.e. au1 = 0.6, Ks1 = -0.21, au2 = 0.8, Ks2 = -0.21, au3 = 1, Ks3 = -0.22, au4 = 1.15 and Ks4 = -0.20. The length of the bars is L =1. Each bar is moddeled with one finite element, which is supported at the left side and pulled at the right side, by imposing the displacement. Note that we considered that both of the nodes are common to all the elements, i.e. the actual vertical distance between the elements is zero. In the left hand side of Figure 5.10 we present the reaction versus Ux Ux Figure 5.10: Reaction versus imposed displacement of structure (left) and structural components (right) Slika 5.10: Krivulje reakcija - vsiljen pomik za celotno konstrukcijo (levo) in za njene posamezne dele (desno) imposed displacement curve for the structure and in the right hand side of the same figure we plot the reaction versus imposed displacement curves for each bar separately. Even though some of the components are in the softening regime the structure still exhibits the hardening like response. The ultimate load of the structure is reached at the moment when the last bar, i.e. Bar 4, reaches its maximum resistance and afterwards the resistance of the structure reduces with the increase of the displacement. 5.5 Concluding remarks and chapter summary In this chapter we have presented the embedded discontinuity concept on a 1D solid finite element. With the introduction of one additional parameter into the standard isoparametric 1d solid finite element we obtained the enriched displacement field 2 Uh(0 = Y, N*(®d* + Mj^, ^ j enrichment uFEM and the space derivative of this field gives us the enriched strain field h2 duh e = ~dx =2-* Ba(i)da + Ga(0a. a=1 "" ^ _/ ea We considered that the virtual strains are interpolated according to 2 e = Ba(C)da + Gaa, a=1 where the following modification ensures the convergence in the spirit of the patch test * 1 f Ga Ga T e I Gadx- We introduced the virtual strains into the principle of virtual work and obtained the virtual work of internal forces as 2 r VA(e) a=1 ^ 5Wnt'(e) = V A(e) daBaodx + A(e) aGaodx. e) a a a ne . JQe _ _ additional standard FEM From the term "additional" we obtained one additional equation per finite element h(e) = A(e) (^j GaOdx + t which ensures that the traction at the discontinuity is in equilibrium with the stress field in the element. We built the rigid plastic cohesive law (see Figure 5.5) by considering the localization criterion T>(t,q) = t - (Ou - q) < 0, and the softening potential = = 1 =2 S(£) = 1 KsC ■ The remaining ingredients of the rigid-plastic response were obtained by defining the plastic dissipation D = tet + qf, and by considering the principle of maximum plastic dissipation min max t,q Y —(t,q,Y) = -D (t,q) + Y 0(t,q) The computational procedure presented in Section 5.3 is split into the local and global phase. In the local phase we determine the updates for the variables related to rigid plasticity by solving the following equation % = tn+1 - - qn+1 (yn+0) = % (fn+o = 0 where we explicitly use the additional equilibrium equation to compute the traction at the discontinuity tn+1 = - GaOn+1dx = ^n+1- jne In the global phase, where the single element contribution to the system of global equations is K(e) K fa Khd K ha (i) ( Ad(e),(i) \ ( Jext,(e) _ fint,(e),(i) Adn+1 | = | J n+1 J n+1 n+^ Aantfw v 0 (= hn+f) we determine the updates for the current iterative values of nodal displacements d(e),(i) = d(e),(i-1) + Ad(e),(i-1) dn+1 = dn+1 + Adn+1 • The static condensation allowed us to form the standard form of the element stiffness matrix K(e)'(i) = K(e)>(i) _ Kfa,(i) (Kha,(i) \-1 Khd,(i) K n+1 = K n+1 K n+1 I ■Kn+1 I K n+1 > thus the global solution procedure is completely the same as in the standard finite element formulation. In the first numerical simulation in Section 5.4 we presented the convergence result for a test with one finite element, while in the second simulation we show the results for a structure composed of four parallel bars. Chapter 6 Failure analysis of metal beams and frames 6.1 Introduction A typical example of a failure (collapse) analysis is the the push-over analysis in earthquake engineering. This is a non-linear static analysis of a building structure, subjected to an equivalent static loading that is pushing a structure over its limit capacity (e.g. [Fajfar et al., 2006]). It has been observed from failure modes, produced by seismic activities and experimental tests, that practical frame structures, composed of columns and beams, fail by exhibiting localized failures in a limited number of critical zones. Those critical zones are usually described as plastic (inelastic) hinges. A usual approach to compute the limit load of a structural frame, or to compute its complete failure, is to model plastic hinges with nonlinear inelastic spring finite elements. Inelastic springs are introduced at predefined critical locations in a mesh of conventional elastic beam finite elements (e.g. [Wilson, 2002]), or, alternatively (e.g. [Powell, 1986]), elastic beam elements with lumped nonlinear spring at both ends are used. In this chapter we present another option, i.e. a nonlinear analysis using the embedded discontinuity finite element method. In the first part of this chapter, we carry on with the developments related to numerical treatment of localized failure in beams in order to study failure of elastoplastic metal frames. To this end, we derive a planar straight stress-resultant beam finite element with the following features: (i) Euler-Bernoulli kinematics, (ii) an elastoplastic stress-resultant constitutive model with isotropic hardening, (iii) a localized softening plastic hinge related to the strong discontinuity in generalized displacements, and (iv) an approximation of the geometrically nonlinear effects by using the von Karman strains for the virtual axial deformations. The derived finite element can be effectively used for the limit load analysis, the push- over analysis and the complete failure analysis of planar metal frames. Localized softening, introduced by embedded discontinuity approach, solves the problem of mesh-dependency. Moreover, the spreading of plasticity over the entire frame and the appearance of the softening plastic hinges in the frame is consistently accounted for in the course of the nonlinear analysis. With respect to the existing embedded discontinuity beam finite elements, see [Ehrlich and Armero, 2005], [Armero and Ehrlich, 2006b], [Armero and Ehrlich, 2004], and [Wackerfuss, 2008], we use more complex material models: stress-resultant elastoplasticity with hardening to describe beam material behavior and stress-resultant rigid-plastic softening to describe material behavior at the discontinuity. The second part of this chapter pertains to a procedure that provides characteristic values of material parameters, used by chosen inelastic models. Those values are the yield and the failure (ultimate resistance) moments of the beam cross-section, the hardening modulus for the stress-resultant beam plasticity, and the softening modulus for the softening plastic hinge. Ideally, one should for any geometry of beam cross-section, any material type and any type of beam stress state seek the appropriate experimental results and fit to them the beam model material parameters with respect to significant quantities (e.g. forces, displacements, energy, dissipation), see e.g. [Kucerova et al., 2009]. In the absence of experimental results for metal beams to make any definitive conclusions, we turn to another approach that belongs under multi-scale label. The material parameters are obtained by numerical simulations on representative part of a beam by using a refined model, which is superior to the beam model in a sense that it is able to describe in more detail the beam response. We focus on rather typical metal frames with thin-walled cross-sections. For this kind of frames, the refined model can be chosen as the nonlinear shell model (e.g. [Brank et al., 1997], [Brank, 2005]). The shell model is superior to the beam model in providing a proper local description of the strain/stress fields and the overall spread of plasticity. It is also capable of describing local buckling of the flanges and the web, which is, in bending dominated conditions, very often the reason for the localized beam failure. Considering the above, the shell model can be seen as the meso-scale model and the beam model as the macro-scale model. The outline of the chapter is as follows. In Section 6.2, we derive an elastoplastic Euler-Bernoulli beam finite element with embedded discontinuity. In Section 6.3, we discuss computation of the beam plasticity parameters and the softening plastic hinge parameters by using the shell model. In Section 6.4, we present details of the computational procedure. Numerical examples are presented in Section 6.5 and concluding remarks in Section 6.6. 6.2 Beam element with embedded discontinuity We consider in this section a planar Euler-Bernoulli beam finite element. The element can represent an elastoplastic bending, including the localized softening effects, which are associated with the strong discontinuity in rotation. The geometrical nonlinearity is approximately taken into account by virtual axial strains of von Karman type, which allows this element to capture the global buckling modes. 6.2.1 Kinematics We consider a straight planar frame member, which middle axis occupies domain Q G R. Spatial discretization of Q leads to Ne1 (Q = [0,L] = U^L^) finite elements. A typical U1,w1,w\ auAe U2,W2,W2 x Xd L^ Figure 6.1: Beam finite element with embedded discontinuity Slika 6.1: Koncni element za nosilce z vgrajeno nezveznostjo 2-node finite element is presented in Figure 6.1. The following notation is used: u are nodal axial displacements, w are nodal transverse displacements, w'i are nodal values of the beam axis rotation (derivative of transverse displacement with respect to the beam axial coordinate x G [0,L(e)]), and i = 1, 2 is node number. In addition to the standard degrees of freedom at the two nodes, we assume strong discontinuity in axial displacement and beam axis rotation a at xd G L(e). We also assume that the domain of the discontinuity influence corresponds to a single element. The axial displacement is thus defined as: uh(x,xd) = Nu(x)u + Mu(x,xd)au, (6.1) where (x) = [1 — x/L(e),x/L(e)], u = [u1,u2]T, and Mu(x,xd) is a function with zero values at the nodes and a unit jump at xd, i.e. Mu(0,xd) = Mu(L(e),xd) = 0 and Mu(x+,xd) = Mu(x-,xd) + 1. Similarly, we can write the transverse displacement as wh(x,xd) = (x)w + (x)w' + Me (x,xd)ae, (6.2) where Nw(x) = [2( —_)3 - 3( )2 + 1, -2( —-))3 + 3( —-))2], w = [W1,W2]t, 'L(e ) —(e) * —(e) ' Nw'(x) = L(e)[(^X))3 - 2(^)2 + (^), (-^)3 - (^)2] x x x x * —(e) * —(e) — (e) W" —(e) w [w1,w2] /iT (6.3) (6.4) and Me (x,xd) is a function with zero values at the nodes and a unit jump of its first derivative at xd, i.e. Me(0,xd) = Me(—(e),xd) = 0 and M&'(x+,xd) = M&'(x-,xd) + 1. The beam axial strain can then be written as: duh e(x,xd) dx Bu(x)u + Gu(x, xd)au + 8xdau (6.5) where Bu(x) = [-1/—(e), 1/—(e)], Gu(x,xd) = d M u/dx, and 8Xd is the Dirac-delta, which appears due to discontinuous nature of axial displacement at xd. We further divide the axial strain into a regular part e and a singular part e. The later can be interpreted as a localized plastic axial strain. The beam curvature is computed as: d 2wh K(x, xd) dx2 Bw (x)w + Bw (x)w' + G°(x, xd)ae + Sxdae, (6.6) where B w (x) Bw' (x) [- — (1-—) — (e) T(e)> —(e) —(e) 6 (1-—)] (1 —(e) )]' 2 3x 2 3x [- (2 - —-j )'- (1 - —-J )]' (6.7) (6.8) and Ge(x,xd) = d2Me/dx2. The curvature k is divided into a regular part K and a singular part K. The later can be interpreted as a localized plastic curvature. The beam strains can be rewritten in a matrix notation as e = e + e, e = Bd + Ga, e = 5xd a, (6.9) (6.10) where e = [e, k]t , e = [e, k]t , e T e = \e, K T and B G Bu 0 0 0 Bw Bw' DIAG [Gu, Ge] , a d = [uT, wT, w'T]T, = [a„, ad]T . (6.11) (6.12) Kinematic description of the element is concluded by derivation of G operator. It may be derived indirectly (i.e. without defining Mu and Md) through requirement that U1,W1 b<- Xd Le) U2 = U1+au W2 = W1 + W1' Le+ag(Le-xd) w1+ae w2' = w{+ae h*->H Figure 6.2: Strain-free mode of the element Slika 6.2: Element v brezdeformacijskem stanju an element has to be able to describe strain-free mode at some non-zero values of aM and , see [Armero and Ehrlich, 2006b]. According to Fig. 6.2, the generalized nodal displacements dhinge = [U1,U2,W1,W2, W', W2]T of such strain-free mode are composed as dhinge drigid + D hinge D hinge 0 1 0 0 0 0 0 0 0 L(e) - xd 0 1 T (6.13) where drigid = [u1,U1,W1,W1 + W'L(e),W', W'] are generalized nodal displacements due to rigid-body motion of a complete beam, and Dhingea are generalized nodal displacements due to rigid-body motion of one part of the beam due to imposed strong discontinuity a = [aM, a^] . If we now set strains (6.9) to zero for dhinge, we have 0 = Bdhinge + Ga = Bdrigid + (G + BDhinge)a. (6.14) Since the above equation should hold for any a, we get the G operator as G = — BDhi hinge , (6.15) which leads to Gu(x,xd) = - Ley, (6.16) n®t \ 1 + 3(1 - L2xey)(1 - L^T) G (x,xd) =---L-. (6.17) The above definition of G matrix concludes kinematic description of the geometrically linear element. Remark 6.1. By using (6.1) and (6.2) to describe strain-free mode of Fig. 6.2, one can also derive interpolation functions Mu and Me. By setting in (6.1) u1 = U1 = 0, u2 = u2 = au, uh = 0 for x < xd, and uh = 0 for x > xd , one can conclude that Mu = H(x — xd) — Nu • [0,1]. Here, H(x — xd) is unit-step function, which is 0 for x < xd and 1 for x > xd. Derivation dMu/dx gives Gu in (6.16). By using similar procedure for bending in (6.2), one can obtain Me = (H(x — xd)) (x — xd)— NW • [0,L(e) — xd^ — Nw' • [0,1]. Derivation d2Mu/dx2 gives Gd in (6.17). In order to account for the geometrically nonlinear effects, and related global buckling, we will use the von Karman axial strain when computing the virtual axial strain. The real axial strain, used for computing the internal forces, will still be assumed as linear, as given in eq. (6.5). The von Karman axial strain is defined as eVK = dX + 2(dwXh)2. The corresponding virtual axial strain is thus: 8E^VK = ^ + dwh ^. (6.18) dx dx dx If we choose to interpolate Suh, wh and 8wh in (6.18) as Suh = Nu(x)Su+ Mu(x,xd)5au, wh = Nw(x)w+ Nw'(x)w' and 5wh = Nw(x)5w+ Nw'(x)Sw', where 5u = [^ui,^u2]T is vector of virtual nodal axial displacements, 5w = [^w1,^w2]T and 5w' = [^w1,^w2]T are vectors of virtual nodal transverse displacements and rotations, and 5au is virtual discontinuity in axial displacement at xd, the chosen interpolations lead to SeVK = Bu (x) 5u + Buw (x) 5w + Bu'w' (x) Sw' + Gu(x, xd)5au + 5Xd5au (6.19) SeVK where dNw ' dNw' (dN w dNw' Bu'w (x) = CdN~, Bu'w (x) = CdN~, C = dN- • w + dNx~ • w'l . (6.20) dx dx dx dx The linear matrix operator B from (6.11) should be thus replaced with the nonlinear matrix operator BVK when computing virtual strains Se = [Se, Sk]t, i.e. Sd + GSa. (6.21) ' Se = SeVK ' Bu B u,W BU,W SK 0 B w bW' "V- z>VK In (6.21) above, we denote with Sd = [SuT,SwT,Sw'T]T the generalized virtual nodal displacements and with Sa = [Sau, Sag]T virtual jumps at xd. Remark 6.2. The tangent stiffness matrix of the beam finite element with von Karman virtual axial strain has symmetric geometric part and non-symmetric material part. The matrix can be symmetrized by using B instead of BVK when computing its material part. Such an approach would lead (for elastic beams) to the element presented in [Wilson, 2002], section 11. In this work we use non-symmetric tangent stiffness matrix. Remark 6.3. (b) If one uses von Karman definition of axial strains for both real and virtual strains, see [Reddy, 2004], section 4-2, the element exhibits severe locking. 6.2.2 Equilibrium equations The weak form of the equilibrium equations (the principle of virtual work) for an element e of a chosen finite element mesh with Ne1 finite elements, can be written as: £nmt,(e) _ ^next,(e) = g. (6.22) By denoting the virtual strains as če = , T, where virtual curvatures čk = + are of the same form as real curvatures k in (6.6), we can write a single element contribution to the virtual work of internal forces as: f L(e) čnmt,(e) = (če)T adx 0 n L(e) r-L(e) = (ByK)T adx + / čaT(GTa + a)dx, (6.23) Jo Jo v _ J v _ standard additional The matrices Band G are defined in (6.21) and (6.15), and a = [N, M]t (6.24) is the vector of beam internal forces that contains axial force N and bending moment M. From the term "standard" in (6.23) we obtain the vector of element internal nodal forces , L(e) f mt,(e) = (B yK)T adx. (6.25) 0 From the virtual work of external forces čnext,(e) we can get the vector of element external nodal forces fext,(e), representing the external load applied to the element. The finite element assembly of vectors fmt,(e) and fext,(e), for all elements of the chosen mesh, leads to a set of global (i.e. mesh related) equations ANi f int,(e) _ fext,(eA = 0, (6.26) where A is the assembling operator. We have only used one part of the right side of equation (6.23) in (6.22) when getting the set of global equations (6.26). The other term in (6.23), denoted as "additional" (since it results from additional enriched kinematics due to embedded discontinuity), will also contribute to the weak form of the equilibrium. However, we will treat this contribution locally element by element. Then, in view of (6.22), the following two equations are obtained for each element of the chosen mesh f-L(e) » /» h(e) h ' N , ' M T I (GT a + 5Xda)dx 0 L(e) ,L(e) TT = GT adx + a|Xd = GT adx + t = 0, Ve G [1,Ne1]. (6.27) jo j0 =t We have defined in (6.27) vector t = a|Xd = [tN,tM]T with components tN and tM that represent axial traction and moment (bending) traction at the discontinuity. By using (6.17) and (6.24), one can obtain the component form of (6.27) , L(e) hN = GuNdx + tN = 0, 0 n L(e) hM = GdMdx + tM = 0, Ve G [1,NeJ]. (6.28) 0 The problem of solving a set of global equations (6.26) together with a set of local (element) equations (6.27) will be further addressed in Section 4. 6.2.3 Constitutive relations We assume that the axial response of the beam material remains always elastic, thus discarding the failure by necking, for example. For the bending behavior of the beam material we choose the following constitutive models: (i) stress-resultant elastoplastic constitutive model with linear isotropic hardening, (ii) stress-resultant rigid-plasticity model with linear softening at the softening plastic hinge. The basic ingredients of the chosen constitutive relations are built on classical plasticity (e.g. [Ibrahimbegovic et al., 1998]) and can be summarized as: • The regular strains e (6.10) can be additively decomposed into elastic part e" and plastic part ep e = ee + ep, e" = [ee,Ke]T, ep = [ep,Kp]T. (6.29) • The axial strain of the beam (6.5) remains always elastic, thus e = e = ee, e =0 ^^ e = 0, au = 0. (6.30) • The free energy for the beam material (before localized softening is activated) is assumed to be the sum of the strain energy function W and the hardening potential " _ _ _ __1 1 _2 *(ee, £) := W (ee) + -(£) = -eeTCee + ^Kh£2, (6.31) where C = DIAG [EA, EI] , E is elastic modulus, A and I are area and moment of inertia of cross-section, £ > 0 is strain-like bending hardening variable, and > 0 is linear bending hardening modulus. • The yield criterion for the beam material is defined in terms of the bending moment. The admissible values of the bending moment and the stress-like bending hardening variable q(£) are governed by the function 0(M,q) = |M|- (My - q) < 0, (6.32) where My > 0 denotes the positive yield moment of the cross-section. Influence of the axial force N on the cross-section yielding is taken into account by defining My and q as functions of N, as shown subsequently. The bending traction tM at the discontinuity is related to the rotation jump as shown in Fig. 6.3 tM = t m (a )• (6.33) tM' i Mu i * s -Mu Ao Figure 6.3: Rigid-plastic cohesive law at discontinuity Slika 6.3: Kohezijski zakon (toga plasticnost) v nezveznosti • Cohesive law (6.33) can be written in terms of localization (failure) criterion that activates softening at discontinuity at and is defined in terms of the bending traction tM and the stress-like softening bending variable q(£) (the later is defined in terms of the bending strain-like softening variable £ > 0) 0(tM,q) = |tMI - (Mu - q) < 0, (6.34) where Mu > My > 0 denotes the positive ultimate (failure) moment of the cross-section. Influence of axial force N on the cross-section failure is taken into account by defining Mu and q as functions of N, as shown below. • The softening potential is assumed to be = = 1 =2 S(£) = ^, (6.35) where < 0 is the linear (bending) softening modulus. The remaining ingredients of the elastoplasticity with hardening can be obtained from the consideration of thermodynamics of associative plasticity and the principle of maximum plastic dissipation (see e.g. [Ibrahimbegovic, 2009], [Lubliner, 1990], [Simo and Hughes, 1998]). In the present beam model the elastoplasticity with hardening happens for a = 0, which leads to e = e and e = 0, see (6.10). By using (6.29) and (6.31) the mechanical dissipation can be written as 0 < D d= aTe - *(r,č) = (a - H)Tee + ^ - , (6.36) where (o) = d (o) /dt and t G [0,T] is pseudo-time. By assuming that the elastic process is non-dissipative (i.e. D = 0), and that the plastic state variables do not change, we obtain from (6.36) d W a = — = Cee N = EAe, M = EI (K - «p) . (6.37) dee We can define the hardening variable q by further considering (6.36) and (6.31) dW 5Š -q = - 7^= - ^ = -KhC. (6.38) de de By replacing (6.37) and (6.38) in (6.36), the plastic dissipation can be obtained as D = aTf + (6.24), (6.30) D^ = MKP + (6.39) The principle of maximum plastic dissipation states that among all the variables (M, q) that satisfy the yield criteria 0 (M, q) < 0, one should choose those that maximize plastic dissipation (at frozen rates Kf and £). This can be written as a constrained optimization problem: min max [XP(M,q,Yy) = -D(M,q) + Y0(M,q)l , (6.40) where 7 > 0 plays the role of Lagrange multiplier. By using (6.39) and (6.32), the last result can provide the evolution equations for internal variables JM = —Jf + YH = 0 ^ K = si«n 0, < < 0, 7< =0, 7< = 0. (6.43) To obtain the remaining ingredients of the rigid-plastic response describing softening at the discontinuity xd, let us isolate the softening plastic hinge. We first define the strain energy function due to softening potential as f = S. The dissipation at xd can be then written as: _ 0 < D d= tMag — f (!) = tMag — dff. (6.44) d! where tM is the discontinuity bending traction given by (6.33). By defining - df dS = = q = — df = — ds = — Ks! = K !, (6.45) d! d! the result in (6.44) can be rewritten as D = D = tM a g + (6.46) The principle of maximum plastic dissipation at the rigid-plastic discontinuity can then be defined as: min max L (tM ,q,7) = —D (tM ,q) + 7<(tM ,q) , (6.47) t m ,q j l j where 7 > 0 is the Lagrange multiplier. By using (6.46) and (6.34), we get from (6.47) above the following evolution equations: =p = dL = d< = -w— = — ag + 7-w— = 0 ag = sign(tM)7, (6.48) dtM dtM =p = dL = — d< = — . . "Tp" = —! + 7^= = 0=^ ! = 7 (6.49) dq dq By observing that sign(tM) = sign(ag) (see (6.33) and Fig. 6.3), it follows from (6.48) that sign(ag)a,g = 1 ^ lag| = 1- (6.50) The Kuhn-Tucker loading/unloading conditions and the consistency condition also apply: Y > 0, 0 < 0, = 0, Y0=0. (6.51) With the above results, we are in position to write the total dissipation of the beam finite element when the element is in the softening regime. Namely, by accounting for the proper definition of strain energy terms for the beam finite element according to L(e)__= ^ = f0 ^dx + the total dissipation in the softening regime can be written as ,L(e) DLt) = jf (^Te - ^(ee, £)) dx + ^mao - *(£)) + Ga - + q£^ dx + (tM a o + f £ ) (6.52) m/ + q£ j dx + yj Go Mdx +tM J a o + • Dp, see (6.39) "-^---' |Ks|l v ' =0, see (6.28) It can be seen from (6.52) that enforcing eq. (6.28) will decouple dissipation in the softening plastic hinge from the dissipation in the rest of the beam. Therefore, eq. (6.28) is further used to compute tM. We conclude description of constitutive relations by defining plastic work of the beam cross-section in the hardening regime W = MK = |M| £"= (My + Kh£) £, (6.53) and plastic work for the beam finite element in the softening regime as W = tMao = |tMI £ = (Mu + Ks£) I• (6.54) 6.3 Computation of beam plasticity material parameters In the previous section, we have built the framework for stress-resultant plasticity for beam finite element with embedded discontinuity. The material parameters that need to be known for chosen material models are: (i) My and Kh for the plasticity with hardening, and (ii) My, and Ks for the softening plastic hinge. In this section we will elaborate on determination of the above parameters. The yield moment My can be determined by considering the uniaxial yield stress of the material ay, the bending resistance modulus of cross-section W, the cross-section area A, and the level of axial force N. One can associate the yield moment of the cross-section with the yielding of the most-stressed material fiber to get My (N) = Way(1 - M). (6.55) The ultimate bending moment Mu can be derived in a closed form by assuming elastic-perfectly-plastic response of material fibers, e.g. [Lubliner, 1990]. However, one may try to determine a better estimate for Mu, which takes into consideration material hardening, as well as possibility of local buckling (e.g. buckling of the flanges and/or the web of the I-beam). This task is addressed in the present work by performing computations with refined finite element model based on geometrically and materially nonlinear shell element, which is able to capture local buckling and gradual spreading of plasticity over the cross-section. The ultimate bending resistance Mu can be obtained by using results of such a shell model computations, as can be moduli Kh and Ks. To obtain desired results, a part of the frame member with a reference length Lref < Ltot (Ltot is the total length of the frame member under consideration) is: (i) modeled with shell finite elements, (ii) subjected to an external axial force N in the first loading step, and (iii) subjected to a varying external bending moment at the end cross-sections in the second loading step, while keeping N fixed, see Fig. 6.4 (b). It is assumed that such a loading pattern would produce approximately constant internal axial force N = NV during the analysis. The computation with shell model takes into account geometrical and material nonlinearity that include: plasticity with hardening and strain-softening, strain-softening regularization, and local buckling effects. The results of shell analysis are cast in terms of diagrams presented in Figs. 6.4 (d) and (g). One can associate the ultimate bending moment Mu with the peak point in the diagram at Fig. 6.4 (d), where applied end moment is plotted versus the end rotation, i.e. Mu (N) = Mruef (N). (6.56) One can also use this point as a border-point between the hardening regime and the softening regime, where the softening can be due to material and/or geometric effects. To determine the values of the beam model hardening and softening parameters, we make an assumption that the plastic work at failure should be equal for both the beam and the shell model. In other words, we want the internal forces of the beam model to produce the same amount of the plastic work as the internal forces of the shell model, when considering the full failure of the part of the frame member of length Lref. Since the plastic work is done in two regimes (hardening and softening), we have to assure that the amount of plastic work in each regime matches for both models, i.e. TT7V 777P>re/ TT7P TT7P,Vef EW (N) = EW (N), EW (N) = EW (N). (6.57) (a) Beam model (b) Refined model M M Lref (c) Moment - curvature/rotation (beam model) Mu I k | = J = y 1 A1 = š = y (e) Moment - hardening variable (f) Moment - softening variable (beam model) mA My (beam model) M Mu * J N = const. J J Figure 6.4: Evaluation of beam material parameters by using results of refined analysis Slika 6.4: Dolocitev parametrov za nosilce z uporabo rezultatov podrobnejse analize M M K, A The plastic work in the hardening regime, EW , and the plastic work in the softening -P, re f regime, EW , are obtained from the shell model analysis, Fig. 6.4 (g). TTrP The plastic work of the beam model in the hardening regime, EW , can be determined by observing that each cross-section in the frame member of length Lref is approximately under the same force-moment state during the hardening regime. Integration of (6.53) allows us to write wp (Lref (t at , - 1 -2 EW= WPdrdl = Lref (My- + -Kh- ). (6.58) Jo Jo 2 In (6.58) above - is the value of hardening variable that corresponds to the bending moment Mu, Fig. 6.4 (e). Since we have assumed linear hardening in the beam model (6.38), we get - Mu — My , , - = K y. (6.59) Kh By using (6.59), (6.58) and (6.57), one can obtain an expression for hardening modulus as (Mu2 (N) - M2 (N))Lref Kh (N) = ( u ( ) Wp,refy ,( )-. (6.60) 2EW f (N) The plastic work of beam model in the softening regime, EW , can be determined by assuming that the softening part of M — (p curve in Fig. 6.4 (d), obtained from the shell model analysis, is produced by a very localized phenomenon (in a single cross-section) related to the local buckling and/or to the localized strain-softening. By using (6.54), one can compute the plastic work in the softening regime for the beam model as =p ft_p ri at t m=o , =\ = 1 -2 EW =1 W dr = J (Mu + Ks-Jd- =2 K(6.61) In (6.61) - is the value of the softening variable that corresponds to the total cross-section failure, Fig. 6.4 (f). Since we have assumed linear softening in the beam model, we obtain !=j^j. <6-62) By using (6.62), (6.61) and (6.57), one can obtain an expression for softening modulus as Mu2 (N) |Ks (N)| = =u,r(e/ ) , Ks < 0. (6.63) 2EW (N) We note, that the choice of the reference length Lref should have very small influence on the values of the searched material parameters. The influence on the value of Kh should be small, since each cross-section is approximately under the same force-moment state during the hardening regime. The influence on the value of should not be too big neither, since the softening effect is localized. However, one should perform large displacement correction of Mu, if the chosen length of Lref enables large deflections, as shown in example 5.4. 6.4 Computational procedure In this section we will present a procedure for solving the set of global (mesh related) and the set of local (element related) nonlinear equations generated by using the stress-resultant plasticity beam finite element with embedded discontinuity presented in section 2. The solution of the set of global nonlinear equations (6.26), along with the set of local nonlinear equations (6.27) (note that (6.27) is reduced to (6.28) due to assumption (6.30)), ought to be computed at discrete pseudo-time values 0, t1, t2,..., tn-1, tn, tn+1,..., T by means of the incremental-iterative scheme. We will consider the solution in a typical pseudo-time incremental step from tn to tn+1. Let us assume that all the variables, related to an element e and its integration points ip = 1, 2, 3 (a 3-point Lobatto integration scheme is used) are given at tn, i.e. given: d^, <'ip, ČT, a£, and xf\ M^e). (6.64) We have also added in (6.64): (i) the yield moment at integration point Myp (which is only true if hardening plasticity has been activated so far) and (ii) position of the discontinuity xde) and the ultimate bending moment M-Ue) (which is only true if softening plastic hinge has been activated so far). We will then iterate in the pseudo-time step in order to compute the converged values of the variables at tn+1, i.e. find: ^ ^ £n^ ^ m+1 and (if not given already) Myp, xde), M^e). (6.65) The moments Myp and M-ie) are computed by using (6.55) and (6.56). Although they depend on axial force N, we keep them fixed once determined. The computation of solution (6.65) is split into two phases: (a) The global (mesh related) phase computes the current iterative values (with (i) as the iteration counter) of nodal generalized displacements at tn+1 while keeping the other variables fixed, i.e. global phase: d^f = dn+f-1) + Ad^f-1^ (6.66) The computation of iterative update Ad^l 1) will be explained below. (b) The local (element and integration point related) phase computes the values of variables K+1, a £n+1 = L . (6.75) In the case of violation of the trial yield criterion (6.74), the values of local variables are updated by backward Euler integration scheme «£+1 = a£ + sign^n+1 > ln++1 = in ) + > (6.76) where = 7yn-+1 (tn+1 — tn). The value of the plastic multiplier is determined from condition 0{e)(tMn+1(a^enn+1(Y^1)),q(1^1(Y^1))) = 0(e)(Yn+) = 0. (6.77) For the linear softening one can determine the plastic multiplier explicitly, whereas for nonlinear softening an iterative solution procedure is needed. Note, that we compute the bending traction in (6.77) as tgn+1 = — / G*(xde),x)M(^«^,<+1)dx. (6.78) Jne The main result of the above described softening plasticity procedure is the new value of softening variable a*e)+1, which influences the stress state of the whole element by giving the new values of the bending moment Mn+1 as M+1 = EI +1) — Knip). (6.79) The updated value of the plastic curvature is Kpv+p1 = Wn>ip. Once the local variables are computed, we turn to the global phase of the iterative loop in order to provide, if so needed, new iterative values of nodal displacements. First, the set of global equilibrium equations (6.26) is checked with newly computed Mn+1 from the local phase ANei r pin,t,(e) A:=1[f n+1 /ext,(e),(i) n+1 ? < tol. (6.80) If the convergence criterion (6.80) is satisfied, we move on to the next pseudo-time incremental step. If the convergence criterion fails, we perform a new iterative sweep within the present pseudo-time incremental step. New iterative values of nodal generalized displacements of the finite element mesh are computed by accounting for each element contribution. A single element contribution can be written as K(e) Kfa Khd Kha (i) n+1 Adn+1 Aae,n+1 f ext,(e) *int,(:),(i) n+1 f n+1 (6.81) where the parts of the element stiffness matrix can be formally written as K (e),(i) n+1 f it,(e) dd(e) K hd,(i) n+1 _ I dhM dd(e) (i) n+1 (i) n+1 K fa,(i) n+1 f it,(e) da (e) (i) n+1 K ha, ( i) n+1 dh (e) da (e) (i) n+1 (6.82) (e),(i) The static condensation of (6.81) allows us to form the element stiffness matrix Kn+1 that contributes to the assembly an=1 fK.::+f)Ad(i) n+1 aN=1 (f r:xt,(e) n+1 -f int,(:),(i) n+1 (6.83) where K (e),(i) n+1 K(:),(i) _ Kfa,(i) I Kha,(i) K n+1 K n+1 \ Kn+1 1 K hd,(i) n+1 . (6.84) Solution of (6.83) gives the values of iterative update Adn+1i), which leads us back to (6.66). 0 9 e 6.5 Examples In this section we illustrate performance of the above derived beam element when analyzing push-over and collapse of steel frames. We also illustrate the procedure, presented in Section 6.3, for computing the beam model plasticity material parameters by using the shell finite element model. The beam model computer code was generated by using symbolic manipulation code AceGen and the examples were computed by using finite element program AceFem, see [Korelc, 2007b] and [Korelc, 2007a]. 6.5.1 Computation of beam plasticity material parameters With this example we illustrate computation of beam plasticity material parameters Mu, and as suggested in Section 6.3. We consider a frame member with an I-cross-section with flange width 6/ = 30 cm, flange thickness t/ = 1.5 cm, web height = 40 cm and web thickness tw = 0.8 cm. The cross-section area is A = 122 cm2 and the bending resistance modulus is I = 43034.2 cm4. We model a part of the frame member of length Lre/ = 2L = 300 cm, which is 7 times the height of the section. This length should be sufficient to capture the local softening effects due to local buckling and/or strain softening. The frame member is made of an elastoplastic material (steel), whose uniaxial response is plotted in Fig. 6.5. Young's modulus is E = 21000 kN/cm2, yield stress is = 24 kN/cm2, ultimate stress is = 36 kN/cm2, yield strain is = /E, strain at ultimate stress is = 0.1 and strain at failure is £/ = 0.12778. The example has been computed with the finite element code ABAQUS [Hobbit et al., 2007] by using shell finite element S4R with 5 integration points through the thickness. Only one half of the considered geometry was discretized, see Figs. 6.6 and 6.7. The symmetry conditions = = = 0 were used in the symmetry plane. The mesh consists of equal squared elements. The free-end cross-section of the model was made rigid by coupling the degrees of freedom of that cross-section. The plasticity models with strain softening are mesh-dependent. In order to minimize that effect, we have adjusted the post peak uniaxial stress-strain relation to fit the mesh size, as suggested in [Ibrahimbegovic, 2009]. According to [Ibrahimbegovic, 2009] the linear softening modulus is computed as 2 = - ^, (6.85) where le is a characteristic dimension of the element (in the present case le = 5 cm is the side-length of the elements) and is the plastic work density (plastic work per unit volume) in the softening regime of the shell material that corresponds to the gray-colored area in Figure 6.5 (in the present case = 0.5 kN/cm2). The strain at failure, adjusted to the mesh, is then / = + J^t = 0.10727. (6.86) E+KF The load was applied in two steps. In the first step we applied a desired level of axial force N at the mid-point MP of the rigid cross-section, see Fig. 6.6. In the second step we applied bending moment M at the point MP and performed nonlinear analysis with the path-following method. Several analyses were carried out with different values J 40 (-Ef,0) , 20 (E y,J y) "" -0.035Ny ' (6.87) Figure 6.10: Approximation of the ultimate bending moment of the cross-section Slika 6.10: Aproksimacija mejnega upogibnega momenta prereza where Mruef'° = MUef (N = 0). We assume that Mr^f (N) can be used to evaluate the ultimate bending moment of the beam model Mu, i.e. Mu (N) = Mruef (N), (6.88) see Fig. 6.10. The values for the ultimate resistance obtained with the shell analyses are marked with dots in Fig. 6.10. The beam model hardening modulus Kh can be evaluated point-wise by using (6.60), (6.55), (6.88) and third column of Table 6.1. We get Kh ranging from 6.26 ■ 106 kN/cm2 to 1.35 ■ 107 kN/cm2 with an average value of 1.06 ■ 107 kN/cm2. Although one could easily find a higher-order function that fits these results, we assume for simplicity that the axial force has no influence on the hardening modulus and adopt Kh(N) = 1.06 ■ 107 kN/cm2. (6.89) The beam model softening modulus Ks can be evaluated point-wise by using (6.63), (6.88) and the last column of Table 6.1. The gray-colored fields in Table 6.1 present ==p,ref unreliable results for EW , since for those cases the shell analysis did not converge. Softening modulus for the first three analyses ranges from -2.85-105 kN/cm2 to -3.97-105 kN/cm2. We assume that the axial force has no influence on softening modulus and adopt the average value Ks(N) = -3.28 ■ 105 kN/cm2. (6.90) In Table 6.2 we make a point-wise comparison between the shell analysis results Muef, -rrrp,ref -—p,ref tt/P U7P EW and EW and the corresponding beam model results Mu, EW and EW , computed by using approximations (6.88), (6.89), (6.90) and expressions (6.58), (6.61) and (6.53). We can see that the error in ultimate bending moment is small, while the error in dissipated plastic work can be quite large. N/Ny Mu [kNm] EWP [kJ] EWP-zW^f- EW [kJ] ==p ==p,ref. I EW -EW I % 1 Mf 1 [%] 1 EWP-ref- 1 [%] 0 550 0.00 50 42 46 12 -0.1 519 0.26 59 20 41 8 -0.2 473 1.54 54 17 34 19 -0.3 426 0.22 48 12 28 ? 0.1 550 1.94 81 14 46 ? 0.2 550 5.07 109 27 46 ? 0.3 550 3.72 134 6 46 ? Table 6.2: Comparison between approximations and shell analyses results Tabela 6.2: Primerjava med aproksimacijami in rezultati analiz z uporabo konCnih elementov za lupine 6.5.2 Push-over of a symmetric frame In this example we present a push-over analysis of a symmetric frame. The geometry is given in Fig. 6.11, where LB = 500 and He = 250 cm. The material and cross-section properties of all frame members are equal. They are the same as those presented in Section 5.1. The vertical load is constant and equals qv = 0.05 kN/cm. The lateral loading is presented in Figure 6.11, where F0 = 1 kN is a concentrated force and A is load multiplier. The mesh consists of eight beam finite elements per each frame member. We performed two analyses, one by using the geometrically linear and the other by using the geometrically nonlinear beam finite elements. The results are presented in Figures 6.12 and 6.13, where utop is horizontal displacement of the top right corner of the frame. In the left part of the Fig. 6.12 we present the total lateral load versus utop curves. The points on those curves mark configurations where the softening hinge was activated in one of the elements of the mesh. We can see that, even though some parts of the frame are failing, the total resistance of the structure is still growing until the maximum load is reached at 1527.3 kN for geometrically linear case, and at 1522.3 kN for geometrically nonlinear case. In the right part of Fig. 6.12 we present the plastic work versus utop displacement curves. The results of the geometrically linear and the geometrically nonlinear elements are completely the same. At the beginning of the analysis there is no energy dissipation since the material response is elastic. The non-dissipative period is followed by a short period with dissipation due to material hardening only, which ends with the first activation A4F0 qv 1 A8F0 qv 1 A8F0 qv 1 A4F0 utop qv qv qv HC A3F0 M A6F0 M A6F0 M A3F0 qv qv qv HC A2F0 M A4F0 M A4F0 M A2F0 qv qv qv HC AF0 M A2F0 M A2F0 M AF0 Hc Lb__LB_^_LB Figure 6.11: Symmetric frame: geometry and loading Slika 6.11: Simetricni okvir: geometrija in obtezba Figure 6.12: Load versus displacement and dissipated energy versus displacement curves for symmetric frame Slika 6.12: Krivulje obtezba - pomik in disipirana energija - pomik pri simetricnem okvirju of softening plastic hinge in one of the beam finite elements. For a while we have a combined hardening and softening energy dissipation, which is finally followed by a period when the structure is dissipating energy only due to softening. On the left part of Fig. Figure 6.13: Symmetric frame: deformed shape and locations of softening plastic hinges at utop ~ 60 cm Slika 6.13: Simetricen okvir: deformirana lega in mesta plasticnih clenkov pri pomiku gornje eteze utop ~ 60 cm 6.13 we present the final deformed configuration of the frame. In the right part of Fig. 6.13 we present locations where the softening plastic hinges appeared during the analysis. 6.5.3 Push-over of an asymmetric frame In this example we analyze an asymmetric frame presented in Fig. 6.14, where LB1 _ 600 cm, LB2 _ 500 cm, LB3 _ 400 cm and He _ 250 cm. All the other geometrical, material and loading parameters are the same as in the previous example. The results are presented in Figs. 6.15 to 6.16. The total lateral load versus utop curves, where utop is horizontal displacement at the top-left corner of the frame, are presented on the left part of Fig. 6.15. The results of geometrically linear and geometrically nonlinear analyses are nearly the same before the ultimate resistance is reached at 1581.8 kN for geometrically linear case and at 1578.4 kN for geometrically nonlinear case. After that point the difference between those two analyses is bigger. The final computed equilibrium configuration for the geometrically nonlinear case is at utop _ 26.52 cm. In the next load step one additional softening plastic hinge is activated, which results in the global failure mechanism. Since our path-following algorithm is only governed by the the increase of utop, we are unable to capture the remaining part of the load-displacement curve. In the right part of Fig. 6.15 we present the dissipated energy versus utop curves. The shapes of the curves are very similar to those from the symmetric frame case. Namely, first A4Fo "top A4F0 A3F0 > 1 ""— A6F0 Sv A3F0 A2F0 a A4F0 ffi A4F0 ffi A2F0 AF0 WW A2F0 Sv A2F0 Sv AF0 L BL L B2_ L B2_ Hc Hc Hc Hc Figure 6.14: Asymmetric frame: geometry and loading Slika 6.14: NesimetriCni okvir: geometrija in obtežba Z -a 5 iooo o G o X "top [cm] 400 - 3 E^ 300 CD G CD -O CD 10 15 20 25 30 35 Hardening Lin. Softening Lin. Total Lin. Hardening Nonlin. Softening Nonlin. Total Nonlin. 10 15 20 25 30 35 "top [cm] Figure 6.15: Load versus displacement and dissipated energy versus displacement curves for asymmetric frame Slika 6.15: Krivulje obtežba - pomik in disipirana energija - pomik pri nesimetricnem okvirju 200 100 0 there is the elastic non-dissipative phase, followed by the pure hardening dissipation phase, followed by the combined hardening and softening dissipation phase and finally the pure softening dissipation phase. In the geometrically nonlinear case we do not have the final pure softening dissipation phase due to activation of global failure mechanism. On the left part of Fig. 6.16 we present the deformed configuration of the frame at utop = 26.52 cm. Locations, where softening plastic hinges were activated at utop ~ 26.52 cm, are presented in the middle part of Fig. 6.16 for the geometrically linear case and in the right part of the same figure for the geometrically nonlinear case. / J 7 7 Figure 6.16: Asymmetric frame: Deformed shape and locations of softening plastic hinges Slika 6.16: Nesimetricen okvir: deformirana lega in mesta plasticnih clenkov 6.5.4 Bending of beam under constant axial force In this example we compare results of the beam model with results obtained by using the shell finite element model from ABAQUS. For the comparison we choose the problem of the bending of the beam of length Lref under a constant axial force, presented in Section 5.1. For that reason, the geometric and material properties are the same as those in the Section 5.1. The beam model analyses are performed with two sets of material parameters, where the first set (SET1) is given by (6.88) to (6.90). In the second set (SET2) we replace the expression (6.88) with mu = MU6-88) - N |Auy | , (6.91) where Mu is the maximum concentrated moment applied at the end cross-section (point MP, see Fig. 6.6), N is the applied axial force (positive when producing tension) and Auy is the relative displacement in the y direction between the point MP and the position of the local failure. The difference between the applied concentrated moment at the point MP (see Fig. 6.6) and Mu thus arises due to large displacements correction. When the yielding and local buckling of the beam are significant and the displacements in the y direction are no longer negligible, the contribution of the axial force N to the bending moment must be taken into account. In this particular case we have N = -0.1Ny, Mu = 521 kNm and Auy = 0.15 m, which leads to M* = 565 kNm. Five beam finite elements are used to model one-half of the beam under consideration, since the symmetry is taken into account. The symmetry conditions at the symmetry plane are u = w = w' = 0. The load was applied in two steps. In the first step the beam was loaded with compressive axial force N = — 0.1Ny. In the second step the moment M was applied at the free-end of the beam. In order to ensure the proper activation of softening in the geometrically linear analyses the ultimate bending moment of the finite element near the symmetry plane was slightly weakened.

re/ TT7P TT7P,ref p ( AT\ __J ( AT\ T-iW ( AT\ _ T-iW Ew (N) = Ew (N), Ew (N) = Ew (N). By assuming that each cross-section in the frame member of length Lref is approximately under the same force-moment state during the hardening regime we obtained the following expression for hardening modulus ) = (M«2 (N) - My2 (N))Lref h ( ) 2Ewp,ref (N) • Contraty, we assumed that all the softening effects are very local and limited to only one cross-section and the softening modulus is then M2 (N) K (N)| = M,(e/ ) , Ks < 0. 2EW (N) The computational procedure presented in Section 6.4 is split into the local and global phase. In the local phase we determine the updates for the hardening variables related to element integration points by solving the following equation 4iP(Mnp+i(d^),Kn^Pi(Yn+i)),q(š:^+i(Yn+i)))=re^i) = 0. The softening variables are determined by considering that the localization criterion is zero in the softening step ^{e)(t^n+i(a^enn+i(Y^i)),q(l^i(Y^i))) = ?(e)(Yne+i) = 0, where we explicitly use the additional equilibrium equation to compute the bending traction at the discontinuity tgn+i = - / G*(xf,x)M(dntf,<-<+i)dx. In the global phase, where the single element contribution to the system of global equations is K(e) Kfa Khd Kha (i) ( Ad(e)'(i) \ ( fext,(e) _ jint,(e),(i) n+i . I — I J n+i J n+i Aa(e)'(i) I \ 0 n+i \ ' V 0 we determine the updates for the current iterative values of nodal displacements d(e),(i) _ d(e),(i-i) , a d(e),(i-i) dn+i = dn+i + Adn+i • The static condensation allowed us to form the standard form of the element stiffness matrix K(e),(i) = K(e),(i) _ Kfa,(i) (Kha,(i)\-i Khd,(i) K n+i = K n+i K n+i I Kn+i I K n+i , thus the global solution procedure is completely the same as in the standard finite element formulation. In the first numerical example in Section 6.5 we presented the shell model computation and by using its results we determined a bilinear approximation function for M (N) = Mref (N) = I Mef,0(L03 + 0.85 N) if N< -0.035NV Mu (N) = Mu (NMreffl y if n >-0.035Ny ' and we also assumed a constant values for hardening and softening modulus Kh(N) = 1.06 ■ 107 kN/cm2, Ks(N) = -3.28 ■ 105 kN/cm2. These values were used in the rest of numerical simulations. The multi-scale procedure proposed in this chapter belongs to the class of weak coupling methods, where we carry out the sequential computations. The results of the shell model computations, accounting for material and geometric localized instability, are stored to be used within the beam model softening response. As presented by numerical simulations, performance of the proposed multi-scale computational approach is very satisfying. One of its main features is that detection and development of the softening plastic hinges in the frame is fully automatic, and spreads gradually in accordance with stress redistribution in the course of the nonlinear analysis. This is in contrast with many standard computational approaches to the limit load, under the push-over and the full collapse analysis of frames, which rely on predefined locations of plastic hinges and the corresponding inelastic deformations. Chapter 7 Failure analysis of 2D solids 7.1 Introduction In this chapter we present a quadrilateral two-dimensional elastoplastic finite element with embedded strong discontinuity. The effective locking-free design of quadrilateral finite element with embedded discontinuity is much more demanding, e.g. [Linder and Armero, 2007], [Manzoli and Shing, 2006], than that of the constant strain triangle. This might be one of the reasons that the constant strain triangle has been used in majority of works related to embedded discontinuity finite element modeling of failure in two-dimensional solids, e.g. [Sancho et al., 2007], [Ibrahimbegovic and Brancherie, 2003], [Brancherie and Ibrahimbegovic, 2008], [Mosler, 2005], [Jirasek and Zimmermann, 2001]. The derived quadrilateral element does not show any locking problems. Its kinematics can model linear jumps in both normal and tangential displacements along the discontinuity line. We use the derived quadrilateral element to describe tensile fracture process in plain concrete two-dimensional solids and to model the crack growth. The same formulation is also suitable for failure analysis of ductile materials where localized shear bands form. The outline of the chapter is as follows. In Section 7.2, we derive an elastoplastic quadrilateral two-dimensional finite element with embedded discontinuity. In Section 7.3 we present details of the computational procedure. Numerical examples are presented in Section 7.4 and concluding remarks and a short summary are given in Section 7.5. 7.2 Family of ED elements for planar problems In this section we present a family of quadrilateral elements in the two-dimensional setting. The finite elements of this kind can represent the elasto-plastic response, including both hardening and the localized softening effects, the latter being associated with the strong discontinuity in displacements. 7.2.1 Kinematics We consider a quadrilateral finite element occupying domain Qe (see Figure 7.1), divided with the discontinuity line re into two sub-domains, and Qe- (Qe = Qe+ U Qe-). Element's geometry is defined by the bilinear mapping £ ^ xh (£ G [-1,1] x [-1,1]; xh G A® Ux4,Uy4 4 Ux1,Uy1 y Ux3,Uy3 3 re An0,Ani ,AmQ,Ami 4 H 3 1 2 X Figure 7.1: Quadrilateral finite element with embedded discontinuity Slika 7.1: Štirivozliscni koncni element z vgrajeno nezveznostjo Qe) with 4 h X" (I) |ne= £ Na (£) xa, xa =[xa,ya]T , £ =[£,n]T , (7.1) a=1 where xa are coordinates of the finite element node a and 1 a 12 3 4 Na (£) = 4 (1 + U) (1 + VaV) , Ca -11 1 -1 . (7.2) 4 na -1 -111 The superscript h is used to denote the discrete approximation of different fields. In Figure 7.1, uxa and uya are the nodal displacements in x and y direction of the node a. In addition to (standard) nodal degrees of freedom we use strong discontinuity interpolation parameters an0 and an1 for the normal direction and am0 and am1 for the tangential direction, along with the unit normal vector n and the unit tangent vector m (see Figure 7.1). The mid-point of the discontinuity line is denoted with xr, see Figure 7.1. We assume that the domain of the discontinuity influence corresponds to a single element and write the displacement field as: 4 [uhx,uh]T = uh(t, re) Na(^)da + a=1 + Mno($, re)ano + Mni($, re)ani + Mmo(£, re)amo + Mmi(£, re)ami, (7.3) < where da = [uxa,uya]T are the nodal values of Uh, the part of the displacement field arising from standard isoparametric interpolation of nodal displacements. Similarly, u^ is the part which arises with the introduction of strong discontinuity and Mn0(£, re), Mn1(£, re), Mm0(£, re) and Mm1(£, re) are the interpolation matrices related to discontinuity parameters that will be derived below. With four discontinuity parameters we model four modes of element separation along re (see Figure 7.2): 1. "n0" - the constant mode of separation in the normal direction, 2. "n1" - the linear mode of separation in the normal direction, 3. "m0" - the constant mode of separation in the tangential direction and 4. "m1" - the linear mode of separation in the tangential direction. We investigate the motion of two rigid bodies, and Qe+, due to the particular separation mode(Figure 7.2). From (7.3) we can obtain the displacements for a mode, mode G (n0,n1,m0,m1), Umode = Ud ,mode + Ua,mode, Ua,mode = Mmodeamode■ (7.4) From this equation we can then determine the interpolation matrix uh _ uh M "'mode u d,mode r\ mode = -=-■ (7.5) amode By examining Figure 7.2 we can determine 7uU!mode and ul^mode for each mode. By using (7.5) we can then derive the interpolation matrices as follows: An0 1 Normal constant mode n0 Normal linear mode n1 Tangential constant mode m0 Tangential linear mode ml 2 2 2 1 Figure 7.2: Different element separation modes Slika 7.2: Razlicni nacini locevanja koncnega elementa Constant normal separation mode - "nO" nan0 for a G He+ d ra,n0 irrh u n0 irh u d,n0 0 otherwise Hr(x)nan0 ; Hr(x) Nada,n0 1 for x G Qe+ O otherwise aene+ M n0 Un0 Uh,n0 an0 Hr(x) - £ Na n aene+ Linear normal separation mode - "nV 0 1 -1 0 xaan1 for a G otherwise 1 — Hr(x) 01 10 x«n1 U d,n 1 M n1 — E Nada,n 1 aene+ Un 0 Uh,n0 ano — Hr(x) 01 10 - E Na aene+ 01 10 Constant tangential separation mode - "mO" mam0 for a G He+ d 'a,m0 wrh U m0 0 otherwise Hr (x)mam0 —h U d,m0 ada,m0 M m0 Y^ Nad, aene+ Uh — Uh Um0 Ud,m0 am0 Hr(x) - V Na) m aene+ Linear tangential separation mode - " m1" d = j (m • xa) mam1 for a G He+ ■xrh _ 0 otherwise U'm1 = Hr(x) (m • x) mam1 Uh ^ Nada U d,m1 = / , N ada,m1 aene+ M m1 Uh — Uh Um1 Ud,m1 am1 Hr(x)m • x - ^ Nam • xa | m aene+ (7.6) (7.7) (7.8) (7.9) (7.10) (7.11) (7.12) (7.13) (7.14) (7.15) (7.16) (7.17) (7.18) (7.19) (7.20) (7.21) Remark 7.1. Motion of the region in the "ml" mode is not a rigid body motion but a stretch. da,n1 0 a Strains We determine the strain field as the symmetric part of the gradient of the displacement 2 [ T field in (7.3) (e = Vsuh = 2 [Vuh + (Vuh)T]) which can then be written in a vector form e duh duh duh duh ^ "'x _y ^ "'x +__y dx ' Sy ' dy dx £b a^a + Gn0 an0 + Gn1an1 + Gm0am0 + Gm1am1, (7.22) (7.23) a=1 where B„ r dNa dx 0 nx 0 mx 0 0 dNa - dy dNa dy dNa dx , Bn = 0 ny 0 my _ ny nx _ my mx G n0 = B a n + £rB„n, aene+ —— ■T^n' Gn0 G n1 £B a aene+ 01 10 Xa + črB „n£p, Gnl Gn G m0 j: B a m + črBram, -' Gm0 aen Gm G m1 HpBmrn - Ba(m • xa)m + fr Bram£r aene+ (7.24) (7.25) (7.26) (7.27) (7.28) Gml In above expressions we have used the folowing derivation rules (see e.g. [Mosler and Bruhns, 2004] and references therein for more details) Vs (Hrn) Vs (Hrm) čr(x) čr (n ® n)s = črBnn , tensor form vector form (m ® n)s = črBmn tensor form vector form to for x G re 0 otherwise ' (7.29) (7.30) (7.31) where Vs = 1 (V + VT). Note that is a coordinate along re, which has 0 value at xr and is positive in the m direction. We further divide the strain field into a regular part e and a singular part e e = e + e, (7.32) 4 e = Bada + Gn0an0 + Gniani + Gm0am0 + Gmiami, (7.33) a=i e = Gn0 an0 + Gniani + Gm0am0 + Gmi ami (7.34) We note that the singular part is just a particular representation of the localized inelastic deformation introduced at the displacement discontinuity. 7.2.2 Equilibrium equations Interpolation of virtual strains is carried out according to 4 c = £ B aa + G nO an0 + GniOtni + G m0am0 + G mi ami. (7.35) a=i The discontinuity parameters can be viewed as additional (incompatible) degrees of freedom of the element. In that sense we use the following expression to compute the interpolation matrices of virtual strains (see e.g. [Ibrahimbegovic and Wilson, 1991]), which ensures the convergence in the spirit of the patch test C 1 f Gmode Gmode I Gmoded^. (7.36) AQe JQe By introducing the last result in (7.36) into (7.25)-(7.28), we obtain c _ 1 r __/r G no = Gno —— G,nodQ —-— B,n n + 5rB,nn, (7.37) AQe JQe AQe " " —v- — — G„0 G nO G ni = Gni - J Gnidtt + StBMt, (7.38) G G n! G nl . _ 1 f _ 1r GmO = GmO --- Gm"dtt - --BnW + Sr^n^l, (7.39) AQe JQe AQe S-V/-' S-v-' G — GmO G mo (Gmi = Gmi - J Gmidtt + 8rBnjnCr . (7.40) G G m1 Gm1 The weak form of the equilibrium equations or the principle of virtual work for an element e of a chosen finite element mesh with Ne1 finite elements, can be written as: OTint,(e) - ^next,(e) = 0. (7.41) By using (7.35) for the definition of the virtual strains, we can write a single element contribution to the virtual work of internal forces as: 5Wnt>(e) = t(e) eTadtt Jne 4 r T = S t(e) / da BTadQ + a=1 Jne standard T T T T ^n1 t(e) / a^G^a + aniG^a + a^G^a + amiGmiadtt, (7.42) Jne V additional where t(e) is the thickness of the element and a =[ax,ay ,aXy ]T (7.43) is the stress vector. From the term "standard" in (7.42) we obtain the vector of element internal nodal forces fint,(e) = If ant'(e)TlT , f Tt,(e) = t(e) I BTadQ. (7.44) ^ ^ Jne From the virtual work of external forces $next'(e) we can get the vector of element external nodal forces fext,(e), representing the consistent external load applied to the element's nodes. We note that the displacement discontinuity parameters do not contribute to the external load vector. The finite element assembly of vectors f mt'(e and fext,(e), for all elements of the chosen mesh, leads to a set of global equilibrium equations A=i (f mt'(e) - f ext'(e)) = 0, (7.45) where A is the assembly operator. Note that we have only used one part of the right side of equation (7.42) in (7.41) to get the set of global equations (7.45). The other term in (7.42), denoted as "additional" (resulting from additional enriched kinematics due to embedded discontinuity), will also contribute to the weak form of the equilibrium. However, we will treat this contribution locally, element by element, only for those elements where the discontinuity was activated. Then, in view of (7.41), the following four equations are obtained for each element of the chosen mesh h (e) n0 t(e) II G^vdn + [ nTBtv dr 'ne ——T — tn t(e) I Gn0vdtt + t(e) I tndr = 0, (7.46) h (e) n1 h(e) hm0 hne hn0 hre hn0 t(e) ^J G Tn1adQ + ^ Cr rn^B^ dr j 'ne -T-T t(e) Gn1adtt + t(e) Crtndr = 0, hn 1 hre hn1 t(e) lj G^vdQ + J mFlBna dr j V Q r — tm ) -t-T t(e) Gm0vdn + t(e) tmdr = o, hre hm0 (7.47) (7.48) h(e) hm1 t(e) \ Gmladn + / Cr mTBTn v dr —tm —T t(e) G m1vdn + t(e) Crtmdr = o, hm 1 hre hm1 (7.49) where tn and tm represent the normal and tangential components of the traction at the discontinuity t [tn, tm] T (7.50) and fQe ^r(o)dn = fre(o)dr was used. We gather the equations (7.46) to (7.49) in a vector form to write the element-based residual as: h (e) = hfi' h"£ + hre hne hn0 r hre 1 hn0 hne hn1 + hre hn1 hQe hm0 hre hm0 hne hm1 hre hm1 (7.51) n e n r e n e m0 n r e e n r e 0 The dedicated solution procedure for a set of global equations in (7.45) together with a set of local (element) equations in (7.51) will be further addressed in Section 7.3. 7.2.3 Constitutive relations In what follows we choose plane stress constitutive relations, although the derived 2D elements with embeded discontinuity could be also used for plane strain problems. We will describe the constitutive behavior of the bulk material \ re with elastoplastic model with isotropic hardening, while the response of the discontinuity re is considered to be rigid-plastic with softening. The basic ingredients of the chosen constitutive relations are built upon the classical plasticity (e.g. [Ibrahimbegovic et al., 1998]) and can be summarized as: • The regular strains e (7.32) can be additively decomposed into elastic part ee and plastic part e^ e = e + e*\ (7.52) • The free energy of the bulk material is assumed to be the sum of the strain energy function W and the hardening potential S 1 tf(r,C) := W (e) + -(£) = -eeTce + s®, where c E 1- v2 1 v 0 v 1 0 0 0 1 — V 2 (7.53) (7.54) E is the Young's modulus of the material, v is the Poisson's ratio of the material and C > 0 is strain-like hardening variable. The yield function is assumed to be based upon the von Mises yield criterion (see [Dujc and Brank, 2008] and references therein for more details), which can be written as 0 (^,q) = ACT - ( 1 -— ) < 0 - \ 2 q a, where A 2a2 2 -1 0 -1 2 0 0 0 6 (7.55) (7.56) whereas q is the stress-like hardening variable conjugate to the equivalent plastic strain C and ay is the uniaxial yield stress. The traction t at a point at the discontinuity re is related to jump in displacements particular point t = t(u). (7.57) u = [un,um]T (see Figure 7.3) at this particular point 1 Figure 7.3: Neighborhood of a point of interest at the discontinuity re Slika 7.3: Okolica obravnavane toCke na nezveznosti • The localization (failure) criterion that activates localized softening at discontinuity re is defined in terms of the traction t and the stress-like softening variable q(1) (the later is defined in terms of the strain-like softening variable 1). % = %(t,q) < 0. (7.58) • The softening potential 1(1), (7.59) defined in terms of the strain-like softening variable 1. The remaining ingredients of the elastoplasticity with hardening,describing bulk material \re, can be obtained from the standard consideration of thermodynamics for associative plasticity based upon the principle of maximum plastic dissipation (e.g. see [Lubliner, 1990], [Simo and Kennedy, 1992], [Ibrahimbegovic, 2009]). By using (7.52) and (7.53) the mechanical dissipation at a point x G \ re can be further written as 0 < D d= aTe - = (a - g)Tf + aTlf - ^, (7.60) where (o) = d (o) /dt denotes the derivative with respect to pseudo-time t G [0,T]. By assuming that the elastic process is non-dissipative (i.e. D = 0), with no change in plastic state variables, we obtain from (7.60) that d^ a = ^ = <7-61> We can also define the stress-like hardening variable q by further considering (7.60) and (7.53) as: 3W dS ^ , q =--=■ =--=■. (7.62) dC dC ( ) By replacing (7.61) and (7.62) in (7.60), the plastic dissipation can be obtained as D = + qC. (7.63) The principle of maximum plastic dissipation states that among all the variables (ct, q) that satisfy the yield criteria 0 (ct, q) < 0, one should choose those that maximize plastic dissipation (at frozen rates eP and C). This can be written as a constrained optimization problem: min max [LP(CT,q, 7) = —D(CT,q) + 70(CT,q)] , (7.64) where 7 > 0 plays the role of Lagrange multiplier. By using (7.63) and (7.55), the last result can provide the evolution equations for internal variables that can be written as: ^ = — ep + 7= 0 f = 72ACT, (7.65) d ct OCT = 4 + # = 0 C = 7-^1 — ^ (7=5) Yr-VCTTACT, (7.66) dq dq ay \ ayJ along with the Kuhn-Tucker loading/unloading conditions and the consistency condition Y > 0, 0 < 0, Y0 = 0, Y0 =0. (7.67) To obtain the remaining ingredients of the rigid-plastic response, describing softening at the discontinuity re, we isolate one point at the discontinuity. We first define the strain energy function due to softening potential in this point as W = S. The dissipation can then be written as: _ 0 < D =f tTU — f(|) = tTU — Zf. (7.68) dC where t is the traction defined in (7.57). Finally, by introducing the traction-like variable - dW dS q = — dW = — ^S, (7.69) dC dC the result in (7.68) can be rewritten as D = D = tTu + qC. (7.70) The principle of maximum plastic dissipation at the rigid-plastic discontinuity can then be defined as: min max t,q y ~L(t,qrf) = -D (t,q) + 70(t,q) , (7.71) where y > 0 is the Lagrange multiplier. By using (7.70) and (7.58), we get from (7.71) above the following evolution equations of displacement dicontinuity: dL — —06 — —06 _ . = -1 + Yf = 0^ i = (7.73) Oq Oq Oq accompanied by the Kuhn-Tucker loading/unloading conditions, and the consistency condition Y > 0, 6 < 0, =0, Y6 = 0. (7.74) 7.3 Computational procedure In this section we present the details of a computational procedure for solving the set of global equilibrium equations along with the local (element based) equations, generated by the finite element with embedded discontinuity as presented in Section 7.2. Figure 7.4: Numerical integration scheme Slika 7.4: Numericna integracijska shema The solution of the set of global nonlinear equations in (7.45) and the local nonlinear equations in (7.51) is obtained at discrete pseudo-time values 0,r1,r2,..., rn-1,rn,rn+1, ... ,T by means of the incremental-iterative scheme. We will consider a single-step scheme providing the solution in a typical pseudo-time increment from Tn to rn+1. Let us assume that all the variables, related to an element e, its bulk integration points bip = 1, 2, 3, 4 and its discontinuity integration points dip =1, 2 are given at Tn, i.e. given: d^ = felT , ^, tn , otf, iT and xrs, xrE, (7.75) where a (e) (e) (e) (e) (e) lT a a a a Lxn0,n> Lxn1,n um0,n umi,n Remark 7.2. 4-point Gauss integration scheme is used for the bulk and 2-point Gauss integration scheme is used for the discontinuity, see Figure 7-4- We also need as indicated in (7.75): (i) the starting point of the discontinuity xrs (which is only true if the discontinuity in one of the neighboring elements ends at the common edge) and (ii) the end point of the discontinuity xrE (which is only true if softening has been activated so far). We will then iterate in the pseudo-time step in order to compute the converged values of the variables at Tn+1, i.e. find: ^ fn+ ^ tnh and (7.76) if not given already and if the softening conditions are met: xrE. The computation of solution (7.76) is split into two phases: (a) The global (mesh related) phase computes the current iterative values (with (i) as the iteration counter) of nodal generalized displacements at Tn+1 while keeping the other variables fixed, i.e. global phase: d^f = dgf-1* + AdJ+f-^. (7.77) The computation of iterative update Ad^!-1'' will be explained further below. (b) The local (element and integration point related) phase computes the values of vari- ables epn+p, iniJ+1, ane+1, in+1 while keeping fixed. The computation procedure depends on weather the softening has been activated in the considered element or not. Therefore, the local computation procedure on the level of a single element can be based either on hardening plasticity procedure or on softening plasticity procedure (excluding each other). In the rest of this section we will first describe in detail the phase (b). This will be followed by the description of the phase (a). The softening plasticity procedure is carried out only in those finite elements where: n • discontinuity has been active at the previous time step Tn, i.e. and are provided and the local softening variables changed, or • discontinuity has not been active in the previous time step but the discontinuity in one of the neighboring elements ends at the common edge, i.e. is provided. In the latter case we first compute the direction of the discontinuity by using the average stress field in the element n n I CTavg (d&1, e^) m m ct avg (d (e) eP 'biP n+1, en (7.78) We determine the direction of the discontinuity in the direction of principal stress if the governing mode of separation is mode I, i.e. n is parallel to the maximum principal stress. If the governing mode of separation is mode II, we set m to be parallel to one of the directions of maximum shear stress. Once the direction of the discontinuity is determined we are able to obtain the end point of the discontinuity xre and the complete description of the finite element with embedded discontinuity geometry (He-, and lr). The main part of the softening plasticity procedure starts with the determination of the trial values of the traction in discontinuity integration points. Since we use the rigid plasticity law at the discontinuity, we are unable to determine the trial values by using the constitutive law. We use instead the local equilibrium equations in (7.51). Since stress in the bulk integration points is well defined we can always obtain the bulk contribution h nutria,! 'n+1 nutria,! 'n0,n+1 n£'tria! 'n1,n+1 nutria,! m0,n+1 nutria,! m1,n+1 CTrfj, en,bip an) (7.79) We then further consider, that the integration at the discontinuity is carried out according to f(Cr)dC £ f (C?p) dip =1 w i^r Cdi 2 , Cr ± lr 2^' wdip = 1, (7.80) where f (Cr) in an arbitrary scalar function, f (Crip) is the value of this function evaluated at the Gauss quadrature point dip =1, 2, Cpip is the coordinate of the quadrature point and wdip is the corresponding weight. By using (7.80), the contribution of the discontinuity to equation (7.51) can then be written as t(e)lr t1 +12 nn Cr tn + Crt a 2 r^n 2 t1 +1 mm Cr tm + Cr tm (7.81) e r 2 r 2 e where t1 = [tn,t]J\[ is the traction in the first discontinuity integration point and t2 [tn,tm]T is the traction in the second discontinuity point. To satisfy the local equilibrium equations we solve the system of four algebraic equations iQe,trial , i re/» 1,trial ,2,trial\ _ pi hn+1 + h (tn+1 , tn+1 ) = 0, (7.82) and determine the four unknown components of the traction in the discontinuity integration points 1,trial tn+1 2,trial tn+1 1,trial 'n,n+1 1,trial m,n+1 t 2,trial t n,n+1 2,trial m,n+1 q /h Qe ,trial 2(hn1,n + 1 hHe ,trial «2 \ ~hn0,n+1 «r) = e («2-Q/i H ,trial -«r )lr j He ,trial j.2 \ ~hm0,n+1 «r) L («2- -«r )lr J rs/1 He ,trial 2(hn1,n + 1 j He ,trial £ 1 \ ~hn0,n+1 «r) = («1-q/7 H ,trial -«2 )lr hHe ,trial «1 \ ~hm0,n+1 «r) L («1- -«2 )lr J (7.83) (7.84) Next we provide the trial values of the failure function at the discontinuity integration points =1,trial = ... , _ =1 _ 0 = 0(tnir1ial,q(£n)) < 0 or 0 ' = (t%r ,q(U)) < 0 . (7.85) =2,trial = ,2,trial =/=2 ? If the trial failure criterion (7.85) is satisfied, i.e. if the yield function is less or equal to zero in at least one discontinuity integration point, the values of softening plasticity local variables remain unchanged =1,trial =2,trial 0 < 0 or 0 < 0 a (e) n+1 a (e) =1 e n+1 =1 =2 en, Cn+1 =2 en' (7.86) In the case of violation of the trial yield criterion (7.85), i.e. the yield function is greater than zero in both discontinuity integration points, the values of discontinuity integration point variables are updated by backward Euler integration scheme r' =1 = =1 + = 1 d0 Un+1 = =n + Y n+1 ^ 1 =2 _ =2 =2 00 Un+1 = Un + T n+1 ~Qt |«r : e 1 = e 1 + =1 en+1 = en + 'n+1 q |«1, f _f ,^2 en+1 = en + 7n+1 q |«2, (7.87) (7.88) where y n+1 _- _2 _- Yn+1(Tn+1 — Tn) and Yn+1 = Yn+1(Tn+1 — Tn). Note that the relation between the jumps in displacements evaluated at the discontinuity integration points and the kinematic parameters related to discontinuity can be obtained as follows a (e)(=1, =2) =1 + =2 =1 =2 =1 + =2 =1 =2 Un + U,n U,n U,n Um + Um Um Um er — er er er T (7.89) n 2 2 _1 r_1 _1 1T —2 r ° o 1 T where u = u„,um and u =2 =2 . By using (7.87)-(7.89), we can thus determine the updated values of discontinuity kinematic parameters (e). = a(e) (U1 (Y an+1 = a _1 _1 _2 _2 \ Un+1 (y n+1 ), Un+1 (Yn+1 ^ . (7.90) The values of the plastic multipliers Yn+1 and Yn+1 are determined with an iteration solution to the system of the following two equations =1 'V , (a(e) fen ^ 1 (= 1 i!(=n ^ \ / ( tn+1 (y n+1,7n+0) ^ (^n+1 (Y1+1))) = / (7n+1,7n+^=o, (7.91) 0 (tn+1 (4+ (yn+1,Yn+0) ^ (L+1 (jn+1))) =0 (Yn+1,Yn+1)= _1 _2 New values of Un+1, Un+1, 1n+1 and 1n+1 are then computed by (7.87), (7.88) and the corresponding is computed by (7.89). Note, that we compute the traction in (7.91)-(7.92) in the same manner as in the case of trial values, i.e. by solving the system of four algebraic equations hne (adh, enbip, at))+hre(tn+1,tn+1) = 0 ^ tn+1, tn+1. (7.93) The end result of the above described softening plasticity procedure are the new values of parameters a[;e+1, which influence the stress state of the whole element by giving the new values of stress in the bulk integration points as * jn+1 ' jv-hd,(i) _ (dhW\ (i) j>^ha,(i) _ ^ (i) ^n+1 UdM/ n+^ n+1 Vd«(e)/ n+1 (7.104) The static condensation in (7.103) above allows us to form the standard form of the -~-(e),(i) element stiffness matrix Kn+1 that contributes to the finite element assembly 4=1 (K^f Ad^) _ (/nx+1(6) — /nn+f,(i)), (7.105) where ^ntf = K^ - K n+f -1 K (7.106) Solution of (7.105) gives the values of iterative update AdJH^, which should be performed as indicated in (7.77). 7.4 Examples In this section we provide the results of a number of numerical simulations that can illustrate a very satisfying performance of the proposed finite element. The code was generated by using symbolic manipulation code AceGen and the examples were computed by using finite element program AceFem, see [Korelc, 2007b] and [Korelc, 2007a]. 7.4.1 Tension test In this example we consider a square block of 20 cm x 20 cm x 0.1 cm subjected to uniaxial tension, see Figure 7.6. The block is made of material with Young's modulus E = 3000 kN/cm2, Poisson's ratio v = 0.2 and the ultimate tensile strength ou = 3 kN/cm2. The softening response is governed by the cohesive law at the discontinuity presented in Figure 7.5. We can write the law in Figure 7.5 in terms of failure criterion Ln (u Ur Figure 7.5: Rigid-plastic cohesive law with linear softening in tension Slika 7.5: Kohezijski zakon (toga plasticnost z linearnim mehcanjem v nategu) $(t, q) = tn - - q) < 0, and the bilinear softening response q = min On, -Ks1 (7.107) (7.108) where Ks = -45 kN/cm3 is the linear softening modulus and 1 = un, see (7.72), (7.73) and (7.107). We model the block with one finite element, which is supported at the left side I Or ux,Rx S o o OJ ux,Rx 20 cm Figure 7.6: Tension test on a square block Slika 7.6: Natezni preizkus kvadratnega bloka and pulled, by imposing displacements, at the right side (see Figure 7.6). Once the tensile strength of the material is reached the discontinuity appears in the direction perpendicular to the maximum principal stress. Note that the behavior of the discontinuity is only defined for Mode I, i.e. the equations for the traction in the tangential direction and the corresponding jumps are not considered in the simulation. In Figure 7.7 we present the ux [cm] Figure 7.7: Reaction force versus imposed displacement curves Slika 7.7: Krivulje reakcija - vsiljen pomik reaction force versus imposed displacement diagram. Note that the results of the present formulation are in complete agreement with the results obtained in [Linder and Armero, 2007] for nonlinear elastic cohesive element, apart the local unloading branch that our rigid-plastic formulation can capture. 7.4.2 Bending test Next we consider the bending test of the block with the same material and geometrical properties as in the previous example. In Figure 7.8 we present the problem definition. 20 cm_^ Pseudo-time Figure 7.8: Geometry of the bending test of the square block and the imposed displacement versus pseudo-time curves Slika 7.8: Geometrija pri upogibu kvadratnega bloka ter krivulje vsiljen pomik - psevdo cas In the first part of loading the displacement at the top uxT and the displacement at the bottom uxB are applied with the same rate. The tensile strength of the material is reached at uxT = uxB = 0.001 cm and at that point the discontinuity in Mode I appears. This is followed by a non-uniform regime of loading with the rate of imposed displacement at the bottom being twice the rate of imposed displacement at the top, see right hand side of Figure 7.8. In Figure 7.9 we present the results of our simulation along with the results obtained in [Linder and Armero, 2007]. The differences in results are hardly noticeable up until the point when the first fibers of the discontinuity fail completely, resulting with tn = 0. After that the results no longer coincide, since the integration scheme in [Linder and Armero, 2007] considered five integration points along the discontinuity and can therefore represent a smoother transition from the softening regime to the complete failure. 7.4.3 Partial tension test In the above numerical examples we have evaluated the performance of the present formulation for the normal opening response. To evaluate the performance of the present UxT [cm] Figure 7.9: Top reaction force versus imposed top displacement curves Slika 7.9: Krivulje reakcija na vrhu - pomik na vrhu formulation for the tangential response, i.e. the Mode II response with linear softening at the discontinuity in its tangential direction, we consider the partial tension test presented in Figure 7.10. The block is of the same dimensions as in the previous sections, with the same value of Young's modulus E = 3000 kN/cm2 but with a Poisson ratio v = 0.0. We consider a pre-existing discontinuity at the center of the block and from the beginning there is no resistance in the discontinuity, i.e. the discontinuity provides no stiffness in the tangential direction and the traction is always zero (tm = 0). The equations involving the normal response are simply left out in this simulation. Again we model the block by using only one finite element. The block is supported at the two nodes on the left hand side and we pull apart the two nodes on the right hand side, thus causing a linear displacement distribution along the height of the Q+ region. In this way the only stress that develops is which is only limited to Q+ region and all the other stress components remain zero throughout the test (axx = = = 0). In Figure 7.11 we present the stresses that develop in the element with the increase of imposed displacements. Again we compare the results of the present work with the results obtained in [Linder and Armero, 2007]. Note that the results presented in Figure 7.11 are in complete agreement. The stress component in the Q- region is unaffected by the imposed displacements and is always equal to zero, while the stress component in the Q+ changes with the increase of imposed displacement according to 0+ = . 7.4.4 Three point bending test In this example we consider a classical benchmark problem of a notched concrete beam under three point bending. In Figure 7.12 we present the geometry of the specimen, a 20 cm Figure 7.10: Partial tension test of the square block Slika 7.10: Delni natezni preizkus kvadratnega bloka u u 2 1.5 ^ 1.0 iT J- present J» present - J-, Armero ---Jjy Armero 0.000 0.002 0.004 0.006 0.008 0.010 Uy [cm] Figure 7.11: Stress at integration points in and Q- region versus imposed displacement curves Slika 7.11: Krivulje napetosti v integracijskih tockah na območjih in v odvisnosti od vsiljenega pomika Ry,U yy 2 cm 10 cm A 99 cm 99 cm A 20 cm Figure 7.12: Three point bending test of a notched concrete beam Slika 7.12: TritoCkovni upogib zarezanega betonskega nosilca 200 cm x 20 cm x 5 cm simply supported concrete beam with a 2 cm x 10 cm x 5 cm notch placed at the bottom of the beam. The beam is loaded by downward displacement imposed at the top in the center. The beam is made of material with Young's modulus E = 3000 kN/cm2, Poisson's ratio v = 0.2 and the ultimate tensile strength au = 3.33 kN/cm2. The softening response is governed by the rigid-plastic cohesive law at the discontinuity presented in Figure 7.13. The law in Figure 7.13 can be also written in terms of failure Figure 7.13: Rigid-plastic cohesive law with exponential softening in tension Slika 7.13: Kohezijski zakon: toga plasticnost z eksponentnim mehcanjem v nategu criterion (7.107) and the exponential softening law _ / _l£u\ q = tfu ( 1 - exp Gf I , (7.109) where Gf = 0.124 • 10-2 kN/cm is the fracture energy. The response of the discontinuity in the tangential direction is not considered in the simulations. In Figure 7.14 we present the two different finite element meshes that were used in simulations. The coarser mesh is made of 530 finite elements and the finer one of 2186 finite elements. On the left hand side of Figure 7.15 we plot the reaction versus imposed displacement diagrams computed for both meshes. The discontinuity starts at the notch when the tensile strength of Figure 7.14: Coarse (top) and fine (bottom) finite element meshes for the three point bending test Slika 7.14: Groba (zgoraj) in fina (spodaj) mreZa koncnih elementov pri tritockovnem upogibu Figure 7.15: Reaction force versus imposed displacement curves and scaled (100 times) deformed meshes Slika 7.15: Krivulje reakcija v odvisnosti od vsiljenega pomika ter povecana (100 krat) deformirana mreza the material is reached and propagates in the direction perpendicular to the maximum principal stress, i.e. in the Mode I fashion. We have encountered a problem when using the above criterion to determine the discontinuity direction, namely the direction of the maximum principal stress at some point suddenly changes for 90 degrees. This causes a problem in convergence in the simulation with the fine mesh and a non-physical response when using the coarse mesh, see left hand side of Figure 7.15. The discontinuity direction problem was also reported in [Mosler and Meschke, 2003] and we direct the reader therein for further discussion. To obtain the solution without the discontinuity direction problem we considered a predetermined direction of the discontinuity, i.e. discontinuity can only propagate perpendicular to the length of the beam. With this modification we were able to obtain with both meshes the results that are within the experimentally established bounds of [Petersson, 1981] and [Rots et al., 1985]. The results of all simulations are given in Figure 7.15 (left). In the center and right hand side of Figure 7.15 we present the deformed configuration (scaled 100 times) of the area near the notch for both course and fine mesh. 7.4.5 Four point bending test In this example we study the four point bending test on a beam with a notch. In Figure 20.3 cm 39.7 cm 12.2 cm 39.7 cm 20.3 cm c CO 0. 3 65.4 cm 65.4 cm Figure 7.16: Four point bending test Slika 7.16: Stiritockovni upogibni test 7.16 we present the specimen geometry along with loading conditions and supports. The specimen is made of material with Young's modulus E = 2880 kN/cm2, Poisson's ratio v = 0.18 and the ultimate tensile strength au = 2.8 kN/cm2. The behavior of the discontinuity is governed by the failure criterion (7.107) and the softening law (7.108), with softening modulus being Ks = -39.2 kN/cm3. In Figure 7.17 we present the mesh Figure 7.17: Finite element mesh for the four point bending test Slika 7.17: Mreza koncnih elementov pri stiritockovnem upogibu that we used in simulations. With respect to the description of the displacements jumps along the discontinuity line we considered three cases: (i) "n0 + m0" - the constant jump in displacements in both normal and tangential direction, (ii) "n0 + n1" - linear jump in displacement in normal direction only and (iii) "n0" - constant jump in displacements in normal direction only. In the mixed mode case ("n0 + m0") we considered a reduced shear stiffnes for the tangential response according to relation 4>(tm) = tm - kmUm = 0, (7.110) where km = 2.88 kN/cm3. The results of the simulations are presented in Figures 7.18 and Figure 7.18: Load versus crack mouth sliding displacement curves and the corresponding crack paths Slika 7.18: Krivulje obtežba v odvisnosti od relativnega zamika na ustju razpoke 7.19. On the left hand side of Figure 7.18 we plot the applied load - crack mouth sliding displacements curves of our simulations along with the envelope of results the we adopted from [Linder and Armero, 2007]. We can see that all the proposed formulations give a good prediction of the limit load of the structure, while only the mixed mode formulation can capture the true softening response of the structure. On the right hand side of Figure 7.18 we plot the crack paths that correspond to curves on the left hand side of the same figure. Figure 7.19: Scaled deformed mesh of the "n0 + m0" (left), "n0 + n1" (middle) and "n0" formulation Slika 7.19: Povecane deformirane mreze za "n0 + m0" formulacijo (levo), "n0 + n1" formulacijo (sredina) in "n0" formulacijo (desno) In Figure 7.19 we present the deformed (scaled 200 times) mesh of the area near the notch for all formulations. We claim that the crack paths presented in Figure 7.18 and the deformed meshes presented in 7.19 are in good agreement with those from [Linder and Armero, 2007]. 7.4.6 Delamination We consider a delamination test shown in Figure 7.20 as presented in [Manzoli and Shing, 2006]. The properties of the material are: Young's modulus E = 50 kN/cm2 and Poisson ratio v = 0.3. We model the interface as an embedded discontinuity, whose properties are determined with the failure criterion (7.107), the ultimate tension stress = 10-1 kN/cm2, and the exponential softening law (7.109), with the fracture energy G/ = 5 ■ 10-3 kN/cm. The simulations are performed by using a coarse and a fine-type 0.05 cm h Ry,Uy interface Ry,Uy 0.45 i Figure 7.20: Delamination test data Slika 7.20: Delaminacija: podatki Figure 7.21: Coarse (top) and fine (bottom) finite element meshes for the delamination test Slika 7.21: Groba (zgoraj) in fina (spodaj) mreža konCnih elementov finite element mesh, as presented in Figure 7.21. The reaction force versus imposed dis- Uy [cm] Figure 7.22: Reaction force versus imposed displacement diagram Slika 7.22: Sila v odvisnosti od vsiljenega pomika placement diagrams are presented in Figure 7.22. One can see, that the results of both coarse and fine mesh are in good agreement with the results obtained in [Manzoli and Shing, 2006] by using a fine mesh and a Q4SH-NU type finite element. Figure 7.23 depicts the deformed meshes that correspond to the imposed displacement uy = 0.2 cm. Note that the deformed meshes are not scaled and one should for a more realistic representation use a geometrically non-linear framework, which is out of the scope of this work. 7.4.7 Elasto-plastic tension test In the last example we consider a tension test of a metal strip. The geometry of the strip is presented in Figure 7.24 and the thickness is 0.055 cm. One of the shorter edges is built-in and the opposite edge is pulled by imposing the displacements as depicted in Figure 7.24. In this example we consider the bulk material as elastoplastic with the following properties: Young's modulus E = 21000 kN/cm2, Poisson's ratio v = 0.3, yield stress ay = 40 kN/cm2, ultimate stress au = 42 kN/cm2, and hardening modulus Kh = 1000 kN/cm2 where we compute the stress-like variable related to isotropic hardening as q = -Khl. (7.111) The softening response is activated once the ultimate shear stress tu = qu = 21 kN/cm2 is reached and then the Mode II propagation of the discontinuity starts. The response of the discontinuity is governed by the cohesive law depicted in Figure 7.25, which one can also represent in the form of failure criterion (t,q) = |tm| - {tu - q) < 0, (7.112) Figure 7.23: Deformed configuration of the coarse (left) and fine (right) mesh Slika 7.23: Deformirana konfiguracija grobe (levo) in fine (desno) mreze S o LO 3.1 cm 7.5 cm <11 <1 Uy,Ry "Čil Figure 7.24: Tension test of a metal strip Slika 7.24: Natezni preizkus metalnega traka tM' Tu k Um X -Tu Figure 7.25: Rigid-plastic cohesive law with linear softening Slika 7.25: Kohezijski zakon: toga plasticnost z linearnim mehcanjem and the stress like softening variable q = min tu, -KSC , (7.113) where Ks = -400 kN/cm3 and £ = |Uun|, see (7.72), (7.73) and (7.113). In our simulations we only considered constant jumps in the tangential direction. We assumed that there is a small imperfection in the metal strip which one can interpret as the starting point of the discontinuity (see Figure 7.24). Several simulations were made with different mesh Figure 7.26: Total reaction force versus imposed displacement curves Slika 7.26: Reakcija v odvisnosti od vsiljenega pomika sizes ranging from 6 to 384 finite elements. The sum of reaction forces at the right edge versus imposed displacement diagrams are presented in Figure 7.26. One can see that the mesh size has very little influence on the results in this particular example. All the curves have three distinguished phases, namely the linear elastic phase, which is followed by the isotropic hardening phase and the final softening phase. Figure 7.27 depicts the Figure 7.27: Discontinuity paths for several mesh sizes Slika 7.27: Lega nezveznosti pri razlicnih mrezah discontinuity paths for different mesh sizes. In Figure 7.28 we present the (scaled 10 times) deformed configurations for different mesh sizes. Figure 7.28: Scaled deformed configurations Slika 7.28: PoveCane deformirane konfiguracije 7.5 Concluding remarks and chapter summary In this chapter we presented the quadrilateral two-dimensional elastoplastic finite element with embedded strong discontinuity. With the introduction of four additional parameters related to discontinuity line into the isoparametric 2D solid finite element we obtained the enriched displacement field 4 r) = Na($)da +Y, Mmoded, re)«mode, mode G (n0,n1 ,m0,m1), a= 1 mode -..-' - where we asociate the aditional parameters with the following modes of separation • "n0" - the constant mode of separation in the normal direction, • "n1" - the linear mode of separation in the normal direction, • "m0" - the constant mode of separation in the tangential direction and • "m1" - the linear mode of separation in the tangential direction. By investigating the rigid body motion of the domains Qe- and separated with the discontinuity line re we determined the interpolation matrices for all additional modes as uh _ uh _ ""mode ""d,mode mode H • amode The enriched strain field is determined as the symmetric part of the gradient of the displacement field 4 e = Bada + Gnoano + Gniani + Gm0am0 + Gm1am1 • a=1 The virtual strains were determined according to 4 e = Bada + Gn0an0 + Gn1an1 + Gm0am0 + Gmlam1, a=1 where we used the following expression Gmode Gmode 1 Gmode •> AQe J Qe in order to ensure the convergence in the spirit of the patch test. We introduced the virtual strains into the principle of virtual work and obtained the virtual work of internal forces as t(e) eTadtt V t(e) I daTbT&dn + a=1 jQe t(e) / an0Gn0a + an1 Gn1a + am0Gm0a + am1Gm1adQ v standard T T T T V additional and from the term "additional" we obtained four additional equations per finite element h(e) = hQe + hre = + ^ = o, r hQe i hn0 r hre i hn0 hQe hn1 + hn1 hQe hm0 hm0 hQe ' m1 hre hm1 which ensure that all the components of the traction at the discontinuity are in equilibrium with the stress field in the element. The basic ingredients of the chosen constitutive relations for the bulk material and the discontinuity are the usual additive decomposition of regular strains into elastic and plastic part e = ee + ep, the free energy of the bulk material considering isotropic hardening ^,0:= W (ee) + Š® = 1 eeT Cee + Š®, the yield criterion for the bulk material 0 (&,q) = ACT - ( 1 ) < 0 - \ 2 q o,. Q e Q e the cohesive law at the discontinuity written in terms of jump in displacements u = t = t(u), which one can also write in terms of localization criterion % = %(t,q) < 0, and the softening potential The remaining ingredients of the bulk elastoplasticity were obtained by computing the plastic dissipation -FTP T-P I -T D = e + q£, and by considering the principle of maximum plastic dissipation min max [ZP(CT,q, 7) = —D(^,q) + 7%(^,q)l . Similarly, we also defined the plastic dissipation at the discontinuity D = D = tTu+qf, and considered the principle of maximum plastic dissipation min max L(t,q,7) = —D(t,q) + 7 %(t,q) which is here limited to discontinuity line only. The computational procedure presented in Section 7.3 is split into the local and global phase. In the local phase we, if not already given, provide the discontinuity geometry n = n (dne)!,)) , m = m (>vg)) , xrE, by considering the average stress field in the element. The updates of the softening variables and the jumps in displacements at the discontinuity line integration points are determined by simultaneously solving the following equations % ^tn+1 (aie-+1 (jn+1)^ = % (^+1^+1)=0 % ^tn+1 (aie-+1 (fn+1,Yn+^) ,q(L+1 (jn+1)^ = % (jn+1,Yn+1)= where we considered the following relation between the jumps in displacements evaluated at the discontinuity integration points and the kinematic parameters related to discontinuity a (e)(U1, U2) =1 + =2 =1 =2 =1 + =2 =1 =2 U„ + U„ U„ — U„ Um + U m Um Um 2 'er - er 2 er - er T and we explicitly used the additional equilibrium equation to compute the traction at the integration points of the discontinuity h (CT(dn+1> ^n P, 0^1^ + h (tn+1, tn+1) 0 t t2 n+1, tn+1' In the local phase we also determine the updates for the hardening variables related to element integration points by solving the following equation 0>n+1(de]f, e^nr1)),q(e^+1(yX 1))) = r (Yn+1) = 0. In the global phase, where the single element contribution to the system of global equations bip,_bip is K(e) Kfa K hd Kha (i) Ad(e)'(i) Adn+1 Aa (e),(i) f e J n+1 ) n+1 \ v"ln+1 we determine the updates for the current iterative values of nodal displacements pext,(e) oint,(e),(i) — J n+1 0 d(e),(i) _ d(e),(i-1) + A .(e),(i-1) dn+1 _ dn+1 + Adn+1 • The static condensation allowed us to form the standard form of the element stiffness matrix K (e),(i) _ K(e),(i) _ Kfa,(i) n+1 K n+1 K n+1 ha,(i) Kn+1 1 hd,(i) K n+1 , thus the global solution procedure is completely the same as in the standard finite element formulation. A finite element with embedded strong discontinuity has been presented and used to model the fracture process in two-dimensional concrete solids, delamination of composite materials and failure of ductile materials. The element has linear interpolations of the displacements jumps (in both normal and tangential directions), which are important for its locking-free response. In order to make the discontinuity propagation algorithm more robust the continuity of the discontinuity line between the elements has been enforced. e Chapter 8 Conclusion The aim of this work is to increase the understanding of structural behavior at the limit loads and at loads that cause structural failure through the development of numerical methods like the finite element method with embedded discontinuity. In the thesis the emphasis was on the treatment of material nonlinearities and localized failure of material. In the research work presented in the thesis we reached the following conclusions: • We have revisited [Ibrahimbegovic et al., 1992], [Ibrahimbegovic and Frey, 1993b], [Ibrahimbegovic and Frey, 1994] and programmed a nonlinear elastic stress-resultant plate finite element for limit load analysis of reinforced concrete plates with Eurocode 2 [Eurocode 2, 2004] recommendations to describe the constitutive behavior of reinforced concrete. The results of the presented approach for the analysis of limit load of reinforced concrete plates are in good agreement with the available experimental results (which are available in the literature) for those plates, where the load is monotonically increased until the failure of the plate. The essence of this approach is that it takes into account the gradual degradation of reinforced concrete due to crushing of the concrete and the yielding of the reinforcement. Although the approach is based on the nonlinear finite element method it is relatively simple and robust. The advantage of this approach, compared to the theory of plastic lines, is the information on the displacements (limit ductility), which may be interesting for studies on limit state of serviceability. • We have derived and programmed the small strain stress resultant elastoplastic and elastoviscoplastic plate finite element, where both isotropic and kinematic hardening were considered. We have also derived an algorithm, that unifies both inelastic formulations within one computational framework. Numerical results of the presented formulation have been compared with the stress formulation from Abaqus [Hobbit et al., 2007] as well as with the stress resultant formulation with a parameter that takes into account gradual spreading of through-the-thickness plasticity. It has been shown that, regarding the accuracy of the limit load computation, the mesh density plays more important role than the type of elastoplastic formulation. • We have revisited [Simo and Kennedy, 1992], [Brank et al., 1997] and programed the inelastic geometrically exact shell finite element, based on stress-resultant plasticity with isotropic and kinematic hardening. We have also derived a computational algorithm for the multi-surface plasticity of the shell stress-resultant formulation. Several numerical examples have been presented, which show a very satisfying performance of the presented approach compared to results from literature. In numerical example presented in Section 4.4.1 we have encountered some problems with convergence of the plasticity loop. We believe this was caused by the relatively large jumps between two equilibrated configurations that occur due to local buckling of the cylinder. This problem could be circumvented by a more sophisticated path following approach. In Section 4.4.3 we present the limit load analysis of rectangular plate. The comparison of the results obtained by the geometrically nonlinear shell formulation with the results of the geometrically linear plate formulation presented in Chapter 3 show that geometrically nonlinear effects can be significant when analyzing thin metal plates. • We have derived and programmed the elastoplastic Euler-Bernoulli beam finite element with embedded discontinuity. We also presented a sequential multi-scale procedure to obtain the beam constitutive model parameters. The multi-scale procedure belongs to the class of weak coupling methods, since we carry out the sequential computations. The results of the shell model computations, accounting for material and geometric localized instability, are stored to be used within the beam model softening response. As presented by numerical simulations, performance of the proposed multi-scale computational approach is very satisfying. One of its main features is that detection and development of the softening plastic hinges in the frame is fully automatic, and spreads gradually in accordance with stress redistribution in the course of the nonlinear analysis. When analyzing frame structures, we noticed that the size of the load step has an influence on the results in the softening stage of the structure. If the step is too big, several discontinuities can appear simultaneously, which leads to convergence problems or to an equilibrated configuration, which would not have been encounter if smaller load steps had been used. In order to circumvent the above mentioned problems, we adapted the size of the load step, so that only one new discontinuity can appear in the frame within one load step. One possibility for the future work in the context of failure analysis of frame structures is a development of a new formulation, where the multi-scale procedure is carried out simultaneously. • We have derived and programmed an elastoplastic quadrilateral two-dimensional finite element with embedded strong discontinuity. The discontinuity kinematics allow to model linear jumps in both normal and tangential displacements along the discontinuity line. In one of the early versions of the derived formulation we allowed the appearance of the discontinuity in any finite element, if the stress level in that element triggered the softening response. This approach often led to convergence problems if several discontinuities appeared at once or to an over stiff response, if the discontinuities in the neighboring elements were not aligned properly. In order to make the discontinuity propagation algorithm more robust, the continuity of the discontinuity line between the elements has been enforced. We performed several numerical simulations where we modeled the crack growth in brittle materials and the delamination of the composite materials. The comparison shows that of our results are in good agreement with those from literature. The proposed formulation can be also used to model the failure of ductile materials where localized shear bands form. The results of our simulations show that the results are independent of the mesh. The convergence of the presented formulation depends on the problem of interest. When analyzing only one finite element, we did not encounter any problems, while the analyses of a more complex structures sometimes demanded a change of the load step (either smaller or larger) in order to obtain a converged configuration. It still remains an opened question, how to make the approach more robust. In numerical simulation presented in Section 7.4.4 we encountered a problem of determining the proper direction of the crack growth. This could be avoided with a formulation, where the direction of the crack growth would be determined by considering the stress state in the patch of the elements instead of only using one element. With some minor modifications of the present formulation one would obtain a formulation that is suitable for analysis of failure in soils. • Finally, we should mention that all finite elements have been generated by the AceGen [Korelc, 2007b] computer program for symbolic expression manipulation and that all the numerical simulations have been preformed in the AceFem [Korelc, 2007a] computer program. These programs have been found to offer a versatile environment for fast production and testing of new finite element formulations. Razširjeni povzetek Najbolj pogosto uporabljeno orodje za racunanje obnasanja konstrukcij je metoda koncnih elementov. To je numericna metoda, s katero resujemo robni problem (t.j. parcialne diferencialne enacbe s pripadajocimi robnimi pogoji), ki opisuje dolocen fizikalni problem. Osnovni gradnik te metode je t.i. koncni element. Z mrezo koncnih elementov iscemo aproksimacijo resitve na obmocju, ki ga obravnavamo. V nasem primeru opisemo z mrezo koncnih elementov celo konstrukcijo ali pa samo del konstrukcije (konstrukcijski element), ki ga analiziramo. Popolna trodimezionalna analiza celotne konstrukcije, pri kateri bi upostevali njeno tocno geometrijo in bi vsebovala podrobne opise vseh geometrijsko in materialno nelinearnih pojavov (npr. plastifikacijo jekla, razpokanje betona, zdrs med betonom in jekleno ojacitvijo, lokalno odpoved materiala, velike pomike in rotacije, lokalni in globalni uklon, ...) je, kljub vedno bolj zmogljivim racunalnikom in vedno boljsim modelom za opis omenjenih nelinearnih pojavov, se vedno neizvedljiva, hkrati pa tudi ni smiselna. Pri modeliranju tako privzamemo doloceno mero fizikalne in matematicne idealizacije. Ponavadi trorazsezno telo (v nasem primeru konstrukcijski element), z upostevanjem dolocenih kinematicnih predpostavk, modeliramo kot dvodimenzionalno (stene, plosce in lupine), oziroma kot enodimenzionalno (nosilci, palice). V ta namen uporabljamo koncne elemente za stene, plosce, lupine, nosilce in palice. Nadalje si modeliranje obnasanja konstrukcij in konstrukcijskih elementov pogosto poenostavimo z uporabo pribliznih ali poenostavljenih konstitutivnih modelov. V disertaciji so izpeljani razlicni koncni elementi, ki so primeri za uporabo pri analizi obnasanja konstrukcije (ali konstrukcijskih elementov) vse do njene popolne porusitve. V ta namen se najprej ukvarjamo z neelasticnimi materialnimi modeli za opis dveh tipicnih gradbenih materialov - jekla in armiranega betona. Nato pa se ukvarjamo se z opisom lokalne odpovedi materiala, ki postopoma vodi k mehcanju konstrukcije, torej k taksnemu odzivu, ko se njeni pomiki povecujejo ob hkratnem zmanjsevanju obtezbe. Za dolocitev mejne obtezbe (t.j. tiste obtezbe, ki jo konstrukcija se lahko prenese, preden se zacne njeno mehcanje), smo tako razvili nelinearno elasticne ter elastoplasticne materialne modele za razlicne konstrukcijske elemente. Skupna lastnost vseh izpeljanih materialnih modelov je v tem, da so doloceni na nivoju rezultant napetosti (t.j. notranjih sil), kar jih naredi robustne in racunsko manj zahtevne. Ko pa se lotimo obravnavati obnasanje konstrukcije vse do njene popolne porusitve, pa uporaba nelinearno elasticnih oziroma elastoplasticnih modelov ne zadostuje vec. Do zacetka porusitve konstrukcije (oziroma njenega dela) ponavadi pride zaradi zelo lokalnih pojavov, ki jih z uporabo materialnih modelov za zvezno snov ni moc primerno opisati. Kot primer si lahko zamislimo vecjo razpoko, ki nastane v betonski steni na mestih natezne odpovedi betona. Širina te razpoke je zelo majhna v primerjavi s preostalimi dimenzijami stene, vendar ima vpliv na obnasanje stene, saj v trdno telo vnese nezveznost. Le-to pa je zelo tezko opisati s standardno metodo koncnih elementov. Z vkljucitvijo singularnih polj, s katerimi opisujemo območja taksnih nezveznosti (t.j. območja lokalne porusitve materiala), v standardno metodo koncnih elementov, pa lahko pridemo do robustnega nacina za opis lokalne odpovedi materiala. Pridemo do t.i. metode koncnih elementov z vkljucenimi nezveznostmi (angl. "embedded discontinuity finite element method"). V doktorski disertaciji smo si zastavili vec nalog povezanih z analizo mejne nosilnosti konstrukcije, oziroma s porusno analizo konstrukcij. Resevanje vsake naloge je prikazano v svojem poglavju. Za prvo nalogo smo si izbrali izpeljavo koncnega elementa za racun mejne nosilnosti armiranobetonskih plosc. V ze obstojeco racunalnisko kodo za koncni element na osnovi Reissner-Mindlinove teorije plosc, ki je predstavljen v [Bohinc et al., 2009], smo vgradili materialni model za armiranobetonske plosce, ki je opisan v [Ibrahimbegovic et al., 1992], [Ibrahimbegovic and Frey, 1993b] ter [Ibrahimbegovic and Frey, 1994], s tem da smo upostevali dolocila Evrokoda 2 [Eurocode 2, 2004] za opis konstitutivnega obnasanja armiranega betona. S tem smo prisli do robustnega in relativno enostavnega modela za materialno nelinearno analizo armiranobetonskih plosc, s katerim lahko ocenimo njihovo mejno nosilnost, kar ima precejsnjo prakticno vrednost. Druga naloga je bila konsistentna izpeljava elastoplasticnega in elastoviskoplasticnega materialnega modela za plosce, ki upostevata tako izotropno kot tudi kinematicno utrjevanje in sta dolocena z rezultantami napetosti. Izpeljali smo nov numericni algoritem, ki hkrati zajema tako elastoplasticno kot tudi elastoviskoplasticno formulacijo. Oba materialna modela smo nato vgradili v ze omenjeni koncni element na osnovi Reissner-Mindlinove teorije plosc. S tem smo si omogocili materialno nelinearno analizo metalnih plosc, ki ima prav tako prakticno vrednost. Tretja naloga je bila izpeljava geometrijsko in materialno nelinearne (elastoplasticne) formulacije za geometrijsko tocne lupine, ki uposteva tako izotropno kot tudi kinematicno utrjevanje in je dolocena z rezultantami nepetosti. Poseben poudarek je bil na izpeljavi tocnega plasticnega algoritma za elastoplasticnost lupin z rezultantami napetosti, ki vsebuje dve funkciji tecenja. V obstojeci geometrijsko tocni koncni element za lupine (glej [Brank et al., 1995] in [Brank and Ibrahimbegovic, 2001]) smo nato vgradili nee- lasticni materialni model s pripadajocimi novimi algoritmi. Prve tri naloge so se torej nanasale na materialno nelinearno analizo plosc in lupin iz kovin in armiranega betona. Modeli, ki omogocajo tovrstne analize imajo precejsnje prakticne vrednosti, vendar ne omogocajo porusne analize, t.j. analize, ki bi omogocala popolno porusitev konstrukcije ali konstrukcijskega elementa. Taksnih analiz pa smo se lotili v naslednjih nalogah. Cetrta naloga se je tako nanasala na porusno analizo ravninskih okvirnih konstrukcij. Ta tema ima prakticno vrednost pri sodobnem projektiranju, se posebej pri potresnem inzenirstvu. Uporabili smo metodo koncnih elementov z vgrajenimi nezveznostmi. Cilj je bil izpeljava novega elastoplasticnega koncnega elementa za ravninske nosilce z vkljuceno nezveznostjo v rotaciji. Koncni element je bil izpeljan na osnovi Euler-Bernoullijeve teorije nosilcev. Ker smo tudi v tem primeru neelasticne materialne modele izpeljali na nivoju rezultant napetosti, smo se srecali s problemom dolocitve potrebnih materialnih parametrov. V ta namen smo uporabili geometrijsko in materialno nelinearne koncne elemente za lupine. Izkazalo se je, da tako konsistentno izpeljani koncni elementi lahko ucinkovito nadomestijo standardne koncne elemente, ki se uporabljajo za porusno analizo okvirjev, pri katerih se vsa nelinearnost veze na vzmeti (ponavadi locirane na obeh koncih koncnega elementa). Zadnja naloga v okviru doktorske disertacije se je nanasala na izpeljavo nove druzine stirivozliscnih ravninskih koncnih elementov, ki upostevajo tako plasticno utrjevanje kot tudi lokalno plasticno mehcanje materiala. Zopet smo uporabili metodo koncnih elementov z vgrajenimi nezveznostmi. Posebnost izpeljanih novih koncnih elementov je v tem, da so skoki v pomikih na mestu nezveznosti linearni tako v smeri normale, kot v smeri tangente na nezveznost. Vecina obstojecih ravninskih koncnih elementov z vkljucenimi nezveznostmi namrec uposteva zgolj konstantne skoke v pomikih znotraj enega elementa. Poleg tega je ponavadi za odziv zvezne snovi upostevan zgolj linearno elasticni materialni model. Z novimi koncnimi elementi je mogoce slediti nastanku in sirjenju razpoke v krhkih dvodimenzijskih telesih ali formiranju lokaliziranega plastificiranega območja v duktilnih dvodimenzijskih telesih. Resitve, ki jih dobimo z novimi koncnimi elementi so neodvisne od mreze koncnih elementov. Analiza mejne obtežbe armiranobetonskih plosc V drugem poglavju se ukvarjamo z materialno nelinearno analizo armiranobetonskih plosc z metodo koncnih elementov. Za armirani beton uporabimo konstitutivni zakon, zapisan na nivoju rezultant napetosti. Ideja je povzeta po [Ibrahimbegovic et al., 1992], [Ibrahimbegovic and Frey, 1993b], [Ibrahimbegovic and Frey, 1994]; spremenimo jo v toliko, da uporabimo priporocila Evrokoda 2 [Eurocode 2, 2004] za opis konstitutivnega modela betona. S tem pristopom lahko identificiramo mejno stanje armiranobetonske plosce pri monotonem povecevanju obtezbe plosce; npr. pri obremenitvi zaradi lastne in stalne obtezbe. Ce pa je obtezevanje plosce taksno, da je prisotno veliko pomembnega razbremenjevanja in ponovnega obremenjevanja, se moramo omenjenemu nacinu racunanja mejne obtezbe odpovedati. Omeniti je potrebno, da je ze pri monotonem povecevanju verjetno, da bo prislo do lokalnega razbremenjevanja zaradi lokalne spremembe togosti plosce, vendar predpostavimo, da so neelasticni efekti zaradi taksnega razbremenjevanja zanemarljivi. Poleg tega predpostavimo tudi, da pomiki plosce niso tako veliki, da bi prisli do izraza membranski efekti, povezani s pojavom osnih (membranskih) sil v plosci. Taksne pojave pri analizi mejnih nosilnosti armiranobetonskih plosc lahko obravnavamo kot sekundarne in zato zanemarljive. Konstitutivni model za armirani beton je razdelejen na dve stanji: na stanje ner-azpokanega betona (stanje I) in na stanje razpokanega betona (stanje II). Razpokano stanje se aktivira, ko pride do porusitve betona v natezno najbolj obremenjenem vlaknu prereza. V stanju I upostevamo linearno elasticne konstitutivne zveze za izotropen material q = CsY, m = Cb k, kjer sta m in q vektorja momentov in precnih sil, Cs in Cs pa sta matriki elasticnih konstant za plosce (glej (2.4)). Podobno kot v stanju I, tudi v stanju II upostevamo linarno elasticen odziv za strizni del q = C sY, medtem ko pri upogibnem delu zanemarimo vpliv Poissonovega kolicnika in upostevamo neodvisna odziva v dveh pravokotnih smereh, dolocenih s kotoma 0 in 0 + 2, v obliki m = m (k) , kjer je K vrednost ukrivljenosti v smeri 0 oz. 0 + 2, m pa pripadajoci moment. Odziv v teh dveh smereh je dolocen s kolicino efektivne armature v teh smereh (slika 2.4) a^ = ^2 ai cos2(0 - ai), 2 = ^ ai cos2(0 + 2 - a), i i kjer je ai kolicina armature v plosci polozene v smeri ai, in z odsekoma linearnim diagramom, ki je dolocen z naslednjimi stanji prereza: • pojav prve razpoke v betonu (tocka A na sliki 2.3), • zacetek plastifikacije armature (tocka B na sliki 2.3), • porusitev betona v tlaku (tocka C na sliki 2.3). Locimo med izotropno armiranimi ploscami, pri katerih je efektivna kolicina armature enaka za vsak kot 0 (razdelek 2.3.1) isotropic . a^ _ const., in anizotropno armiranimi ploscami, pri katerih efektivna kolicina armature ni konstantna a^(0) _ const.. Pri slednjem obravnavamo dve možnosti (razdelek 2.3.2): s fiksno smerjo razpoke ter z rotirajočo smerjo razpoke. Glede na gornje kriterije dolocimo kot 0. Uporabljen pristop za racun mejne obtezbe armiranobetonskih plosc se izkaze kot zelo zadovoljiv. Rezultati analiz se dobro ujemajo z razpolozljivimi eksperimentalnimi rezultati (ki so na voljo v strokovni literaturi) za tiste plosce, kjer se je obtezba monotono povecevala vse do porusitve (razdelek 2.4). Bistvo nase analize je, da uposteva postopno degradacijo armiranega betona zaradi razpokanja betona in zaradi plastifikacije armature. Ceprav analiza temelji na nelinearni metodi koncnih elementov, je razmeroma preprosta in robustna in zato tudi uporabna v projektantski praksi. Neelastična analiza metalnih plosc V tretjem poglavju izpeljemo materialni model za elastoplasticne plosce, ki uposteva tako izotropno kot tudi kinematicno utrjevanje. Model je dolocen z rezultantami nepetosti (notranjimi silami). Elastoplasticno formulacijo nadgradimo v elastoviskoplasticno tipa Perzyna, glej npr. [Kojic and Bathe, 2005] in [Kleiber and Kowalczyk, 1996]. Oba nee-lasticna modela sta izpeljana ob predpostavki hipnega elasticnega odziva in principa maksimalne plasticne disipacije (za plasticnost) oz. kazenske oblike principa maksimalne plasticne disipacije (za viskoplasticnost), glej npr. [Ibrahimbegovic et al., 1998]. Osnovni gradniki tega konstitutivnega modela so: • Razcep celotnih deformacij £ na elasticni £e in plastični £p (viskoplastični £vp) del _ _e i „.p,vp £ _ £ + £ . • Predpostavka, da ima deformacijska energija naslednjo obliko, ki vsebuje tudi clena zaradi izotropnega in kinematicnega utrjevanja 4 (£e,e, K) _ 1 £e,T C£e + S (e) + 1 Q Hkin) XT D K, kjer je C matrika elastičnih konstant za plosce (glej (3.9)), £ parameter, ki doloca izotropno utrjevanje, S je funkcional povezan z izotropnim utrjevanjem, Hkin konstanta kinematičnega utrjevanja, D je matrika konstant za kinematično utrjevanje (glej (3.10)) in K je vektor "zaostalih" deformacij (t.i. "back-strain"). Funkcija tecenja 0, dolocena z rezultantami napetosti ct, q in a 0 (ct, q, a) = (ct + a)T A (ct + a) - (1 - —\ =0, kjer je q napetosti podoben parameter povezan z izotropnim utrjevanjem, a vector napetosti podobnih kolicin povezanih s kinematicnim utrjevanjem, A je matrika konstant, s katerimi opisemo izotropno plasticnost pri ploscah (glej (3.12)), pa je napetost na meji tecenja pri enoosnem preizkusu. • Princip maksimalne plasticne disipacije Dp, ki ga pri plasticnosti opisemo z iskanjem minimuma naslednjega funkcionala Lp (ct, q, a, 7) = -Dp (ct, q, a) + Y0 (ct, q, a), kjer je 7 Lagrangev mnozitelj. • Kazenska oblika principa maksimalne plasticne disipacije, ki jo v primeru visko-plasticnosti zapisemo v obliki iskanja minimuma naslednjega funkcionala Lvp (ct, q, a) = -Dvp (ct, q, a) + 1 g (0 (ct, q, a)), n kjer je n parameter viskoznosti in g(0)=( 002 ce0, y vrj [0 ce 0 < 0 ' kazenski funkcional. Resevanje se izvede v casovnih korakih. Nove vrednosti plasticnih spremenljivk na koncu tipicnega casovnega koraka [tn,tn+1] dolocimo z resitvijo ene skalarne enacbe 0n+1((CTra+1 + an+1) (Yn+1) , qn+1 (Yn+1)) = 0n+1 (Yn+1) = 0 ^ Yn+1, kjer je Yn+1 plasticni mnozitelj, Yn+1 pa njegova skonvergirana vrednost na koncu koraka. V primeru viskoplasticne analize zgornjo enacbo zamenjamo z n _ - A^Yn+1 + 0n+1 (Yn+1) = 0 ^ kjer je At = tn+1 — tn. V kolikor v viskoplasticni formulaciji postavimo viskozni parameter na n = 0, dobimo plasticno formulacijo in tako z enim racunskim algoritmom zajamemo obe neelasticni formulaciji. Rezultate nasih formulacij primerjamo (glej razdelek 3.4) z rezultati formulacije, ki je dolocena na nivoju napetosti [Hobbit et al., 2007], kakor tudi z rezultati formulacije, dolocene na nivoju rezultant napetosti s parametrom a, ki uposteva postopno sirjenje plastificiranega območja v smeri debeline plosce. Pokazalo se je, da ima gostota mreze koncnih elementov vecji vpliv na nivo limitne obtezbe, kot pa izbira same formulacije. NeelastiCna analiza metalnih lupin V cetrtem poglavju se ukvarjamo z elastoplasticno analizo tankih metalnih lupin z metodo koncnih elementov. Vse konstitutivne zveze so napisane na nivoju rezultat napetosti, s cemer se stevilo integracijskih tock po debelini lupine zmanjsa na vsega eno, kar naredi analizo mnogo hitrejso v primerjavi s pristopom, ki je dolocen z napetostmi, glej npr. [Brank et al., 1997]. Podoben pristop je bil z uporabo nekoliko drugacnega koncnega elementa in z drugacnim utrjevanjem ze prikazan v [Simo and Kennedy, 1992]. Osnovni gradniki konstitutivnega modela so: • Razcep celotnih deformacij £ na elasticni ee in plasticni £p del £ = £e + £p. Predpostavka, da ima deformacijska energija naslednjo obliko, ki vsebuje tudi clena zaradi izotropnega in kinematicnega utrjevanja 4 (£e,e ,*) = 2 £e'T C£e + S ) + 2 d K, kjer je C matrika elasticnih konstant za lupine (glej (4.31)), parameter, ki doloca izotropno utrjevanje, S je funkcional povezan z izotropnim utrjevanjem, D 3 Hkin I8 je matrika konstant za kinematicno utrjevanje, Hkin je konstanta kinematicnega utrjevanja, I8 je enotska matrika velikosti 8 x 8 in x je vektor "zaostalih" deformacij. Funkcija tecenja 0, dolocena z rezultantami napetosti a, ki jo omejujeta dve ploskvi (slika 4.1) 0, (a, q, a) = (a + a)T A, (a + a) — — p =1, 2, kjer je q napetosti podoben parameter povezan z izotropnim utrjevanjem, a vector napetosti podobnih koliCin povezanih s kinematiCnim utrjevanjem, A^ je matrika konstant, s katerimi opiSemo izotropno plastiCnost pri lupinah (glej (4.33)), ay pa je napetost na meji teCenja pri enoosnem preizkusu. PlastiCna disipacija Dp = aTep + qC1 + aTk > 0. • PrinCip maksimalne plastiCne disipaCije, ki ga opisemo z iskanjem minimuma naslednjega funkCionala 2 Lp (a,q, a, 71,72) = -Dp (a,q, a) + ^ 7^ (&,q, a), m=1 kjer sta Y in Y2 Lagrangeva mnozitelja. Notranje spremenljivke raCunamo na dva naCina (glej razdelek 4.3.2). Pri prvem naCinu resujemo le enaCbe povezane s funktijo teCenja RP (l1,n+1,l2,n+1) 01 (Y1,n+1, Y2,n+1) 02 (71,n+1,72,n+1) 0, ki jih zapisemo v vektor Rp. Z resitvijo za plastiCna mnozitelja 71,n+1 in 72,n+1 nato doloCimo nove vrednosti plastiCnih spremenljivk. Pri drugem naCinu pa v vektorju Rp nastopajo tudi evoluCijske enaCbe RP (£I+1,CI+1l1,n+1,Y2,n+1) P I P I -En+1 + En + 2^ M=1 lp.!+1V V,n+1 Cn +1 + Ci + E \i=1 2 01,n+1 02,n+1 0, tako da tudi plastiCne spremenljivke doloCimo neposredno. Ko imamo opravka z veCploskovnimi funktijami teCenja, nimamo vedno informaCije o tem, katera ploskev bo v skonvergiranem stanju aktivna (v zgornjih enaCbah smo predstavili le moznost, ko sta obe ploskvi aktivni). Tako pri dveh ploskvah loCimo med tremi moznimi sistemi enaCb R12 01,n+1 02,n+1 0, R1 01,n+1 Y2,n+1 0, R2 Y1,n+1 02,n+1 0, s katerimi pokrijemo vse mozne konCne izide (R12 obe ploskvi sta aktivni, R le prva ploskev je aktivna in R le druga ploskev je aktivna). Predstavljeni sta dve proCeduri za doloCitev pravilne resitve (glej razdelek 4.3.2). Pri prvi proCeduri poisCemo resitve za vse tri mozne sisteme enaCb in na konCu izberemo tisto resitev, ki zadosCa Kuhn-TuCkerjevim pogojem obremenjevanja in razbremenjevanja (7^,n+1 > 0, 0^,n+1 < 0). Pri drugi proceduri pa iscemo resitev le enega sistema enacb, ki pa se lahko med postopkom spreminja, in sicer upostevamo, da morajo biti pri vseh aktivnih ploskvah plasticni mnozitelji vedno nenegativni (7^,n+1 > 0). Predstavljenih je vec numericnih primerov (glej razdelek 4.4), ki kazejo dobro ujemanje nase formulacije z rezultati iz literature. Predstavitev koncepta vključene nezveznosti V petem poglavju je na primeru koncnega elementa za palico predstavljen koncept vkljuce-ne nezveznosti (angl. "embedded discontinuity"). Obravnavamo natezni test idealizirane palice, ki je predstavljen na sliki 5.1. Odziv palice je linearno elasticen dokler ni dosezena nosilnost materiala, nato pa se nosilnost palice linearno manjsa z vecanjem pomika, glej sliko 5.1 desno spodaj. Na zgornjem levem delu slike 5.1 je prikazana kljucna lastnost porusitve, in sicer, da pride do zmanjsanja nosilnosti zaradi zelo lokalnih efektov, ki so skoncentrirani v okolici najsibkejse tocke palice (mi smo izbrali sredino palice). Na sliki 5.1 levo spodaj je prikazana porusitev zaradi t.i. "necking efekta", ki je znacilen pri natezni porusitvi metalnih palic. Omeniti je potrebno, da je razen sredine idealne palice, kjer se pojavijo nepovratne (plasticne) deformacije, celotna palica v elasticnem obmocju. Tako je celoten odziv palice v mehcanju odvisen le od obnasanje sredinske tocke palice. Na zgornjem desnem delu slike 5.1 je predstavljena razporeditev pomikov vzdolz osi palice, ko nosilnost palice se ni v celoti izcrpana. Skok v pomikih u si lahko razlagamo tudi kot lokalno plasticno deformacijo na mestu nezveznosti. Izpeljave koncnega elementa, ki je predstavljen v petem poglavju, sledi naslednjemu postopku. V standarni izoparametricni koncni element za palico dodamo parameter a, s katerim opisemo kinematiko nezveznosti ter dobimo obogateno polje vzdolznih pomikov kjer s h oznacimo diskretno aproksimirano kolicino, £ G [-1,1] je brezdimenzionalna lokalna koordinata, a = 1, 2 je indeks vozlisca, N je interpolacijska funkcija (glej (5.2)), d« je vozliscni pomik in Ma je interpolacijska funkcija povezana z nezveznostjo (glej (5.5)). Z «MkE oznacimo izraz, ki izhaja iz standardne metode koncnih elementov. Z odvajanjem polja pomikov po koordinati x (glej (5.1) in sliko 5.3) pridemo do polja deformacij 2 uh(£) = £ N«(£K + Ma(£)a, / obogatitev "mke kjer je Ba deformacijski operator, ki ga povezujemo z vozliscnimi pomiki (glej (5.7)), Ga = Ga + Ga je deformacijski operator, ki ga povezujemo s kinematiko nezveznosti in ga razdelimo na zvezni Ga in nezvezni Ga del (glej (5.7)), td je del deformacij, ki izhaja iz vozlisCnih pomikov, ta pa del deformacij, ki izhaja iz nezveznosti. Upostevamo, da se virtualne deformacije znotraj koncnega elementa interpolirajo kot 2 e = Ba(Qda + Gaa, a=1 kjer je da virtualni vozliscni pomik, a je virtualni parameter nezveznosti ter * 1 f Ga Ga Tel Gadxj operator, s katerim zagotovimo konvergenco elementa v smislu patch testa. Ce v princip virtualnega dela notranjih sil vstavim nastavek za virtualne deformacije, dobimo OTmt'(e) = ^ A(e) [ daBaodx + A(e) [ aGaadx, a=1 v _, v-S/-' v standardna MKE kjer je A(e) povrsina precnega prereza elementa e, obmocje elementa in a napetost v elementu. Iz izraza "dodatno" tako dobimo eno dodatno ravnotezno enacbo h(e) = A(e) ( f Gaadx + t) =0, ki zagotavlja, da je napetost na nezveznosti t v ravnotezju z napetostnim stanjem v koncnem elementu. Togi plasticni kohezijski zakon, predstavljen na sliki 5.5, gradimo na osnovi kriterija zacetka nezveznosti ^(t,q) = t - (au - q) < 0, ter potenciala, povezanega z mehcanjem materiala = = 1 =2 ) = ^K.4 , kjer je au nosilnost materiala, q napetosti podobna kolicina povezana z mehcanjem, - je parameter, ki opisuje mehcanje materiala, Ks pa je modul linearnega mehcanja. Preostale sestavine kohezijskega zakona dolocimo z upostevanjem plasticne disipacije D = ta + qf, n e ter principa maksimalne plasticne disipacije, ki ga zapisemo kot min max t,q y ~L(t,q,Yf) = —D (t,q)+ 7 0(t,q) kjer je 7 Lagrangev mnozitelj. Resevanje se izvede v casovnih korakih, racunski postopek pa je razdeljen na lokalni in na globalni del. V lokalnem delu na koncu casovnega koraka [Tn, Tn+1] dolocimo popravke k spremenljivkam mehcanja z resitvijo ene skalarne enacbe 0 = tn+1 (7n+0 — ( — — qn+1 (jn+1)) = 0 (jn+1) = 0, kjer je 7n+1 mnozitelj povezan s togo plasticnostjo na nezveznosti. Da dolocimo napetost na nezveznosti, eksplicitno upostevamo dodatno ravnotezno enacbo tn+1 = — Gaan+1dx = a n+l, Jne kjer s a 0 pa linearni modul utrjevanja. • FunkCija teCenja za zvezno snov 0(M,q) = |M|- (My - q) < 0, kjer je My > 0 upogibni moment na meji teCenja, q pa upogibnemu momentu podobna koliCina povezana z linearnim utrjevanjem materiala. • Kriterij zaCetka nezveznosti 0(tM,q) = |tM| - (Mu - q) < 0, kjer je Mu > 0 upogibna nosilnost prereza, q pa upogibnemu momentu podobna koliCina povezana z mehCanjem. Potencial, povezan z mehcanjem = = 1 =2 =(-) = 1 Ks- , kjer je Ks < 0 linearni modul mehCanja, - pa koliCina, ki opise mehCanje. Preostale sestavine za elastoplastiCni odziv materiala doloCimo z upoštevanjem plastiCne disipacije D = Mlf + q-, ter prinCipa maksimalne plastiCne disipatije, ki ga zapisemo kot min max [LP(M,q,Yf) = -D(M,q) + Yf$(M,q)] , M,q kjer je y Lagrangev mnozitelj. Podobno definiramo tudi plastiCno disipaCijo na nezveznosti D = D = t M a 0 + q-, ter tudi tu upostevamo prinCip maksimalne plastiCne disipaCije min max tM,q Y L (tM, q, f) = -D (tM, q) + Y0(tM, q) kjer je Y Lagrangev mnozitelj. Materialno nelinearni odziv konCnega elementa za nosilCe, je doloCen s parametri: My, Kh, Mu in Ks. Upogibni moment, ki oznanja zaCetek materialno nelinearnega odziva, je doloCen z momentom, ki povzroCi plastifikaCijo najbolj obremenjenaga vlakna v preCnem prerezu in je odvisen od nivoja osne sile My (N ) = Way (1 - M), AOy kjer je W upogibna odpornost prereza. Da dobimo ustrezne vrednosti za preostale tri parametre, se obrnemo na modeliranje dela konstrukdje s podrobnejsim modelom, ki v nasem primeru sloni na geometrijsko in materialno nelinearnih konCnih elementih za lupine. Maksimalna upogibna nosilnost nosilCa Mu je tako doloCena z maksimalno nosilnostjo, ki smo jo zabelezili pri analizi z lupinastimi konCnimi elementi MUef (N) Mu (N) = Mruef (N). Ta vrednost nam tudi pomaga, da v rezultatih lupinastega modela loCimo med rezimom utrjevanja in rezimom mehCanja, do katerega lahko pride zaradi odpovedi materiala (mehCanja) ali pa lokalne geometrijske nestabilnosti (uklona). Nadalje predpostavimo, da mora biti plastiCno delo (glej (6.53) in (6.54)), ki ga opravi model nosilCa, enako plastiCnemu delu, ki ga opravi model lupine. Ker pa plastiCno delo poteka v dveh rezimih (utrjevanje in mehCanje), zagotovimo, da je tudi v posameznem rezimu plastiCno delo pri obeh modelih enako 777P TTTP.r ef TTrP -—ef EW (N) = EW (N), EW (N) = EW (N) , -P -P ,ref kjer je EW plastiCno delo pri utrjevanju za model z nosilti, EW plastiCno delo pri utrjevanju za model z lupinami, EW plastiCno delo pri mehCanju za model z nosilti ter ==P,ref ^ EW plastiCno delo pri mehCanju za model z lupinami. Ce upostevamo, da so vsi preCni prerezi v vzorcu dolzine Lref pod (priblizno) enakim napetostnim stanjem, dobimo naslednji izraz za modul utrjevanja (Mu2 (N) - My2 (N))Lref Kh (N ) = --p re f-• 2EWP, f (N) Nasprotno pa v primeru mehCanja upostevamo, da se vse zgodi le v enem prerezu in tako dobim modul mehCanja kot M2 (N) |K (N)| = MU, ) , K < 0. 2EW (N) Resevanje se izvede v Casovnih korakih, raCunski postopek pa je razdeljen na lokalni in na globalni del. V lokalnem delu na konCu tipiCnega Casovnega koraka [tn,tn+1] doloCimo nove vrednosti spremenljivk povezanih z utrjevanjem materiala z resitvijo enaCbe ^(M^dSf^(ti+i )),q(čn+i(Yn+i)))=r (Yi+i) = 0, kjer je ip = 1, 2, 3 indeks integraCijske toCke, i steveC globalne iteraCijske sheme, 7n+i pa plastiCni mnozitelj povezan z utrjevanjem. Spremenljivke, povezane z mehCanjem, pa doloCimo z resitvijo naslednje enaCbe nti(YSi)))=0(e)(Yie+i) = 0, =(e) kjer je 7n+i plastiCni mnozitelj povezan z mehCanjem. Upogibni moment na nezveznosti doloCimo z uporabo dodatne ravnotezne enaCbe ,(e) _ I ( (e) )M(d(e)'(i) "P'iP (e) )d tM,n+i = — / G (Xd ,x)M (dn+i , Kn , a0,n+i)dx. Jne V globalni fazi doloCimo popravke Ad^H i) za trenutne vrednosti vozlisCnih pomikov d(e),(i) = d(e),(i-i) + Ad(e),(i-i) dn+i = dn+i + Adn+i • Prispevek enega končnega elementa k sistemu globalnih enačb je K(e) Kfa Khd K^a (i) ( Ad(e),(i) \ / j?ext,(e) ~int,(e),(i) I+U = 1 f n+1 - f n+1 , Aa(e)'W n+1 V Aae,n+1 kjer so K(e), Kfa, Khd in Kha podmatrike tangentne matrike končnega elementa (glej (6.82)), fI+'1(e)'(i) so notranje vozlisčne sile (glej (6.25)), fn+1(e) pa zunanja vozlisčna obtežba. S pomočjo statične kondenzacije pridemo do tangentne matrike končnega elementa K(e)'(i) = K(e),(i) _ Kfa,(i) (Kha,(i)\-1 Khd,(i) K n+1 = K n+1 K n+1 l Kn+1 I K n+1 , ki je povsem enake oblike kot v standardni metodi končnih elementov, zato je globalni sistem resevanja nespremenjen. V prvem numeričnem primeru v razdelku 6.5 je prikazan izračun s končnimi elementi za lupine in z rezultati tega izračuna smo določili bilinearno aproksimačijo za upogibno nosilnost M (N) = Mref (N) = I M^a-OS + O.85NNy) if N< -0.035Ny Mu (N) = Mu (NMref,0 y if N >-0.035Ny ' medtem ko smo za modula utrjevanja in mehčanja upostevali kar konstantni vrednosti Kh(N) = 1.06 ■ 107 kN/čm2, Ks(N) = -3.28 ■ 105 kN/čm2, kjer z Mref označimo aproksimirano funkčijo na osnovi analize z lupinami, M^f'° je upogibna nosilnost prereza določena pri analizi z lupinami pri nični vrednosti osne sile in Ny = Aay. Materialni parametri določeni pri analizi z lupinami so bili nato uporabljeni v vseh nadaljnih numeričnih primerih z uporabo končnih elementov za nosilče. Predstavljen pristop spada v kategorijo simulačij na več nivojih (angl. "multi-sčale"), kjer najprej izvedemo izračun na podrobnejsem modelu in dobljene rezultate nato nesemo v bolj robustni makro model. Rezultati analize z lupinami, kjer so upostevane tako materialne kot tudi geometrijske nelinearnosti, so tako shranjeni in upostevani v modelu za nosilče. Ena glavnih lastnosti tega pristopa je avtomatsko zaznavanje in razvoj plastičnih členkov (nezveznosti), ki se pojavijo postopoma, v skladu s prerazporeditvijo obremenitev tekom nelinearne analize. Pri mnogih standardnih postopkih za račun t.i. "push-over" analiz je namreč potrebno vnaprej določiti kritična mesta ter njihov odziv. Porusna analiza z ravninskimi končnimi elementi V sedmem poglavju obravnavamo druzino stirivozlisčnih elastoplastičnih ravninskih elementov z vgrajeno nezveznostjo. Izdelava učinkovitih stirivozlisčnih končnih elementov z vgrajeno nezveznostjo, ki niso podvrzeni blokiranju, je mnogo bolj zahtevna (glej npr. [Linder and Armero, 2007] in [Manzoli and Shing, 2006]), kot pa v primeru trikotnih elementov, ki so sposobni opisati zgolj konstantno napetostno stanje. Morda je tu potrebno iskati vzrok, da je bila do sedaj vecina raziskav povezanih z vključevanjem nezveznosti v ravninske končne elemente izvedena ravno s trikotnimi elementi, glej npr. [Sancho et al., 2007], [Ibrahimbegovic and Brancherie, 2003], [Brancherie and Ibrahimbegovic, 2008], [Mosler, 2005], [Jirasek and Zimmermann, 2001]. Kinematika discontinuitete omogoca linearna skoke tako v smeri normale kot tudi v smeri tangente nezveznosti. S predstavljenim koncnim elementom lahko analiziramo natezno porusitev ravninskih betonskih vzorcev s pripadajocim sirjenjem razpoke. Ista formulacija se lahko uporabi tudi za analizo porusitve duktilnih materialov, pri katerih pride do pojava striznih pasov in tudi za analizo delaminacije pri kompozitnih materialih. Izpeljava koncnega elementa, ki je predstavljen v sedmem poglavju, sledi naslednjemu postopku. Z vpeljavo stirih dodatnih parametrov (an0, an1, am0 in am1), povezanih s kine-matiko nezveznosti, v izoparametricni ravninski koncni element, pridemo do obogatenega polja pomikov 4 vh(£, re) = 52 Na($)da +Y, Mmode(i, re)amode, mode G (n0, nI, m0, ml), a=1 mode -..-' -..-' kjer povezujeme dodatne parametere z eno od naslednjih oblik obnasanja nezveznosti • "n0" - konstaten skok pomikov v smeri normale nezveznosti, • "nl" - linearen skok pomikov v smeri normale nezveznosti, • "m0" - konstaten skok pomikov v smeri tangente nezveznosti, • "ml" - linearen skok pomikov v smeri tangente nezveznosti. S h oznacimo diskretno aproksimirano kolicino, £ = [C,v]T G [-1,1] x [-1,1] je brezdi-menzionalana koordinata znotraj elementa, z re oznacimo nezveznost znotraj elementa e, a = 1, 2, 3, 4 je indeks vozlisca, Na je interpolacijska funkcija (glej (7.2)), da = [ux,uy]T so vozliscni pomiki, Mmode je interpolacijska matrika povezana s parametri nezveznosti, z v/h oznacimo del polja pomikov, ki izhaja iz interpolacije vozliscnih vrednosti, z v^ pa del, ki izhaja iz parametrov nezveznosti. S preucevanjem togih pomikov obmocji Qe- in Qe+ (Qe = Qe- U He+), ki sta loceni z nezveznostjo re (glej sliko 7.2 ), dolocimo interpolacijske matrike povezane s parametri nezveznosti kot vh _ vh _ ""mode ""d,mode mode H , amode kjer sta u}mnode in u}^ mode doloCena v (7.4), amode pa je vrednost togega pomika na nezveznosti. Polje obogatenih deformatij dobimo s simetriCnim gradientom polja pomikov e = Bada + Gn0an0 + Gn1an1 + Gm0am0 + G m1am1, a=1 kjer je Ba deformatijski operator povezan z vozlisCnimi pomiki (glej (7.24)), Gn0, Gn1, Gm0 in Gm1 pa so deformaCijski operatorji, ki so povezani s kinematiko nezveznosti (glej (7.25)-(7.28)). Upostevamo, da se virtualne deformaCije znotraj konCnega elementa inter-polirajo kot 4 e = £ B aa + Gn0 an0 + Gn1an1 + Gm0am0 + Gm1am1-, a=1 kjer so da vozlisCni virtualni pomiki, an0, an1, am0 in am1 so virtualni parametri nezveznosti ter Gn0, Gn1, Gm0 in Gm1 operatorji, ki jih doloCimo kot 1 G mode G mode Ac Gmode , J Qe da zagotovimo konvergenCo v smislu t.i. "patCh" testa. Ce v printip virtualnega dela vstavimo izraz za virtualne deformaCije, dobimo OTnt'(e) = t(e) / eTadtt 4 (e) V t(e) f daTbTadn + a=1 Jce standardna MKE T T T T t(e) / an0Gn0a + an1 Gn1a + am0Gm0a + 0 pa parameter povezan z utrjevanjem materiala. FunkCija teCenja za zvezno snov 0 (u,q) = uTAu - - < 0, kjer je A matrika konstant, s katerimi opisemo izotropno plastiCnost pri ravninskih problemih (glej (7.56)), ay je napetost na meji teCenja pri enoosnem preizkusu, q pa napetosti podoben parameter povezan z izotropnim utrjevanjem. • Kohezijski zakon na nezveznosti zapisan s skoki v pomikih na nezveznosti u = t = t(u), ki ga lahko zapisemo tudi s kriterijem zaCetka nezveznosti 0 = 0(t,q) < 0, in potenCialom povezanim z mehCanjem materiala i(l), kjer je £ koliCina, ki opise mehCanje. Preostale sestavine za elastoplastiCni odziv materiala doloCimo z upostevanjem plastiCne disipaCije D = u e + q£, ter prinCipa maksimalne plastiCne disipaCije, ki ga zapisemo kot min max rZP(u,q,'y) = —D(u,q) + 70(u,q)! , kjer je 7 Lagrangev mnozitelj. Podobno definiramo tudi plastično disipačijo na nezveznosti d = D = tTt + qf, ter tudi tu upostevamo prinčip maksimalne plastične disipačije min max t,q Y L,(t,qn) = -D (t, q) + 70(t,q) kjer je 7 Lagrangev mnozitelj. Resevanje se izvede v časovnih korakih, računski postopek pa je razdeljen na lokalni in na globalni del. V lokalnem delu na konču tipičnega časovnega koraka [rn,rn+1j, če ni ze določeno, najprej, z upostevanjem povprečnega napetostnega stanja v končnem elementu, določimo geometrijo nezveznosti n n a avg (d(e) —P,bip I+1, —n ) m m a avg (d (e) —p,bip I+1, —n ) xrE, kjer je n normala nezveznosti, m tangenta nezveznosti, aavg povprečna napetost v elementu, določena na osnovi napetostnega stanja v integračijskih točkah bip = 1, 2, 3, 4, ter xrE končna točka nezveznosti v elementu (glej sliko 7.4). Nove vrednosti spremenljivk za mehčanje materiala ter skokov v pomikih v dveh integračijskih točkah nezveznosti, določimo s hkratnim resevanjem naslednjih enačb 0 (tn+1 (aie+1 (ji+1 ,q(L+1 (jn+1))) = 0 (ji+1,Yn+1)=° 0 (tn+1 (aie+1 (fi+1 ^(in+1 (jn+1))) =0 (Yn+1,YI+^=° kjer z indeksom 1 opisemo količino, ki se nanasa na prvo integračijsko točko, z indeksom 2 —1 —2 opisemo količino, ki se nanasa na drugo integračijsko točko, 7n+1 in 7n+1 pa sta plastična mnozitelja povezana z mehčanjem. Upostevali smo naslednje zveze med skoki v pomikih v integračijskih točkah nezveznosti ter kinematičnimi parametri nezveznosti a (e)(U1, U2) =1 + =2 =1 =2 =1 + =2 =1 =2 Un + Un U, — U„ Um + U m Um Um 2 er er 2 er er T kjer sta eh in eT koordinati integračijskih točk na nezveznosti (glej (7.80)). Pri izračunu napetosti v integračijskih točkah nezveznosti smo ekspličitno uporabili dodatne ravnotezne enačbe h' (a(dI+1, -Ibip, + h (tI+1, tI+1) 0 t1 , ^n+1, t n+1 V lokalni fazi določimo tudi nove vrednosti spremenljivk v integračijskih točkah elementa, ki so povezane z utrjevanjem materiala, z resitvijo naslednje enačbe 0 (a II+ 1(d+(i), i))rq(tp+ 1(yII+1))) = 0Uv(i7+1) = 0, e 2 kjer je i steveC globalne iteraCjske sheme, 7^+1 pa plastiCni mnozitelj povezan z utrjevanjem. V globalni fazi doloCimo popravke Adn+1 1) za trenutne vrednosti vozlisCnih pomikov d(e),(i) = d(e),(i-1) + A d(e),(i-1) dn+1 = dn+1 + Adn+1 • Prispevek enega konCnega elementa k sistemu globalnih enaCb je K(e) K/a K hd ^^ha (i) / Ad(e),(i) \ / J?ext,(e) -int,(e),(i) Adn+1 \ = l Jn+1 - Jn+1 4 Aa(e)'(i) = l 0 n+1 V Aan+1 ' \ 0 kjer so K(e), K/a, Khd in Kha podmatrike tangentne matrike konCnega elementa (glej (7.104)), Jr+f'(i) so notranje vozlisCne sile (glej (7.44)), /n+i(e) pa zunanja vozlisCna obtezba. S pomoCjo statiCne kondenzaCije pridemo do tangentne matrike konCnega elementa (e),(i) _ (e),(i) ~/a,(i) /ha,(i)\-1 ^hd,(i) v , _ l^(e),(i) _ |>-/a,(iw |>-ha,(i) \ K n+1 = K n+1 K n+1 l K n+1 J K n+1 , ki je povsem enake oblike kot v standardni metodi konCnih elementov in je s tem globalni sistem resevanja nespremenjen. V razdelku 7.4 so prikazani numeriCni primeri, kjer analiziramo dvodimenzionalne betonske vzorce, delaminaCijo kompozitnih materialov ter porusitev duktilnih materialov. Element omogoCa opis linearnih skokov pomikov tako v smeri tangente kot tudi normale. Daje postal postopek sirjenja razpoke bolj robusten, smo povezali nezveznosti med dvema sosednjima elementoma. Zaključki Namen dela je poglobitev znanja o obnasanju konstruktij v okoliti mejne nosilnosti, ter pri obremenitvah, ki povzroCijo porusitev konstruktije, z uporabo numeriCnih metod, kot je na primer metoda konCnih elementov z vkljuCenimi nezveznostmi. V disertaCiji smo se osredotoCali predvsem na modeliranje materialnih nelinearnosti in na modeliranje lokalne porusitve materiala v kontekstu omenjene numeriCne metode. V raziskovalnem delu, ki je predstavljeno, smo izvedli naslednja dela in prisli do naslednjih ugotovitev: • Izpeljali in sprogramirali smo nelinearno elastiCen konCni element za raCun mejne obtezbe plosC, ki je definiran z rezultantami, njegovo delovanje smo pa preverili z raCunskimi primeri. Upostevali priporoCila Evrokoda 2 [EuroCode 2, 2004] za opis konstitutivnih zvez armiranega betona. Rezultati raCunskih primerov se dobro ujemajo z razpolozljivimi eksperimentalnimi rezultati (ki so na voljo v strokovni literaturi) za tiste plosCe, kjer se je obtezba monotono poveCevala vse do porusitve. Bistvo uporabljene analize je, da uposteva postopno degradacijo armiranega betona zaradi razpokanja betona in zaradi plastifikacije armature. Ceprav analiza temelji na nelinearni metodi koncnih elementov, je razmeroma preprosta in robustna. Prednost prikazanega pristopa, glede na precej uporabljano teorijo plasticnih linij, je informacija o velikosti pomikov pri dosezeni mejni nosilnosti, ki je lahko zanimiva za studij mejnega stanja uporabnosti. • Izpeljali smo elastoplasticno in elastoviskoplasticno formulacijo za plosce in ju, z novim algoritmom, ki hkrati zajame obe neelasticni formulaciji, vgradili v koncni element za Reissner-Mindlinove plosce. Obe formulaciji sta definirani z rezultantami napetosti in upostevata tako izotropno kot tudi kinematicno utrjevanje materiala. Koncni element smo preverili z racunskimi primeri, kjer smo rezultate nase formulacije primerjali s formulacijo, ki je dolocena na nivoju napetosti in s formulacijo z rezultantami napetosti, ki uposteva postopno plastifikacijo v smeri debeline plosce. Ugotovili smo, da ima gostota mreze koncnih elementov vecji vpliv na nivo mejne obtezbe, kot pa izbira same formulacije. • Izpeljali smo geometrijsko in materialno nelinearno (elastoplasticno) formulacijo za geometrijsko tocne lupine, ki vsebuje tako izotropno kot tudi kinematicno utrjevanje in je dolocena z rezultantami nepetosti. Razvili smo tudi algoritem za elastoplasticne lupine z rezultantami napetosti, ki vsebuje dve funkciji tecenja. V obstojeci geometrijsko tocni koncni element za lupine smo vgradili neelasticni materialni model s pripadajocimi novimi algoritmi in ga preverili z racunskimi primeri. Primerjava rezultatov nase formulacije z rezultati iz literature kaze dobro ujemanje. Pri racunskem primeru v razdelku 4.4.1, smo imeli nekaj tezav s konvergenco plasticne zanke, ki najverjetneje nastane zaradi velikih skokov med dvemi ravnoteznimi konfiguracijami, zaradi pojava lokalnih uklonov. Z uporabo boljsega pristopa za sledenje obtezne poti bi se problemom s konvergenco izognili. V razdelku 4.4.3 smo na primeru jeklene plosce primerjali elastoplasticno formulacijo za lupine z elastoplasticno formulacijo za plosce, ki je predstavljena v tretjem poglavju. Ugotovili smo, da so lahko geometrijsko nelinearni efekti pri analizi tankih plosc zelo veliki. • Izpeljali in sprogramirali smo elastoplasticno formulacijo z vkljucenimi nezveznostmi za analizo ravninskih nosilcev. Formulacija je izpeljana na osnovi Euler-Bernoullijeve teorije nosilcev in vsebuje tako utrjevanje materiala kot tudi njegovo lokalno mehcanje. Izpeljan je tudi postopek zaporednega racuna na vec nivojih, s katerim dolocimo konstitutivne parametre za nosilce. Z rezultati analize, ki je izvedena s koncnimi elementi za lupine, zajamemo materialne in geometrijske nelinearnosti, vkljucno z lokalnimi efekti. Te rezultate shranimo in jih uporabimo pri analizi s konCnimi elementi za nosilCe. Z numeriCnimi testi smo ugotovili dobro delovanje predstavljenega pristopa. Ena glavnih lastnosti tega pristopa je avtomatsko zaznavanje in razvoj plastiCnih Clenkov (nezveznosti), ki se pojavljajo postopoma v skladu s prerazporeditvijo obremenitev tekom nelinearne analize. Pri analizi ravninskih okvirjev se je izkazalo, da lahko velikost obteznega koraka vpliva na resitev v obmoCju mehCanja konstrukCije. Ce je korak prevelik, se lahko nezveznosti naenkrat pojavijo na veC mestih, kar lahko privede do problemov s konvergenCo ali pa je poslediCa prevelikih korakov neka nova ravnotezna lega, ki je drugaCna od tiste, ki bi jo dobili z manjsim korakom. Tako smo pri raCunskih primerih uporabili tako velikost koraka, da je v enem obteznem koraku nastala kvečjemu ena nova nezveznost. Naslednji korak pri analizi okvirjev bi lahko bil razvoj nove formulaCije, kjer bi veCnivojska analiza potekala istoCasno. • Izpeljali smo elastoplastiCno formulaCijo z vgrajenimi nezveznostmi za analizo ravninskih problemov. FormulaCijo smo vgradili v izoparametriCni stirivozlisCni konCni element in ga preverili z raCunskimi testi. Pri opisu kinematike nezveznosti smo upostevali linearne skoke v pomikih, tako v smeri tangente kot tudi normale nezveznosti. V prvih razliCiCah formulaCije se je lahko nezveznost pojavila v kateremkoli konCneme elementu, Ce je tako narekovalo napetostno stanje v njem. Tak pristop je privedel do problemov pri konvergenCi ali pa je zaradi razliCno orientiranih nezveznosti v sosednjih elementih prislo do pretogega odziva. Da je postal postopek sirjenja razpoke bolj robusten, smo povezali nezveznost med dvema sosednjima elementoma. Primerjava rezultatov nase formulaCije pri analizi porusitve betonskih vzorcev ter analizi delaminaCije kompozitnih materialov kaze dobro ujemanje z rezultati iz literature. Na primeru analize elastoplastiCnega vzorca z duktilno porusitvijo smo z razliCnimi gostotami mrez pokazali, da je rezultat neodvisen od mreze konCnih elementov. KonvergenCa predstavljenega pristopa je zelo odvisna od obravnavanega problema. Analize z enim konCnim elementom delujejo brez tezav, medtem ko je pri analizi bolj komplidranih konstruktij vCasih potrebna intervendja in je potrebno spremeniti dolzino obteznega koraka, da pridemo do uravnotezene konfiguraCije. Zaenkrat obstaja se odprto vprasanje, kako narediti raCunski postopek se bolj robusten. Pri raCunskem primeru v razdelku 7.4.4 smo imeli tezave z doloCitvijo pravilne smeri sirjenja nezveznosti. Temu problemu bi se lahko izognili z razvojem formulaCije, ki bi za kriterij smeri sirjenja nezveznosti, upostevala napetostno stanje v sirsi okoliCi elementa in ne le v enem elementu. Z manjsimi spremembami ze obstojeCe formulaCije, pa bi dobili orodje, ki bi bilo primerno za analizo drsin pri zemljinah. • Vsi, v tej disertaciji predstavljeni koncni elementi, so bili pripravljeni v programskem okolju AceGen [Korelc, 2007b]. Vsa testiranja koncnih elementov so bila izvedena v programskem okolju AceFem [Korelc, 2007a]. Obe programski okolji sta se izkazali za vsestranski orodji, ki omogocata hitro "proizvodnjo" in testiranje novih formulacij koncnih elementov. Bibliography Armero, F., Ehrlich, D. 2004. An analysis of strain localization and wave propagation in plastic models of beams at failure. Comput. Meth. Appl. Mech. Eng., 193: 3129-3171. Armero, F., Ehrlich, D. 2006a. Finite element methods for the multi-scale modeling of softening hinge lines in plates at failure. Comput. Meth. Appl. Mech. Eng., 195: 1283-1324. Armero, F., Ehrlich, D. 2006b. Numerical modeling of softening hinges in thin Euler-Bernoulli beams. Computers and Structures, 84: 641-656. Auricchio, F., Taylor, R.L. 1994. A generalized elastoplastic plate theory and its algorithmic implementation. International Journal for Numerical Methods in Engineering, 37: 2583-2608. Basar, Y., Itskov, M. 1999. Constitutive model and finite element formulation for large strain elasto-plastic analysis of shells. Computational Mechanics, 23: 466-481. Bazant, Z.P., Oh, BH. 1983. Crack band theory for fracture of concrete. Materials and Structures, 16: 155-177. Bazant, Z.P., Belytschko, T., Chang, T.P. 1984. Continuum model for strain softening. Journal of Engineering and Mechanics 110, 12: 1666-1691. Babuska, I., Melenk, J.M. 1997. The partition of unity method. International Journal for Numerical Methods in Engineering, 40: 727-758. Bathe, K.J., Dvorkin, E. 1985. A four-node plate bending element based on Mindlin-Reissner plate theory and a mixed interpolation. Int J Numer Meth Engng, 21: 367-383. Bittencourt, T.N., Wawrzynek, P.A., Ingraffea, A.R., Sousa, J.L. 1996. Quasi-automatic simulation of crack propagation for 2D LEFM problems. Engineering Fracture Mechanics, 55 :321-334. Bohinc, U., Ibrahimbegovic, A. 2005. Robustni koncni elementi za plosce. Zbornik Kuhljevi dnevi 2005, 33-40. Bohinc, U., Ibrahimbegovic, A., Brank, B. 2009. Model adaptivity for finite element analysis of thin and thick plates based on equilibrated boundary stress resultants. Engineering Computations, 26: 69-99. Brancherie, D., Ibrahimbegovic, A. 2008. Novel anisotropic continuum-discrete damage model capable of representing localized failure of massive structures. Part I: theoretic formulation and numerical implementation. Engineering Computations 26, 1/2: 100-27. BranCherie, D. 2003. Modeles Continus et disCrets pour les problemes de loCalisation et de rupture fragile et/ou duCtile. DoCtoral Thesis, CaChan, ENS CaChan: 180 p. Brank, B. 1994. Velike deformaCije lupin pri nelinearnih materialnih modelih. Doktorska disertaCija. Ljubljana, Univerza v Ljubljani, Fakulteta za gradbenistvo in geodezijo: 201 str. Brank, B. 2005. Nonlinear shell models with seven kinematiC parameters. Comput. Meth. Appl. MeCh. Eng., 194: 2336-2362. Brank, B., IbrahimbegoviC, A. 2001. On the relation between different parametriza-tions of finite rotations for shells. Engineering Computations, 18: 950-973. Brank, B., PeriC, D., DamjaniC, F.B. 1995. On implementation of a nonlinear four node shell finite element for thin multilayered elastiC shells. Comput. meCh., 16: 341-359. Brank, B., PeriC, D., DamjaniC, F.B. 1997. On large deformation of thin elasto-plastiC shells: Implementation of a finite rotation model for quadrilateral shell elements. International Journal for NumeriCal Methods in Engineering, 40: 689-726. CamaCho, G.T., Ortiz, M. 1996. Computational modeling of impaCt damage in brittle materials. International Journal of Solids and StruCtures, 33: 2899-2938. Carter, B.J., Wawrzynek, P.A., Ingraffea, A.R. 2000. Automated 3-D CraCk growth simulation. International Journal for NumeriCal Methods in Engineering, 47: 229253. Coleman, B.D., Hodgon, M.L. 1985. On shear bands in duCtile materials. Archive for Rational MeChaniCs and Analysis, 90: 219-247. Crisfield, M.A., Peng, X. 1992. Effitient nonlinear shell formulations with large rotations and plastitity. In: D.R.J. Owen et al. Computational plastitity: models, software and appliCations, Part 1, Pineridge Press, Swansea, p. 1979-1997. Crisfield, M.A. 1981. Finite element analysis for Combined material and geometriC nonlinearities. In: W. WundeliCh (ed.), Nonlinear Finite Element Analysis in StruCtural MeCaniCs. Springer-Verlag, p. 325-338. Darvall, P.L., Mendis, P.A. 1985. ElastiC-plastiC-softening analysis of plane frames. Journal of StruCtural Engineering, 111: 871-888. de Borst, R., Sluys, L.J. 1991. LoCalization in a Cosserat Continuum under statiC and dynamiC loading Conditions. Computer Methods in Applied MeChaniCs and Engineering, 90: 805-827. de Borst, R., Remmers, J.J.C., Needleman, A., Abellan, M.A., 2004. DisCrete vs smeared CraCk models for ConCrete failure: bridging the gap. International Journal for NumeriCal and AnalytiCal Methods in GeomeChaniCs, 28: 583-607. DujC, J., Brank, B. 2006. Limit load analysis of steel plates. Gradbeni vestnik, 55: 288-299. DujC, J., Brank, B. 2008. On stress resultant plastitity and visCoplastitity for metal plates. Finite Elements in Analysis and Design, 44: 174-185. Dujč, J., Brank, B., Ibrahimbegovič, A. 2009. Multi-sčale čomputational model for failure analysis of metal frames that inčludes softening and ločal bučkling, Comput. Methods Appl. Mečh. Engrg. 199, 21: 1371-1385. Dujč, J., Brank, B., Ibrahimbegovič, A., Brančherie, B. 2010. An embedded čračk model for failure analysis of čončrete solids, Computers and Cončrete, submitted January 2010. Dvorkin, E., Cuitino, A., Goia, G. 1990. Finite elements with displačement interpolated embedded ločalization lines insensitive to mesh size and distortions. International Journal for Numeričal Methods in Engineering, 30: 541-564. Ehrličh, D., Armero, F. 2005. Finite element methods for the analysis of softening plastič hinges in beams and frames. Comput. Mečh., 35: 237-264. Euročode 2: Design of čončrete stručtures - Part 1-1: General rules and rules for buildings, 2004. Euročode 3: Design of shell stručtures - Part 1-6: General rules - Supplementary rules for the shell stručtures, 2007. Fajfar, P., Dolsek, M., Marusič D., Stratan, A. 2006. Pre- and post-test mathemat-ičal modelling of a plan-asymmetrič reinforčed čončrete frame building. Earthquake Engng. and Stručt. Dyn., 35: 1359-1379. Hobbit, Karlsson, Sorensen. 2007. Abaqus manuals. Hillerborg, A., Modeer, M., Petersson, P.E. 1976. Analysis of čračk formation and čračk growth in čončrete by means of fračture mečhaničs and finite elements. Cement and Cončrete Researčh, 6: 773-782. Ibrahimbegovič, A. 1993. Quadrilateral finite elements for analysis of thičk and thin plates. Computer Methods in Applied Mečhaničs and Engineering, 110: 195-209. Ibrahimbegovič, A. 2006. Mečanique non lineare des solides deformables: formulation theorique et implatantion elements finis. Paris. Hermes Sčienče-Lavoisier: 558 p. Ibrahimbegovič, A. 2009. Nonlinear Solid Mečhaničs: Theoretičal Formulations and Finite Element Solution Methods. Dordrečht. Springer: 574 p. Ibrahimbegovič, A., Brančherie, D. 2003. Combined hardening and softening čon-stitutive model of plastičity: prečursor to shear slip line failure. Comput. Mečh., 31: 88-100. Ibrahimbegovič, A., Frey, F. 1993. An effičient implementation of stress resultant plastičity in analysis of Reissner-Mindlin plates. International Journal for Numeričal Methods in Engineering, 36:303-320. Ibrahimbegovič, A., Frey, F. 1993. Stress resultant finite element analysis of reinforčed čončrete plates. Engineering Computations, 10: 15-30. Ibrahimbegovič, A., Frey, F. 1994. An effičient approačh to servičeability analysis and ultimate load design of reinforčed čončrete plates. Computational modelling of čončrete stručtures (H. Mang, N. Bičanič,R. de Borst, uredniki), Pineridge Press, 875-884. IbrahimbegoviC, A., Melnyk, S. 2007. Embedded disContinuity finite element method for modelling of loCalized failure in heterogeneous materials with struCtured mesh: an alternative to extended finite element method. Comput. MeCh., 40: 149-155. IbrahimbegoviC, A., Wilson, E.L. 1991. A modified method of inCompatible modes. Commun. Appl. Numer. Meth., 7: 187-194. IbrahimbegoviC, A., Frey, F., Sarf, J.L. 1992. Limit load analysis of plates with partiCular referenCe to steel and reinforced ConCrete. Engineering Modeling 5, 3-4: 75-82. IbrahimbegoviC, A., Gharzeddine, F., Chorfi, L. 1998. ClassiCal plastitity and vis-CoplastiCity models reformulated: TheoretiCal basis and numeriCal implementation. International Journal for NumeriCal Methods inEngineering, 42: 1499-1535. Ingraffea, A.R., Saouma, V. 1985. NumeriCal modelling of disCrete CraCk propagation in reinforCed and plain ConCrete. In FraCture meChaniCs of ConCrete: StruCtural appliCation and numeriCal CalCulation, eds. Sih GC, Tomasso D, 171-225. Jirasek, M., Zimmermann, T. 2001. Embedded CraCk model: I. BasiC formulation. International journal for numeriCal methods in engineering, 50: 1269-1290. Jirasek, M. 1997. AnalytiCal and numeriCal solutions for frames with softening hinges. ASCE J. Eng. MeCh. Din., 123: 8-14. Jirasek, M. 2000. Comparative study on finite elements with embedded disContinuities. Comput. Methods Appl. MeCh. Engrg., 188: 307-330. Khan, A.S., Huang, S. 1995. Continuum Theory of Plastitity. New York. John Wiley: 440 p. Kleiber, M., KowalCzyk, P. 1996. Sensitivity analysis in planestress elasto-plastitity and elasto-visCoplastiCity. Computer Methods in Applied MeChaniCs and Engineering, 137: 395-409. KoeChlin, P., Potapov, S. 2007. Global Constitutive Model for Reinforced ConCrete Plates. Journal of engineering meChaniCs 133, 3: 257-266. KojiC, M., Bathe, K.J. 2005. InelastiC Analysis of Solids and StruCtures. Berlin. Springer: 414 p. KorelC, J. 1997. AutomatiC generation of finite element Code by simultaneous optimization of expressions. Theor Comput Sti, 187: 231-248. KorelC, J. 2002. Multi-language and multi-environment generation of nonlinear finite element Codes, Engineering with Computers 18, 4: 312-327. KorelC, J. 2007a. ACeFem. http://www.fgg.uni-lj.si/SymeCh, 2007. KorelC, J. 2007b. ACeGen. http://www.fgg.uni-lj.si/SymeCh, 2007. KuCerova, A., BranCherie, D., IbrahimbegoviC, A., Zeman, J., Bittnar, Z. 2009. Novel anisotropiC Continuum-disCrete damage model Capable of representing loCalized failure of massive struCtures. Part II: identifiCation from tests under heterogeneous stress field. Engineering Computations, 26: 128-144. Linder, C. 2007. New finite elements with embedded strong discontinuities for the modeling of failure in solids. Doctoral Thesis. Berkeley. University of California: 321 p. Linder, C., Armero, F. 2007. Finite elements with embedded strong discontinuities for the modeling of failure in solids. Int. J. Numer. Meth. Engng., 72: 1391-1433. Lubliner, J. 1990. Plasticity Theory. New York. MacMillian: 495 p. Manzoli, O.L., Shing, P.B. 2006. A general technique to embed non-uniform discontinuities into standard solid finite elements. Computers and Structures, 84: 742-757. Marusich, T.D., Ortiz, M. 1995. Modelling and simulation of high-speed machining. International Journal for Numerical Methods in Engineering, 38: 3675-3694. Melenk, J.M., Babuska, I. 1996. The partition of unity finite element method: basic theory and applications. Computer Methods in Applied Mechanics and Engineering, 139: 289-314. Mosler, J., Bruhns, O.T. 2004. A 3D anisotropic elastoplastic-damage model using discontinuous displacement fields. Int. J. Numer. Meth. Engng., 60: 923-948. Mosler, J., Meschke, G. 2003. 3D modelling of strong discontinuities in elastoplastic solids: fixed and rotating localization formulations. Int. J. Numer. Meth. Engng, 57: 1553-1576. Mosler, J. 2005. On advanced solution strategies to overcome locking effects in strong discontinuity approaches. Int. J. Numer. Meth. Engng., 63: 1313-1341. Moy, S.J. 1996. Plastic methods for stell and concrete structures. London. MacMil-lian: 271 p. Needleman, A. 1987. A continuum model for void nucleation by inclusion debonding. Journal of Applied Mechanics, 54: 525-531. Ngo, D., Scordelis, A.C. 1967. Finite element analysis of reinforced concrete beams. Journal of the American Concrete Institute, 64: 152-163. Nielsen, M.P. 1984. Limit analysis and concrete plasticity. Englewood Cliffs. Prentice-Hall: 420 p. Oliver, J. 1995. Continuum modelling of strong discontinuities in solid mechanics using damage models. Computational Mechanics, 17: 49-61. Oliver, J. 1996. Modelling strong discontinuities in solid mechanics via strain softening constitutive equations. Parts 1: Fundamentals. International Journal for Numerical Methods in Engineering, 39: 3575-3600. Oliver, J. 1996. Modelling strong discontinuities in solid mechanics via strain softening constitutive equations. Part 2: Numerical simulation. International Journal for Numerical Methods in Engineering, 39: 3601-3623. Oliver, J., Simo, J.C. 1994. Modelling strong discontinuities by means of strain softening constitutive equations. In Proceedings EURO-C Computer Modelling of Concrete Structures, eds. Mang H, Bicanic N, de Borst R, 363-372. Ortiz, M., Pandolfi, A. 1999. Finite deformation irreversible Cohesive elements for three-dimensional CraCk-propagation analysis. International Journal for NumeriCal Methods in Engineering, 44: 1267-1282. Ortiz, M., Quigley, J.J. 1991. Adaptive mesh refinement in strain loCalization problems. Computer Methods in Applied MeChaniCs and Engineering, 90: 781-804. Ortiz, M., Leroy, Y., Needleman, A. 1987. A finite element method for loCalized failure analysis. Computer Methods in Applied MeChaniCs and Engineering, 61: 189-214. Park, R., Gamble, W. L. 2000. ReinforCed ConCrete slabs. Wiley. Petersson, P.E. 1981. CraCk growth and development of fraCture zones in plain ConCrete and similar materials. Report No. TVBM-1006. Lund, Sweden, University of Lund, Division of Building Materials. Powell, G.H. 1986. Nonlinear struCtural analysis by Computer Code INSA, UC Berkeley reports SEMM 86-15. Press, W.H., Teukolsky, S.A., Vetterling, W.T., Flannery, B.P. 1992. NumeriCal ReCipes in Fortran. 2nd edition, Cambridge University Press. RadosavljeviC, Z., BajiC, D. 1990. Armirani beton knjiga 3, Elementi armiranobetonskih konstrukCija. Beograd. Gradevinska knjiga: 431 str. Rashid, R. 1968. Analysis of prestressed ConCrete pressure vessels. Nudear Engineering and Design 7, 4: 334-355. Reddy, J.N. 2004. An introduCtion to nonlinear finite element analysis. Oxford. Oxford University Press: 463 p. Rots, J.G., Nauta, P., Kusters, G.M.A., Blaauwendraad, J. 1985. Smreared CraCk approaCh and fraCture loCalization in ConCrete. Heron 30, 1: 1-48. SanCho, J.M., Planas, J., Cendon, D.A., Reyes, E., Galvez, J.C. 2007. An embedded CraCk model for finite element analysis of ConCrete struCture. Engineering FraCture MeChaniCs, 74: 75-86. Save, M.A., Massonet, C.E. 1972. PlastiC analysis and design of plates, shells and disks. Amsterdam. Nort-Holland: 478 p. SawCzuk, A. 1989. MeChaniCs and Plastitity of StruCtures. ChiChester. Ellis Hor-wood: 203p. Shi, G., Voyiadjis, G.Z. 1992. A simple non-layered finite element for the elasto-plastiC analysis of shear flexible plates. International Journal for NumeriCal Methods in Engineering, 33: 85-100. Simo, J.C., Hughes, T.J.R. 1998. Computational Inelastitity. New York. Springer: 392 p. Simo, J.C., Kennedy, J.G. 1992. On a stress resultant geometriCally exaCt shell model. Part V. Nonlinear plastitity: formulation and integration algorithms. Computer Methods in Applied MeChaniCs and Engineering, 96: 133-171. Simo, J.C., Oliver, J. 1994. A new approačh to the analysis and simulation of strong disčontinuities. In Fračture and Damage in Quasibrittle Stručtures, eds. Bazant ZP, Bittnar Z, Jirasek M, Mazars J, 2-6: 25-39, Boundary Row, London. Simo, J.C., Oliver, J., Armero, F. 1993. An analysis of strong disčontinuities in-dučed by strain-softening in rate independent inelastič solids. Computational Me-čhaničs,12: 277-296. Skallerud, B., Myklebust, L.I., Haugen, B. 2001. Nonlinear response of shell stručtures: effečts of plastičity modelling and large rotations. Thin-Walled Stručtures, 39: 463-482. Tvergaard, V. 1990. Effečt of fibre debonding in a whisker-reinforčed metal. Material Sčienče and Engineering A, 125: 203-213. Voyiadjis, G.Z., Woelke, P. 2006. General non-linear finite element analysis of thičk plates and shells. International Journal of Solidsand Stručtures, 43: 2209-2242. Wačkerfuss, J. 2008. Effičient finite element formulation for the analysis of ločalized failure in beam stručtures. Int. J. Numer. Meth. Engng., 73: 1217-1250. Wilson, E.L. 2002. Three dimensional statič and dynamič analysis of stručtures. Berkley. Computers and Stručtures, Inč.: 423 p. Xu, X.P., Needleman, A. 1994. Numeričal simulations of fast čračk growth in brittle solids. Journal of the Mečhaničs and Physičs of Solids, 42: 1397-1434. Zahlten, W. 1993. A čontribution to the Physičally and Geometričally Nonlinear Computer Analysis of General Reinforčed Cončrete Shells. Ruhr - Universitaet Bočhrum, 212 - 215. Zeng, Q., Combesčure, A., Arnaudeau, F. 2001. An effičient plastičity algorithm for shell elements appličation to metal forming simulation. Computers and Stručtures, 79: 1525-1540.