Numerical determination of the material properties of porous dust cakes
aa r X i v : . [ a s t r o - ph ] F e b Astronomy&Astrophysicsmanuscript no. draft c (cid:13)
ESO 2018October 22, 2018
Numerical determination of the material properties of porous dustcakes.
Paszun D. ⋆ and Dominik C. Astronomical Institute ”Anton Pannekoek”, University of Amsterdam, Kruislaan 403 1098 SJ AmsterdamReceived < date > / Accepted < date > ABSTRACT
The formation of planetesimals requires the growth of dust particles through collisions. Micron-sized particles must grow by manyorders of magnitude in mass. In order to understand and model the processes during this growth, the mechanical properties, and theinteraction cross sections of aggregates with surrounding gas must be well understood. Recent advances in experimental (laboratory)studies now provide the background for pushing numerical aggregate models onto a new level. We present the calibration of apreviously tested model of aggregate dynamics. We use plastic deformation of surface asperities as the physical model to bringcritical velocities for sticking into accordance with experimental results. The modified code is then used to compute compressionstrength and the velocity of sound in the aggregate at di ff erent densities. We compare these predictions with experimental results andconclude that the new code is capable of studying the properties of small aggregates. Key words.
Dust – Aggregates – Interplanetary Dust
1. Introduction
It is commonly accepted that planets form in protoplanetarydisks. Giant planets are believed to be produced via one of twocompeting scenarios. The first one is a gravitational instabil-ity (Boss 2007). Clumps of gas (intermixed with dust) may be-come gravitationally bound and contract forming giant planets.This model is hotly debated (i.e. Pickett et al. 2000; Boley et al.2006, 2007; Mayer et al. 2002, 2004, 2007). The second, com-monly accepted scenario is known as the core accretion model(Pollack et al. 1996). In this mechanism the solid core of a planetforms first and gas is accreted onto that core later. The core itselfis believed to form by collisional accumulation of planetesimals.Terrestial planets, due to their similarity to the giant planets’cores, must form out of planetesimals with help of gravitationalinteractions. The formation of planetesimals is also a subjectof debate. The scenario initially proposed by Goldreich & Ward(1973) assumed a laminar structure of the disk. When the dustsublayer reaches a critical thickness, a gravitational instabil-ity may form planetesimals from the dense and gravitationallybound concentrations of dust. However, this mechanism requiresan extremely laminar nebula. The shearing motion of the dustlayer causes a Kelvin-Helmholtz instability and therefore limitsthe thickness of the sublayer. The gravitational instability of thedust layer is then prevented. This was proved by Weidenschilling(1984), Cuzzi et al. (1993) and Johansen et al. (2006a).The core accretion model requires to grow micron sizeddust grains into km sized planetesimals. These are over 27 or-ders of magnitude in mass. Radial drift alone may move dustparticles inwards onto the central star before they reach big-ger size. The drift can move meter sized particles all the wayin within only 100 orbits (Weidenschilling 1977). Relative ve-
Send o ff print requests to : D. Paszun, e-mail: [email protected] ⋆ Present address:
Anton Pannekoek Institiute, University ofAmsterdam, Kruislaan 403, 1098SJ Amsterdam The Netherlands locities of large bodies set another barrier against accumulationof large particles (Weidenschilling & Cuzzi 1993; Wurm et al.2005a; Ormel & Cuzzi 2007). A recent study by Johansen et al.(2006b) showed however, that planetesimals might be producedby a gravitational collapse in the presence of turbulence. Metersized boulders concentrated in high pressure turbulent eddies inthe disk may form gravitationally bound clumps. Even if thisprocess really can be made to work, dust particles must first growover 18 orders of magnitude in mass to be able to form that denseaccumulations. This growth must happen through the collisionalsticking of dust particles.The initially small, micron sized, dust grains (further referredto also as monomers ) collide and stick to each other due to at-tractive surface forces (Johnson et al. 1971). The fine dust isvery well coupled to the gas, meaning that the particles movecollectively. The relative velocities, due to the Brownian mo-tion (Brown 1828), are very small. These conditions lead toa quasi-monodisperse growth that preferentially collides parti-cles of similar size. The aggregates (further referred to alsoas agglomerates or particles) formed in this case are frac-tal (Kempf et al. 1999; Blum et al. 2000; Krause & Blum 2004;Paszun & Dominik 2006).The Brownian growth produces aggregates with very openstructure and the fractal dimension D f = . ff ected by turbulence, sedimentation andradial drift (Weidenschilling 1977, 1980, 1984). The relative ve-locities thus increase and collisions occur preferentially between D. Paszun & C. Dominik: Properties of porous aggregates agglomerates of di ff erent sizes (Weidenschilling & Cuzzi 1993).When the collision energy becomes higher than some thresholdrestructuring energy, aggregates are compacted and their frac-tal dimension approaches D f = ff ectively and therelative velocities increase further. Dominik & Tielens (1997)and Blum & Wurm (2000) showed that the outcome of collisionscan be categorized in terms of collision energy. Perfect stickingwithout restructuring occurs when the collision energy is lowerthan the threshold rolling energy, which is the energy needed toroll a single particle over an angle of π/
2. When the collisionenergy reaches this limit, the aggregates start to compress uponimpact. The maximum compaction occurs when the collision en-ergy equals the rolling energy times the number of contacts inthe aggregates. The particles start to lose monomers, when theimpact energy reaches the number of contacts in the aggregatetimes the energy needed to break one contact. Catastrophic dis-ruption occurs when the energy reaches several times the num-ber of contacts times the breaking energy. Therefore, as particlesgrow and decouple from the gas, eventually relative velocitiesare reached that will lead to shattering of the colliding aggre-gates (Ormel & Cuzzi 2007).This general picture has several missing elements. Althoughnumerous experiments have been performed in di ff erent sizeand velocity regimes (Krause & Blum 2004; Blum & Wurm2000; Wurm et al. 2005a,b, for a review see Blum & Wurm,ARAA,2008 in press), our understanding of the processes in-volved in the growth of aggregates is still incomplete. The re-structuring mechanism for instance is understood only qualita-tively. The degree of the collisional compaction still remains amystery. This, however, is crucial for the growth of the metersized aggregates. The sticking e ffi ciency as a function of the par-ticle density along with the density evolution must be determinedin order to understand quantitatively the growth of meter sizedboulders.The low strength of the aggregates due to the fractal struc-ture leads to fragmentation in the case of fast impacts. The dis-tribution of fragment sizes must be also determined as a func-tion of the collision energy. This makes it possible to keeptrack of the realistic size distribution in the disk. Small par-ticles, if not replenished in the disk, are very quickly sweptup by larger grains (Dullemond & Dominik 2005). Thus thesmall particles should be supplied to the disk by some pro-cesses. Dominik & Dullemond (2008) show that the infall ofgrains from the parent cloud is rather unrealistic and requiresfine tuning of parameters to reproduce observational data. Thuscollisional fragmentation is the most likely mechanism that canexplain the population of small grains in protoplanetary disks.However, long before fractal aggregates approach velocities highenough to cause fragmentation, they undergo collisional com-pression (Blum & Wurm 2000; Weidenschilling & Cuzzi 1993).Thus the only scenario leading to fragmentation of fractals iscollision with much larger non–fractal particle.The current understanding of processes like sticking, bounc-ing and fragmentation of aggregates is poor. The understandingof this processes on micro scales may fill these gaps and allowextrapolations to the larger sizes. This may be resolved by twoapproaches. Experiments may be performed in a laboratory toprovide useful data. However, this way is available only for ag-gregates below centimeter in size. Larger particles are currentlyinaccessible for experiments. The second, theoretical approach,understanding the material properties of porous matter reach ag-gregates of meter size and beyond. Thus it would be ideal to provide theoretical predictions for centimeter sized and largeraggregates.Sirono (2004) developed a model capable to simulate big– meter sized and larger aggregates, using Smoothed ParticlesHydrodynamics (SPH). In this case one particle in the model isan aggregate characterized by material properties like, compres-sive strength, tensile strength, density and sound speed. In orderto obtain some of these properties he fitted power-law functionsto experimental data of compression and tensile strength.This method was then also used by Sch¨afer et al. (2007).Further studies require, however, a more precise determinationof the material properties of porous bodies. These may be ob-tained in laboratory experiments or in computer simulations aspresented in this work.More realistic Solar Nebula dust analogs were used in ex-periments by Blum et al. (2006). They investigated aggregatesmade of micron sized dust particles. Moreover di ff erent aggre-gates used in these experiments consisted of spherical or irreg-ularly shaped monomers. The compression and tensile strengthcurves for this dust cakes were determined.Although experiments give a quantitative description of pro-cesses in cm sized aggregates, it is di ffi cult to access for small,a few microns size aggregates. Aggregate dynamics models(Dominik & Tielens 1997; Dominik & N¨ubold 2002) are idealfor simulations of these small scale structures. This method spa-tially resolves single monomers, which is certainly needed to un-derstand physics of bigger – meter sized bodies. Until now themain drawback of this aggregate dynamics model was missingquantitative agreement with experiments, even though the qual-itative agreement has already been established (Blum & Wurm2000). Thus a calibration of this model is required.Another aggregate dynamics model was recently presentedby Wada et al. (2007). They made use of potential energies in or-der to derive forces acting between grains in di ff erent degrees offreedom. In this case they present just a 2D case, but their resultsare in agreement with findings of N–body model presented ear-lier by Dominik & Tielens (1997). Wada et al. (ApJ,2007; sub-mitted) shows the first results of compression and disruption of3 – dimensional aggregates in head – on collisions. Althoughthe model is qualitatively in agreement with experiments, as isstudied by Dominik & Tielens (1997), the quantitative mismatchis still present. Wada et al. (2007) do not show any solution orworkaround to the quantitative disagreement between theory andlaboratory experiments. The work presented here addresses thisissue and provides mechanisms that fit our model to the empiri-cal data.In this paper we present the calibration of the aggregatesdynamics model developed by Dominik & Tielens (1997) andDominik & N¨ubold (2002). We fit the code using experimentaldata and further test it. Modifications that are implemented inthe code are presented together with a few possible applicationof the model.
2. The model
In order to study the agglomeration mechanisms involved in thegrowth of planetesimals, we use the SAND code (Soft AggregateNumerical Dynamics) developed by Dominik & N¨ubold (2002).It is a N–body model of a system of spherical particles interact-ing via surface forces.Two monomers feel the attractive force only when they arein contact. The surfaces of the particles deform and form a con-tact area. When the particles are pulled outwards, increasingthe relative distance, the area decreases and the monomers are . Paszun & C. Dominik: Properties of porous aggregates 3 pulled back to the equilibrium position by the surface forces.The compression of the system on the other hand leads to in-crease of the contact area and repulsive force pushing the parti-cles apart. A detailed description of the surface forces was pro-vided by Johnson et al. (1971) (further referred to as JKR). Theinfluence of adhesion forces on the contact between particleswas also studied by Derjaguin et al. (1975). The model is ableto treat long range magnetic forces (Dominik & N¨ubold 2002),however these are not the subject of our current study.There are two main processes that govern the events dur-ing a collision. The first one is breaking a contact between twomonomers. This dissipates part of the energy and weakens or de-stroys aggregates participating in the collision. The second pro-cess is a rolling motion of a monomer over another grain. Thisalso dissipates energy and causes a restructuring of an aggre-gate. This restructuring may be attributed to both compression,which results in strengthening the aggregate, or decompression,i.e. weakening the aggregate’s structure.The JKR theory predicts a critical force that is needed toseparate two particles. This prediction was tested for the case ofmicron size Silica spheres by Heim et al. (1999). Two monomerswere pulled o ff each other using an Atomic Force Microscope.The measured pull-o ff force was in agreement with the force pre-dicted by the JKR theory. The pull-o ff force is F c = πγ R , (1)where γ is a surface energy and R = (cid:16) R + R (cid:17) − is a reduced ra-dius of the spheres in contact. Thus the results of the experimentsconfirm the theoretical predictions.A similar experiment was performed in order to determinethe horizontal forces acting on the particles in contact (first stud-ied theoretically byDominik & Tielens (1997)). The horizontaldisplacement of the contact zone causes a torque acting againstthe displacement. The resulting rolling friction was measured inlaboratory experiments, again using an AFM (Heim et al. 1999). Rolling friction is one of the most important energy dissipationchannels in the restructuring of aggregates (Dominik & Tielens1997). It is thus of a great importance to treat it in a correct way.The theoretical derivation by Dominik & Tielens (1997) showsthat the energy associated with initiating rolling is expressed as e roll = πγξ , (2)where ξ crit is a critical displacement at which the inelastic behav-ior occurs and energy is dissipated. ξ crit was initially assumed tobe of the order of inter – atomic distances (Dominik & Tielens1997). The rolling motion causes shift of the contact area, im-plying that at one end the contact has to be broken in order toform it at the other side of the contact area. The experimentallydetermined friction showed that the critical displacement mustbe larger. The determined value was ξ crit = . ξ crit =
10 Å. This howeverwas applied to di ff erent material than the one investigated in thelab. We followed Heim et al. (1999) and applied larger value of ξ crit =
20 Å, which is approximately 10 inter atomic distances.
At low velocities the dominating mechanism that dissipatesenergy is rolling. At higher impact speeds, however, otherchannels become more important. Chokshi et al. (1993) andDominik & Tielens (1997) calculated how much of the initialenergy is dissipated in the collision between two particles. Thisgives the maximum energy at which the particles can stick in ahead–on collision. The critical velocity is given by v crit = . γ / E ∗ / R / ρ / , (3)where γ , R and ρ are surface energy, reduced radius and massdensity, respectively. E ∗− = − ν E + − ν E is a reduced elastic-ity modulus with ν i and E i being Poisson’s ratio and Young’smodulus, respectively, of grain i . When eq. 3 is applied to a R = . µ m silica grain, impacting flat silica surface, it givesthe critical velocity v crit = .
18 ms − . This was again tested inexperiments. Poppe et al. (2000) showed that such a particle canstick to the surface at significantly higher velocities, of the or-der of 1 ms − . The measured velocities were 1 . ± . − for0 . µ m grains and 1 . ± . − for 0 . µ m grains. They usedslightly di ff erent definition of the critical velocity. The criticalvelocity was defined as the velocity at which the probabilities ofsticking and bouncing are equal 50%. They measured the stick-ing velocity for di ff erent materials, shapes and sizes of particles.The resulting critical velocity was shown to depend on size ofthe monomer as a power law with index of about 0 .
53. Althoughthe theoretically derived slope R − / of the power law is di ff er-ent from the empirical data, they are still consistent within errorbars. Moreover the power law was fitted to only two data points.We will assume that the theoretical dependence on the ra-dius is correct and consistent with the experiment. The discrep-ancy between the absolute values of the critical velocity howeverhave to be revised, doing so is absolutely essential for meaning-ful results.In order to better understand the mechanisms involved in thesticking of the monomers we refer to experiments by Poppe et al.(2000) and Dahneke (1975). They investigated the bouncing ofmicron sized spheres o ff a flat silica surface. In the former exper-iment the impacting particle was made of silica, while the latterone used soft polystyrene grains.Dahneke (1975) determined experimentally the critical ve-locity to be about 1 ms − . The ratio of the rebound velocity tothe incident speed (coe ffi cient of restitution) decreases in thiscase with decreasing impact velocity. At very large velocities ofabout 20 ms − Dahneke (1975) noticed a decrease of the resti-tution coe ffi cient with increasing impact speed. He related thatbehavior to plastic deformation.In the experiment by Poppe et al. (2000) the first positivevalue of the coe ffi cient of restitution indicates a critical veloc-ity of the order of ∼ − .In order to scale the model to the experimental data we in-troduce a mechanism that dissipates part of the initial energy atthe moment of first contact. The plastic deformation of small as-perities on the surface of the grain is an easy way to dissipate theenergy and increase the critical velocity. Below we will estimatethe central pressure in order to compare that to the yield strengthof 104 MPa (Callister 2000) for silica . Chokshi et al. (1993) onthe other hand argue that the yield strength for very small bodies Since the yield strength of brittle material like silica is not defined,flexular strength is given. D. Paszun & C. Dominik: Properties of porous aggregates is of the order of 0 . p ( r , a ) = F c π a aa (1 − (cid:18) ra (cid:19) ) / − F c π a aa ! − / (1 − (cid:18) ra (cid:19) ) − / . (4) a in this equation is the equilibrium contact radius given as a = πγ R E ∗ ! / . (5)Now we need to estimate the maximum contact radius that isreached during the collision. The radius of the contact area in-creases with increasing normal load. The force applied to thesphere during the collision is F = d ( mv ) / dt , where d ( mv ) is themomentum of the impacting particle and dt ≈ − s is the col-lision time. With this force applied, the contact radius is a = R E ∗ ( F + πγ R + q (6 πγ R ) + πγ RF ) ! / , (6)as derived by Johnson et al. (1971). Thus in the case of the ve-locity v = / s, the central pressure is of the order of 1 GPa.This pressure however is exactly in between the two values ofthe yield strength specified above, making the plastic deforma-tion possible. The lower limit of the yield strength (104 MPa) isreached at radii r / a < . E asp = YV asp , (7)where Y is the yield strength of the deforming material and V asp is the volume of the asperities that are flattened during the colli-sion. To calculate the dissipated energy we need to get the vol-ume of asperities. We follow Tsai et al. (1990) again and calcu-late the maximum contact area given by a max = R E ∗ ! [ F eq + πγ R + q (6 πγ R ) + πγ RF eq ] / (8)(Chokshi et al. 1993) and the equivalent impact force given byTsai et al. (1990) as F eq = . / E / R / E ∗ / . (9)Now we can get the volume of the deformed material by multi-plying the contact area by the volume of a single asperity, andthe number of asperities per surface area V asp = π r N asp ∗ π a , (10)where r asp is the radius of a single asperity. Here we assume thebumps to be hemispheres distributed homogeneously over thesurface of a particle or a target. In our case we use a parameterthat describes the total volume plastically deformed per unit ofthe surface area. This way we have only one parameter that de-fined how e ffi ciently the energy is dissipated. Note that in realitythe pressure in the contact area is compressive in the center andtensile on the edge. Thus our parameter may be slightly higherthan what we present. Fig. 1.
Critical velocity as a function of a grain size. Squaresindicate the model without additional dissipation process, dia-monds show present model with the energy dissipation by plasticdeformation of surface roughness and triangles show experimen-tally determined values by Poppe et al. (2000) with error bars.Fig. 1 shows the critical velocities determined experimen-tally (Poppe et al. 2000) for two di ff erent grain sizes. Also thecritical velocity obtained using our model is plotted showing thatwe successfully fitted our model and we can reproduce the ex-perimental results.In our model the volume of the asperities deformed uponcollision is very small. When we assume surface roughness tobe hemispheres with radii r asp = ff erence is that the energy regime has beenshifted to the level measured in experiments.In order to avoid confusion, we must stress here again thatwe do not claim that the plastic dissipation takes place duringthe collisions of silica spheres. We simply implement additionalscaling parameter in order to fit our model to the experimen-tal data. The discrepancy between the empirical data and theoryshould be further investigated. We also think that plastic defor-mation might need somewhat more attention, because of the veryhigh pressure that is present in the center of the contact area andrelatively low strength of the material considered in this work. In parallel to the energy dissipation due to contact breaking, en-ergy can be dissipated via other channels in the vertical degree offreedom (along the line connecting centers of two grains in con-tact). Two grains held together by the surface force vibrate rela-tive to each other (Dominik & Tielens 1995). Ideally the oscilla-tion is frictionless and no dissipation occurs. In reality, however,the vibration causes an oscillation in the size of the contact area.When decreasing the contact size, part of the energy is dissipatedby breaking the connections at the edges of the contact area. Thisultimately leads to a cooling of the aggregate and damps the vi-brations. If this mechanism is not taken into account, successiveslow collisions may heat up the aggregate and ultimately lead to“evaporation” of monomers from the aggregate surface. . Paszun & C. Dominik: Properties of porous aggregates 5
In fact we observed this phenomenon in our simulations oflinear chains. Successive collisions of grains with an aggregatecaused an increase of the amplitude of the oscillation and eventu-ally lead to breaking of several contacts. All individual collisionvelocities were far below the sticking velocity, proving that thismechanism may produce artificial results.To resolve this problem we introduced a weak dampingforce, which was intended to slowly dissipate the vibrational en-ergy. The force is acting in the vertical direction in respect to thecontact area. The damping force is expressed as F damp = const v z , (11)where v z is the vertical component of the relative velocity, andconst is arbitrarily chosen to damp exponentially 95 percent ofthe vibration energy within ∼
100 oscillation periods. By tuningthe constant to this value we made sure that the damping forceinfluence is not significant within short timescales of a few firstvibration periods, where the sticking or bouncing event is deter-mined. In the first period the dissipation takes place and if theparticle looses not enough energy it will bounce. If the dampingforce is too strong, it may remove enough energy to allow stick-ing. We made sure, that the main energy dissipation process thatsets the critical velocity is the plastic deformation mechanismand the sticking velocity is not a ff ected by the presence of thedamping force.To show the energy leakage due to the damping force wecalculated the total energy in a vibrating pair of monomers. Theparticles were displaced from the equilibrium position and thekinetic and potential energies were calculated. In the case of africtionless oscillation the potential energy is E p = Z δ F ( δ ′ )d δ ′ , (12)where δ is a displacement such that the distance between centersof the two grains is r + r − δ . The vertical force F ( a ) = F c aa ! − aa ! / ! (13)depends on the contact size a and equilibrium contact radius a .The integral (12) may be then changed into E p = Z a F ( a ′ )d a ′ . (14)To get the energy we solve eq. 14 by changing the variable from δ to a using the following relation δ = a R aa ! − aa ! / ! . (15)The potential energy is then given as E p = F c a R " aa ! − aa ! / + aa ! . (16)The kinetic energy is simpler. We just add kinetic energies of allmonomers together E k = X i m i v i , (17)where m i and v i are mass and absolute velocity of the i -th grain.The total energy E tot in the case of a frictionless oscillationis plotted in fig.2 a. To better see variations in the kinetic energy we shifted it to the level of the potential energy. The shift isequal to the potential energy of the system in the equilibrium E p = R δ F ( δ ′ )d δ ′ . The total energy in this case is conserved andthe amplitude is constant. When we enable the damping forcethe energy starts to leak. The potential and kinetic energies maybe calculated using the same formulas. The energy dissipation ispresented in fig. 2 b.The energy loss due to the damping force is very weak andremoves about 95% of the total energy within approximately 100vibration periods. Thus other processes like rolling, and break-ing dominate the energy dissipation. In particular the critical ve-locity for sticking is not influenced at all by this damping force,because the assumed plastic deformation of asperities dissipatesnearly the entire collision energy within the first vibration pe-riod.
3. Applications
The presented model is a very good tool that can be used tostudy the aggregation of dust particles. However, the modifi-cations that we introduced in the previous section need furtherverification. Although the critical velocity is well fitted to theexperiments, it is very useful to further test the model againstlaboratory experiments. When the model is successfully verified,it can be applied confidently in the study of aggregates dynam-ics. Compaction and fragmentation of dust aggregates are pro-cesses that require further investigation (Dullemond & Dominik2005). Our model is perfectly suited to this task. Dynamical pro-cesses, involving small sub-mm sized aggregates, can be studiedin depth using our tool. We can provide detailed understandingof micro-physics that governs compaction and fragmentation inthis scale. Phenomena involving aggregates of mm and largersizes require an entirely di ff erent approach. They must be treatedwith di ff erent methods i.e. Smoothed Particles Hydrodynamics(SPH). However, to do this, material properties, such as com-pressive strength and sound speed, are needed. Below we presentsimulations which can be directly compared to experiments andprovide great test of our model. We simulate compression of adust cake and determine the compressive strength of such cake.We also determine the sound speed in such a porous aggregates. Blum & Schr¨apler (2004) measured the compressive strength ofdust cakes in the laboratory. The sample was prepared by ran-dom ballistic deposition (RBD). Single monomers were shotfrom one direction at low velocity (hit–and–stick) and grow thedust agglomerate. The resulting “cake” was about 2 cm in di-ameter and similar in height. The finished cake was later placedbetween two flat surfaces. The load was applied onto the uppersurface and caused it to move towards the lower surface com-pressing the dust sample. In order to simulate compression of adust cake we need to prepare adequate setup.
In the experiment by Blum & Schr¨apler (2004) the applied pres-sure resulted in compression of the dust cake. Measurement ofthe cake volume resulted in determination of a filing factor φ .Our setup was organized in a similar manner to the experi-mental one. Our code, however, can handle only spherical par-ticles. Thus instead of two planes we used two very big “wall”grains, each with diameter of 2 · − cm. The sample was shaped D. Paszun & C. Dominik: Properties of porous aggregates
Fig. 2.
The total energy is redistributed between the kinetic energy (dashed line) and the potential energy (dotted line). At themaximum and minimum separation the potential energy is maximum, while in the equilibrium position the kinetic energy is thehighest. The kinetic energy was shifted down for better overview (see text). No damping case – a, and the energy leak form thesystem due to the damping force – b.as a cylinder and the distance between two compressing “wall”grains was adjusted to exactly fit the micro dust cake.Our dust cake was grown via particle cluster aggregationmethod (PCA). The target aggregate (initially a monomer) wasrandomly oriented before each collision with a grain approach-ing at a random impact parameter. Any successful hit resultedin perfect sticking without any restructuring. The new target wasthen randomly oriented, and a new grain was shot at it again. Theresulting aggregate was then shaped into a cylinder by removingall grains outside of the desired contour.The size of the cake was 13 . µ m in diameter and it was13 . µ m high. When filled by 291 monomers with radius of0 . µ m the filling factor of the cake was φ = . φ = . ± .
01. The monomer size in that experiment was alsoof similar size r = . µ m.The compression was done by moving one of the large grain-surfaces at constant velocity. While it was approaching, the sec-ond grain was fixed in its position and was unable to move evenupon extreme pressures. The sample was thus compressed byone surface against the other one. The initial setup is presentedin fig. 3a. Two large grains on both sides of the aggregate arethe back wall, fixed plane in the right and the approaching, com-pressing plane in the left. The aggregate is placed in betweenand the particles can escape sideways increasing this way thediameter of the cake.In order to simulate quasi static compression we fixed thevelocity of the compressing “wall” grain to 0 .
05 m / s. This isover an order of magnitude lower than a critical velocity andmuch lower than the sound speed in this medium (see section3.2). Thus the assumption that we are in a quasi static regime isreasonable. The dominant acceleration of a single monomer, inthis case, is due to surface forces.For the purpose of this compression simulation, we disablethe net force acting onto the compressing “wall” particle, butwe save the record of this force for each time step. Thus theapproaching surface cannot be stopped and moves with constantspeed. For each time step the net force is stored and later used todetermine a pressure. At each successive step the dust cake becomes slightly morecompressed. The degree of compression can be related to theforce that was needed to get the cake into this state. The small size of the dust cake used in this numerical experi-ment causes a low number of contacts of the compressed aggre-gate with the approaching surface. Consequently, the net forceapplied to the compressing surface was strongly variable. Everynew contact formed between the dust cake and the incoming sur-face resulted in a sudden decrease of the force. The new con-nection is initially stretched and thus the surface is attracted.Similarly, oscillations of monomers at the surface cause addi-tional variation. This makes it more di ffi cult to uniquely deter-mine the compression force. To overcome this problem, for teneach successive time steps we choose the one with the highestforce, because ultimately this is the force required to compressthe cake.The pressure was calculated as P = FS , where F is the normalload and S is the cross-section area.Fig. 3 shows the initial setup and the results of compressionwith increasing pressure. The lowest pressure cannot restructurethe aggregate. The height of the dust cake is almost unchanged.A higher pressure of 2 kPa compresses the cake. The monomersin the cake’s surface are pushed down into the cake. At a pres-sure of 1 · Pa the cake is compressed significantly, causinga horizontal flow of the particles and thus an increase of thecake diameter. The evolution of the dust cake cross-section isshown in fig. 4. The relative increase is a ff ected by the size ofthe dust cake. What may be a negligible boundary e ff ect in alarge macroscopic aggregate, here it has a large impact onto theentire dust cake. When the diameter of the cake increases just bya diameter of a single monomer, the area increases by 40 %. Inthe experiment by Blum & Schr¨apler (2004) a cake of about 2cm diameter was compressed and the final increase of the cross-section was measured to be larger by a factor of 1 . . Paszun & C. Dominik: Properties of porous aggregates 7 Fig. 3.
The setup of the experiment. The dust cake in the center is compressed with di ff erent pressure. Initial arrangement (a), resultsof compression at 2 · Pa (b),2 · Pa (c), 5 · Pa (d), 1 · Pa (e).
Fig. 4.
The ratio of the cross-section of the dust cake to the initialcross-section as a function of the applied compressive pressure.The data presented here is from a single experiment. Each pointis plotted with a 40% uncertainty in determination of the cross-section (see text).mination of the dust cake cross-section. The initially cylindricalshape changed upon compression into an irregular profile.To determine the volume filing factor, we used only the innerpart of the cake. A cylindrical volume enclosing initially the dustaggregate, and used to determine the filling factor, was decreasedin height only. In this way we reduced the uncertainty that arosefrom the boundary e ff ects. The initially porous aggregate is de-formed and can expand. Particles are pushed into the cake, fill-ing voids. Thus the filling factor must increase as an e ff ect ofcompression. Fig. 5 shows our results together with the data ob-tained in the laboratory experiments (Blum & Schr¨apler 2004).Low pressures are unable to a ff ect the aggregate. However, theboundary e ff ects also cause problems. Initially the “wall” parti-cle connects to the cake by only one monomer. This causes verystrong variations in pressure. We show the compression curvestarting in the point when the “wall” particle has made 5 con-tacts with the cake. At this point, however, a few grains are al-ready pressed into the dust cake, causing an increase of the vol-ume filling factor. Thus our compression curve in fig. 5 is shiftedupwards. The filling factor is also overestimated by using bigspherical grains as the compressing surfaces. The height of thedust cake we used is calculated as D − R compress , where D is dis-tance between centers of the two big grains and R compress is theirradius. Thus the volume occupied by the dust cake is actually Fig. 5.
The volume filling factor vs normal pressure. Thesolid line indicates the results of the laboratory experiment byBlum & Schr¨apler (2004). The dashed area is the error to thatdata. Squares show results of our single simulations. Error barsare due to di ffi culty in determination of the volume of the finalaggregate.slightly larger. This e ff ect is initially less important, as it con-tributes only about 2 % error. However, it gets more relevant atlarger pressure, where the cross-section of the dust cake is largerand the volume smaller.Error estimation was done using the central parts of the dustcake, where boundary e ff ects are smaller. The filling factor wasdetermined by calculating the volume of monomers enclosed ina cylindrical volume. The error was then estimated to be the dif-ference between the filling factors determined for two di ff erentcylinders. The first one had a radius of 6 µ m, while the secondone was one monomer radius smaller.The volume filing factor is initially constant, until the pres-sure reaches the value of about 5 · Pa. The filling factorincreases until it reaches the maximum compression of about φ = .
35 and remains constant. The most interesting result is theresemblance of our findings to the lab experiments. The onsetof compression fits the experimental data of Blum & Schr¨apler(2004). Our compression curve follows the laboratory datatightly ending up at a similar value of the filling factor. The smalldi ff erences between two curves are most likely due to small sizeof our simulated dust cake. D. Paszun & C. Dominik: Properties of porous aggregates
One of the very important properties in the porous material isthe sound speed. It must be lower than the bulk sound speed be-cause the mass is being moved by a force acting on a very smallcontact area. The collision of two grains with supersonic veloc-ities can result in complete disruption of those bodies. We willfirst derive an analytical formula for the sound speed. For this weapply the JKR theory and assume that the signal transported isa very small perturbation. The main assumption is that the twomonomers in contact behave like a perfect spring which is thecase for small amplitudes. However, the contact forces are asym-metric in respect to the equilibrium position. Thus compression of two particles results into di ff erent forces than stretching it tothe same displacement. We can therefore expect that large am-plitudes lead to modified sound speeds. In the following sectionswe present an analytical approach and later we show our simula-tions for both the simplified case of a linear chain of monomersand the general case of a non fractal aggregate. Every two monomers that are in contact are held together bysurface forces acting in the contact area (Johnson et al. 1971).Mutual attraction inevitably leads to a vibrational spring-likemotion, which in linear approximation can be written as F = − k ( δ − δ ) , (18)where δ , δ , k and F are displacements, equilibrium displace-ment, spring constant and restoring force, respectively. Note thatthe displacement δ is defined in a way that it increases, when thedistance between two monomers is decreasing. In our case wewant to determine the spring constant, which is necessary to findthe sound speed. Johnson et al. (1971) showed that the surfaceforce in the contact area is related to the contact area radius a by F = F c aa ! − aa ! / ! , (19)We also know the relation between the displacement and contactradius δ = / δ c (2 a − a − / a − / a / ) , (20)where δ c = / a / R is a critical displacement for breaking thecontact. With these equations we can determine the spring con-stant and later the sound speed. First we di ff erentiate eq. 19 and20 at a = a to get dFda (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) a = a = F c a (21)andd δ d a (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) a = a = δ a , (22)where δ = a R is an equilibrium displacement. Now we can writed F d δ = d F d a d a d δ = d F d a (cid:18) d δ d a (cid:19) − (23)and substituting eq. 21 and 22 we getd F =
65 (9 πγ R E ⋆ ) / d δ. (24) Table 1.
Material properties used in this work for silica. E ∗ [dyn / cm ] γ [erg / cm ] ρ [g / cm ]2 . · . . Since the restoring force is exactly opposite, we may add theminus sign here and see that the spring constant k is k =
65 (9 πγ R E ⋆ ) / . (25)With the spring constant of the oscillating system of twospheres in contact we can proceed to the calculation of the soundvelocity. For a spring the velocity of sound is given by c s = s kL ρ l , (26)where L is the total length of the spring (or set of the springs) and ρ l is a mass of the spring per unit length. We can then determinethe velocity in a linear chain of n grains in contact, each r inradius. The length of such a string of grains is L = n (2 r − δ ) + δ − r . (27)Since the spring constant for a set of springs is simply k L = k / n ,the sound velocity can be written as c s = s k L L ρ l = s kL / nnm = Ln s k / π r ρ . (28)Substituting eqs. 25 and 27 we get c s = − n ! (2 r − δ ) s (9 πγ R E ⋆ ) / / π r ρ / . (29)Now we can try to see how the sound velocity in such systemdepends on the size of a single monomer. c s ∝ r vt r / r = r − / . (30)This means that to double the speed of sound we need to use 64times smaller grains. Therefore, the sound speed in this case is avery weak function of monomer size.The sound speed of an infinitely long chain of 0 . µ m silicagrains is c s =
513 m / s. If we apply our findings to a chain of50, 0 . c s =
503 m / s. We can compare it now to the previousresults obtained in research of granular medium composed ofmm sized and bigger grains. In all our simulations and analyticalcalculations we used material properties as specified in tab. 1.Hasco¨et et al. (1999) developed a model of macroscopicgrains to study the propagation of sound in a granular chain.The centimeter sized grains they use do not interact in this casevia attractive surface forces. The signal is transported due tothe Hertzian stress that arises as an e ff ect of overlap of grainsin contact. They also apply the spring theory but with an arbi-trarily chosen spring constant k . For values of k in the rangebetween k = N / m and k = N / m, and mass den-sity of 1 . × kg / m , the sound velocity they derive is in arange of 300 m / s to 3000 m / s. Similar results were obtainedby Mouraille et al. (2006). A sound speed c s ≈
200 m / s wasfound for closely packed grains. In this case monomer radiuswas 1 mm, density 2 × kg / m , and the spring constant k = N / m. . Paszun & C. Dominik: Properties of porous aggregates 9 In order to verify our findings we performed a numerical exper-iment. We prepared a linear chain of 50 monomers. The firstone in line was slightly displaced from it’s equilibrium posi-tion. When the simulation started the grain started to move inorder to reach it’s equilibrium and the second grain was dis-turbed. Such motion was propagating until it finally reached thelast grain. The total distance traveled by the density wave was99( r − δ ). For 50 silica monomers with radii r = . µ mthe travel time took t = .
16 10 − s, which results in the soundspeed of c s =
512 m / s. We can now compare the value with thetheoretically derived sound speed c s =
503 m / s. The di ff erencebetween the theoretical velocity and the one that was obtainednumerically is only about 3%. Moreover in the simulation thefirst grain is displaced by a finite distance meaning that the as-sumption of low displacements may not be entirely correct. Thedisplacement of about 0 . δ is relatively large. We performeda series of simulations with stronger perturbations and the dataobtained shows that the sound velocity increases as the pertur-bation strength increases. Since a linear chain of monomers is a special case, we willdiscuss now the more general agglomerates of irregular shape.The sound speed in such a system is a ff ected by several things.Firstly, the path length that a signal has to travel is not a straight Fig. 6.
Sketch of a sound wave propagation in a linear aggregate(a) and in a porous aggregate (b).line and thus with a given speed, longer distances result in lower e ff ective sound speed. For a small RBD aggregate (approxi-mately 300 monomers), the e ff ective distance was larger by onlya factor of 1 .
5. Thus structure of an aggregate has a limited im-pact on the sound velocity by increasing the path length. Thesecond factor that may have an e ff ect is the tangential force. Ina linear aggregate, the signal is passed forward due to the pres-ence of the vertical force, and the tangential force has no e ff ect.However, in irregular aggregates grains interact also via rollingand sliding. This might lower the contribution of the verticalforce and in this way change the sound speed. Fig. 6 shows asketch of how a sound wave propagates in di ff erent aggregates.Indeed linear aggregates involve only stronger, vertical forcesand thus the signal travels faster. Irregular particles also makeuse of the tangential, rolling force.In order to estimate the sound speed due to the rolling fric-tion, we use again the spring theory. Dominik & Tielens (1995)gives the recipe for the rolling friction to be F roll = πγξ, (31) where the ξ is a displacement from the equilibrium position. Thisforce can be then expressed as F roll = k ξ, (32)with k = πγ – a spring constant. We then apply this in eq. 28 toget c s = − n ! (2 r − δ ) s πγ / π r ρ . (33)The dependence of the sound speed on a grain size, c s ∝ r − / , isin this case much stronger than in the case of the vertical force( c s ∝ r − / ).When we apply eq. 33 to the 1 . µ m silica grains we getsound speed c s = . / s. This is a factor of 30 lower thanthe velocity derived for the linear chain. This suggests that thesound velocity in a porous aggregate might be significantly dif-ferent from the speed derived for a linear chain of grains.Sirono (2004) derived compressive and tensile strengths de-pendence on density, based on experimental data. That relationswere later used to develop an SPH model of large, mm sized andlarger aggregates collisions. The sound speed in an aggregatemade of 0 . µ m ice monomers was calculated to be c s = p E /ρ .He used values for the Young’s modulus and the density to be E = × Pa and ρ = . / cm , respectively. The resultingsound speed is then c s = . / s.When we apply our formula for the sound speed to 100aligned 0 . µ m ice grains, we get sound speed c s =
885 m / s.This is almost an order of magnitude higher. We have to keep inmind however that the speed in a linear chain may be consider-ably di ff erent to the one in a real porous aggregate. In the rollingcase, the theoretical sound speed is also high with c s =
250 m / s.Teiser (2007) performed a lab experiment, where he mea-sured the sound speed in a RBD aggregate with filling factor φ = .
15 and 1 . µ m silica monomers. He hit the dust cake frombelow and measured the response at the surface with a force sen-sor. The measured velocity was c s = ± / s, very much con-sistent with our results when the signal is assumed to be trans-ferred through the rolling degree of freedom.In order to numerically derive the sound speed and furthertest our model we performed a simulation of an RBD dust cake.The cake was being pushed from one side and we determined theresponse time of di ff erent particles in the aggregate. We com-bined the positions of the particles with their response time,which results in an average sound speed in the dust cake. Fig 7ashows the result of our simulations. For each monomer in thedust cake we plotted a corresponding sound speed. Initially largespread in the data shows that inside the cake very linear struc-tures are present. This leads to a few particles with very largesound speed. At larger distances, however, the range of soundspeeds is much lower and shows that a signal is transportedmainly via rolling degree of freedom. The average sound speedthat is determined at the far end of the dust cake is only a factorof about 1 . ff erent forces to the cake and later deter-mined mean sound speeds in the cake for both cases. When astronger force of 5 × − N was applied to the cake the soundspeed reached c s ≈
60 m / s, while with the lower force of 5 × − N the sound velocity was c s ≈
20 m / s. For one monomer the lim-iting sound velocity, as calculated for linear chain of monomerswas overcome. This, however, happened in the case of larger Fig. 7.
Sound velocity in m / s versus distance from the plane thathits the aggregate. Squares indicate applied force of 5 10 − Nand triangles 5 10 − N. The left panel (a) shows results of soundspeed determination in a RBD aggregate, while the right one(b) in a CCP aggregate. The inset image shows the CCP ag-gregate placed in between two planes. The lines indicate soundvelocities determined theoretically using vertical force in a lin-ear chain of monomers (dashed line) and using rolling friction(solid line). The Dashed area limited by dotted lines indicates thesound speed range determined experimentally by Teiser (2007)force, and thus the low displacement approximation was mostlikely violated leading to larger sound speeds.The simulations show that the sound propagation is domi-nated by the rolling friction. In order to verify this we run anothersimulation to determine the sound speed in the aggregate withmonomers arranged according to the cubic close packing (CCP)because this arrangement disables any rolling motion. The insetin fig. 7b shows the aggregate. We applied a force on one sideof the agglomerate and determined the response time of di ff er-ent particles. Data shown in fig. 7b shows that indeed closelypacked aggregates are characterized by much higher velocitiesthan the porous ones. The monomers in each layer received thesignal at a di ff erent time because the compressing planes weresimulated by big grains and thus central particles were hit firstand they forwarded the signal. Thus particles at the edges of theaggregate responded later than the ones in the center of eachlayer. The velocity c s ≈
600 m / s shows that indeed sound speedin a porous medium strongly depends on porosity because ofthe forces involved in the transport of the signal and the longerpath for the signal to travel in more porous aggregates. Usinga constant sound speed in SPH simulations is bound to lead tospectacularly wrong results.
4. Conclusions
We have modified the original code by Dominik & Tielens(1997) and Dominik & N¨ubold (2002). We have shown that anumerical N-particle method for studying the properties of ag-gregates can be calibrated to experimental results by includingthe flattening of small asperities on the surface of the grains, andby using critical displacements for rolling of grains consistentwith measured values. The new code reproduces the measuredcritical velocities for sticking very well. We would like to em-phasize that there is currently no proof for plastic deformation actually happening at the grain surfaces, but as a model it worksvery well.We went on to measure the compressive strength of the ag-gregates and compared the results with experiments. We get verysimilar results with a small o ff set in porosity which is mostlikely due to the relatively small aggregates used in our simu-lations. Finally, we computed the sound velocity in an adhesion-dominated porous material and showed that this leads to veryinteresting results. Three very di ff erent velocities play a role insuch a medium, and the di ff erent velocities are separate by fac-tors of the order of 10. The fastest speed is the bulk materialsound speed which only plays a role after a body has been moltenand re-hardened. The second speed, a factor of ∼
10 lower is thespeed at which a signal is transported in a longitudinal wavein a linear chain. This speed applies either in a perfectly linearchain, or in an aggregate that has been compacted su ffi ciently sothat rolling of grains is no longer possible. Finally, the slowestspeed is the one transmitted by rolling forces in a non-straightchain of grains. A small decrease stems from the longer paththe sound has to take in a porous aggregate. By far the largestfraction stems from the weak forces in the rolling degree of free-dom. Experiments show that this is indeed the dominant speedin porous aggregates. However, our results show that the soundspeed should be a very steep function of density once a signif-icant number of monomers has 3 or more contacts with theirneighbors. It is to be expected that SPH approaches to model theproperties of dust aggregates (Sirono 2004; Sch¨afer et al. 2007)will fail strongly if these e ff ects are not taken into account prop-erly.In summary, we conclude that we do now have a workingmodel of dust aggregates that can be applied in parameter stud-ies. Acknowledgements.
This work was supported by the Nederlandse Organisatievoor Wetenschapelijk Onderzoek, Grant 614.000.309. We thank J¨urgen Blumfor useful discussions and hospitality during several visits. We also acknowledgea financial support of Leids Kerkhoven-Bosscha Fonds.
References
Blum, J. & Schr¨apler, R. 2004, Physical Review Letters, 93, 115503Blum, J., Schr¨apler, R., Davidsson, B. J. R., & Trigo-Rodr´ıguez, J. M. 2006,ApJ, 652, 1768Blum, J. & Wurm, G. 2000, Icarus, 143, 138Blum, J., Wurm, G., Kempf, S., et al. 2000, Physical Review Letters, 85, 2426Boley, A. C., Hartquist, T. W., Durisen, R. H., & Michael, S. 2007, ApJ, 656,L89Boley, A. C., Mej´ıa, A. C., Durisen, R. H., et al. 2006, ApJ, 651, 517Boss, A. P. 2007, ApJ, 661, L73Brown, R. 1828, Philosphical Magazine, Vol. 4, p. 161, 4, 161Callister, Jr., W. D. 2000, Materials Science and Engineering: An Introduction,5th Edition (Materials Science and Engineering: An Introduction, 5th Edition,by William D. Callister, Jr., pp. 871. ISBN 0-471-32013-7. Wiley-VCH ,January 2006.)Chokshi, A., Tielens, A. G. G. M., & Hollenbach, D. 1993, ApJ, 407, 806Cuzzi, J. N., Dobrovolskis, A. R., & Champney, J. M. 1993, Icarus, 106, 102Dahneke, B. 1975, Journal of Colloid and Interface Science, 51, 58Derjaguin, B. V., Muller, V. M., & Toporov, Y. P. 1975, Journal of Colloid andInterface Science, 53, 314Dominik, C. & Dullemond, C. P. 2008, in pressDominik, C. & N¨ubold, H. 2002, Icarus, 157, 173Dominik, C. & Tielens, A. G. G. M. 1995, Philosphical Magazine A, 72, 783Dominik, C. & Tielens, A. G. G. M. 1996, Philosphical Magazine A, 73, 1279Dominik, C. & Tielens, A. G. G. M. 1997, ApJ, 480, 647Dullemond, C. P. & Dominik, C. 2005, A&A, 434, 971Goldreich, P. & Ward, W. R. 1973, ApJ, 183, 1051Hasco¨et, E., Herrmann, H. J., & Loreto, V. 1999, Phys. Rev. E, 59, 3202Heim, L., Blum, J., Preuss, M., & Butt, H. 1999, Physical Review Letters, 83,3328Johansen, A., Henning, T., & Klahr, H. 2006a, ApJ, 643, 1219 . Paszun & C. Dominik: Properties of porous aggregates 11. Paszun & C. Dominik: Properties of porous aggregates 11