Numerical Simulations of Highly Porous Dust Aggregates in the Low-Velocity Collision Regime

5.1 MB·English

by  Ralf J. Geretshauser

Astronomy&Astrophysicsmanuscriptno.technical˙paper˙v10 (cid:13)c ESO2010
January18,2010
Numerical Simulations of Highly Porous Dust Aggregates in the
Low-Velocity Collision Regime
Implementation and Calibration of a Smooth Particle Hydrodynamics Code
R.J.Geretshauser1,R.Speith1,C.Gu¨ttler2,M.Krause2,andJ.Blum2
1 Institut fu¨r Astronomie und Astrophysik, Abteilung Computational Physics, Eberhard Karls Universita¨t Tu¨bingen, Auf der
Morgenstelle10,D-72076Tu¨bingen,Germany
0 e-mail:[email protected]
1 2 Institut fu¨r Geophysik und extraterrestrische Physik, Technische Universita¨t zu Braunschweig, Mendelssohnstr. 3, D-38106
0 Braunschweig,Germany
2

n Preprintonlineversion:January11,2010;acceptedbyAstronomy&AstrophysicsDecember22,2009
a
J ABSTRACT

1
1 Context.A highly favoured mechanism of planetesimal formation is collisional growth. Single dust grains hit each other due to
relativevelocitiescausedbygasflowsintheprotoplanetarydisc,whichtheyfollow.TheystickduetovanderWaalsforcesandform
] fluffyaggregatesuptocentimetresize.Themechanismoffurthergrowthisunclearsincetheoutcomeofaggregatecollisionsinthe
P relevantvelocityandsizeregimecannotbeinvestigatedinthelaboratoryunderprotoplanetarydiscconditions.Realisticstatisticsof
E theresultofdustaggregatecollisionsbeyonddecimetresizeismissingforadeeperunderstandingofplanetarygrowth.
. Aims.Joiningexperimentalandnumericaleffortswewanttocalibrateandvalidateacomputerprogramthatiscapableofacorrect
h
simulation of the macroscopic behaviour of highly porous dust aggregates. After testing its numerical limitations thoroughly we
p
willcheck theprogramespeciallyforarealisticreproduction ofthecompaction, bouncing andfragmentationbehaviour. Thiswill
-
demonstrate the validity of our code, which will finally be utilised to simulate dust aggregate collisions and to close the gap of
o
fragmentationstatisticsinfuturework.
r
t Methods.Weadoptthesmoothparticlehydrodynamics(SPH)numericalschemewithextensionsforthesimulationofsolidbodies
s
and a modified version of the Sirono porosity model. Experimentally measured macroscopic material properties of SiO dust are
a 2
[ implemented.Bysimulatingthreedifferentsetupswecalibrateandtestforthecompressivestrengthrelation(compactionexperiment)
andthebulkmodulus(bouncingandfragmentationexperiments).Datafromexperimentsandsimulationswillbecompareddirectly.

1 Results.SPHhasalreadyproventobeasuitabletooltosimulatecollisionsatrather highvelocities.Inthisworkwedemonstrate
v thatitsareaofapplicationcannotonlybeextendedtolow-velocityexperimentsandcollisions.Itcanalsobeusedtosimulatethe
7 behaviour ofhighlyporous objectsinthisvelocityregimetoaveryhighaccuracy. A correctreproduction ofdensitystructuresin
1 thecompactionexperiment,ofthecoefficientofrestitutioninthebouncingexperimentandofthefragmentmassdistributioninthe
6 fragmentation experiment show thevalidity and consistency of our code for the simulationof theelastic and plasticproperties of
1 thesimulateddustaggregates.TheresultofthiscalibrationprocessisanSPHcodethatcanbeutilisedtoinvestigatethecollisional
. outcomeofporousdustinthelow-velocityregime.
1
0 Key words. hydrodynamics - methods: laboratory - methods: numerical - planets and satellites: formation - planetary systems:
0 formation-planetarysystems:protoplanetarydisks
1
:
v
i 1. Introduction (Heimetal. 1999). In this process, which has been investi-
X
gated experimentally (Blumetal. 2000; Blum&Wurm 2000;
r In gaseouscircumstellar discs, the potentialbirthplaceof plan- Krause&Blum 2004) and numerically (Dominik&Tielens
a
etary systems, dust grains smaller than a micrometre grow to
1997; Paszun&Dominik 2006, 2008, 2009; Wadaetal. 2007,
kilometre-sizedplanetesimals,whichthemselvesproceedtoter-
2008, 2009), they form fluffy aggregateswith a high degreeof
restrialplanetsandcoresofgiantplanetsbygravity-drivenrun-
porosity.
away accretion. Depending on their size the dust grains and
aggregates perform motions in the disc and relative to each Due to restructuring the aggregates gain a higher mass to
other. Brownian motion, radial drift, vertical settling and tur- surfaceratioandreachhighervelocities.Blum&Wurm(2000,
bulent mixing cause mutual collisions (Weidenschilling 1977; 2008)andWadaetal.(2008)showedthatcollisionsamongthem
Weidenschilling&Cuzzi1993). lead to fragmentation and mass loss. Depending on the model
Since real protoplanetary dust particles are not available of the protoplanetary disc, which provides the kinetic colli-
for experiments in the laboratory much of the following sion parameters, this means that direct growth ends at aggre-
work has been carried out with dust analogues such as SiO gate sizes of a few centimetres. However, Wurmetal. (2005)
2
(Blum&Wurm 2008). Theoreticalmodels also refer to micro- andTeiser&Wurm(2009)havedemonstratedinlaboratoryex-
scopic and macroscopic properties of these materials. Initially periments in the centimetre regime and with low-porosity dust
thedustgrainshitandstickoncontactbyvanderWaalsforces thattheprojectilecanstickpartiallytothetargetatvelocitiesof
2 R.J.Geretshauseretal.:NumericalSimulationsofHighlyPorousDustAggregates
morethan20ms−1.Thus,collisionalgrowthbeyondcentimetre 2009b) and utilised to understand the formation of an aster-
sizeseemstobepossibleandtheexactoutcomeofthefragment oid family (Jutzietal. 2009a). However, pumice is a mate-
distributioniscrucialfortheunderstandingofthegrowthmech- rial whose strength parameters decrease when it is compacted
anism. (crushed). Additionally, the underlyingthermodynamicallyen-
Numerical models that try to combine elaborate pro- hancedporositymodelisdesignedtodescribeimpactsofsome
toplanetary disc physics with the dust coagulation prob- kms−1. Thus, it is perfectly suitable to simulate high-velocity
lem(Weidenschillingetal.1997;Dullemond&Dominik2005; collisionsofporousrock-likematerial.
Braueretal.2008;Zsom&Dullemond2008;Ormeletal.2007,
2009)havetomakeassumptionsabouttheoutcomeofcollisions
In contrast, collisions between pre-planetesimals occur at
betweendustyobjectsforallsizesandrelativevelocities.Since
relative velocities of some tens of ms−1 and compressive,
data for these are hardly available, in the most basic versions
shear and tensile strengths are increasing with increasing den-
ofthesemodelsperfectsticking,inmoreelaborateonespower-
sity (Blum&Schra¨pler 2004; Blumetal. 2006; Gu¨ttleretal.
law fragment distributions from experiments and observations
2009b).
(Mathisetal.1977;Davis&Ryan1990;Blum&Mu¨nch1993)
areassumed.Theresultsofalikesimulationshighlydependson
the assumed fragmentation kernel. However, the given experi- Scha¨feretal. (2007) have used an SPH code based on the
mentalreferenceshavebeenmeasuredonlyforsmallaggregate porosity model by Sirono (2004) to simulate collisions be-
sizes.Theinfluenceofinitialparameterssuchastherotationof tween porous ice in the ms−1 regime. They found that a suit-
the objects or their porosity have not been taken into account able choice of relations for the material parameters can pro-
forthesizeregimesbeyondcentimetresize,althoughtheymight duce sticking, bouncing or fragmentation of the colliding ob-
playanimportantrole(Sirono2004;Ormeletal.2007). jects. Therefore,theystressedtheimportanceofcalibratingthe
A new approach to model the growth of protoplanetary material parametersof porousmatter with laboratorymeasure-
dustaggregateswasrecentlydevelopedbyGu¨ttleretal.(2009a) ments.Numericalmolecular-dynamicssimulationsareaboutto
andZsometal.(2009)whodirectlyimplementedtheresultsof use asufficientnumberofmonomerstobeclose enoughto the
dustaggregationexperimentsintoaMonte-Carlogrowthmodel. continuumlimitandtoprovidetherequiredmaterialparameters
They found that bouncing of protoplanetary dust aggregates (e.g.Paszun&Dominik2008,reproducingexperimentalresults
playsamajorrolefortheirevolutionasitisabletoinhibitfur- of Blum&Schra¨pler 2004). They represent an important sup-
thergrowthandchangestheiraerodynamicproperties.Although portto thedifficultexperimentaldeterminationofthesequanti-
their model relies on the most comprehensivedatabase of dust ties.InGu¨ttleretal.(2009b)wehavemeasuredthecompressive
aggregate properties and their collisional behaviour, they were strengthrelationforsphericalSiO dustaggregatesforthestatic
2
unable to make directpredictionsfor any arbitraryset of colli- caseinthelaboratoryandhavegivenaprescriptionhowtoapply
sionparametersandwerethusobligedtoperformextrapolations thistothedynamiccaseusingacompactioncalibrationexperi-
overordersofmagnitude.Someoftheseextrapolationsarebased mentand2D simulations.We havepointedoutrelevantbench-
onphysicalmodelswhichneedtobesupportedbyfurtherexper- markfeatures.Itwasshownthatthecodeisinprincipalcapable
imentsandwherethisisnotpossiblebysophisticatednumerical of simulating not only fragmentationbut also bouncing,which
modelssuchastheonepresentedhere. cannotbeseeninmolecular-dynamicscodessofar.
Sticking,bouncing,andfragmentationaretheimportantcol-
lision outcomes which need to be implemented into a coagu-
InthisworkwepresentourSPHcodewithitstechnicalde-
lation code and affect the results of these models. Thus, it is
tails (Sect. 2), experimental reference (Sect. 3), and numerical
not only importantto correctly implementthe exact thresholds
properties(Sect. 4).On the basisof the compactioncalibration
betweenthese regimesbutalso detailsconcerningtheoutcome
simulationwedemonstratethattheresultsconvergeforincreas-
suchasthefragmentsizedistribution,thecompactioninbounc-
ing spatial resolutionand choose a sufficient numericalresolu-
ingcollisions,andmaybeeventheshapeoftheaggregatesafter
tion(Sect.4.2). We investigatethedifferencesbetween2D and
they merged in a sticking collision. Due to the lack of impor-
3Dnumericalsetup(Sect.4.3)andtherebyimprovesomedraw-
tantinputinformationthenecessityforasystematicstudyofall
backs in Gu¨ttleretal. (2009b). Artificial viscosity will be pre-
relevantcollisionparametersarisesandwillbeaddressedinthis
sentedasstabilisingtoolforvariousproblemsinthesimulations
work.
anditsinfluenceonthephysicalresultswillbepointedout(Sect.
Becauseofrestrictionsinsizeandrealisticenvironmentpa-
4.4).
rameters this task cannot be achieved in the laboratory alone.
Therehasbeenalotofworklatelyonmodellingthebehaviour
of dust aggregates on the basis of molecular interactions be- Mostprominentlywecontinuethecalibrationprocessstarted
tween the monomers (Paszun&Dominik 2006, 2008, 2009; inGu¨ttleretal.(2009b)utilisingtwofurthercalibrationexperi-
Wadaetal.2007,2008,2009).However,simulatingdustaggre- mentsforbouncing(Sect.5.2)andfragmentation(Sect.5.3).In
gates with a model based on macroscopic material properties theendwepossessacollectionofmaterialparameters,whichis
such as density, porosity, bulk and shear moduli and compres- consistentforallbenchmarkexperiments.Finally,theSPHcode
sive,tensileandshearstrengthsremainsanopenfieldsincethese has gained enough reliability to be used to enhance our infor-
quantities are rarely available. The advantage of this approach mationabouttheunderlyingphysicsofdustaggregatecollisions
overthemoleculardynamicsmethod,whichiscomputationally beyondthecentimetreregime.Infutureworkitwillbeapplied
limited to a few ten thousandmonomers,is the accessibility of to generatea catalogue of pre-planetesimalcollisionsand their
aggregatesizesbeyondthecentimetreregime. outcomeregardingallrelevantparametersforplanetformation.
Recently, Jutzietal. (2008) have implemented a porosity Jointlywithexperimentsandcoagulationmodelsthiscanbeim-
model into the smooth particle hydrodynamics (SPH) code plementedintoprotoplanetarygrowthsimulations(Gu¨ttleretal.
by Benz&Asphaug (1994). It was calibrated for pumice 2009a;Zsometal.2009)toenhancetheirreliabilityandpredic-
material using high-velocity impact experiments (Jutzietal. tivepower.
R.J.Geretshauseretal.:NumericalSimulationsofHighlyPorousDustAggregates 3
2. PhysicalModelandNumericalMethod The conservation of momentum is ensured by the second
equation
2.1.SmoothParticleHydrodynamics
dvα 1∂σαβ
The numericalLagrangianparticle methodsmooth particle hy- = . (3)
drodynamics (SPH) was originally introduced by Lucy (1977) dt ρ ∂xβ
andGingold&Monaghan(1977)tomodelcompressiblehydro-
InSPHformulation,themomentumequationreads
dynamic flows in astrophysicalapplications. Later, the method
has been extended (Libersky&Petschek 1990), and improved dvα σαβ σαβ ∂Wij(h)
extensively (e.g. Liberskyetal. 1993; Benz&Asphaug 1994; i = m i + j . (4)
Relaanstdiclesan&dLpilbaesrtisckyde1fo9r9m6a;tiLonibseorsfksyoelitdalm. a1t9er9i7a)ls.toA sciommuplarete- dt Xj j ρ2i ρ2j  ∂xβi
Due to the symmetry in the interaction terms, conservation of
hensive description of SPH and its extensions can be found in
momentum is ensured by construction. Additionally we ap-
Monaghan(2005).
plythestandardSPH artificialviscosity(Monaghan&Gingold
IntheSPHscheme,continuoussolidobjectsarediscretized
1983). This is essential in particular for stability at interfaces
intointeractingmasspackages,socalled“particles”.Thesepar-
with highly varying densities. The influence of artificial vis-
ticlesformanaturalframeofreferenceforanydeformationand
cosity on our simulation results is investigated thoroughly in
fragmentationthatthesolidbodymayundergo.Allspatialfield
Sec.4.4.
quantitiesoftheobjectareapproximatedontotheparticleposi-
The third equation, the energy equation, is not used in our
tions x by a discretizedconvolutionwith a kernelfunctionW.
i
model. Hence, we assume that kinetic energy is mainly con-
ThekernelWdependsonparticledistance|x −x|andhascom-
j i
vertedintodeformationenergyandenergydissipatedbyviscous
pact support, determined by the smoothing length h . We are
effectsisconvertedintoheatandradiatedaway.Thestresstensor
using the standard cubic spline kernel (Monaghan&Lattanzio
σcanbesplitintoapartrepresentingthepurehydrostaticpres-
1985)butnormalisedsuchthatitsmaximumextensionisequal
sure p and a traceless part for the shear stresses, the so called
toonesmoothinglengthh.
deviatoricstresstensorSαβ.Hence,
We apply a constant smoothing length. This also allows
to model fragmentation of solid objects in a simple way. σαβ =−pδαβ+Sαβ. (5)
FragmentationoccurswhensomeSPHparticleswithinthebody
lose contact with their adjacent particles. Two fragments are Anydeformationof a solid bodyleadsto a developmentof in-
completelyseparatedassoonastheirrespectivesubsetsofparti- ternalstressesinaspecificallymaterialdependentmanner.The
cleshavereachedadistanceofmorethan2hsothattheirkernels relation between deformationand stresses is not taken into ac-
donotoverlapanymore. countwithintheregularequationsoffluiddynamics.Therefore,
Time evolutionof theSPH particlesiscomputedaccording theyareinsufficienttodescribeaperfectlyelasticbodyandhave
totheLagrangianformoftheequationsofcontinuummechan- tobeextended.Themissingrelationsaretheconstitutiveequa-
ics, while transferring spatial derivatives by partial integration tions,whichdependonthestraintensor
ontotheanalyticallygivenkernelW.
1 ∂x′α ∂x′β
ǫαβ = + . (6)
2 ∂xβ ∂xα!
2.2.Continuummechanics
Itrepresentsthelocaldeformationofthebody.Theprimedco-
Asystemofthreepartialdifferentialequationsformstheframe-
ordinatesdenotethepositionsofthedeformedbody.
workofcontinuummechanics.Ascommonlyknowntheyfollow
FollowingHooke’slaw a proportionalrelation betweende-
fromtheconstraintsofconservationofmass,momentumanden-
formation is assumed involving the material dependent shear
ergy.Thefirst,accountingfortheconservationofmass,iscalled
modulusµ=µ(ρ),whichdependsitselfonthedensity:
continuityequation
1
dρ ∂vα Sαβ ∝2µ ǫαβ− δαβǫγγ . (7)
+ρ =0. (1) d !
dt ∂xα
However,thisis onlythe constitutiveequationforthe traceless
Followingtheusualnotationρandvα denotedensityandveloc-
shearpart.Forthehydrostaticpartofthestresstensorweadopt
ity, respectively. Greek indices run from 1 to d, the dimension
amodificationoftheMurnaghanequationofstatewhichispart
oftheproblem.Einsteinsummingnotationholdsthroughoutthe
oftheSirono(2004)porositymodel:
entirepaper.IncontrasttotheusualSPHscheme,wheretheSPH
densityρiiscalculateddirectlyfromtheparticledistribution,we p(ρ)= K(ρ′)(ρ/ρ′ −1), (8)
0 0
solvethecontinuityequationaccordingto
where ρ′ is the so called reference density, the density of the
0
dρi =−ρ mj(vα−vα)∂Wij(h) (2) materialatzeroexternalstress,andK(ρ)isthebulkmodulus.
dt i ρ j i ∂xα The density dependence of the bulk K(ρ) and shear µ(ρ)
Xj j i moduliismodelledbyapowerlaw
(e.g.Randles&Libersky1996). Here the sum runsoverallin- K(ρ)=2µ(ρ)= K (ρ/ρ)γ. (9)
teractionpartners j of particlei, m is the particlemass ofpar- 0 i
j
ticle j,andWij denotesthekernelfortheparticularinteraction. AlthoughaccordingtoSirono(2004) ρ isthe initialdensityof
i
Althoughthisapproachismoreexpensive,asitrequirestosolve thematerialatthe beginningofthesimulation,we,in contrast,
anadditionalordinarydifferentialequationforeachparticle,itis want to ensure that our dust material possesses the same bulk
morestableforhighdensitycontrastsanditavoidsartifactsdue modulus K(ρ) even for simulations with different initial densi-
tosmoothingatboundariesandinterfaces. ties. In this work the dust material has two different densities
4 R.J.Geretshauseretal.:NumericalSimulationsofHighlyPorousDustAggregates
at the beginningof the bouncing(Sect. 5.2) and fragmentation 2.3.Porosityandplasticity
(Sect.5.3)calibrationsetup.AccordingtoSirono(2004)thema-
FollowingtheSirono(2004)model,theporosityΦismodelled
terials should feature two differentρ. As a consequence K(ρ),
i bythedensityoftheporousmaterialρandtheconstantmatrix
andinparticularK ,dependsontheinitialsetup.Sincewewant
0 densityρ
to compare and validate K by using two different setups we s
0
have to fix a unique ρi for all simulations. We choose ρi such ρ
thatK(ρ)= K isthebulkmodulusofthegenericuncompressed Φ=1− =1−φ, (17)
i 0 ρ
s
dustmaterialthatisproducedbytherandomballisticdeposition
(RBD)method(Blum&Schra¨pler2004). Inthefollowingwewillusethefillingfactorφ=ρ/ρ .
s
The time evolution of the pressure is directly given by the Wemodelplasticitybyreducinginnerstressesgivenbyσαβ.
time evolution of the density (Eq. 1) and since the pressure Forthisreasonweneedconstitutiverelationsdescribingthebe-
is a scalar quantity it is intrinsically invariant under rotation. haviour of the material during plastic deformation.These rela-
However, in order to gain a frame invariantformulation of the tions are specific for each material and have to be determined
timeevolutionofthedeviatoricstresstensor,i.e.thestressrate, empirically. Particularly for highly porous materials it is ex-
correctiontermshavetobeadded.Averycommonformulation tremely difficult to acquire them. Therefore, it is an advanta-
forSPH(seee.g.Benz&Asphaug1994;Scha¨feretal.2007)is geousfeatureofourmodelthatitisbasedonmeasurementsby
theJaumannrateform Blum&Schra¨pler (2004), Blumetal. (2006) and Gu¨ttleretal.
dSαβ 1 (2009b).
=2µ ǫ˙αβ− δαβǫ˙γγ +SαγRγβ+SβγRγα, (10)
The main idea of the adopted plasticity model is to reduce
dt d !
inner stress once the material exceeds a certain plasticity cri-
whereRαβistherotationratetensor,definedby terion. In the elastic case, described by Eqns. 7 and 8, inner
1 ∂vα ∂vβ stressesgrowlinearlywithdeformation.Hence,thematerialre-
Rαβ = − (11) turnstoitsoriginalshapeatvanishingexternalforces.Reducing
2 ∂xβ ∂xα!
inner stresses, i.e. deviating from the elastic deformation path,
andǫ˙αβ denotesthestrainratetensor reducestheinternalabilityofthematerialtorestoreitsoriginal
1 ∂vα ∂vβ shape.Therefore,bystressreductiondeformationbecomesper-
ǫ˙αβ = + . (12)
2 ∂xβ ∂xα! manent, i.e. plastic. Following and expanding the approach by
Sirono(2004), we treatplasticity forthe purehydrostaticpres-
Todeterminerotationrateandstrainratetensorandthustheevo- sure pandthedeviatoricstresstensorSαβseparately.
lutionofthestresstensorEq.10inSPHrepresentation,theSPH
For the deviatoric stress tensor we follow the approach by
velocityderivatives∂vα/∂xβhavetobecalculated.Thestandard Benz&Asphaug(1995)andScha¨feretal.(2007)assumingthat
i i
SPH formulation,however,doesnot conserveangularmomen- ourmaterialis isotropic,whichmakesthe von Mises yieldcri-
tum due to the discretization error by particle disorder, which terion applicable. This criterion is characterised by the shear
leads to a rotational instability and in particular inhibits mod- strengthY, which in our modelis a compositeof the compres-
ellingrigidrotationofsolidbodies.Toavoidthiserror,weapply
sive andtensile strengths:Y(φ) = Σ(φ)|T(φ)|. Thesuitability
acorrectiontensorCγβ(Scha¨feretal.2007)accordingto
ofthischoicewasalreadydemonstraptedinGu¨ttleretal.(2009b).
∂vα m ∂Wij Since Y(φ) is a scalar we have to derivea scalar quantityfrom
∂xiβi =Xj ρjj(vαj −vαi)Xγ ∂xγi Cγβ, (13) SseαcβonfodrirrereadsouncisboleficnovmarpiaanratbJilit=y.SWαβeSdαoβ.tFhiinsabllyy,ctahlceurleadtiuncgtioitns
2
wherethecorrectiontensorCγβistheinverseof ofthedeviatoricstressisimplementedinthefollowingway
m ∂Wij
j(xα−xα) , (14) Sαβ → fSαβ, where f =min Y2(φ)/3J ,1 . (18)
ρ i j ∂xγ 2
Xj j i h i
The hydrostaticpressure is limited by the tensile strengthT(φ)
thatis
for p<0andbythecompressivestrengthΣ(φ)for p>0:
m ∂Wij
j(xα−xα) Cγβ =δαβ. (15)
Xj ρj j i Xγ ∂xγi p(φ)= Σ(φ) φ>φ+c . (19)
This approach leads by construction to first order consistency (T(φ) φ<φ−c
wheretheerrorsduetoparticledisordercanceloutandthecon- For φ− ≤ φ ≤ φ+ the material is in the elastic regime and Eq.
servationofangularmomentumisensured.Onlythisallowsthat 8isapcplied.φ− acndφ+ denotethefillingfactorswheretheelas-
rigidrotationcanbesimulatedcorrectly. c c
ticpathintersectsthetensilestrengthandcompressivestrength,
Finally,thesoundspeedofthematerialisgivenby
respectively(seeFig.1).
c (ρ′)= K(ρ′)/ρ′. (16) The pressure reduction process is implemented such that
s 0 0 0 at any time step p is computed with Eq. 8. If for a given φ
q
Together with Eq. 9 this relation shows that the soundspeed is p(φ) > Σ(φ) and φ > φ+ the pressure p(φ) is reduced to Σ(φ).
c
a strong function of density. This behaviour has been seen in The deformation becomes irreversible once the new reference
moleculardynamics simulations by Paszun&Dominik (2008), density ρ′ is computed through Eq. 8 and the elastic path is
0
butthereisnodataavailablefromlaboratorymeasurements. shiftedtowardshigherdensities.Herebyalsothelimitingfilling
Up to this point the set of equations describes a perfectly factors φ− and φ+ are set anew. In principalthere are two pos-
c c
elastic solid body. Additionally, the material simulated in this sible implementations for this. (1) Plasticity becomes effective
workareSiO dustaggregates,whichfeaturesahighdegreeof immediatelyandρ′ iscomputedwhenever p > Σ.(2)Plasticity
2 0
porosity and, thus, plasticity. The modificationsaccounting for becomeseffectiveafterpressuredecrease,whichisequivalentto
thesefeaturesaredescribedinthenextsection. φ<φ+.Wetestedbothimplementations.Forourunderstanding
c
R.J.Geretshauseretal.:NumericalSimulationsofHighlyPorousDustAggregates 5
p
Σ(φ) Since we do not solve the energy equation thermodynami-
callyenhancedfeatureslikeanyphasetransitionsuchasmelting
andfreezingcannotbesimulated.Thisschemealsodoesnotfea-
tureadamagemodel.Especiallyinthesectionaboutfragmenta-
c o pmlapsrteics sio n tbtiheoentare(kSseuenlctitin.nt5go.3far)cactghomeuenpntrteydseeits,ntacrlietbhouoftuioagnnhy.thfleaywmsiinghtthheamveatienrfliaulecnacnenoont
n
c o 3. Experimentalreference
sti ati
φ−c φ′0 φ+c elaorm φ′1 φ 3.1.Materialparameters
ef The material used for the calibration experiments are highly-
d porous dust aggregates as described by Blum&Schra¨pler
plastic tension T(φ) (2004),consistingofsphericalSiO2 sphereswithadiameterof
1.5 µm. For these well defined dust aggregates it is possible
to reproducibly measure macroscopic material parameters like
tensilestrength,compressivestrength,and,potentially,alsothe
Fig.1. In the modified Sirono porosity model the regime of shear strength as needed for the SPH porosity model (see sec-
elasticdeformationislimitedbythecompressivestrengthΣ(φ), tion2.3).Thetensilestrengthforthismaterialwasmeasuredfor
whichrepresentsthetransitionthresholdtoplasticcompression highlyporousandcompactedaggregates(φ = 0.15 .. 0.66)by
for p>Σ(φ),andthetensilestrengthT(φ),whichrepresentsthe Blum&Schra¨pler(2004).Thesemeasurementssupportalinear
transition threshold to plastic tension or rupture for p < T(φ). dependencebetweenthetensilestrengthandthenumberofcon-
Returningfromoneoftheplasticregimestovanishingexternal tactspermonomer(increasingwithincreasingφ),whichyields
pressure via an elastic path leads to a φ′ that differs from the thetensilestrengthas
1
initialφ′.Hence,thematerialhasbeendeformedirreversibly.
0 T(φ)=− 102.8+1.48φ Pa. (20)
(cid:16) (cid:17)
The compressive strength was measured in the experimental
counterpart of this paper (Gu¨ttleretal. 2009b) with an experi-
possibility(1)isclosertotheunderlyingphysicalprocess.Inad-
mental setup to determinethe static omni-directionalcompres-
ditionitprovedtobemorestable.Accordingtothebenchmark
sion(ODC),wherebythesampleisenclosedfromallsidesand
parameterstheresultswereequivalent.
the pressure is constant within the sample. We found that the
For the tensile regime, i.e. for φ < φ−, we do not adopt
c compressivestrengthcurvecanbewelldescribedbytheanalytic
the damage and damage restoration modelpresented in Sirono
function
(2004).Thisdamagemodelforbrittlematerialsuchasrocksor
pumicehasbeendevelopedforSPHbyBenz&Asphaug(1994, φ −φ ∆·ln10
1995) and was recently used by Jutzietal. (2008, 2009a,b). Σ(φ)= pm· φ2 −φ1 −1! , (21)
It is assumed that a material contains flaws, which are acti- 2
vated and developunder tensile loading (Grady&Kipp 1980). where the free parameters were measured to be φ = 0.12,
1
Scha¨feretal. (2007) did not find the model applicable in their φ = 0.58,∆=0.58,and p =13kPa.However,theseparame-
2 m
simulationsofporousicebecauseitincludescompressivedam- tersweremeasuredwithastaticsetupandweexpectadifferent
age effects. Brittle material like pumice and rocks tend to dis- strengthforadynamiccollision.Theparametersφ andφ deter-
1 2
integrate when they are compressed: they crush. Whereas in minetherangeofthevolumefillingfactor,whereφ isdefined
1
our material, highly porousSiO dust, tensile and compressive bytheuncompressedmaterialandφ correspondstothedensest
2 2
strengthincrease with compression.Thisis due to the fact that packing.These parametersare expected to be the same for the
themonomersareabletoformnewbondswhentheygetincon- dynamiccompression,whilep and∆aretreatedasfreeparam-
m
tact.Therefore,weadoptthesameapproachasinthecompres- eters, which we will calibrate in the forthcomingsections. The
sive regime and reducethe pressure p(φ) to T(φ) once the ten- lastimportantmaterialparametertodescribethedustaggregates
silestrengthisexceeded.Finally,thematerialcanrupturedueto is the elasticity. We can determine the Young’s modulus from
plastic flow. However, material that is plastically stretched can measurementandsimulationofthesoundspeed(Blum&Wurm
becompressedagainuptoitsfullstrength.Bychoosingthisap- 2008; Paszun&Dominik 2008), which is c = 30 m s−1. From
proacha “damagerestorationmodel”is implementedin a very this approachthe bulk modulusfor the uncompressedmaterial
naturalway. is K = ρc2 = 300kPawithρ = 300kgm−3. However,other
0 i i
Finally, a remark has to be made about energy.Apart from plausiblecalculations(Weidlingetal.2009)indicateanelastic-
energy dissipation by numerical and artificial viscosity we as- ityofrather1kPa.Wewillthereforealsovarythisparameter.
sume intrinsic energy conservation. We suppose that heat pro-
duction in the investigated physical processes is negligible.
3.2.Calibrationexperiment
Therefore,our modelis limited to a velocity regime below the
soundspeedc ofthedustmaterial(≈30m/s,seeBlum&Wurm As a setup for an easy and well-definedcalibrationexperiment
s
2008; Paszun&Dominik 2008). By choosing the approach of (see Gu¨ttleretal. 2009b), we chose a glass beadwith a diame-
modelling plasticity via stress reduction, we assume that most terof1to3mm,whichimpactsintothedustaggregatematerial
ofthe energyisdissipated byplasticdeformation,since there- with a velocity between 0.1 and 1 m s−1 under vacuum condi-
ductionofinternalstressesaccountsforthereductionofinternal tions(pressure0.1mbar).Wewereabletomeasurethedeceler-
energy. ationcurve,stoppingtime,andintrusiondepthoftheglassbead
6 R.J.Geretshauseretal.:NumericalSimulationsofHighlyPorousDustAggregates
(for various velocities and projectile diameters) and the com-
pactionofthedustundertheglassbead(fora1.1mmprojectile
withavelocityof0.65ms−1).Theseresultswillserveforcali-
bratingandtestingtheSPHcode.
For the measurementof the decelerationcurve,we used an
elongatedepoxyprojectileinsteadoftheglassbead.Thebottom
shape and the mass resembled the glass bead, while the lower
density and the therefore longer extension made it possible to
observe the projectile during the intrusion. The projectile was
observed by a high-speed camera (12,000 frames per second)
and the position of the upperedge was followed with an accu-
racyof≈ 3 µm.We foundthat–independentlyofvelocityand
projectile diameter – the intrusioncurve can be described by a
sinecurve
π t
h(t)=−D·sin ; t≤T . (22)
s
2T !
s
Here,DandT denotetheintrusiondepthandthestoppingtime
s
of the projectile. The stopping time is defined as the time be- Fig.2.Cumulativemassdistributionofthefragmentsafteradis-
tweenthefirstcontactandthedeepestintrusionatzerovelocity. ruptive collision, which can be described by a power law. The
In Sect. 5.1.2we will use the normalisedform of the intrusion divergence for low masses is due to depletion of small aggre-
curvewithD=T =1. gatesbecauseofthecameraresolution.
s
Fortheintrusiondepthwefoundgoodagreementwithalin-
earbehaviourof
A further experiment deals with the disruptive fragmenta-
D= 8.3·10−4 m2 s · mv , (23) tionofdustaggregates(fordetailsseeGu¨ttleretal.2009a,b).In
kg ! A thiscase,adustaggregatewithadiameterof0.57mmconsist-
ing of 1.5 µm spherical SiO dust with a volume filling factor
where v is the impact velocity of the projectile and m and 2
of φ = 0.35 collides with a solid glass target at a velocity of
A = πr2 are the mass and cross-sectional area of the projec-
8.4ms−1(alsoseeFig.20inGu¨ttleretal.2009b).Theprojectile
tile,respectively.WefoundthestoppingtimeT tobeindepen-
s fragments and the projected sizes of these fragments are mea-
dent of the velocity and only depending on the projectile size
suredwith ahigh-speedcamerawitharesolutionof16µmper
(3.0±0.1 ms for 1 mm projectiles and 6.2±0.1 ms for 3 mm
pixel.Asthe massmeasurementisrestrictedtothe2Dimages,
projectiles).
theprojectedareaofeachfragmentisaveragedoverasequence
Thecompactionofdustunderneaththeimpactedglassbead
of images where it is clearly separated from others. From this
wasmeasuredbyx-raymicro-tomography.Theglassbeaddiam-
projectedarea, the fragmentmasses are calculatedwith the as-
eter and velocityfor these experimentscorrespondto the com-
sumptionsof a spheric shape and an unchangedvolume filling
pactioncalibrationsetupdescribedinTable1.Thedustsample
factor. Fig. 2 shows the mass distribution in a cumulativeplot.
withembeddedglassbeadwaspositionedontoarotatablesam-
For the larger masses, which are not depleted due to the finite
plecarrierbetweenanx-raysourceandthedetector.Duringthe
camera resolution, we find good agreement with a power-law
rotationaround360◦,400transmissionimagesweretaken,from
distribution
which we computed a 3D density reconstructionwith a spatial
resolutionof21µm.Theaccordingresultsofthedensityrecon- m α
n(m)dm= dm, (24)
structioncanbefoundinGu¨ttleretal.(2009b),wherewefound
µ!
that roughly one sphere volume under the glass bead is com-
pressed to a volume filling factor of φ ≈ 0.23, while the sur- wheremisthenormalisedmass(fragmentmassdividedbypro-
rounding volume is nearly unaffected with an original volume jectilemass)andµ = 0.22isameasureforthestrengthoffrag-
filling factor of φ ≈ 0.15. In this work, we will focus on the mentation, being defined as the mass of the largest fragment
verticaldensityprofilethroughthecentreofthesphereandthe divided by the mass of the projectile. The exponent α has a
compressedmaterial(seesection4). value of 0.67. A similar distribution was already described by
Blum&Mu¨nch (1993) for aggregate-aggregate fragmentation
ofZrSiO aggregateswithacomparableporosity.
4
3.3.Furtherbenchmarkexperiments
In this section, we will present two further experiments which
4. NumericalIssues
will be used for the validation of the SPH code in sections 5.2
and 5.3. Heißelmannetal. (2007) performed low-velocity col- Before we performthe calibration process, some numericalis-
lisions (v = 0.4 m s−1) between cubic-shaped, approx. 5 mm- sues have to be resolved. For instance, it is unfeasible to sim-
sized aggregates of the material as described in Sect. 3.1 and ulate the dustsample, into which the glass bead is droppingin
found bouncing whereby approx. 95% of the energy was dis- thecompactioncalibrationexperimentpresentedinsection3,as
sipatedinacentralcollision.Detailedinvestigationofthecom- a whole. Itis also unfeasibleto carryoutallnecessary compu-
pactioninthesecollisions(Weidlingetal.2009)revealedsignif- tationsin 3D.Therefore,we willonlysimulate partofthedust
icantcompactionoftheaggregates(fromφ = 0.15toφ = 0.37) sampleandexploreatwhichsizespuriousboundaryeffectswill
afterapprox.1,000collisions.Theenergyneededforthiscom- emerge.Most of the calibration process has been conductedin
pactionisconsistentwiththeenergydissipationasmeasuredby 2D, butthe differencesbetween2D and 3D results will be dis-
Heißelmannetal.(2007). cussedandquantified.
R.J.Geretshauseretal.:NumericalSimulationsofHighlyPorousDustAggregates 7
sphere PhysicalQuantity Symbol Value Unit
particles
Glassbead
Bulkdensity(∗) ρ 2540 kgm−3
0
dust boundary Bulkmodulus(∗) K 5×109 Pa
0
particles particles Murnaghanexponent(∗) n 4 -
v Radius r 0.55×10−3 m
0
Impactvelocity v 0.65 ms−1
0
r
sample Dustsample
Initialdensity ρ 300 kgm−3
i
Bulkdensity ρ 2000 kgm−3
s
Referencedensity ρ′ 300 kgm−3
0
Initialfillingfactor φ 0.15 -
i
vertical Bulkmodulus K 3×105 Pa
0
density profile ODCmeanpressure pm 1300 Pa
ODCmax.fillingfactor φ 0.58 -
2
ODCmin.fillingfactor φ 0.12 -
Fig.3. Compaction calibration setup in 2D or, respectively, 1
ODCslope ∆ 0.58 -
cross-sectionof3Dcompactioncalibrationsetup.Aglasssphere
Table 1. Numerical parameters for the compaction calibration
impactsintoadustsample(radiusr )withinitialvelocityv .
sample 0
setup. ODC stands for omni-directional compression relation
Thedustsampleissurroundedbyboundaryparticlesatitsbot-
(Eq.21).Quantitiesmarkedby(*)representtheparametersfor
tom.Theiraccelerationissettozeroateverytimestep.Thever-
sandstoneinMelosh(1989)whichweadoptforglasshere.
ticaldensityprofileatmaximumintrusionservesascalibration
feature.
In this contextwe make use of 2D simulations in cartesian Ford <3.3mmwefindspuriousdensitypeaksatthelower
sample
coordinatesalthoughthesymmetryoftheproblemsuggeststhe boundaries(D≈−1.4mmandD≈−2.2mm).
usage of cylindrical coordinates. However, the SPH scheme in In order to reduce the computation time we simulated the
cylindricalorpolarcoordinatesbattleswiththeproblemofasin- dust sample as a semicircle with the same radius variation as
gularityattheoriginofthekernelfunction.Thereexistonlyfew above. The resulting density profiles are shown in Fig. 4 (bot-
attemptstosolvethisissue(e.g.Omangetal.2006),buttheyare tom).Incontrasttothecorrespondingsimulationswiththebox-
stillunderdevelopmentandrequirehighimplementationefforts. shaped samples we find for r ≤ 1.375mm an increased
sample
Sinceinourcase2Dsimulationsonlyserveasindicatorforcal- maximum filling factor and a slightly reduced intrusion depth.
ibrationand 3D simulations are aimed at, we stick to cartesian Due to the greater amount of volume lateral to the intrusion
coordinates. channel,materialcanbepushedasidemoreeasilythaninsidethe
TheglassbeadissimulatedwiththeMurnaghanequationof narrow boundaries of the semicircle. Therefore, a higher frac-
state tion of the material is compressed to higher filling factors. For
K ρ n rsample > 3.3mm the spurious boundary effects become negli-
p(ρ)= 0 −1 , (25) gible within the compaction calibration setup and the density
(cid:18) n (cid:19)" ρ0! # structure shows no significant difference for box-shaped and
followingthe usual laws of continuummechanicsas presented semicircleshapeddustsamples.
in section 2.2. The compaction calibration setup is initialised Hence, all computations of section 4 are conducted on the
with the numerical parameters shown in Table 1, unless stated basisofasemicirclein2Dorahemispherein3Dwitharadius
otherwise in the text. Our tests showed that the maximum in- ofr =3.3mm.
sphere
trusiondepthandthedensityprofilearethecalibrationparame- In all cases the dust sample is bordered by a few layers of
terswhicharemostsensitivetochangesofthenumericalsetup. boundaryparticles. The acceleration of these particles is set to
Density profiles (e.g. Figs. 4 and 5) display the filling factor φ zeroateachintegrationstep,simulatingreflectingboundarycon-
alonga linethroughthecentreofthe sphereandperpendicular ditions.Apartfromthat,i.e.intermsofequationofstate,theyare
to the bottomof the dustsample (see Fig. 3). Theoriginin the treatedlikedustparticles.Wealsotesteddampingboundarycon-
diagramsrepresentsthesurfaceoftheunprocesseddustsample. ditionsbysimulatingtwolayersofboundaries.Theouterlayer
Theglassspherehasbeenremovedinthesefigures. was treated as described above, the inner (sufficiently large)
layerwassimulatedwithahighartificialα-viscosity.Sincethere
was no significant difference in the outcome we fix all bound-
4.1.Computationaldomainandboundaryconditions
ariesintheaforementionedwayandtreatthemasreflecting.
In2Dsimulationswetestedtheeffectofchangingsizeandshape
ofthedustsample.Initiallytheparticlesweresetonatriangular
4.2.ResolutionandConvergence
latticewithalatticeconstantof25µm.Tobegeometricallycon-
sistentwiththecylindricexperimentalsetup,firstly,weutilised Jutzi (2008) found in his studies of a basalt sphere impacting
a box (width 8 mm) and varied its depth: 1.375 mm, 2.2 mm, intoaporoustargetthattheoutcomeofhissimulationsstrongly
3.3 mm, and 5.5 mm. This is equivalent to 2.5×,4×,6×,and dependsontheresolution.Withacalibrationsetupsimilartothe
10 × r . Comparing the density profiles (Fig. 4, top), two one used in this paper Geretshauser (2006) confirms that also
sphere
features are remarkable:(1) The maximum filling factor at the in simulations of the type presented in this work a strong res-
top of the dustsample (φ ≈ 0.27at D ≈ −0.6mm)and the in- olution dependence is present. He has found that the intrusion
trusiondepthDisnearlythesameforalldustsamplesizes.(2) depthofthe glassbeadcan bedoubledbydoublingthe resolu-
8 R.J.Geretshauseretal.:NumericalSimulationsofHighlyPorousDustAggregates
0.4 0.4
d = 1.375 mm l = 100 µm
dsample = 2.2 mm lc = 50 µm
0.35 dsample = 3.3 mm 0.35 lc = 25 µm
dsample = 5.5 mm lc = 12.5 µm
0.3 sample 0.3 c
φ φ
or 0.25 or 0.25
t t
c c
fa 0.2 fa 0.2
g g
n n
i 0.15 i 0.15
l l
l l
i i
f 0.1 f 0.1
0.05 0.05
0 0
-5 -4 -3 -2 -1 0 -2.5 -2 -1.5 -1 -0.5 0
sample height [mm] intrusion depth [mm]
0.4 0.4
r = 1.1 mm l = 100 µm
rsample = 1.375 mm lc = 50 µm
0.35 rsample = 2.2 mm 0.35 lc = 25 µm
sample c
r = 3.3 mm
0.3 sample 0.3
φ rsample = 5.5 mm φ
or 0.25 or 0.25
t t
c c
fa 0.2 fa 0.2
g g
n n
i 0.15 i 0.15
l l
l l
i i
f 0.1 f 0.1
0.05 0.05
0 0
-5 -4 -3 -2 -1 0 -2.5 -2 -1.5 -1 -0.5 0
sample height [mm] sample height [mm]
Fig.4. Vertical density profile at maximum intrusion for the Fig.5.Convergencestudy ofthe verticaldensityprofilefor2D
compactioncalibrationsetupanddifferentshapesofthe2Ddust (top)and3D(bottom)compactioncalibrationsetups.Thefilling
sample(boxandsemicircle).Depthd ofan8mmwidebox factorincreasetowardsthe surfaceofthe dustsampleaccounts
sample
(top)andradiusr ofthesemicircle(bottom)werevaried.In fortheglassbeadwhichisnotremovedinthisplot.Simulations
sample
bothcasesspuriousboundaryeffectsappearford <3.3mm wereperformedfordifferentspatialresolutions.Allcurvesshow
sample
andr <3.3mm,respectively. acharacteristicfillingfactorminimumbetweensphereanddust
sample
sample and a characteristic filling factor maximum indicating
thedustsamplesurface.
tion.SincethecalibrationexperimentspresentedinGu¨ttleretal.
(2009b) are extremely sensitive even to minor changes of the
setup,theconvergencepropertiesofporositymodelandunderly- rialanddustmaterialhavetobeseparatedbyartificialviscosity
ingSPHmethodwillbeinvestigatedcarefullyinthisparagraph. due to stability reasons. This will be discussed below. (2) The
Additionally,wewillstudythedifferencesintheoutcomeof2D volumeofthesphererepresentsanareaofextremelyhighden-
and3Dsetup. sity and pressure with respect to the dust sample. This area is
Forthe2Dconvergencestudyparticleswereinitiallyputona smoothedoutduetotheSPHmethod.Thewidthofthesmooth-
triangularlatticeagain.Thelatticeconstantsl were100,50,25, ing is givenby the smoothinglength.Althougha clear conver-
c
and12.5µmforthecompactioncalibrationsetup.Thesmooth- gence behaviour is visible in Fig. 5 for both the 2D and the
inglengthhwaskeptconstantrelativetol ataratioof5.6 × l . 3D case, a moreunique convergencecriterion has to be found.
c c
The maximum numberof interaction partnerswas I ≈ 180, Forthispurposewechoosethemaximumintrusiondepthwhich
max
theaverageI ≈100andtheminimumI ≈30. provedtobeverysensitivetoresolutionchanges.Theshapeof
av min
Inthe3Dconvergencestudyweareusingacubiclatticewith thefillingfactorprofileofferstwochoicestodeterminetheintru-
edgelengthsl =100,50,and25µm.Thelatterwassimulated siondepth:(1)Thefillingfactorminimumwhichrepresentsthe
c
with3.7millionSPHparticles,whichrepresentthelimitofour middle betweensphere and dust sample and (2) the filling fac-
computationalresources.Wefixedh = 3.75 × l whichyielded tormaximum(peak)ofthedustmateriallefttothegapbetween
c
I ≈ 370,I ≈ 240,and I ≈ 70.Theresultsarepresented sphereanddustsample.
max av min
inFig.5.IncontrasttotheplotsinFig.4theglasssphereisnot Fig. 6 shows the results for both cases in 2D (top) and 3D
removedhere.Comingfromtherightsideoftheplot,thefilling (bottom).The error bars aroundthe minimumvalues represent
factorrapidlydecreasesfroma highvalue outsidethe scopeof thesmoothinglengthandgiveanindicationofthemaximumer-
theplotindicatingthesphere.Thefillingfactorreachesitsmin- ror. The position of the density peak remains almost constant,
imumatanartificialgapbetweensphereandsurfaceofthedust convergingto ≈ −0.9mm(2D) and ≈ −0.65mm(3D), respec-
sample. The width of this gap is about one smoothing length tively, for higher resolutions. The position of the density min-
h. The existence of the gap has two reasons: (1) Sphere mate- imum at low resolutions significantly differs from the position
R.J.Geretshauseretal.:NumericalSimulationsofHighlyPorousDustAggregates 9
0 0.4
2D density min h = 0.05 mm
2D density peak 0.35 h = 0.075 mm
] -0.2 h = 0.1 mm
m
h = 0.125 mm
m 0.3
n [ -0.4 φor 0.25 hh == 00..1157 5m mmm
o t
i c
us -0.6 fa 0.2
tr g
n n
max. i -0-.18 filli 0 .01.51
0.05
-1.2
0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0
-1 -2.5 -2 -1.5 -1 -0.5
0
3D density min sample height [mm]
3D density peak
] -0.2 Fig.7. Convergence study for the density profile using the 2D
m
m setupandvaryingthesmoothinglengthh.Throughthisvariation
n [ -0.4 thenumberofinteractionpartnersisvariedaccordingtoTable2.
o Theglassbeadhasbeenremovedinthisplot.Forh≤0.075mm
i
us -0.6 clearsignsofinstabilitiesarevisible.Forh≥0.1mmthefilling
tr factorhasthesamevalueanditspositionremainsconstant.The
n
x. i -0.8 smoothing of the boundaryof the dust sample is increased for
a increasingh.
m -1
Table 2. Parametersfor the convergencestudy regardinginter-
-1.2
actionnumbers
0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09
resolution [µm-1]
h h/l I I I T N
c min av max comp steps
Fig.6. Convergencestudy of the maximum intrusion depth for 0.050mm 2 3 13 25 16.2h 132401
0.075mm 3 10 30 55 14.8h 96912
the 2D (top) and 3D (bottom) compaction calibration setups.
0.100mm 4 16 53 92 19.8h 71175
Filled symbols represent the position of the filling factor peak
0.125mm 5 23 82 142 19.0h 56347
ofthedustmaterial,whereasemptysymbolsdenotetheposition
0.150mm 6 32 116 205 21.6h 46782
oftheminimumfillingfactoratthegapbetweenglassbeadand
0.175mm 7 43 158 274 24.3h 39980
dust material. The values are derived from the density profiles
in Fig. 5. The smoothing length is indicated by the error bars.
Whilethepeakpositionremainsalmostconstantat≈ −0.9mm lengthshcanbefoundinTable2.Additionally,wemeasuredthe
(2D) and ≈ −0.65mm (3D) with increasing spatial resolution, computationtimeT thatthesimulationstookon4coresofa
comp
thepositionofthe fillingfactorminimumquicklyconvergesto clusterwithIntelXenonQuad-Coreprocessors(2.66GHz)fora
thesamevalue.Thisisduetotheartificialseparationofdustand simulatedtimeof5msandthenumberofintegrationstepsN
step
glassmaterials,whichisintheorderofasmoothinglength. ofouradaptiveRunge-KuttaCash-Karpintegrator.
Comparing the density profiles in Fig. 7 (where the glass
bead has been removed) instabilities in the form of filling fac-
ofthedensitypeak,butconvergesquicklytothesameintrusion tor fluctuations due to insufficient interaction numbers appear
depthwithhigherresolutions.However,thedifferencesbetween for smoothing lengths h ≤ 0.075mm, i.e. for I ≤ 30. For
av
the extrema remain well within one smoothing length. This is h ≥ 0.1mm the density profile maintains essentially the same
duetothefactthatsphereanddustsampleareseparatedbyabout shape:Thepositionandheightofthefillingfactorpeakremains
one smoothing length. Comparing2D and 3D convergencethe nearlythesameandφsmoothlydropsto≈0.18towardsthebot-
3Dcaseseemstoconvergemorequickly. tomofthedustsample.Onlythesharpedgeatthetopofthedust
Duetothefindingsofthisstudywechooseaspatialresolu- sampleissmoothedoutoverawiderrangeduetotheincreased
tionofl = 25µmforfurthersimulationsin2D.Inthe3Dcase smoothinglength.
c
l = 50µmissufficient,butl ≤ 50µmisdesirableifitisfeasi- Table2showsthatthenumberofintegrationstepsN isde-
c c step
ble. Afterdefiningsuitablevaluesforthespatialresolutionwe creasingwithincreasinginteractionnumbers.Thisisduetothe
nowturntothenumericalresolution,whichfortheSPHscheme factthatelastic wavesinside the dustsample are smoothedout
isgivenbythenumberofinteractionpartnersofeachsinglepar- overawiderrangecausingtheadaptiveintegratortoincreasethe
ticle.Fortheinvestigationofthisfeatureweperformedastudy durationofatimestep,sincedensityfluctuationsdonothaveto
utilisingthe2Dcompactioncalibrationsetupwithaspatialres- beresolvedassharplyasatsmallersmoothing.Asexpected,the
olution of l = 25µm and varied the ratio between smoothing computationtime T is generallyincreasingwith increasing
c comp
length and lattice constant h/l from 2 to 7 in steps of one. numberofinteractions.Therearetwoexceptions:h=0.075mm
c
h/l determinestheinitialnumberofinteractionpartnersthatis and h = 0.125mm. Here, the decrease of N overcompen-
c steps
smoothedover.Theresultingmaximum,average,andminimum satestheincreaseofinteractionsleadingtoadecreaseofT .
comp
interactionsI ,I ,andI andthecorrespondingsmoothing Hence,aratioh/l ≈5yieldsthenecessaryaccuracyandanac-
max av min c
10 R.J.Geretshauseretal.:NumericalSimulationsofHighlyPorousDustAggregates
ceptable amount of computation time. This study also justifies
0
the choice of h/l = 5.6 in Gu¨ttleretal. (2009b) and we will 2D peak corr
c
sticktothisratiothroughoutthispaper. 3D peak
] -0.2 2D density peak
m
Accordingtothesefindings,for3Dsimulationstheoretically m
an average interaction number of Ia3v/2 ≈ 750 would be needed n [ -0.4
to achieve the same numerical resolution. However, alike sim- o
i
ulations are unfeasible and our choice of Iav ≈ 240 in 3D is us -0.6
equivalenttoIav ≈40in2Dwhichshouldprovidesufficientand ntr
reliableaccuracy. x. i -0.8
a
m -1
4.3.Geometricaldifference–2Dand3Dsetups
-1.2
0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09
AsonecaneasilyseeinFig.62Dand3Dsimulationshavesig-
nificantly different convergence values for the intrusion depth. resolution [µm-1]
Thisdeviationiscausedbythegeometricaldifferenceofthe2D
and3Dsetup.The2Dsetup(glasscircleimpactsintodustsemi-
Fig.8. Verification of the 2D-3D correction factor. Filled sym-
circle) represents a slice through a glass cylinder and a semi-
bols denote the position of the filling factor peak of the dust
cylindricaldustsample,whichimpliesaninfiniteexpansioninto
material in Fig. 5. Triangles represent 3D and squares 2D val-
thethirdspatialdirection.Whereasthe3Dsetuprepresentsareal
ues.Theconversionfrom2Dto3Dintrusiondepthutilisingthe
spheredroppingintoa“bowl”ofdust.Therelationfortheintru-
correctionfactorfromEq.26andEq.23duetothegeometrical
siondepthfoundbyGu¨ttleretal.(2009b)(seeSect.3)contains differenceisindicatedbythelinewithoutsymbols.The3Dval-
ageometricaldependence:D ∝ mv/A,whereDistheintrusion
uesareinverygoodagreementwith theveryroughtheoretical
depth,m themassoftheimpactingglassbead,vits impactve-
prediction.Theyliewellwithintheerrors.
locity, and A its cross-section. m is a mass per unit length.
2D
Gu¨ttleretal. (2009b) already exploited this relation in order to
determinearoughcorrectionfactorbetween2Dand3Dsimual- 4.4.ArtificialViscosity
tionalsetups:
Since artificial viscosity plays an eminent role for the stability
ofSPHsimulationsweinvestigateitsinfluenceontheoutcome
of our compaction calibration setup. Only artificial α-viscosity
m v 4πr3ρ·v 8 πr2ρ·v 8 m v is applied in ourtest case, but in a threefoldway: (1) It damps
3D = 3 = = 2D . (26)
A πr2 3π 2r 3π A high oscillation modes of the glass bead caused by the stiff
3D 2D
Murnaghanequationofstate(EOS).Therebyitenlargesthetime
step of ouradaptiveintegratorandsaves computationtime. (2)
Hence,the2Dintrusiondepthhastobecorrectedbyafactorof Itisused toprovidethedustmaterialwith a basic stability.(3)
≈ 8 inordertodeterminethe3Dintrusiondepth.Thecompar- ItseparatestheareasofMurnaghanEOSanddustEOSandpre-
3π
ison is shown in Fig. 8. Again we choose the 2D and 3D data ventsasocalled“cannonballinstability”.Forallthreecasesalso
gainedintheconvergencestudyforthepeakfillingfactorvalues theinfluenceofβ-artificialviscosity hasbeentested,butitsin-
showninFig.6.Fig.8showstheoriginal2Ddata,thecorrected fluenceonallbenchmarkparameterswasnegligible.
2D data,andthe according3D data(with errorbarsdisplaying (1) The choice of the α-viscosity of the glass bead proved
the smoothing length). The 3D values nicely follow the rough to be very uncritical. We choose the canonical value α = 1.0.
correction and remain well within the maximum error. This There was no influence on the physical benchmarkparameters
comparisonjustifiesthecorrectionoftheresultsinGu¨ttleretal. for all α values, except α = 0 which produces an instability.
(2009b). With the aid ofthis – nowverified– correctionfactor Values for α > 1.0 show no significant effect on the damping
the2Ddataoftheintrusiondepthcanbeconvertedinto3Ddata andtheinfluenceof0.1 < α < 1.0onitisnottoohigh,butstill
andallcalibrationtestsinvolvingtheintrusiondepthcanbecar- observable.Hence,westicktothecanonicalvalue.
riedoutin2D,whichsavesasignificantamountofcomputation (2)Sirono(2004)appliesnoartificialviscositytohisporous
time. Comparingthe verticaldensityprofilesof2D and 3D se- ice material because of its spurious dissipative properties. Our
tupsinFig.5resolvesalsoanotherissuethatremainedunsolved findings, shown in Fig. 9, confirm this assessment regarding
inGu¨ttleretal.(2009b).Accordingtotheexperimentaldatathe the dust material. Within our 2D compaction calibration setup
fillingfactordropstoavalueofφ≈0.16within≈0.6mmaway (l = 25µm,h/l = 5.6)wevaryαfrom0to 2andobserveits
c c
fromthebottompointoftheglassbead.Forhigh-resolution2D influence on the density profile (Fig. 9, top) and the maximum
simulations (Fig. 5, top) the filling factor does not drop to this intrusionrepresentedbythefillingfactorpeakofthedustmate-
value within the whole dust sample. However, the 3D simula- rial(Fig.9,bottom).
tions (Fig. 5, bottom) show that this effect again is due to the The position of the filling factor peak ranges from ≈
difference of 2D and 3D geometry. With the 3D setup the fill- −0.92mmwithα = 0.0to≈ −0.62mmatα = 2.0.Thisclearly
ing factor drops to φ ≈ 0.16 within ≈ 0.9mm. All deviations demonstrates the dissipative feature of the α-viscosity, since a
from experimental findings due to this effect, in particular the smalleramountofkineticenergyoftheglassbeadistransformed
presenceofalargeamountofvolumewithφ≤0.2inthecumu- intoplasticdeformationwithhigherα.Theresidualenergymust
latedvolumeoverfillingfactordiagram(Fig.15inGu¨ttleretal. have been dissipated. However, the α-viscosity-intrusioncurve
2009b),caninprincipalberemovedbyswitchingto3Dsimula- seems to saturate at a value of ≈ −0.6mm. The decrease of
tions.Theydonotrepresentafundamentalerrorintheporosity the maximum intrusion can also be seen in the density profile
model. (Fig.9,top).Whiletheprofilemaintainsnearlythesameshape,