EDP Sciences
Free Access
Issue
A&A
Volume 541, May 2012
Article Number A59
Number of page(s) 11
Section Planets and planetary systems
DOI https://doi.org/10.1051/0004-6361/201218855
Published online 27 April 2012

© 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, hit-and-stick 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 sweep-up 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 meter-sized 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 shear-strength 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 cross-section 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 104. 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. Particle-particle 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: ri denotes the radius of a monomer i, γi the surface energy per area, Ei the Young’s modulus, νi the Poisson’s ratio, and Gi 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.

thumbnail Fig. 1

The four types of particle interaction: compression/adhesion a), Rolling b), Sliding c), and Twisting d). The left panel depicts particle-particle interactions while on the right the corresponding particle-wall version is displayed.

Open with DEXTER

Table 1

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/m2. 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/m2 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 x1 and x2, 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 a0 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 Fc = 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 vstick at which particles stick on head-on impacts in JKR-theory (4)For micron-sized SiO2 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 JKR-theory. 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 micron-sized SiO2 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 x1 and x2 that are in contact with each other. The contact normal nc is (5)The relative velocity vrel in normal direction of the contact is then given by (6)where v1, v2 denote the velocity of particle 1 and 2, respectively. The viscous damping force Fdamp 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 n1 and n2 be the corresponding contact pointers describing the contact. The rolling displacement ξ, sliding displacement ζ, and twisting displacement φ are then defined by: where t0 denotes the time when the contact has been established.

For the rolling motion, the forces Fr and torques Mr 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 kr, ks, and kt 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 inner-atomic 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 SiO2 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).

thumbnail 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. Particle-Wall 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 particle-particle 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 rwall → ∞.

2.2.1. Compression

The force of the particle-wall interaction in normal direction is very similar to the case of the particle-particle 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 nw, 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 particle-particle interaction. In the limit of rwall → ∞ the reduced radius R equals the radius r of a monomer:

2.2.2. Rolling

Keeping in mind that R = r, the torque Mr acting on the particle caused by rolling along the surface of the wall is given by (20)where kr,wall is equivalent to the rolling constant kr given in Eq. (17) taking the different reduced radius of the particle-wall interaction into account. nw 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 nw. 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 particle-particle interaction. Taking r2 → ∞ 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 n2 of the “wall particle” is equivalent to the normal vector nw 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 nc → nw and n2 → nw. 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 + rn1, 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 kt,wall = kt. Here we assume that particles may twist elastically around the same angle for both the particle-particle and particle-wall 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 box-shaped 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 Fw acting upon the top wall is calculated by summarizing the forces Fi exerted on the wall by grains that are currently in contact with it, where only the component in normal direction nw 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, Fw may change considerably from one integration step to the next due to the normal oscillations of the particle-wall 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 Vmat = 4/3πr3N denotes the volume of N particles of radius r and Vtot is the volume that the aggregate currently occupies. Calculating the filling factor is trivial in the case of omni-directional compression (see Fig. 2, left panel) as Vtot = Ah for a box with base size A and current height h.

However, in the uni-directional 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 Vtot 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 Aproj and the current height h of the aggregate, which equals the distance between the top and the bottom wall. Aproj 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 Leap-Frog 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 Us(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 ms = 2.5 (see Sect. 5.1) the tangential oscillation period decreases to T = 7.83   ns, which limits our timestep to 0.3   ns.

thumbnail 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 omni-directional compression experiment. (Figure taken from Güttler et al. 2009.)

Open with DEXTER

5. Results

Calibration experiments using a continuum SPH-code 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 quasi-static 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 quasi-static compression the speed vwall 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 vwall, the runtime of the simulation constitutes a lower limit of vwall. Paszun & Dominik (2008) considered vwall = 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 vwall = 1   cm/s. To ensure it is a reasonable choice we checked lower velocities down to vwall = 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  × 103 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 pm = 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, pm = 5.6   kPa, and Δ = 0.33.

thumbnail 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. Quasi-static case

In Fig. 4, we compare the experimental fitting curve to results from our simulations of box-shaped 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 pull-off 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 kr and ks (see Eqs. (17) and (18)) with correction factors that we further refer to as rolling/sliding modifiers mr and ms. 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.

thumbnail 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 mr = 8 and ms = 2.5, we obtain the red-dotted 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 quasi-static 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 RBD-generated 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 quasi-static case remain small. For higher velocities we observe an increasing deviation from the quasi-static curve. Thus, the results from Güttler et al. (2009) cannot be applied for higher velocities.

thumbnail 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 box-shaped samples are composed of about 40   000 particles and feature a base length and height of  ≈60   μm. Compared to the quasi-static 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.

thumbnail Fig. 7

Dynamic omnidirectional compression with vwall between 1 − 5   m   s-1.

Open with DEXTER

thumbnail Fig. 8

Dynamic omnidirectional compression with vwall 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 vp 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 vp = 6.98 ± 0.16   m   s-1.

thumbnail 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 continuum-simulations 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 pm and Δ serve as fitting parameters. Thus, we obtain values of pm 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 quasi-static case for low velocities and the dynamic case for higher velocities. In the beginning, pm 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 vcrit ≈ 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 pm on the compression speed v is shown for values of v ≤ 1   m   s-1. Using the ansatz pm(v) = av2 + bv + c we obtain (32)where the compression speed v is given in m s-1.

In Fig. 11 the dependence of pm 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 pm(v) = avb + 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)

thumbnail Fig. 10

Dependence of the fit parameter pm of Eq. (31) on the compression speed v in the quasi static regime of v ≤ 1   m   s-1.

Open with DEXTER

thumbnail Fig. 11

Dependence of the fit parameter pm of Eq. (31) on the compression speed v in the dynamic regime of v ≥ 1   m   s-1.

Open with DEXTER

thumbnail 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 vwall = 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 mr = 8 and ms = 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 104   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.

thumbnail 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

thumbnail 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 (kr and ks) for rolling and sliding. Indeed, the higher stiffness can be accommodated by an increase of mr = 8 and ms = 2.5 (for kr and ks, 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

All Tables

Table 1

Material parameters.

All Figures

thumbnail Fig. 1

The four types of particle interaction: compression/adhesion a), Rolling b), Sliding c), and Twisting d). The left panel depicts particle-particle interactions while on the right the corresponding particle-wall version is displayed.

Open with DEXTER
In the text
thumbnail 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
thumbnail 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 omni-directional compression experiment. (Figure taken from Güttler et al. 2009.)

Open with DEXTER
In the text
thumbnail 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
thumbnail 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
thumbnail Fig. 6

Compressive strength for low velocity omnidirectional compression.

Open with DEXTER
In the text
thumbnail Fig. 7

Dynamic omnidirectional compression with vwall between 1 − 5   m   s-1.

Open with DEXTER
In the text
thumbnail Fig. 8

Dynamic omnidirectional compression with vwall between 7−15   m   s-1.

Open with DEXTER
In the text
thumbnail 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
thumbnail Fig. 10

Dependence of the fit parameter pm 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
thumbnail Fig. 11

Dependence of the fit parameter pm of Eq. (31) on the compression speed v in the dynamic regime of v ≥ 1   m   s-1.

Open with DEXTER
In the text
thumbnail Fig. 12

Dependence of the fit parameter Δ of Eq. (31) on the compression speed.

Open with DEXTER
In the text
thumbnail 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
thumbnail 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

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.