Downloaded from orbit.dtu.dk on: Jan 14, 2019 Development of a user element in ABAQUS for modelling of cohesive laws in composite structures Feih, Stefanie Publication date: 2006 Document Version Publisher's PDF, also known as Version of record Link back to DTU Orbit Citation (APA): Feih, S. (2006). Development of a user element in ABAQUS for modelling of cohesive laws in composite structures. Denmark. Forskningscenter Risoe. Risoe-R, No. 1501(EN) General rights Copyright and moral rights for the publications made accessible in the public portal are retained by the authors and/or other copyright owners and it is a condition of accessing publications that users recognise and abide by the legal requirements associated with these rights.  Users may download and print one copy of any publication from the public portal for the purpose of private study or research.  You may not further distribute the material or use it for any profit-making activity or commercial gain  You may freely distribute the URL identifying the publication in the public portal If you believe that this document breaches copyright please contact us providing details, and we will remove access to the work immediately and investigate your claim. Risø-R-1501(EN) Development of a user element in ABAQUS for modelling of cohesive laws in composite structures Stefanie Feih Risø National Laboratory Roskilde Denmark January 2005 Author: Stefanie Feih Risø-R-1501(EN) Title: Development of a user element in ABAQUS for modelling of January 2005 cohesive laws in composite structures Department: AFM Abstract (max. 2000 char.): ISSN 0106-2840 ISBN 87-550-3410-1 The influence of different fibre sizings on the strength and fracture toughness of composites was studied by investigating the characteristics of fibre cross-over bridging in DCB specimens loaded with pure bending moments. These tests result in bridging laws, which are obtained by simultaneous measurements of the crack growth resistance and the end opening of the notch. The advantage of this method is that these bridging laws represent Contract no.: material laws independent of the specimen geometry. However, the adaption of the experimentally determined shape to a numerically valid model shape is not straight forward, and most existing Group's own reg. no.: publications consider theoretical and therefore simpler softening 16150043-00 shapes. In this article, bridging laws were implemented into an interface element in the UEL user subroutine in the finite element code ABAQUS. Comparison with different experimental data Sponsorship: points for crack opening, crack length and crack shape show the sensitivity of these results to the assumed bridging law shape. It is furthermore shown that the numerical predictions can be used to Cover : improve the bridging law fit. One shape with one adjustable parameter then fits all experimental data sets. Pages: 52 Tables: 3 Figures: 26 References: 15 Risø National Laboratory Information Service Department P.O.Box 49 DK-4000 Roskilde Denmark Telephone +45 46774004 [email protected] Print: Pitney Bowes Management Services Denmark A/S, Fax +45 46774013 2005 www.risoe.dk Contents 1 Introduction 5 2 ModeIbridginglawmeasurement 5 3 Experimentalresults 6 4 Implementationoftheuserelement 8 4.1 Backgroundandbasicequations 9 4.2 Unloading 12 4.3 Statevariablesfortheanalysis 13 4.4 Numericalintegration 13 5 Numericalmodel 14 5.1 Symmetrichalfmodelversusfullmodel:boundaryconditions 14 5.2 Applicationofpuremomentbendingbydisplacementcontrol 14 5.3 Bridginglawadjustmentsfornumericalmodelling 17 6 Numericalresults 20 6.1 Comparisonofnumericalintegrationprocedures 20 6.2 Numericalbridgingresultswith(cid:2)nitestressvalue 21 6.3 Numericalbridgingresultswithpower-lawincrease 22 6.4 Adjustmentofbridginglaw 23 7 Summary 29 A Implementationoftheplaneinterfaceelement 31 B Quadraticshapefunctionsandderivatives 33 C Integrationpointsandweights 35 D ABAQUScodingforthequadraticlineelement 37 E Cohesiveelementveri(cid:2)cation 49 3 Acknowledgements I would like to thank Bent F. Srensen for readily available input and suggestions regarding cohesivezone models. A special thanks is going to Lars Pilgaard Mikkelsen for all the help anddiscussionsabout(cid:2)niteelementanalysisandABAQUS. 4 1 Introduction Forglass(cid:2)brecomposites,theinterfacialpropertiesarecontrolledbythesizing,whichisap- plied to the glass (cid:2)bres during manufacture.For the same matrix system, a changeofsizing results in changes of these properties, thereby in(cid:3)uencingthe mechanical properties such as strengthandfracturetoughness.Theconceptofstrengthisusedforcharacterisingcrackinitia- tionincompositedesign,whilefracturetoughnessdeterminescrackgrowthanddamagedevel- opment.In modeI crackgrowthin unidirectional(cid:2)bercomposites,(cid:2)bre cross-overbridging occursduringcrackingalongthe(cid:2)berdirection.Thisfailuremodeplaysanimportantroledur- ingdelaminationof(cid:2)brecompositesandsplittingcracksaroundholesandnotches.The(cid:2)bre bridgingzone must be modelledas a discrete mechanismon its own;failure is not just con- trolledbythecrackingatthecracktip.Thefailureprocesscanbedescribedbyabridginglaw, whichde(cid:2)nestherelationshipbetweenthecrackopeningdisplacementandthelocalbridging tractionsresultingfromthebridgingligaments.Cohesivelawsweremeasuredexperimentally in previouswork.This report derivesthe necessary basics and equationsto implementthese lawsintothecommercial(cid:2)niteelementcodeABAQUSwithacohesiveuserelement.Different numericaladjustmentsofthebridginglawarediscussedindetail.Crackaspects,suchascrack openingshapeandthein(cid:3)uenceofbridginglawparameters,arestudiedbasedonthenumerical results.Itisfurthermoreofinteresttoidentifytheexperimentalmeasurementswhichshowthe highestsensitivitywithrespecttothebridginglawshape. 2 Mode I bridging law measurement The approachfor the measurementsof bridginglaws is based on the application of the path independentJ integral[1],andhasbeenusedrecentlytodeterminethebridgingcharacteristics of unidirectional carbon (cid:2)bre/ epoxy composites [2] and glass (cid:2)bre composites [3]. A sym- metricDCBspecimenisloadedwithpurebendingmomentsM (Figure1)underpuremodeI. Thisspecimenisoneofthefewpracticalspecimengeometries,forwhichtheglobalJ integral (i.e.theintegralevaluatedaroundtheexternalboundariesofthespecimen)canbedetermined analytically[1]: M2 J =12(1(cid:0)(cid:23)13(cid:23)31)b2H3E11 (1) E11 is theYoung’smodulusreferringtothematerialdirections,(cid:23)13 and(cid:23)31 arethemajor andminorPoisson’sratio,bisthewidthandH thebeamheight. M x 2 H+Du* 2H 2 x 1 M Figure1.DCBspecimenwithpurebendingmoment Nowconsiderthespecimenhavingacrackwithbridging(cid:2)bresacrossthecrackfacesnear thetip.Theclosurestress (cid:27) (x2-direction)canbeassumedtodependonlyonthelocalcrack 5 opening (cid:14), i.e. the crack grows in pure mode I. The bridging law (cid:27) = (cid:27)((cid:14)) is then taken as identicalateachpointalongthebridgingzone.Since(cid:2)breswillfailwhenloadedsuf(cid:2)ciently, weassumetheexistenceofacharacteristiccrackopening(cid:14)0,beyondwhichtheclosuretraction vanishes. Shrinkingthe pathofthe J integraltothe crackfaces andaroundthecracktip [4] gives (cid:14)(cid:3) J = (cid:27)((cid:14))d(cid:14)+J ; (2) tip Z0 whereJ istheJ integralevaluatedaroundthecracktip(duringcrackingJ isequaltothe tip tip fractureenergyofthetip,J0).Theintegralistheenergydissipationinthebridgingzoneand (cid:14)(cid:3) istheend-openingofthebridgingzoneatthenotchroot. Connectionis made to the overallR-curvesas follows.By de(cid:2)nitionJ is the valueof J R duringcrackgrowth.Initially,thecrackisunbridged.Thus,byEq.(2),crackgrowthinitiates when JR = Jtip = J0. As the crack grows,JR increases in accordancewith Eq. (2).When (cid:3) theendopeningofthebridgingzone(cid:14) reaches(cid:14)0,theoverallR-curveattainsitssteadystate valueJ . ss ThebridginglawcanbedeterminedbydifferentiatingEq.(2)[4]. @J (cid:3) R (cid:27)((cid:14) )= (3) @(cid:14) Theappliedmomentandtheendopeningofthebridgingzone(cid:1)u(cid:3)2arerecorded.Assumingthat (cid:14)(cid:3) (cid:1)u(cid:3)2,where(cid:1)u(cid:3)2 isthenotchopeningmeasuredattheneutralaxisoftheDCBspecimen (cid:25) (seeFigure1),thebridginglawcanbedetermined.Thisapproachmodelsthebridgingzoneas a discrete mechanismonits own.Contraryto crackgrowthresistance curves(R-curves),the bridginglawcanbeconsideredamaterialpropertyanddoesnotdependonspecimensize[2]. The test above has been modi(cid:2)ed with two different bending moments to result in mixed modetesting[5].Inthiscase,mixedbridginglawscanbemeasured.Thishasnotbeenunder- takenforthecurrentmaterialselection. 3 Experimental results Recently,wehave,bytheuseofaJintegralbasedapproach,measuredthebridginglawsunder modeIfractureduringtransversesplittingofunidirectionalglass-(cid:2)ber/epoxyandglass-(cid:2)ber/ polyestercompositeswithdifferentinterfacecharacteristics[3].Withincreasingappliedmo- ment,crackpropagationtookplace.Fibrecross-overbridgingdevelopedinthezonebetween thenotchandthecracktip. J iscalculatedaccordingtoEq.(1).Thespecimenwidthbwas5mmwithabeamheightof R H=8mm.Assumingthattheunidirectionalcompositeistransverselyisotropic,thefollowing elasticcompositedatawereappliedforEq.(1)aspreviouslymeasured:E11;epoxy=41.5GPa, E33= 9.2 GPa, E11;polyester = 42 GPa, E33;polyester= 10 GPa and (cid:23)13=0.3 (assumption). The analyticalfunction (cid:14)(cid:3) 12 (cid:3) JR((cid:14) )=J0+(cid:1)Jss (cid:14)0 (4) (cid:18) (cid:19) wasfoundto(cid:2)tallexperimentaldatacurvesofcrackgrowthresistanceversuscrackopening well, resulting in curve (cid:2)ts as shown in Figure 2. J0 is the initial value of the experimental curve and equal to the fracture energy of the tip during crack growth, while (cid:1)J , which is ss equalto(Jss J0),istheincreaseincrackgrowthresistance.SłrensenandJacobsen[2]found (cid:0) thatthesamefunction(cid:2)tthedataofcarbon(cid:2)brecompositesystemswell. The experimentalvalues for the bridginglaws are given in Table 1. The starting value J0 indicatesthepointofcrackgrowthinitiationandcaneasilybedeterminedduringtheexperi- 2 ment.Thehighestvalueof345J/m wasobservedforthesizingB/epoxysystem.Thecrack 6 6000 2m] J/ e [ c n 4000 a st si e r h wt o 2000 gr A / epoxy k B / epoxy c a A / polyester r C B / polyester 0 0 2 4 6 8 Crack opening [mm] Figure2.Comparisonofcrackgrowthresistance Table1.Experimentalvaluesforthebridginglawfordifferentcompositesystems Compositesystem J0 [J/m2] (cid:1)Jss[J/m2] (cid:14)0 [mm] sizingA/epoxy 321 40 4000 1000 2.0 0.2 (cid:6) (cid:6) (cid:6) sizingB/epoxy 345 30 3700 500 2.0 0.2 (cid:6) (cid:6) (cid:6) sizingA/polyester 150 20 3800 5.5 (cid:6) sizingB/polyester 120 30 >4100 >5.0 (cid:6) 2 initiationvalueissigni(cid:2)cantlylowerforthesizingB/polyestersystemwith120J/m ,which isalsorelatedtoasigni(cid:2)cantlylowertransversestrengthofnotmorethanhalfthestrengthof theothercomposites[3].Infact,therelationbetweentransversestrength(measuredwiththe transverse bending test [3]) and crack initiation for the different composite systems is (cid:2)tted well by a linear relationship as seen in Figure 3. This veri(cid:2)es the assumption that the crack initiationvaluefortheDCB test is controlledbythestrengthofthe(cid:2)ber-matrixbond.It can furthermorebeseen thatbothaxesshowaboutthesameratiolowstrengthandhighstrength ofnearlyafactor3;however,thetransversebendingtestresultsinmuchlowerstandarddevi- ations.ThehigherstandarddeviationfortheDCB fractureinitiationis mostlikelyexplained by specimen-to-specimendifferences of the manufacturednotch, which in(cid:3)uences the crack initiationprocess. Theendopeningvalue(cid:14)0 attheonsetofsteady-statecrackingwasdeterminedtobe2mm fortheepoxysystems.ForthesizingB/polyestersystem,steady-statecrackingcouldnotbe determinedwith thepresentspecimens,as the(cid:2)brescontinuedto bridgethe wholelengthof thecrackafterthemaximummeasurablenotchopeningof5mmwasobtained.Sincenoupper boundwasfoundfor(cid:1)J ,thisbridgingbehaviourwastermed’in(cid:2)nitetoughening’. ss Differentiatingequation(4)resultsinthebridginglaw (cid:0)1 (cid:1)J (cid:14) 2 (cid:27)((cid:14))= ss ; for 0<(cid:14) <(cid:14)0; (5) 2(cid:14)0 (cid:14)0 (cid:18) (cid:19) where(cid:1)J istheincreaseincrackgrowthresistanceduetobridging(fromzerotosteadystate ss bridging),and(cid:14)0 isthecrackopeningwherethebridgingstressvanishes. Thebridginglawsforthedifferent(cid:2)bresystemsarecomparedinFigure4.Thebridginglaw can be considereda material property[2,6] and is in an accessible form for implementation into(cid:2)niteelementcodes. 7 a] 100 P M h [ gt 75 n e r st g n 50 ndi R2= 0.87 e b e s 25 r e v s n a r T 0 0 100 200 300 400 Crack initiation J [J/m2] 0 Figure3.RelationshipbetweencrackinitiationvalueJ0andthetransversecompositestrength forthedifferentcompositesystems 3 A / epoxy B / epoxy a] A / polyester P M B / polyester s [ 2 s, s e r st g n 1 gi d ri B 0 0 2 4 6 8 Crack opening [mm] Figure4.Comparisonofresultingbridginglaws 4 Implementation of the user element There are a variety of possible methods for implementing cohesive laws within commercial (cid:2)nite element programs. The most versatile is the development and programming of cohe- sive elements[7(cid:150)10].Theseelements arein most cases de(cid:2)nedwith zero thickness andpre- scribestressesbasedontherelativedisplacementofthenodesoftheelement.Similarworkhas alsobeenundertakenwithspringelements(force-openingrelation),althoughinthiscasethere mightbesimpli(cid:2)cationsrequiredwhencalculatingtheequivalentnodalspringforcesfromthe surrounding elements. The procedure is not straight forward when springs are connected to elementswithnon-linearshapefunctions,suchas8-nodedelements[11]. 8 4.1 Background and basic equations Figure5showstheinterfaceelementsin2-Dand3-D.Theinterfaceelementismadeupoftwo quadraticlineelements(a),ortwoquadratic2-Dplaneelements(b).Theseelementsconnect thefacesofadjacentelementsduringthefractureprocess.Theimplementationisbasedon[10]. Quadraticelements are chosen as the cantileverbeam mostly deformsunder bending,which isbestmodelledbyquadraticsolidelements.Thenodesoftheinterfaceelementneedto(cid:2)tto theseelements.ThenodenumberingischosenaccordingtoABAQUSconventionsofquadratic 2-Dand3Dsolids.Theelementscanalsobederivedinlinearformbysubstitutingthequadratic shapefunctionswithlinearonesintheappropriateequations.Consequently,thenodenumbers anddegreeoffreedomsfortheelementswilldecrease. 12 15 11 h 4 7 3 16 14 4 6 5 x 8 13 10 9 6 x y 1 3 2 y 1 5 2 x x z Figure5.(a)Quadraticlineinterfaceelement(3nodepairs)and(b)Quadraticplaneinterface element(8nodepairs) Thetwosurfacesoftheinterfaceelementinitiallylietogetherintheunstresseddeformation state(zerothickness)andseparateastheadjacentelementsdeform.Therelativedisplacements of the element faces create normal and shear displacements, which in turn generate element stresses depending on the constitutive equations (stress - opening relations) of the material. Theconstitutiverelationshipwasderivedintheexperimentalsection,andisindependentofthe elementformulation. Theimplementationofageneralinterfaceelementisexplainedinthefollowing.Detailsare givenforthespecialcaseofthequadraticlineelementfor2-Dsimulations,whiletheequivalent formulationsforthequadraticplaneelementfor3-DsimulationsareprovidedinAppendixA. Thelineinterfaceelementhas12(2x6)degreesoffreedom.Thevector(12 1)ofthenodal (cid:2) displacementsintheglobalcoordinatesystemisgivenas: d = d1 d1 d2 d2 d6 d6 T (6) N x y x y (cid:1)(cid:1)(cid:1) x y Theorderfo(cid:0)llowstypicalABAQUS’con(cid:1)ventions,andthisisconsideredinthederivedformu- lationsbelow. Theopeningoftheinterfaceelementisde(cid:2)nedasthedifferenceindisplacementsbetween thetopandbottomnodes: (cid:1)u= u top u bot; (7) f g (cid:0)f g thereby leading to the following de(cid:2)nition of the interface opening (cid:1)u in terms of nodal N displacementsofpairednodes: (cid:1)uN =(cid:8)dN =[ I6(cid:2)6 I6(cid:2)6]dN (8) (cid:0) j whereI6(cid:2)6 denotesaunitymatrixwith6rowsandcolumns.uN isa6 1vector. (cid:2) From the nodal positions, the crack opening is interpolated to the integration points with the help of standard shape functions. Let N ((cid:24)) be the shape functions for the node pair i i (i = 1;2;3),where(cid:24) stands forthe local coordinateof the element with 1 (cid:24) 1. The (cid:0) (cid:20) (cid:20) 9