Research Article Combined Numerical-Statistical Analyses of Damage and Failure of 2D and 3D ...

13 Pages·2015·4.5 MB·English

Hindawi Publishing Corporation
Mathematical Problems in Engineering
Volume 2015, Article ID 702563, 12 pages
http://dx.doi.org/10.1155/2015/702563
Research Article
Combined Numerical-Statistical Analyses of Damage and
Failure of 2D and 3D Mesoscale Heterogeneous Concrete
XiaofengWangandAndreyP.Jivkov
SchoolofMechanical,AerospaceandCivilEngineering,theUniversityofManchester,ManchesterM139PL,UK
CorrespondenceshouldbeaddressedtoXiaofengWang;[email protected]
Received6May2015;Accepted16July2015
AcademicEditor:EugeneKogan
Copyright©2015X.WangandA.P.Jivkov.ThisisanopenaccessarticledistributedundertheCreativeCommonsAttribution
License,whichpermitsunrestricteduse,distribution,andreproductioninanymedium,providedtheoriginalworkisproperly
cited.
Generationandpackingalgorithmsaredevelopedtocreatemodelsofmesoscaleheterogeneousconcretewithrandomlydistributed
elliptical/polygonal aggregates and circular/elliptical voids in two dimensions (2D) or ellipsoidal/polyhedral aggregates and
spherical/ellipsoidal voids in three dimensions (3D). The generation process is based on the Monte Carlo simulation method
whereintheaggregatesandvoidsaregeneratedfromprescribeddistributionsoftheirsize,shape,andvolumefraction.Acombined
numerical-statisticalmethodisproposedtoinvestigatedamageandfailureofmesoscaleheterogeneousconcrete:thegeometrical
modelsarefirstgeneratedandmeshedautomatically,simulatedbyusingcohesivezonemodel,andthenresultsarestatistically
analysed.Zero-thicknesscohesiveelementswithdifferenttraction-separationlawswithinthemortar,withintheaggregates,and
attheinterfacesbetweenthesephasesarepreinsertedinsidesolidelementmeshestorepresentpotentialcracks.Theproposed
methodologyprovidesaneffectiveandefficienttoolfordamageandfailureanalysisofmesoscaleheterogeneousconcrete,anda
comprehensivestudywasconductedforboth2Dand3Dconcreteinthispaper.
1.Introduction An extensive literature review recently carried out by
Wang et al. [3] shows that numerical image processing
Concrete is the most widely used construction material in
technique and parameterization modelling technique are
the world due to its good strength and durability relative
two main approaches to generate mesostructure models of
to its cost. Numerical simulations coupled with theory and
concrete.Althoughthemeshgeneratedfromimagesisclosest
experiment are considered to be an important tool for
to nature, it is presently very costly and time consuming
examiningthemechanicalbehaviourthroughcomputational
to generate sufficient scanned images and convert them to
materials science. Wittmann [1] proposed three levels of
−6 −3 meshformeaningfulstatisticalanalyses(e.g.,[4–6]).Inthe
observation: microlevel (10 m), mesolevel (10 m), and
0 parameterization modelling technique, both direct (e.g., [3,
macrolevel(10 m).Themicrolevelrepresentsthemostbasic
7–9])andindirectapproaches(e.g.,[10–13])tocharacterize
level, in which the internal structure of cement paste is
materialrandomheterogeneityhavebeenused.Theindirect
consideredandwherephysicalandchemicalforcesdominate.
approach is able to generate a large number of samples
Inthemesolevel,concreteisusuallyrepresentedasamaterial
made up of coarse aggregates, mortar with fine aggregates with ease but the key multiphasic parameters which could
and cement paste embedded inside, and interfaces between significantly influence the mechanical behaviour cannot be
coarseaggregatesandmortar.Atthemacroscale,concreteis considered.Thedirectapproachcantakeintoaccountmost
treatedasahomogenousandcontinuousmaterialinwhich of the key parameters such as shape, size, gradation and
the internal structure is not considered [2]. As concrete is spatial distribution of voids and aggregates, phase volume
generally used in large-sized structures and its dependence fractions,andaggregate-mortarinterfaces.Soamongallthe
of mechanical behaviours on mesostructures is significant, approaches,itseemsthatthedirectapproachisparticularly
mesoscalemodellingofconcretebecomesessentialandwas suitable for statistical analysis of concrete samples and was
investigatedinthisstudy. usedinthisstudy.
2 MathematicalProblemsinEngineering
Mortar Voids
Aggregates
(a) Glass ball concrete (b) Rolled aggregate concrete (c) Crushed aggregate concrete
Figure1:Tomographiccross-sectionalviewofthreeconcretesvaryingbycoarseaggregateshape(after[14]).
Identificationandgenerationofunitcellgeometryarea Table1:Aggregatesizedistribution[28].
vital step in the mesoscale analysis of concrete. Both shape
Sievesize Totalpercentage Totalpercentagepassing
andsizeofaggregateshavesignificantinfluencesonthestress
(mm) retained(%) (%)
distribution, crack initiation, and damage accumulation up
tothemacroscopicfailurewithintheconcretematerial[16]. 12.70 0 100
Mostoftheexistingalgorithmsformesostructuregeneration 9.50 39 61
of concrete only consider random aggregates but neglect 4.75 90 10
voids[2,7,17–19].However,theX-rayCTimages[4,14,20]
2.36 98.6 1.4
clearlyshowthatvoidsexistinconcreteatthisscaleandhave
severely adverse effects on the specimen strength [3, 9]. So
automatic generation of morphological details of materials
usingthecombinednumerical-statisticalmethod,leadingto
withrandomlydistributeddifferentshapeofaggregatesand
valuable statistical results that may help improve designs of
voidsposesnewchallenges.Inthispaper,acomputationally
concretematerials.
efficient and versatile in-house program of heterogeneous
material generator (HMG) has been developed considering
both2Dand3Drandomaggregatesandvoidsbasedonthe 2.2Dand3DHeterogeneousConcrete
prescribedparameters.
Severalnumericalmodelsforcrackpropagationatmeso- 2.1. Aggregate and Void Size Distribution. The size distribu-
scopic level have been used to study the heterogeneity of tionofaggregatesinconcreteisoftendescribedbytheFuller
concrete. Continuum based finite element models are the curve [16], which is discretised into a certain number of
mainapproachesemployedintheliterature[21–24],primar- segments. The aggregate size distributions found in Hirsch
ily based on cohesive zone model [25]. Other alternative [28]andsummarisedinTable1areusedinthisstudy.Inthis
approaches have also been developed, such as discrete ele- study, concrete is treated as a material consisting of coarse
mentmodel(DEM)[26]andlatticemodel[11,12].Thebiggest aggregates,mortarincorporatingfineaggregatesdissolvedin
disadvantage of DEM and lattice model is the difficulty it,voids,andinterfacesbetweenaggregatesandmortar.Here
in determining the equivalent model parameters, which is only coarse aggregates are explicitly modelled as mesoscale
relativelystraightforwardforcontinuumbasedfiniteelement features.Fornormalstrengthconcretecastwithamouldin
models[27].Thenumericalmodelusedinthisworkforcrack thelaboratory,coarseaggregatesusuallyrepresent30–50%of
propagation at mesoscopic level is based on the cohesive theconcretevolume.
zone model [3]. It is becoming more and more popular in From X-ray CT images, it was found that gravel aggre-
modelling crack propagation due to its simple formation gates have circular/elliptical shapes in 2D and spheri-
and easy implementation in the form of cohesive interface cal/ellipsoidal shapes in 3D, while crushed aggregates have
elements(CIEs). polygonal shapes in 2D and polyhedral shapes in 3D [14].
WiththepreprocessingfunctionalityinANSYS[29]and A number of voids can be clearly observed in concrete at
explicit solver in ABAQUS [30], nonlinear modelling of mesoscale from CT images (see Figure1). Therefore, voids
multiple crack propagation with a considerable number of with size ranging from 2 to 4mm as presented in [31] are
concrete samples has been performedby high performance included in the mesostructure of concrete in this study.
computinginthisstudy.Statistically,analysisofmechanical Coarseaggregatesareconsideredtohaveelliptical/polygonal
behaviourofboth2Dand3Dspecimenshasbeenconducted shapes in 2D and ellipsoidal/polyhedral shapes in 3D, and
MathematicalProblemsinEngineering 3
Loop over all the previous generated inclusions until prescribed volume fraction is achieved
No
Hierarchy I: radial distance check Save and go for next inclusion
Yes
Hierarchy II: check if the centres of any previous inclusions Yes
Exit loop
filtered in Hierarchy I lie inside the new inclusion
No
Hierarchy III: check if any node on new inclusion lies inside Yes
Exit loop
of inclusions filtered in Hierarchy II
No
Save and go for next inclusion
End loop
Figure2:Flowcharttohierarchicalintersectionandoverlappingchecking.
(a) Ellipticalaggregatesandcir- (b) Ellipticalaggregatesandellip- (c) Ellipsoidalaggregatesandspheri- (d) Ellipsoidal aggregates and ellip-
cularvoids ticalvoids calvoids soidalvoids
Figure3:Numericallygenerated2Dand3Dmodelsforconcretewithgravelaggregates(𝑃 =30%,𝑃 =2%).
agg void
the shape of voids is assumed to be circle/ellipse in 2D or 2.3. Mesostructure Models. Elliptical aggregates and circu-
sphere/ellipsoidin3D. lar/ellipticalvoidsin2Dorellipsoidalaggregatesandspheri-
cal/ellipsoidalvoidsin3Dwereusedforconcretewithgravel
2.2. Specimen Generation. The basic idea is to create the aggregates.Aseriesofconcretesampleswithdimensionsof
aggregates and voids in the concrete in a repeated manner, 50mm × 50mm in 2D and 50mm × 50mm × 50mm in
until the target area/volume is filled with aggregates and 3Daregenerated.Figure3showsthenumericallygenerated
voids. The generation starts with the information input modelswithgravelaggregates(aggregatesingreenandvoids
process and is followed by taking and placing voids within in red). Aggregate gradation in Table1, aggregate volume
size range; the aggregates within the grading segments are fraction𝑃 = 30%,voidcontent𝑃 = 2%,aspectratios
agg void
produced last. The “input” process reads the controlling forelliptical/ellipsoidalaggregatesandvoids(theratioofthe
parameters for generating a heterogeneous concrete; the major axis to the minor axis) 𝑅1 = 𝑅2 = [1,2], minimum
“taking” process generates an individual void or aggregate
spacebetweentheedgeofanaggregateandtheboundaryof
inaccordancewiththerandomsizeandshapedescriptions. theconcretespecimen𝛾1 =0.2mm,andminimumgapwidth
The“placing”processpositionstheaggregatesandvoidsinto
thepredefinedarea/volumeinarandommanner,subjectto between any two aggregates 𝛾2 = 0.2mm are adopted here
[31].
theprescribedphysicalconstraints.Therearethreephysical
Polygonalaggregatesandcircular/ellipticalvoidsin2Dor
conditions to be satisfied simultaneously: (1) the whole
inclusion(voidoraggregate)mustbewithintheconcrete,(2) polyhedral aggregates and spherical/ellipsoidal voids in 3D
thereisnooverlappingofinclusions,and(3)thereisnointer- areconsideredforconcretewithcrushedaggregates.Aseries
sectionbetweenanytwoinclusions.Athree-levelhierarchical ofconcretesampleswiththesamedimensionsastheprevious
methodisproposedtoreducethecomputationalcostandis gravel aggregate samples are generated. Figure4 shows the
outlinedinFigure2.Theprocedurewasprogrammedusing numericallygeneratedmodelswithcrushedaggregates.The
MATLAB[32]inthisstudy,andthedetailedalgorithmsand minimumandmaximumedgesofthepolygons/polyhedrons
procedureingeneratingmesostructurescanbefoundinour aresetto𝑁 = 8,𝑁 = 16.Alltheotherparametersare
min max
previouswork[3,33]. thesameasthevaluesusedforgravelaggregates.
4 MathematicalProblemsinEngineering
(a) Polygonalaggregatesandcir- (b) Polygonal aggregates and (c) Polyhedralaggregatesandspheri- (d) Polyhedral aggregates and ellip-
cularvoids ellipticalvoids calvoids soidalvoids
Figure4:Numericallygenerated2Dand3Dmodelsforconcretewithcrushedaggregates(𝑃 =30%,𝑃 =2%).
agg void
3.CombinedNumerical-StatisticalMethod in Figure5. A bilinear cohesive crack model is used here
topredictdiscretecrackpropagationinconcretespecimens
3.1. Description of the Method. A combined numerical- undertensionloading.Alinearascendingbranchisaddedin
statistical method is proposed in this paper to study the thesofteningcurvetomodeltheinitiallyuncrackedmaterial.
material behaviour of concrete in a statistical sense. The The separation displacement is difficult to derive from
detailedprocedureisasfollows: experiments, so cohesive fracture energy and cohesive
strengthareusuallyused.Amongthem,onlytwoparameters
(1)Generate a model with prescribed variables, for
areindependent,andthefractureenergycanbeobtainedas
example,samplesize,aggregatevolumefraction,void
content,andaggregateshape. 𝛿
𝑓 1
𝐺=∫ 𝑡(𝛿)𝑑𝛿= 𝑡 𝛿 , (3)
(2)Performafiniteelementsimulationofthesamplefor 2 0 𝑓
0
givenboundaryconditions.
where 𝐺 is the cohesive fracture energy, 𝑡0 is the cohesive
(3)Compute the mean value, standard deviation, and strength,and𝛿𝑓istheseparationdisplacement.
coefficientofvariation(CoV)ofeffectivepropertyfor
Inthe3Dcohesivezonemodel,itisassumedthatthere
theconsideredsamplesize. existanormaltraction𝑡𝑛andtwotangentialtractions(shear
(4)Repeatsteps(1)to(3)forsufficientnumberofrandom cohesion)𝑡𝑠and𝑡𝑛acrossthecracksurfaces,throughmech-
samplestomeettherequiredprecision,andconduct anisms such as material bonding, aggregate interlocking,
statisticalanalysis. and surface friction in the fracture process zone. Figures
6(a) and 6(b) illustrate typical linear softening curves for
Thisprocedureisautomatedbyrunningabatchfileinthis 𝑡𝑛 − 𝛿𝑛 and 𝑡𝑠(𝑡𝑡) − 𝛿𝑠(𝛿𝑡), where 𝛿𝑛𝑓 and 𝛿𝑠𝑓(𝛿𝑡𝑓) are the
study.
critical relative displacements when the tractions diminish
Resultsfromallrealisationsareevaluatedstatistically.The
for normal traction and tangential traction components,
standarddeviation𝑠withinaseriesof𝑛samplesisgivenby
respectively.Theunloadingpathsarealsoindicated.
𝑛 It is worth noting that the initial tensile stiffness
𝑠2 = 𝑛−11∑𝑖=1(𝑥𝑖−𝑥)2, (1) 𝑘ac𝑛0rosshsothueldcboehehsiigvheeinntoeurgfahcteosebnefsourreetdhiesptelancseimleesntrtecnogntthin𝑡𝑛u0itiys
reached,butnottoohightocausenumericalinstabilitydue
where𝑥 = (1/𝑛)∑𝑛𝑖=1𝑥𝑖 istheseriesaverageresultand𝑥𝑖 is to ill-conditioned global stiffness matrix. Reasonable initial
theresultfromsample𝑖. shear stiffness, 𝑘𝑠0 and 𝑘𝑡0, is also needed before the shear
To compare results obtained with different sample sizes strength, 𝑡𝑠0 or 𝑡𝑡0, is reached. The effects of initial stiffness
quantitatively,weusethecoefficientofvariationgivenby on computational results have been investigated previously
[37, 38]. The following relation is suggested in [38] as a
𝑠
𝜀= . (2) guidelineforinitialstiffnessselection:
𝑥
𝑐(1−V)
Thisexpressesthefluctuationofmeasuredpropertyrelative 𝑘𝑛0 =𝑘𝑠0 =𝑘𝑡0 = 𝑏(1+V)(1−2V)𝐸, (4)
toitsaveragevalue.
where𝐸andVareYoung’smodulusandPoisson’sratio,𝑏is
3.2. Cohesive Zone Model in ABAQUS. The cohesive zone the characteristic size of elements, and 𝑐 is taken as 10∼100
model developed by Barenblatt [34, 35] and Dugdale [36] fromtheexperiencein[38].
enables the simulation of the energy dissipation process in The areas under the curves in Figures 6(a) and 6(b)
thefractureprocesszone(FPZ)duringfracture,asillustrated (calculated by (3)) represent the normal facture energy
MathematicalProblemsinEngineering 5
Aggregate
Concrete
Real crack Fracture process zone (FPZ)
Traction distribution along FPZ
t
0
Figure5:Characterizationofcohesivezonemodel(after[15]).
t
n
A
tn0 ts(tt)
A
t (t )
s0 t0
G (G )
k (k ) sf tf
s0 t0
k
n0
G k(k)
nf s t
−𝛿 (−𝛿 )
sf tf
𝛿 (𝛿 )
kn s0 t0 𝛿sf(𝛿tf)𝛿s(𝛿t)
𝛿 𝛿 𝛿
n0 nf n
t (t )
s0 t0
Loading Loading
Unloading Unloading
(a) 𝑡𝑛−𝛿𝑛curveinnormaldirection (b) 𝑡𝑠(𝑡𝑡)−𝛿𝑠(𝛿𝑡)curveintangentialdirection
Figure6:Bilinearsofteninglawsforcohesiveelements.
𝐺𝑛𝑓 and twice the tangential fracture energy, 𝐺𝑠𝑓 and 𝐺𝑡𝑓, displacements𝛿𝑚 combiningtheeffectsof𝛿𝑛,𝛿𝑠,and𝛿𝑡 can
respectively: beobtainedas
𝛿 =√⟨𝛿 ⟩2+𝛿2+𝛿2, (6)
𝑚 𝑛 𝑠 𝑡
𝛿
𝑛𝑓 1
𝐺 =∫ 𝑡 (𝛿 )𝑑𝛿 = 𝑡 𝛿 ,
𝑛𝑓 𝑛 𝑛 𝑛 2 𝑛0 𝑛𝑓 where⟨⟩istheMacaulaybracketand
0
𝐺𝑠𝑓 =∫0𝛿𝑠𝑓𝑡𝑠(𝛿𝑠)𝑑𝛿𝑠 = 21𝑡𝑠0𝛿𝑠𝑓, (5) ⟨𝛿𝑛⟩={{{0𝛿,𝑛, 𝛿𝛿𝑛𝑛 ≥<00 ((tceonmsiporne)ssion). (7)
𝛿
𝑡𝑓 1
𝐺 =∫ 𝑡 (𝛿)𝑑𝛿 = 𝑡 𝛿 .
𝑡𝑓 𝑠 𝑠 𝑠 2 𝑠0 𝑡𝑓 Thedamageevolutionlawisgivenby
0
𝛿 (𝛿 −𝛿 )
𝐷= 𝑚𝑓 𝑚,max 𝑚0 ,
(8)
Cohesive elements in ABAQUS based on the cohesive 𝛿 (𝛿 −𝛿 )
𝑚,max 𝑚𝑓 𝑚0
zonemodelareusedhere.Thedamageischaracterisedbya
scalarindex𝐷representingtheoveralldamageofthecrack where𝛿𝑚,maxisthemaximumeffectiverelativedisplacement
caused by all physical mechanisms. The effective relative attainedduringtheloadinghistory.𝛿𝑚0and𝛿𝑚𝑓areeffective
6 MathematicalProblemsinEngineering
relative displacements at damage initiation and complete mesh generation for a large number of samples required
failure,respectively.Itisobviousthatthedamagevariable𝐷 forstatisticalanalysis.Followingthemeshingprocedure,the
evolvesmonotonicallyfrom0to1uponfurtherloadingafter generatedmeshwillautomaticallyhavesharednodesatthe
theinitiationofdamage. interfacesbetweentwoelements.Iftheinterfacessurround-
The damage initiation and evolution will degrade the ing the elements are to be modelled explicitly as potential
unloadingandreloadingstiffnesscoefficients𝑘𝑛,𝑘𝑠,and𝑘𝑡, microcracksources,aduplicatesetofnodeswillberequired
whichcanbecalculatedas at the interface locations. Here 4-node or 6-node zero in-
planethicknessCIEsarepreinsertedintotheexistingelement
𝑘𝑛 =(1−𝐷)𝑘𝑛0, interfaces in 2D or 3D. This is conducted by a purpose
𝑘𝑠 =(1−𝐷)𝑘𝑠0, (9) winrtietrtfeanceaenldemfleenxtisbilnesienr-tihoonu(sCeIcEoImN)p.uThterrepersoegtsraomfCcIoEhsewsiivthe
𝑘𝑡 =(1−𝐷)𝑘𝑡0. differenttraction-separationsofteninglawsareinserted(see
Figure7), namely, CIE-AGG inside the aggregates (grey in
Thetractionsarealsoaffectedbythedamageaccordingto Figures7(b)and7(d)),CIE-MORinsidethemortar(greenin
Figures7(b)and7(d)),andCIE-INTontheaggregate-mortar
{(1−𝐷)𝑡𝑛, 𝑡𝑛 ≥0 interfaces (yellow in Figures 7(b) and 7(d)). The element
𝑡𝑛 ={ and node numbers are denoted by 𝐸 and 𝑁, respectively.
{𝑡𝑛, 𝑡𝑛 <0, Thedetailednumberingofelementsandnodesintheinitial
(10)
mesh and the mesh with inserted cohesive elements shows
𝑡𝑠 =(1−𝐷)𝑡𝑠, theinsertionprocedurewiththenewnodesgeneratedatthe
samepositionsandCIEsbetweenthecontinuumelements.
𝑡𝑠 =(1−𝐷)𝑡𝑡,
where 𝑡𝑛 and 𝑡𝑠 are the traction components predicted by 4.NumericalSimulationandResults
the elastic traction-displacement behaviour for the current
separationwithoutdamage. 4.1. Description of the Numerical Model. 2D concrete spec-
In this study, damage in the cohesive zone model is imens with elliptical aggregates and circular voids and 3D
assumed to initiate when a quadratic interaction function specimens with ellipsoidal aggregates and spherical voids
involvingthenominalstressratiosreachesavalueofone (𝑃agg =30%,𝑃void =2%,𝑅1 =𝑅2 =[1,2],𝛾1 =𝛾2 =0.2mm)
were modelled, loaded under uniaxial tension. 25mm ×
⟨𝑡 ⟩ 2 𝑡 2 𝑡 2 25mmconcretesquarein2Dand25mm × 25mm × 25mm
{ 𝑛 } +{ 𝑠 } +{ 𝑡 } =1. (11) concretecubicin3Dareusedhere.Alltheothermesoscale
𝑡 𝑡 𝑡
𝑛0 𝑠0 𝑡0 parameters were fixed to be the same for 2D and 3D
models.Figure8showsthegeometry,boundary,andloading
Thematerialproperties,suchasdensity,Young’smodu-
conditionsofthenumericalmodel.Horizontaldisplacements
lus,Poisson’sratio,tensilestrength,andfractureenergy,are
were prescribed to all nodes on the left and right surfaces
set forcontinuumelementsofaggregates andmortar,three
of the specimen, with a value equal to zero for the left
different interface cohesive elements. The material hetero-
surface, and a uniformly distributed displacement on the
geneityisinvestigatedbyconsideringdifferentphasesinthe
rightsurface.Verticaldisplacementsforthesamenodesare
concretespecimenwithcorrespondingmaterialproperties.
leftfree,exceptforthenodesattheleftlowercornerofthe
Inthefractureprocesszonefora2Dcase,tractionsexist
in the normal direction 𝑡𝑛 and shear direction 𝑡𝑠 across the specimen,whicharefixedtopreventrigidbodytranslation.
Adisplacementcontrolledloadingschemewasusedandall
crackinterface,andthecorrespondingrelativedisplacements
are 𝛿𝑛 and 𝛿𝑠. So the 2D constitutive law could be easily tlohaedainnagltyismesew𝑡=ere0.e0n0d5esd[3a]t.Caodnisspidlaecreinmgetnhte𝑑ba=lan0c.1embemtwaenend
deducedfrom3Dcharacterizationdescribedabovebytaking
accuracyandefficiency,theelementlength1mmwasusedfor
oneofthesheartractionsanddisplacementsout.
allthemeshesinthisstudy[3].
Generally, aggregates have much higher strength than
3.3. Cohesive Interface Elements Insertion. In this study, all
mortar and interfaces in normal concrete. In this section,
finiteelementmeshingisperformedwiththepreprocessing
no cracks are allowed to initiate inside the aggregates by
functionality in commercial FE package ANSYS [29]. To
assuming elastic behaviour without damage in CIE-AGG.
meshthemesoscopicstructureofconcrete,differentmaterial
The linear tension/shear softening laws described above
parts(mortarandaggregates)shouldmaintaincontinuityon
were used to model CIEs with quadratic nominal stress
theirsurfaces.Hence,thefiniteelementboundariesarecoin-
cidentwithmaterialsurfacesandtherearenomaterialdis- initiation criterion and linear damage evolution criterion.
continuitieswithintheelements.Themortarandaggregates Similarmaterialpropertiesextractedfrom[39]wereusedin
aremeshedwithtriangularelements(plainstress)in2Dand thisstudy,whicharelistedinTable2.Thefractureproperties
hexahedralelementsin3Dsothatmorerealisticcrackpaths relatedtotheshearbehaviourwereassumedtobethesame
canbeobtained.ThespeciallydevelopedANSYSparametric asthenormalonesduetothelackofexperimentaldata.Both
designlanguage(APDL)programsincombinationwiththe initialtensileandshearstiffnessforallthecohesiveelements
6
ANSYSbatchprocessingprovideapowerfultoolofautomatic weresetto10 MPa/mmaftertrialanderror.
MathematicalProblemsinEngineering 7
N12 N13 N11 N44 N45 N46N47N48 N43N42
E14 E16 E14 E16
E15 E34 E35 E15 E36 E33
E13 E12 E13 N29N30N31 N36N37N35E12
N8 E8 N9 N26N28N27 E24 N33N34N32
E31 E8 E32
N7 E6 E9 N10 E11 E7 N6 N25 E18 E17 N22
N24 E6 E22 E9 N39N38N41 E11 E23 E7 N21
N40
E10
N23 N20
N2 N4
E5 E4 E30 E19 E10 E20 E29
E2 N6N7 N8 N16N17N15
E1 E3 E5 N5N3N4 E21 N12N13N14E4
E2
N1 N3 N5 E26 E1 E25 E27 E3 E28
Mortar
N2 N1 N9N10N11 N18N19
Aggregate
CIE-MOR Mortar
CIE-AGG Aggregate
CIE-INT
(a) 2Dinitialmesh (b) 2Dmeshwithzero-thicknessCIEs
N1 N2N3
N1 N4
E5 E7
E6
N2 E1 N4 E4 N7 N5 E1 N8 N11 E4 N16
N9N10
E2 E3 E2
E3
N3 N5 N6 N6 N15
N7 N12N13 N14
Mortar
CIE-MOR Mortar
Aggregate
CIE-AGG Aggregate
CIE-INT
(c) 3Dinitialmesh (d) 3Dmeshwithzero-thicknessCIEs
Figure7:Insertingdifferentcohesiveelementstotheinitialmesh.
Table2:Materialproperties.
mYo(oMduuPnluag)s’s𝐸 Proaitsisoon]’s (10−9Dteonnsnitey/m𝜌m3) E𝑘l𝑛as(tMicPsati/ffmnmes)s sCtr(oeMnhgePstahi)ve𝑡𝑛 e(FnNrear/cgmtyum𝐺re)𝐹
Aggregate 70000 0.2 2.8 — — —
Mortar 25000 0.2 2.0 — — —
CIE-AGG — — 2.8 106 — —
CIE-MOR — — 2.5 106 4 0.06
CIE-INT — — 2.0 106 2 0.03
8 MathematicalProblemsinEngineering
25
mm
25mm
m
m
5
2
m
m
5
2
m
m
25
(a) 2D (b) 3D
Figure8:Specimendimensions,loading,andboundaryconditions.
4 4
A-3D
3 3
A-2D
Pa) - 3m)
Mean stress (M 2 Toughness (kJ/ 2
B-2D/3D
1 1
0 0
0.00 0.02 0.04 0.06 0.08 0.00 0.02 0.04 0.06 0.08
Displacement (mm) Displacement (mm)
2D 2D
3D 3D
(a) Stress-displacementcurves (b) Toughness-displacementcurves
Figure9:Comparisonofstress-displacementandtoughness-displacementcurves.
ABAQUS/explicit with small time increments (typically Fiftyrandomsamplesweremodelledinthisstudytoensure
Δ𝑡 = 3 × 10−8s) was utilised to solve the highly nonlinear thattheresultsarestatisticallyconverged.
equationsystemsbythehighperformancecomputing(HPC) In Figure9, resulting mean stress-displacement and
clusterattheUniversityofManchestercomputationalshared toughness-displacement curves, together with the variation
facility(CSF).Atypicalrunofsimulationswith50samplesin range for both 2D and 3D, are depicted. It can be seen
thisstudytakesabout6hoursbyparallelcomputationwith thatthepredictedstress-displacementcurvesarequalitatively
32CPUs. similar with a clear peak and sharp initial postpeak drop
followed by a long shallow tail. Numerical mean stress for
4.2. Tensile Behaviour. Due to the statistical nature of every curve is obtained by summing the nodal reaction
mesostructuremodels,anextensiveseriesofnumericalsimu- forces in the constraints and dividing by the specimen
lationswouldbenecessarytocapturetherangeofbehaviours. cross section area. The mean peak stress predicted by 2D
MathematicalProblemsinEngineering 9
0.5 6 80
B-3D
5
Pa) 0.4 60
M
on ( 0.3 B-2D %) 4 %)
eviati V-A ( 3 40 V-B (
d d 0.2 Co Co
ar 2
and A-2D 20
St 0.1 1
B-3D
0.0 0 0
0.00 0.02 0.04 0.06 0.08 0 10 20 30 40 50
Displacement (mm) Sample number
2D 2D-A 2D-B
3D 3D-A 3D-B
Figure10:Comparisonofstandarddeviation-displacementcurves. Figure11:InfluenceofsamplenumberontheCoVofstress.
and 3D modelling is 2.65MPa and 3.49MPa, respectively, testingofdifferentsamples.ItcanbeobservedthatCOVof
representing an increase of 24.1%. The mean toughness peakstressfordifferentvoidandaggregatedistributionislow
(energy absorbed per unit volume up to the break point, (4.2%for2Dmodellingand1.5%for3Dmodelling),which
i.e., the area under the stress-strain curve up to the break demonstrates that the peak stress is relatively insensitive to
3
point,measuredinJ/m )predictedin2Dand3Dmodelling thevoidandaggregatedistribution.ThehighCoVofstress
3 3
is 2.47KJ/m and 3.15KJ/m , respectively, representing an at point B further demonstrates the necessity to conduct
increase of 21.6%; Figure10 shows the standard deviation- a statistical analysis as the results vary greatly for different
displacement curves for both 2D and 3D modelling. The aggregateandvoiddistributioninsofteningrange.
standard deviation of peak stress is 0.11MPa and 0.05MPa, Figure12showsthestatisticaldistributionofpeakstress,
respectively,representingadecreaseof54.5%.Theseincreases
togetherwiththeprobabilitydensityandthebestfitGaussian
in peakstress ordecreases in standarddeviationaredue to
ProbabilityDensityFunction(PDF).Theprobabilitydensity
theconstrainteffectsinthethicknessdirectionin3Dandthe
curvecanbeusedtocalculatestructuralreliabilityorfailure
lesssmoothcrackingsurfacesfrom3Dmodelling.However,
reliabilityagainstgivenexternalloadingsandmaterialprop-
the standard deviation of 3D modelling at softening range
ertiesforstructuraldesign.
islargerthanthatin2Dmodellingwhichindicatesthat3D
modelling is more complicated and dependent on different
samples. This is due to the 3D heterogonous distribution 4.3.CrackPatterns. Thecomplexmesoscalecrackpropaga-
of aggregates and voids. It is clear that the third dimension tion is realistically simulated using the proposed method,
cannotbeneglectedduetothenatureoffractureprocess.The andtypicalcrackpatternsforboth2Dand3Datfailureare
importanceofthethicknesseffecthasalsobeenpointedout
showninFigure13.Toclearlyvisualisethefracturesurfaces,
experimentallybyVanMierandSchlangen[40].
models with damaged cohesive elements, models without
The mean curve, mean value, and standard deviation
cohesiveelements,andfailedinterfacesonlyareusedfor2D
of the stress and toughness shown in Figures 9 and 10 are
models(seeFigures13(a)–13(c)),andmodelswithdamaged
extractedfromthestress-displacementcurvesfor50samples
cohesive elements, models without cohesive elements, and
with different aggregate and void distribution. In order to
morphology of failed surface are used for 3D models (see
evaluate the effect of sample number on simulation results
Figures13(d)–13(f)).ThecracksshowninFigures13(a),13(c),
asrequiredfortheproposedcombinednumerical-statistical
and 13(d) are represented by red CIEs with damage index
method,twospecialloadingpoints,namely,pointAatpeak
𝐷 ≥ 0.95(𝐷 = 1meanscompletefailure).Amagnification
stress and point B at displacement 0.03mm where largest
factorof20and200wasusedforthedeformedmodelswith
standarddeviationisobserved(denotedinFigures9(a)and
10), are investigated. The influence of sample number on damagedCIEsandthosewithoutdamagedCIEs,respectively.
the coefficient of variation (CoV) of stress at points A and Thenoticeabledifferenceof2Dand3Dmodellingdiscussed
B is illustrated in Figure11. It can be seen that 50 samples above could also be summarized in that the 3D modeling
are enough to get convergent values of CoV of stress with predicts more realistic fracture surfaces in the thickness
astablefluctuation.Thisisanindicationthatthecombined directionwhereas2Dmodelingonlyassumesplanefracture
numerical-statistical method with 50 random realisations problems. This result confirms the importance of including
offers a good balance between general applicability and thethirddimensionintotheanalysis.
10 MathematicalProblemsinEngineering
0.25 0.30
Mean:2.65MPa, standard deviation:0.11MPa Mean:3.49MPa, standard deviation:0.05MPa
0.25
0.20
y y 0.20
nsit 0.15 nsit
e e
d d
y y 0.15
abilit 0.10 abilit
b b
o o 0.10
Pr Pr
0.05
0.05
0.00 0.00
2.3 2.4 2.5 2.6 2.7 2.8 2.9 3.35 3.40 3.45 3.50 3.55 3.60
Peak stress (MPa) Peak stress (MPa)
(a) 2D (b) 3D
Figure12:StatisticaldistributionofpeakstressfromMonteCarlosimulation.
(a) 2DmodelwithdamagedCIEs (b) 2DmodelswithoutCIEs (c) 2Dmodelswithfailedinterfaces
(d) 3DmodelwithdamagedCIEs (e) 3DmodelwithoutCIEs (f) 3Dmorphologyoffailedsurface
Figure13:Failureofconcreteundertension.