Development of a user element in ABAQUS for modelling of cohesive laws in composite structures

54 Pages·2017·1.77 MB·English

Downloaded from orbit.dtu.dk on: Jan 12, 2023

Development of a user element in ABAQUS for modelling of cohesive laws in
composite structures
Feih, S.
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. Risø National Laboratory. 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 [phone]
[email protected]
Print: Pitney Bowes Management Services Denmark A/S,
Fax [phone]
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