Compression behavior of porous dust agglomerates
^{1} Institut für Astronomie and Astrophysik, Eberhard Karls Universität Tübingen, Auf der Morgenstelle 10, 72076 Tübingen, Germany
email: alexs@tat.physik.unituebingen.de
^{2} Physikalisches Institut, Eberhard Karls Universität Tübingen, Auf der Morgenstelle 14, 72076 Tübingen, Germany
Received: 20 January 2012
Accepted: 11 March 2012
Context. The early planetesimal growth proceeds through a sequence of sticking collisions of dust agglomerates. Very uncertain is still the relative velocity regime in which growth rather than destruction can take place. The outcome of a collision depends on the bulk properties of the porous dust agglomerates.
Aims. Continuum models of dust agglomerates require a set of material parameters that are often difficult to obtain from laboratory experiments. Here, we aim at determining those parameters from ab initio molecular dynamics simulations. Our goal is to improve on the existing model that describe the interaction of individual monomers.
Methods. We use a molecular dynamics approach featuring a detailed microphysical model of the interaction of spherical grains. The model includes normal forces, rolling, twisting and sliding between the dust grains. We present a new treatment of wallparticle interaction that allows us to perform customized simulations that directly correspond to laboratory experiments.
Results. We find that the existing interaction model by Dominik & Tielens leads to a too soft compressive strength behavior for uni and omnidirectional compression. Upon making the rolling and sliding coefficients stiffer we find excellent agreement in both cases. Additionally, we find that the compressive strength curve depends on the velocity with which the sample is compressed.
Conclusions. The modified interaction strengths between two individual dust grains will lead to a different behavior of the whole dust agglomerate. This will influences the sticking probabilities and hence the growth of planetesimals. The new parameter set might possibly lead to an enhanced sticking as more energy can be stored in the system before breakup.
Key words: planets and satellites: formation / methods: numerical / protoplanetary disks
© ESO, 2012
1. Introduction
Unraveling the question of planetesimal formation is a crucial issue for the core accretion model of planet formation proposed by Pollack et al. (1996). To this day it remains unclear how dust and ice particles can grow several orders in magnitude from micron to kilometer sized objects. In the beginning micron sized dust grains (further referred to also as monomers or particles) may grow by low velocity, hitandstick collisions resulting in highly porous fractal aggregates (Kempf et al. 1999; Blum et al. 2000). As the aggregates grow larger their motion increasingly decouples from the surrounding gas leading to higher collision velocities (Weidenschilling 1977). With increasing impact energy colliding aggregates get compacted (Blum & Wurm 2000; Suyama et al. 2008; Wada et al. 2008; Paszun & Dominik 2009). At some point relative velocities become large enough that fragmentation is supposed to set in (e.g. Blum & Wurm 2008), which may limit the collisional growth of aggregates. Apart from fragmentation, radial drift and bouncing (e.g. Langkowski et al. 2008; Weidling et al. 2009; Güttler et al. 2010) may hamper the growth of planetesimals (Zsom et al. 2010). Despite these obstacles, experimental evidence indicates possible growth in high speed collisions (a few 10 m s^{1}) through sticking and reaccretion of ejecta (Wurm et al. 2005; Teiser & Wurm 2009) or the sweepup of smaller particles (Windmark et al. 2012). To understand the possible growth regimes requires material properties and simulations beyond the current data base.
In the context of planetesimal formation, collisions of micron sized aggregates have been studied theoretically using a molecular dynamics approach (e.g. Dominik & Tielens 1997; Wada et al. 2007). Here, the motions of all grains that make up the aggregate are followed individually, considering suitable interaction forces. However, billions of such grains would be necessary to model metersized objects. For computational reasons it is therefore necessary in such cases to use continuum models to simulate collisions between macroscopic aggregates. Smooth Particle Hydrodynamics (SPH) constitutes such an approach (Sirono 2004; Schäfer et al. 2007; Geretshauser et al. 2010). However, the outcome of such simulations depends strongly on the underlying porosity model (Güttler et al. 2009). To simulate collisions of large dust boulders with SPH, several parameters describing the behavior of the material such as the compressive, tensile, and shearstrength must be known in advance.
To this day, only few attempts have been made to measure the required material parameter of porous dust aggregates in laboratory experiments. Blum & Schräpler (2004) studied the case of unidirectional compression, where a sample aggregate is compacted between two plates. While the bottom plate was fixed a certain load was applied to the top plate. To calculate the pressure the crosssection of the compressed aggregate had to be determined. Their samples were produced by random ballistic deposition (RBD) and featured an initial filling factor of 0.15. The maximum filling factor they could reach was limited to approximately 0.33 as grains began to flow outwards as the pressure in the center increased.
Later on, Güttler et al. (2009) used a similar approach to study the omnidirectional compression. Their sample was put into a solid box and the load was applied by a movable piston on the top (see Fig. 3). As the grains could not escape the box they reached a significantly higher filling factor of roughly 0.58 upon compression. One advantage of this approach lies in the elimination of any uncertainty in the determination of the filling factor since the volume currently occupied by the sample is unambiguous.
The compressive strength of porous dust aggregates was determined numerically by Paszun & Dominik (2008) using a molecular dynamics approach based on an interaction model by Dominik & Tielens (1995, 1996). Yet, they only modeled the unidirectional compression and their simulations were limited to very few particles ( ≈ 300). Using similar material and model parameters the model has been applied to explore the growth regimes of dust agglomerates under protoplanetary conditions (Suyama et al. 2008; Wada et al. 2009, 2011; Ormel et al. 2009).
In this work we step back again and present a method to obtain the continuum material parameter of porous particle agglomerates, in particular the compressive strength curve, from ab initio simulations using a molecular dynamics approach. Our approach is based on the model by Dominik & Tielens (1995, 1996) using extensions by Wada et al. (2007). To test the applicability of the model we perform customized simulations of both, omni and unidirectional compression, with a much greater number of particles on the order of 10^{4}. As we will see, modifications of the model are required to properly reproduce the experimental results from Güttler et al. (2009). Having calibrated our model to the case of the slow compression as measured by Güttler et al. (2009) we will subsequently present new results on how the compressive behavior changes with increasing compression speed.
First, we briefly summarize the underlying physical model and present our extensions in Sect. 2. Our numerical setup is explained in Sect. 3. In Sect. 5, we first describe the calibration of our model. Afterwards, we present our results of studying the dynamic compression and provide simple analytical approximations describing the dependence of the compressive strength on the compression speed.
2. Interaction model
The agglomerates used in our simulations are composed of spherical monomers of equal size. To describe the interaction between individual monomers we adopt the physical model that has been presented by Dominik & Tielens (1997). In this model elastic grains may establish adhesive contacts when touching each other. Upon deformation of these contacts kinetic energy is dissipated. The forces and torques acting upon the monomers can be derived from corresponding potentials (Wada et al. 2007).
2.1. Particleparticle interaction
In Dominik & Tielens (1997) the interaction between two spherical particles is divided into four types (see Fig. 1):

1.
compression/adhesion;

2.
rolling;

3.
sliding;

4.
twisting.
In accordance with Wada et al. (2007) the following notation will be used in this work: r_{i} denotes the radius of a monomer i, γ_{i} the surface energy per area, E_{i} the Young’s modulus, ν_{i} the Poisson’s ratio, and G_{i} the shear modulus. Furthermore, the reduced radius R is given by
and Note that for simplification all monomers in our simulations feature the same properties.
Fig. 1 The four types of particle interaction: compression/adhesion a), Rolling b), Sliding c), and Twisting d). The left panel depicts particleparticle interactions while on the right the corresponding particlewall version is displayed. 

Open with DEXTER 
Material parameters.
2.1.1. Material parameters
The basis of our simulations are the material parameters of silicate as summarized in Table 1. With the exception of the surface energy γ these values comply with Paszun & Dominik (2008), where they used γ = 25 mJ/m^{2}. These parameters are also in reasonable agreement with experimental data as quoted by Blum & Schräpler (2004); Güttler et al. (2009). In this work we use a slightly lower value of γ = 20 mJ/m^{2} which agrees with recent measurements of Gundlach et al. (2011).
2.1.2. Compression/adhesion
The particle interaction in the normal direction has been developed by Johnson et al. (1971) (often referred to as JKR theory). They extend the Hertzian theory by taking adhesion due to surface forces into account.
To very briefly summarize their model let us first note that the compression length δ of two monomers, which are in contact with each other and located at x_{1} and x_{2}, is defined by (1)where ∥ u ∥ denotes the norm of the vector u. The radius a of the corresponding contact area can be obtained by (2)where δ_{0} and a_{0} denote the equilibrium compression length and contact radius, respectively, where the repulsive and attractive normal forces are equal. It applies The force acting upon the monomers is (3)where F_{c} = 3πγR is the force needed to break the contact.
A new contact is established if two freely moving particles touch each other, which means δ = 0. However, once a contact has formed it will not break until the compression length exceeds a certain threshold δ_{c} = (9/16)^{1/3}δ_{0}. Thus, contacts may be stretched out a little bit before they finally break.
2.1.3. Sticking velocity
On contact formation or breaking there is a jump in the potential corresponding to the JKR force (see Wada et al. 2007), which reflects the dissipation of kinetic energy upon contact formation or breaking. Chokshi et al. (1993) proposed that this energy dissipation is caused by the excitation of elastic waves. From this energy dissipation one can calculate the maximum velocity v_{stick} at which particles stick on headon impacts in JKRtheory (4)For micronsized SiO_{2} grains, we obtain velocities on the order of 0.1 m s^{1}. However, laboratory experiments by Poppe et al. (2000) yield a significantly higher sticking velocity of roughly 1.2 m s^{1} for similar sized particles. To overcome this discrepancy, Paszun & Dominik (2008) introduced another damping mechanism that dissipated additional kinetic energy upon the contact creation and thus increased the sticking velocity to match the laboratory experiments.
Since the sticking velocity has to be increased by approximately one order of magnitude a significant amount of kinetic energy has to be dissipated aside from JKRtheory. Dissipating the kinetic impact energy all at once leads to a significant change of the relative velocity between the colliding particles. As they may be in contact with other particles as well, significantly modifying the velocity of one particle during one integration step may introduce numerical hazards. For low collision velocities (the compression velocity used by Paszun & Dominik 2008, was 0.05 m s^{1}) the additional damping is low and therefore this problem does not arise. However, since in this work we also want to study the dynamic compression behavior at high velocities we choose a different approach to adjust the sticking velocity as explained in the next section.
2.1.4. Normal oscillations
The normal force may be both attractive or repulsive depending on the current compression length δ. When two particles come too close to each other they are repelled, whereas they get attracted to each other due to adhesive surface forces while the contact is stretched out. This leads to oscillations of two monomers in the normal direction of a contact (from now on referred to as normal oscillations). In reality, one would expect that these oscillations are eventually damped away. For instance Brilliantov et al. (2007) proposed a viscoelastic damping mechanism.
From a numerical point of view these normal oscillations constitute a major nuisance. Not only may they artificially heat up aggregates (Paszun & Dominik 2008) but most importantly they need to be properly resolved in time. For micronsized SiO_{2} grains the typical timescale of these oscillations is on the order of 10 ns. Therefore our integration timestep is limited to ≈0.1−0.3 ns.
Paszun & Dominik (2008) tackle the problem of artificial heating by introducing an additional, weak damping mechanism which has hardly any influence on the sticking velocity but eventually damps away normal oscillations. In this work we follow a similar approach. As before we consider two monomers located at x_{1} and x_{2} that are in contact with each other. The contact normal n_{c} is (5)The relative velocity v_{rel} in normal direction of the contact is then given by (6)where v_{1}, v_{2} denote the velocity of particle 1 and 2, respectively. The viscous damping force F_{damp} is (7)where κ is a damping constant determining the strength of the damping. We use κ = 1 × 10^{6} g/s.
2.1.5. Rolling, sliding, and twisting
The forces and torques resulting from the tangential motion of the contact area have been formulated by Dominik & Tielens (1995, 1996). So called “contact pointers” (Dominik & Nübold 2002) provide a convenient way to track the evolution of the contact area. They are unit vectors that initially point from the center of the particle to the center of the contact area. Due to the rotation of the particles their orientation changes over time (Dominik & Nübold 2002). Contact pointers are used to define the rolling, sliding, and twisting displacement which quantify how much rolling, sliding, or twisting of a contact has occurred.
Using these displacements, Wada et al. (2007) derived the forces and torques which agree with Dominik & Tielens (1995, 1996), from corresponding potentials. All three types of interaction remain elastically as long as the displacements do not exceed certain thresholds (from now on referred to as critical displacements).
Let n_{1} and n_{2} be the corresponding contact pointers describing the contact. The rolling displacement ξ, sliding displacement ζ, and twisting displacement φ are then defined by: where t_{0} denotes the time when the contact has been established.
For the rolling motion, the forces F_{r} and torques M_{r} acting upon particle 1 due to being in contact with particle 2 are For the sliding interaction, the forces and toques are given by And for the twisting interaction it applies The constants k_{r}, k_{s}, and k_{t} are given by
2.1.6. Inelastic interaction
As already mentioned before, inelastic motion sets in when the displacements exceed a critical displacement ξ_{crit}, ζ_{crit}, or φ_{crit}. Physically this corresponds to energy being dissipated upon relocation of the contact area. The critical sliding and twisting displacements have been derived theoretically (Dominik & Tielens 1996) as The value of the critical rolling displacement ξ_{crit} is still debated. At first, Dominik & Tielens (1997) assumed ξ_{crit} close to inneratomic distances and chose ξ_{crit} = 0.2 nm. However, subsequent laboratory experiments by Heim et al. (1999) indicate a much higher value of ξ_{crit} = 3.2 nm for spherical SiO_{2} grains of 1.9 μm in diameter. In this work we follow Paszun & Dominik (2008) and set ξ_{crit} = 2 nm for 1.2 μm sized grains.
When a displacement exceeds its critical value it will be restored to the elastic limit. If for instance the rolling displacement exceeds its critical value ∥ξ∥ > ξ_{crit}, the contact pointers will be corrected , and thus ξ → ξ^{c} such that ∥ ξ^{c} ∥ = ξ_{crit}. This modification of the contact pointers reflects a change of the corresponding potential energy. Therefore, we can keep track of the amount of dissipated energy. For further details on how the inelastic motion is applied we refer to Wada et al. (2007).
Fig. 2 Setup of the numerical simulations to obtain material parameter for porous dust agglomerates. Left: omnidirectional compression, while the top wall is moving downwards at constant speed, the sample is enclosed in a box with fixed walls on all sides. Right: unidirectional compression, as the sample is getting compressed between two walls, particles can leave the initial volume to the sides. Particles colored red are actually in contact with the walls. 

Open with DEXTER 
2.2. ParticleWall interaction
To study the compression of a sample aggregate it is necessary to constrain the motion of the monomers using suitable boundary conditions. In the corresponding experiments the sample has been put in a solid box with a movable piston on top (see Güttler et al. 2009, Fig. 2). In our simulations we model the experimental setup by putting the sample aggregate into a box of fixed walls. During the simulation the top wall is moving downwards with constant speed in order to compress the sample (see Fig. 2).
We assume that the monomers may interact with a wall in a similar way as they interact with other grains. In accordance with the particleparticle interaction we derive the corresponding forces and torques following the approach presented by Wada et al. (2007). For this purpose we assume that the wall can be described as a very huge particle in the limit r_{wall} → ∞.
2.2.1. Compression
The force of the particlewall interaction in normal direction is very similar to the case of the particleparticle interaction. The compression length δ is given by
where r is the radius of the monomer and d denotes the distance between the surface of the wall and the center of the grain. Given a point p located on the surface of the wall and the surface normal n_{w}, we can easily obtain d by
where x denotes the position of the particle. The force in normal direction can then be calculated using Eq. (3). Note that here the reduced radius R is different to the case of particleparticle interaction. In the limit of r_{wall} → ∞ the reduced radius R equals the radius r of a monomer:
2.2.2. Rolling
Keeping in mind that R = r, the torque M_{r} acting on the particle caused by rolling along the surface of the wall is given by (20)where k_{r,wall} is equivalent to the rolling constant k_{r} given in Eq. (17) taking the different reduced radius of the particlewall interaction into account. n_{w} denotes the surface normal of the wall. Note that it is important on which side of the wall the particle is located with respect to the direction of n_{w}. In the simulation we must either ensure that the particles remain on the positively oriented side all the time or check on which side of the wall the particle is located and correct the sign of Eq. (20) if necessary.
To model the motion of a particle which is rolling inelastically over a wall we use a similar approach as for the inelastic particleparticle interaction. Taking r_{2} → ∞ into account, the correction of the contact pointers is then given by where α is a correction factor derived in detail in Wada et al. (2007) and Δξ is given by (23)The contact pointer n_{2} of the “wall particle” is equivalent to the normal vector n_{w} of the wall and is not modified during the inelastic rolling motion.
2.2.3. Sliding
To describe the sliding motion of a particle on a wall it is not suitable to start with Eq. (9) and assume that n_{c} → n_{w} and n_{2} → n_{w}. Therefore we choose a different approach that takes into account how far the contact has slided from the location where it has initially formed.
Let p denote the position, where the surface of the particle touched the wall when establishing the contact. For any later time, the center of the contact area is given by x + rn_{1}, where x is the current position of the particle. We define
The sliding displacement is then given by (24)To model the inelastic wall sliding we modify the initial center of the contact area p. If ∥ ζ ∥ > ζ_{crit}, we apply the correction (25)
2.2.4. Twisting
The torque caused by the twisting motion is calculated in the same way as if two particles are in contact. Starting with the twisting displacement φ given in Eq. (10) we obtain (26)under the assumption that the wall does not rotate. The torque is then given by (27)where k_{t,wall} = k_{t}. Here we assume that particles may twist elastically around the same angle for both the particleparticle and particlewall interaction. According to test simulations where we measured the relative importance of the different wall interaction dissipation channels, inelastic wall twisting is of only minor importance.
3. Setup
In order to calibrate the model and to compare simulations in detail with experimental results we focus here on two different numerical experiments that follow closely the experimental setup. Specifically, we will deal with agglomerates enclosed in a box and aggregates confined between two plates as depicted in Fig. 2.
3.1. Sample generation
In accordance with the mentioned laboratory experiments by Blum & Schräpler (2004); Güttler et al. (2009), our samples are produced by random ballistic deposition (RBD), which means that single grains are successively poured down on the existing aggregate impacting from the same direction. In order to prevent any restructuring upon impact the impact velocity of a monomer hitting the sample is kept very low. The resulting samples feature a filling factor φ between 0.12 − 0.15. Filling factors of 0.12 − 0.14 result from the fluffier surface and are therefore only observed for small samples. As the size of the samples increases this surface effect becomes negligible and the filling factor converges to 0.15 which agrees well with the dust cakes used in the corresponding laboratory experiments and with numerical studies by Watson et al. (1997).
3.2. Measurements
3.2.1. Pressure
In the numerical (and laboratory) experiments a boxshaped sample is enclosed between walls that constrain the motion of the grains (see Fig. 2). Then, the top wall is being moved downwards at a constant speed until the filling factor exceeds a certain threshold φ_{crit}. Typically, φ_{crit} is set to 0.7. As the top wall is moving downwards the volume of the sample decreases and it is increasingly more compressed. The total force F_{w} acting upon the top wall is calculated by summarizing the forces F_{i} exerted on the wall by grains that are currently in contact with it, where only the component in normal direction n_{w} to the wall is taken into account
The pressure P is then given by
where A denotes the base size of the box.
If there are only a few particles in contact with the wall, F_{w} may change considerably from one integration step to the next due to the normal oscillations of the particlewall contact (see Sect. 2.1.4). Since these vibrations occur on a timescale of nanoseconds whereas the compression timescale is typically orders of magnitudes higher, it is reasonable to average over several integration steps (covering a few normal particle oscillations) to reduce the noise in the pressure determination. Typically, we averaged over 100 integration steps in this work.
3.2.2. Filling factor
The volume filling factor φ is defined as (28)where V_{mat} = 4/3πr^{3}N denotes the volume of N particles of radius r and V_{tot} is the volume that the aggregate currently occupies. Calculating the filling factor is trivial in the case of omnidirectional compression (see Fig. 2, left panel) as V_{tot} = Ah for a box with base size A and current height h.
However, in the unidirectional case (Fig. 2, right panel) there are no side walls containing the aggregate, and particles will leave the initial volume. They flow to the sides as the top wall is moving downwards. This complicates the determination of V_{tot} in the numerical as well as experimental setup. In the following we assume that the volume the aggregate is currently occupying is given by its projected cross section A_{proj} and the current height h of the aggregate, which equals the distance between the top and the bottom wall. A_{proj} is obtained by projecting the aggregate in the plane of the top/bottom wall.
3.3. Previous work
A first attempt to determine the compressive strength numerically has been presented by Paszun & Dominik (2008). With the exception of the damping of the normal interaction (see Sects. 2.1.3 and 2.1.4) they used the same particle interaction model that we use here. However, instead of flat walls they used two huge particles with radii much bigger than the dust grains to model the boundary conditions. The sample was put between the “wall particles” and the upper particle was being moved downwards at constant speed while the force acting on this “wall particle” was stored for later analysis. Since the motion of the grains has only been constrained in the vertical direction, grains could dodge to the sides as the pressure increased. Thus, an increase of the initial cross section of the aggregate was observed during the compression that led to a significant uncertainty in the determination of the volume the aggregate occupied at a certain point of time. Since the number of monomers was also limited to very small numbers ( ≈ 300) it is questionable if the results hold in the continuum limit.
Paszun & Dominik (2008) confined themselves to the case of unidirectional compression. To our knowledge, the case of omnidirectional compression has not been simulated so far for this material. The compressive strength was only determined for a compression speed of 0.05 m s^{1}. A possible dependency of the compression behavior on the speed of the compression has not been studied before.
4. Numerical method
To integrate the equations of motion we use a second order, symplectic LeapFrog scheme. The main reason behind this choice lies in the fact that the forces and torques have to be determined only once during each integration step. Likewise, the evolution of the contact pointers and the twisting displacement (see Eq. (10)) is calculated with second order accuracy.
The timestep is limited by the oscillations in direction of the contact normal caused by the normal force (see Sect. 2.1.4) as well as the oscillations in the tangential plane of the contact caused by the sliding force. As already mentioned the normal oscillation period is on the order of 10 ns.
The period of the oscillations in the tangential plane can be obtained in the following way: Wada et al. (2007) derive the sliding force from the corresponding sliding potential U_{s}(29)We can get an estimate of the oscillation period T by (30)where ω denotes the corresponding angular frequency of the oscillation and m is the mass of a monomer. For the material parameters given in Table 1 we obtain T = 12.8 ns. Applying the sliding modifier m_{s} = 2.5 (see Sect. 5.1) the tangential oscillation period decreases to T = 7.83 ns, which limits our timestep to 0.3 ns.
Fig. 3 The compressive strength (filling factor φ versus pressure p) as obtained from experiments. The dark shaded region with the solid line fit refer to the omnidirectional compression experiment performed by Güttler et al. (2009) whereas the light shaded region and the dashed line fit were obtained for unidirectional compression by Blum & Schräpler (2004). The small image on the top left depicts the experimental setup of the omnidirectional compression experiment. (Figure taken from Güttler et al. 2009.) 

Open with DEXTER 
5. Results
Calibration experiments using a continuum SPHcode indicate that the compressive strength of a porous dust aggregate depends on how fast the compaction takes place (Güttler et al. 2009; Geretshauser et al. 2010). So far, in laboratory experiments the compressive strength could only be determined for a slow quasistatic compression process, where the compressed aggregate had been given sufficient time for relaxation (in the following referred to as static compression) (Güttler et al. 2009).
The static compression provides us with the possibility to check how well our model is able to describe the compression behavior of porous dust aggregates. In the first step we will therefore use the case of the omnidirectional static compression to calibrate our molecular dynamics model. Afterwards we will increase the speed of the top wall. At a sufficiently high speed the relaxation of the aggregate will not be possible any longer (from now on referred to as dynamic compression).
5.1. Calibration of our model
To model the quasistatic compression the speed v_{wall} at which the top wall is moving downwards should be as low as possible. As the time required to reach a certain filling factor is inversely proportional to v_{wall}, the runtime of the simulation constitutes a lower limit of v_{wall}. Paszun & Dominik (2008) considered v_{wall} = 5 cm/s to be low enough to model static compression. However, we observe that the curves still change when using even lower velocities (see Fig. 6 below). To model the case of static compression we set v_{wall} = 1 cm/s. To ensure it is a reasonable choice we checked lower velocities down to v_{wall} = 0.2 m s^{1} but observed only a tiny deviation of the resulting curves.
Depending on the number of particles we use, the diameter of our sample aggregates is roughly 40 − 60 μm which is about × 10^{3} times smaller than the samples used in the laboratory experiments. While the sample is getting compacted the particles on the edges of the sample must overcome the sliding resistance of the side wall. Owing to the small diameter of the sample this has a significant influence on the resulting force on the top wall. To mitigate this effect we reduce the strength of the rolling, sliding, and twisting interaction between particles and the side walls by a factor of 1000.
The results of the corresponding laboratory experiment are shown in Fig. 3. The solid black curve is a good fit to the omnidirectional experimental data (see Fig. 4 and dark shaded area in Fig. 3) and is given by (31)Using φ_{1} = 0.15 and φ_{2} = 0.58 we obtain p_{m} = 16.667 kPa and Δ = 0.562. The black curve shown in Fig. 4 depicts this fit.
The results of the unidirectional experiments, also displayed in Fig. 3, have been fitted by Blum & Schräpler (2004) using a similar curve with φ_{1} = 0.15, φ_{2} = 0.33, p_{m} = 5.6 kPa, and Δ = 0.33.
Fig. 4 Compressive strength for omnidirectional compression of cubical aggregates. The result of the unaltered interaction model is compared to our improved model. The black line represents a fit to experimental data obtained from Güttler et al. (2009). 

Open with DEXTER 
5.2. Omnidirectional compression
5.2.1. Quasistatic case
In Fig. 4, we compare the experimental fitting curve to results from our simulations of boxshaped samples featuring approximately 11 000 particles and an edge length of about 40 μm. All curves have been obtained from averaging the results of six independent runs with samples having statistically equal bulk properties.
Obviously, simulations using the Dominik and Tielens model do not reproduce the experimental data very well. For low pressures the blue dashed curve in Fig. 4 lies significantly above the solid black one, i.e. applying the same pressure the sample has been compressed to a higher filling factor in the simulations than compared to the experiments. In order to make it harder to compress an aggregate we tried to increase the stiffness of the aggregates by modifying the particle interaction.
To our knowledge the equations describing the rolling and sliding interaction have not been experimentally tested yet, whereas the pulloff force of the normal interaction has been measured using atomic force microscopy (Heim et al. 1999). Thus, we vary the strength of the rolling and sliding interaction. For this purpose we simply multiply the constants k_{r} and k_{s} (see Eqs. (17) and (18)) with correction factors that we further refer to as rolling/sliding modifiers m_{r} and m_{s}. This constitutes a straightforward approach to increase the stiffness of monomer chains. In fact, we also modified the strength of the twisting interaction but found that it had very little to no impact on the compressive curve. Therefore we do not alter twisting in this work.
Fig. 5 Snapshots of the vertical profile of the aggregate’s filling factor averaged horizontally. Shown are results at 5 evolutionary times using 3 different speeds. The left (solid) curves represent the initial state at t = 0. As the top wall moves slowly downwards the aggregate is compacted almost homogeneously as is indicated by the nearly vertical curves. At the top and bottom of the box wall effects produce a slight deviation. 

Open with DEXTER 
Choosing m_{r} = 8 and m_{s} = 2.5, we obtain the reddotted curve in Fig. 4. All in all, our modified interaction model is able to reproduce experimental results much better than the original version. In particular, for pressures below 100 kPa we observe a very good agreement with experimental data. However, we observe a deviation for pressures above 300 kPa. In our simulations a pressure of more than 1 MPa is required to further compress aggregates when a filling factor of about Φ = 0.52 is reached.
In the quasistatic case the aggregate is given sufficient time to restructure and counteract the pressure exerted on the walls. Thus, we expect the filling factor increases homogeneously in the sample. In Fig. 5 the vertical profile of the filling factor is plotted for different stages of the compression process. To determine the filling factor profile, the sample is split vertically into equidistant intervals with the length of one particle diameter. Then, the average filling factor is calculated for each interval. Note that the fractal structure resulting in a filling factor of φ = 0.15 in the bulk part of RBDgenerated aggregates is not present at the bottom of the aggregate. Therefore the filling factor there exceeds the average value of φ = 0.15.
During the compression process several snapshots of the filling factor profile have been taken. As expected the filling factor increases almost homogeneously for slow compression speeds. Keep in mind that it requires higher pressure to compact particles close to a wall since neighbouring particles have to be pushed away. Therefore the filling factor of the particle layers close to the top or bottom wall is lower compared to the rest of the aggregate for highly compacted aggregates. This effect causes the crescent shaped curves observed for highly compacted aggregates in Fig. 5.
The corresponding compressive strength curve for these low velocities is shown in Fig. 6. For compression speeds of 5 cm/s and below the differences to the quasistatic case remain small. For higher velocities we observe an increasing deviation from the quasistatic curve. Thus, the results from Güttler et al. (2009) cannot be applied for higher velocities.
Fig. 6 Compressive strength for low velocity omnidirectional compression. 

Open with DEXTER 
5.2.2. Dynamic case
In a second step, we studied the influence of large compression speeds. If the compression speed exceeds around 1 m s^{1} the compression behavior changes considerably. The required simulation time is inversely proportional to the compression speed. Therefore we use a higher number of particles for compression experiments with wall speeds above 1 m s^{1}. The boxshaped samples are composed of about 40 000 particles and feature a base length and height of ≈60 μm. Compared to the quasistatic case the shape of the curves changes drastically, see Figs. 7 and 8. Instead of a smooth transition, three distinguished regimes emerge: the filling factor does not increase until a certain critical pressure is reached. Then, only a small additional pressure is required to compact the aggregate. When the aggregate is close to its final compaction the pressure again increases sharply.
Fig. 7 Dynamic omnidirectional compression with v_{wall} between 1 − 5 m s^{1}. 

Open with DEXTER 
Fig. 8 Dynamic omnidirectional compression with v_{wall} between 7−15 m s^{1}. 

Open with DEXTER 
This can be easily explained by looking at Fig. 9. When the compression speed exceeds a value of 1 m s^{1} the aggregate is compacted inhomogeneously. The compression occurs too fast to allow the propagation of the top pressure through the entire sample. We clearly see the emergence of a very dense layer right beneath the top wall. This compact layer propagates downwards at the speed of the wall similar to a snowplough clearing freshly fallen snow. While pushing the dense layer downwards the pressure remains constant. After it reaches the bottom of the sample the pressure required to compress the sample a little bit further increases drastically. The sharp spikes shown in Figs. 7 and 8 result from the highly compacted layer reaching the bottom of the sample. The density wave is reflected from the bottom causing heavy fluctuations of the pressure on the top and bottom wall.
By comparing the filling factor profile during the compression to the initial one we can determine at which speed v_{p} the compaction is propagating downwards trough the sample. For this purpose we measure the height where the initial and current filling factor profile deviate from each other and use the time that has passed since the start of the simulation to calculate the speed. Averaging over six different samples we obtain v_{p} = 6.98 ± 0.16 m s^{1}.
Fig. 9 Snapshots of the vertical profile of the filling factor. As the top wall moves downwards rapidly the compaction of the lower parts of the aggregate is lagging behind. The color indicates the compression speed whereas the line type indicates the position of the top wall. The solid lines show the filling factor profile of the initial uncompressed sample. 

Open with DEXTER 
To provide continuumsimulations with a simple recipe for the dynamic compressive strength we performed simulations using compression speeds up to 25 m s^{1}. A few examples are shown in Figs. 7 and 8. For every compression speed we determine a fit curve similar to Eq. (31), where p_{m} and Δ serve as fitting parameters. Thus, we obtain values of p_{m} and Δ for different compression speeds. In the last step we determine for each parameter an analytic approximation describing the dependency on the compression speed.
We observe different behavior between the the quasistatic case for low velocities and the dynamic case for higher velocities. In the beginning, p_{m} decreases which means that the sample can be compressed more easily. This can be explained by the fact, that the aggregate is given less time to restructure and counteract the external pressure exerted on it by the wall. However, this effect will be reversed when the compression speed exceeds a critical value of v_{crit} ≈ 0.9 m s^{1}. In the dynamic regime, it gets considerably harder to compress the sample with increasing velocity of the wall. Therefore it is helpful to distinguish between the two regimes.
In Fig. 10 the dependence of p_{m} on the compression speed v is shown for values of v ≤ 1 m s^{1}. Using the ansatz p_{m}(v) = av^{2} + bv + c we obtain (32)where the compression speed v is given in m s^{1}.
In Fig. 11 the dependence of p_{m} on the compression speed v is shown for values of v ≥ 1 m s^{1}. To find a simple analytic approximation we choose a power law of the form p_{m}(v) = av^{b} + c and we obtain (33)The fitted parabola describes the data points well. As the exponent of 1.93 is close to 2, we determine a second fit where the exponent was set to 2 and get (34)Similarly, for the parameter Δ (see Fig. 12) we obtain (35)
Fig. 10 Dependence of the fit parameter p_{m} of Eq. (31) on the compression speed v in the quasi static regime of v ≤ 1 m s^{1}. 

Open with DEXTER 
Fig. 11 Dependence of the fit parameter p_{m} of Eq. (31) on the compression speed v in the dynamic regime of v ≥ 1 m s^{1}. 

Open with DEXTER 
Fig. 12 Dependence of the fit parameter Δ of Eq. (31) on the compression speed. 

Open with DEXTER 
5.3. Unidirectional compression
Additionally, we simulated the unidirectional compression of cylindrical samples of different sizes using the non modified model. The results are shown in Fig. 13 where again each curve is obtained by averaging the results for six different samples of equal size. To compare our results to Paszun & Dominik (2008) the speed of the top wall was set to v_{wall} = 5 cm/s. Apparently there is a noticeable discrepancy between our simulations and the laboratory results. As in the case of omnidirectional compression, the pressure required to reach a certain filling factor is significantly lower in our simulations.
As we can see in Fig. 13 the deviation for pressures above 10 kPa becomes more apparent if we increase the size of the samples. As Paszun & Dominik (2008) compressed very small samples using only about 300 particles this may be the reason why their results showed better agreement with laboratory results for higher pressures. However, their compressive strength curve was also shifted in the same direction as in this work.
Afterwards we tested the modified model with the same m_{r} = 8 and m_{s} = 2.5 as found above. As it is shown in Fig. 14 (red dotted curve), the modified model agrees very well with the experimental results for pressures below 10^{4} Pa. However, we still end up with considerably higher filling factors. A possible explanation is given by the fact that the size of the samples used in the laboratory experiments was on the order of centimeters and thus about 2000 times larger than our samples. Blum & Schräpler (2004) measured that the projected cross section increased by a factor of 1.6 during the compression process. To reach similar pressures in our simulations the sample has to be compressed until its height is only about two times the diameter of a single monomer where the cross section increased roughly by a factor of 4. Due to the larger diameter of the laboratory samples it is harder for monomers in the center to flow in the outward direction.
Fig. 13 Compressive strength for unidirectional compression of cylindrical aggregates of different size using the original Dominik & Tielens model. The black line has been obtained from experiments by Blum & Schräpler (2004). 

Open with DEXTER 
Fig. 14 Unidirectional compression of cylindrical aggregates using the modified model with 18 000 particles. The black line has been obtained from experiments by Blum & Schräpler (2004). 

Open with DEXTER 
6. Conclusions
We have performed molecular dynamics simulations to study the compressive strength of dust agglomerates which plays an important role in determining the outcome of mutual collisions. Using a special setup for the simulations we were able to compare our results in detail to the outcome of of dedicated laboratory experiments.
Our simulations using the frequently applied interaction model by Dominik & Tielens (1997) indicate that real aggregates composed of micron sized silicate grains feature a greater stiffness. Since the primary bulk properties of material used for the individual monomers are known experimentally very well (see Table 1), we decided to vary the force constants (k_{r} and ks) for rolling and sliding. Indeed, the higher stiffness can be accommodated by an increase of m_{r} = 8 and m_{s} = 2.5 (for k_{r} and k_{s}, respectively) in comparison to the quoted values in Eqs. (17) and (18). After modifying the interaction model as described in Sect. 5.1, we have been able to reproduce the experimental results much better, and found very good agreement for both, unidirectional and omnidirectional compression. This work reveals the importance of the rolling and sliding interaction for the restructuring of aggregates. As these interactions currently lack experimental testing we feel it desirable to study in particular the rolling and sliding of micron sized grains in laboratory experiments.
We have also studied the influence of the wall speed on the compression behavior. If an aggregate is compressed slowly the filling factor increases homogeneously and the pressure needed to further compact the aggregate increases with increasing filling factor. For higher compression velocities a compacted layer emerges underneath the moving wall, similar to the shovel of a snow plow, when pushing away snow. Once this layer has formed the pressure remains nearly constant until the layer has reached the bottom of the sample. The transition from the static towards the dynamic case occurs at compression speeds on the order of 1 m s^{1}. Since impact velocities of typical collisions of planetesimals lie within the range of 1 − 10 m s^{1} the dynamic compression behavior must be taken into account when simulating such collisions.
To determine the impact of the rescaling of the rolling and sliding forces on the very early phases in the planetesimal formation process, we plan to perform detailed collision simulations for a wide set of collision parameter. This will allow us to calculate ab initio the division between bouncing, sticking and fragmentation. This has recently been under experimental and theoretical scrutiny (Zsom et al. 2010; Geretshauser et al. 2011). It is hard to estimate the consequences of the new (stiffer) parameter set on the outcome of agglomerate collisions. We suspect that this can possibly lead to an enhanced sticking as more energy can be stored in the system before it breaks. Using a much larger particle number than previously will allow us to determine more accurately important bulk parameters such as the sound speed.
Acknowledgments
A. Seizinger acknowledges the support through the German Research Foundation (DFG) grant KL 650/16. Very helpful discussions with Carsten Güttler and Jürgen Blum from the experimental laboratory astrophysics group at Braunschweig, as well as with Carsten Dominik, Hidekazu Tanaka and Ralf Geretshauser are gratefully acknowledged. The authors acknowledge support through DFG grant KL 650/7 within the collaborative research group FOR 759 The formation of planets.
References
 Blum, J., & Wurm, G. 2000, Icarus, 143, 138 [NASA ADS] [CrossRef] [Google Scholar]
 Blum, J., & Wurm, G. 2008, ARA&A, 46, 21 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Blum, J., & Schräpler, R. 2004, Phys. Rev. Lett., 93, 115503 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Blum, J., Wurm, G., Kempf, S., et al. 2000, Phys. Rev. Lett., 85, 2426 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Brilliantov, N. V., Albers, N., Spahn, F., & Pöschel, T. 2007, Phys. Rev. E, 76, 051302 [NASA ADS] [CrossRef] [Google Scholar]
 Chokshi, A., Tielens, A. G. G. M., & Hollenbach, D. 1993, ApJ, 407, 806 [NASA ADS] [CrossRef] [Google Scholar]
 Dominik, C., & Nübold, H. 2002, Icarus, 157, 173 [NASA ADS] [CrossRef] [Google Scholar]
 Dominik, C., & Tielens, A. G. G. M. 1995, Philos. Mag. A, 72, 3, 783 [NASA ADS] [CrossRef] [Google Scholar]
 Dominik, C., & Tielens, A. G. G. M. 1996, Philos. Mag. A, 73, 5, 1279 [NASA ADS] [CrossRef] [Google Scholar]
 Dominik, C., & Tielens, A. G. G. M. 1997, ApJ, 480, 647 [NASA ADS] [CrossRef] [Google Scholar]
 Geretshauser, R. J., Speith, R., Güttler, C., Krause, M., & Blum, J. 2010, A&A, 513, A58 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Geretshauser, R. J., Meru, F., Speith, R., & Kley, W. 2011, A&A, 531, A166 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Gundlach, B., Kilias, S., Beitz, E., & Blum, J. 2011, Icarus, 214, 717 [NASA ADS] [CrossRef] [Google Scholar]
 Güttler, C., Krause, M., Geretshauser, R. J., Speith, R., & Blum, J. 2009, ApJ, 701, 130 [NASA ADS] [CrossRef] [Google Scholar]
 Güttler, C., Blum, J., Zsom, A., Ormel, C. W., & Dullemond, C. P. 2010, A&A, 513, A56 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Heim, L.O., Blum, J., Preuss, M., & Butt, H.J. 1999, Phys. Rev. Lett., 83, 3328 [NASA ADS] [CrossRef] [Google Scholar]
 Johnson, K., Kendall, K., & Roberts, A. 1971, Proc. Roy. Soc., 324, 301 [NASA ADS] [CrossRef] [Google Scholar]
 Kempf, S., Pfalzner, S., & Henning, T. K. 1999, Icarus, 141, 388 [NASA ADS] [CrossRef] [Google Scholar]
 Langkowski, D., Teiser, J., & Blum, J. 2008, ApJ, 675, 764 [NASA ADS] [CrossRef] [Google Scholar]
 Ormel, C. W., Paszun, D., Dominik, C., & Tielens, A. G. G. M. 2009, A&A, 502, 845 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Paszun, D., & Dominik, C. 2008, A&A, 484, 859 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Paszun, D., & Dominik, C. 2009, A&A, 507, 1023 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Pollack, J. B., Hubickyj, O., Bodenheimer, P., et al. 1996, Icarus, 124, 62 [NASA ADS] [CrossRef] [Google Scholar]
 Poppe, T., Blum, J., & Henning, T. 2000, ApJ, 533, 454 [NASA ADS] [CrossRef] [Google Scholar]
 Schäfer, C., Speith, R., & Kley, W. 2007, A&A, 470, 733 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Sirono, S.I. 2004, Icarus, 167, 431 [NASA ADS] [CrossRef] [Google Scholar]
 Suyama, T., Wada, K., & Tanaka, H. 2008, ApJ, 684, 1310 [NASA ADS] [CrossRef] [Google Scholar]
 Teiser, J., & Wurm, G. 2009, MNRAS, 393, 1584 [NASA ADS] [CrossRef] [Google Scholar]
 Wada, K., Tanaka, H., Suyama, T., Kimura, H., & Yamamoto, T. 2007, ApJ, 661, 320 [NASA ADS] [CrossRef] [Google Scholar]
 Wada, K., Tanaka, H., Suyama, T., Kimura, H., & Yamamoto, T. 2008, ApJ, 677, 1296 [NASA ADS] [CrossRef] [Google Scholar]
 Wada, K., Tanaka, H., Suyama, T., Kimura, H., & Yamamoto, T. 2009, ApJ, 702, 1490 [NASA ADS] [CrossRef] [Google Scholar]
 Wada, K., Tanaka, H., Suyama, T., Kimura, H., & Yamamoto, T. 2011, ApJ, 737, 36 [NASA ADS] [CrossRef] [Google Scholar]
 Watson, P. K., Mizes, H., Castellanos, A., & Pérez, A. T. 1997, Powders & Grains, 109 [Google Scholar]
 Weidenschilling, S. J. 1977, MNRAS, 180, 57 [NASA ADS] [CrossRef] [Google Scholar]
 Weidling, R., Güttler, C., Blum, J., & Brauer, F. 2009, ApJ, 696, 2036 [NASA ADS] [CrossRef] [Google Scholar]
 Windmark, F., Birnstiel, T., Güttler, C., et al. 2012, A&A, 540, A73 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Wurm, G., Paraskov, G., & Krauss, O. 2005, Icarus, 178, 253 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Zsom, A., Ormel, C. W., Güttler, C., Blum, J., & Dullemond, C. P. 2010, A&A, 513, A57 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
All Tables
All Figures
Fig. 1 The four types of particle interaction: compression/adhesion a), Rolling b), Sliding c), and Twisting d). The left panel depicts particleparticle interactions while on the right the corresponding particlewall version is displayed. 

Open with DEXTER  
In the text 
Fig. 2 Setup of the numerical simulations to obtain material parameter for porous dust agglomerates. Left: omnidirectional compression, while the top wall is moving downwards at constant speed, the sample is enclosed in a box with fixed walls on all sides. Right: unidirectional compression, as the sample is getting compressed between two walls, particles can leave the initial volume to the sides. Particles colored red are actually in contact with the walls. 

Open with DEXTER  
In the text 
Fig. 3 The compressive strength (filling factor φ versus pressure p) as obtained from experiments. The dark shaded region with the solid line fit refer to the omnidirectional compression experiment performed by Güttler et al. (2009) whereas the light shaded region and the dashed line fit were obtained for unidirectional compression by Blum & Schräpler (2004). The small image on the top left depicts the experimental setup of the omnidirectional compression experiment. (Figure taken from Güttler et al. 2009.) 

Open with DEXTER  
In the text 
Fig. 4 Compressive strength for omnidirectional compression of cubical aggregates. The result of the unaltered interaction model is compared to our improved model. The black line represents a fit to experimental data obtained from Güttler et al. (2009). 

Open with DEXTER  
In the text 
Fig. 5 Snapshots of the vertical profile of the aggregate’s filling factor averaged horizontally. Shown are results at 5 evolutionary times using 3 different speeds. The left (solid) curves represent the initial state at t = 0. As the top wall moves slowly downwards the aggregate is compacted almost homogeneously as is indicated by the nearly vertical curves. At the top and bottom of the box wall effects produce a slight deviation. 

Open with DEXTER  
In the text 
Fig. 6 Compressive strength for low velocity omnidirectional compression. 

Open with DEXTER  
In the text 
Fig. 7 Dynamic omnidirectional compression with v_{wall} between 1 − 5 m s^{1}. 

Open with DEXTER  
In the text 
Fig. 8 Dynamic omnidirectional compression with v_{wall} between 7−15 m s^{1}. 

Open with DEXTER  
In the text 
Fig. 9 Snapshots of the vertical profile of the filling factor. As the top wall moves downwards rapidly the compaction of the lower parts of the aggregate is lagging behind. The color indicates the compression speed whereas the line type indicates the position of the top wall. The solid lines show the filling factor profile of the initial uncompressed sample. 

Open with DEXTER  
In the text 
Fig. 10 Dependence of the fit parameter p_{m} of Eq. (31) on the compression speed v in the quasi static regime of v ≤ 1 m s^{1}. 

Open with DEXTER  
In the text 
Fig. 11 Dependence of the fit parameter p_{m} of Eq. (31) on the compression speed v in the dynamic regime of v ≥ 1 m s^{1}. 

Open with DEXTER  
In the text 
Fig. 12 Dependence of the fit parameter Δ of Eq. (31) on the compression speed. 

Open with DEXTER  
In the text 
Fig. 13 Compressive strength for unidirectional compression of cylindrical aggregates of different size using the original Dominik & Tielens model. The black line has been obtained from experiments by Blum & Schräpler (2004). 

Open with DEXTER  
In the text 
Fig. 14 Unidirectional compression of cylindrical aggregates using the modified model with 18 000 particles. The black line has been obtained from experiments by Blum & Schräpler (2004). 

Open with DEXTER  
In the text 