EDP Sciences
Highlight
Free Access
Issue
A&A
Volume 554, June 2013
Article Number A4
Number of page(s) 12
Section Planets and planetary systems
DOI https://doi.org/10.1051/0004-6361/201321325
Published online 28 May 2013

© ESO, 2013

1. Introduction

Planetesimal formation is a key issue in the study of how planets form in protoplanetary disks (Hayashi et al. 1985; Weidenschilling & Cuzzi 1993). However, the collisional growth of the dust from submicron-sized dust to kilometer-sized planetesimals is still unknown.

In the growth process, one of the most important but unresolved problems is to determine the internal structural evolution of dust aggregates. The internal structure of dust is important in planetesimal formation because the dynamics of dust aggregates in protoplanetary disks are determined by coupling of gas and dust, in other words, the size and internal density of dust aggregates. In the early stage of dust coagulation in protoplanetary disks, the collision energy of the aggregates is too low to cause collisional compression (Blum 2004; Ormel et al. 2007; Zsom et al. 2010, 2011; Okuzumi et al. 2012). As a result, the internal mass density ρ decreases to ρ < 1.0   g   cm-3.

Both theoretical and experimental studies have shown that mutual collisions lead dust aggregates to have their fractal dimension D ~ 2, which is so-called ballistic cluster-cluster aggregation (BCCA; Smirnov 1990; Meakin 1991; Kempf et al. 1999; Blum & Wurm 2000; Krause & Blum 2004; Paszun & Dominik 2006). The dust aggregates would be gradually compacted or disrupted in coagulation because of the increase in impact energy. This compaction has been investigated with numerical N-body simulations that consider particle-particle interactions (Dominik & Tielens 1997; Wada et al. 2007, 2008, 2009; Suyama et al. 2008, 2012; Paszun & Dominik 2008, 2009; Seizinger et al. 2012).

Table 1

Material parameters in our simulation.

In most previous studies investigating dust growth in protoplanetary disks, dust grains have been assumed to have constant internal mass density for simplicity (Nakagawa et al. 1981; Tanaka et al. 2005; Brauer et al. 2008; Birnstiel et al. 2010). However, dust porosity evolves during dust growth in real protoplanetary disks. In recent dust coagulation calculations, porosity evolution has been considered to be based on experimental and theoretical results (Ormel et al. 2007; Okuzumi et al. 2009, 2012; Zsom et al. 2011). These results also suggest that ρ decreases as ρ ≪ 0.1 g cm-3.

In the most recent work, the dominant coagulation mode has been shown to be similar-size collisions of dust aggregates though dust grains have size distribution (Okuzumi et al. 2012). As a result, their fractal dimension is approximately equal to two, and their internal mass density ρ has been shown to become 10-5 g cm-3 (equivalent to be the filling factor φ = 10-5 for ice particles with a density of 1.0   g   cm-3). Such fluffy dust aggregates are believed to become planetesimals. Since comets in our solar system, which would be remnants of planetesimals, have their internal mass density of  ~0.1 g cm-3 (A’Hearn 2011), dust aggregates must be compressed from ρ ≪ 0.1 g cm-3 to ρ ~ 0.1 g cm-3 in protoplanetary disks.

Compression at dust aggregate collisions has been investigated in previous studies. When collisional impact energy exceeds the critical energy, dust aggregates are compacted by their collision (e.g., Dominik & Tielens 1997; Suyama et al. 2008; Wada et al. 2007, 2008, 2009). However, the collisional compression is not effective at compressing dust aggregates (Okuzumi et al. 2012).

One of the other compression mechanisms in protoplanetary disks is static compression by disk gas or self-gravity. The static compressive strength of dust aggregates has been investigated both experimentally and numerically (Paszun & Dominik 2008; Güttler et al. 2009; Seizinger et al. 2012). However, the compressive strength has examined only relatively compact aggregates with ρ ≳ 0.1 g  cm-3 because their initial aggregates are ballistic particle-cluster aggregation (BPCA) clusters. Because ρ decreases to ρ ≪ 0.1 g cm-3, at least in the early stage of dust growth, we need to reveal the static compressive strength with ρ ≪ 0.1 g cm-3.

In this work, we investigate the static compression of highly porous aggregates with ρ < 0.1 g cm-3 by means of numerical simulations and an analytical approach. It is challenging to perform numerical simulations of the static and uniform compression of highly porous aggregates. Because such porous aggregates have low sound speed, we have to compress them at a much slower velocity than in the case of compact aggregates, as is shown in our simulations. Such a slow compression of the fluffy aggregates costs much computational time.

In previous numerical studies of static compression, a dust aggregate is compressed by a wall moving in one direction (Paszun & Dominik 2008; Seizinger et al. 2012). However, this method has disadvantages when reproducing uniform and isotropic compression. There are also side walls that do not move. These side walls also obstruct the tangential motion of monomers in contact with the walls, causing artificial stress on the aggregate, which restructures them. Moreover, since they measure the pressure with the force on the moving wall, the side walls may affect the pressure measurement. In the present work, we develop a new method reproducing static compression. Instead of the walls, we adopt periodic boundary conditions and the boundaries get closer to each other. With these slowly moving periodic boundaries, the aggregate is compressed uniformly and naturally. The periodic boundary condition also enables us to represent a much larger aggregate than inside the computational region. This saves on computational time remarkably.

This paper is organized as follows. We describe the model of our numerical simulations in Sect. 2. We show the results of our simulations and find the compressive strength in Sect. 3. We confirm the obtained compressive strength formula analytically in Sect. 4, and present our conclusion in Sect. 5.

2. Simulation setting

We performed three-dimensional numerical simulations of the compression of a dust aggregate consisting of a number of spherical monomers. As the initial aggregate, we adopted a BCCA cluster. In this method, we solved interactions between all monomers in contact in each time step. Interactions between monomers in contact are formulated by Dominik & Tielens (1997) and reformulated by using the potential energies by Wada et al. (2007). We used the interaction model proposed by Wada et al. (2007) in this work. We briefly summarize the particle interaction model and material constants (see Wada et al. 2007, for details). Moreover, we describe the additional damping force in normal direction and the simulation setting in this section. In our simulations, the aggregate is gradually compressed by its copies over the moving periodic boundaries. This is an appropriate method of simulating uniform and isotropic compression. We also describe the boundary condition in this section. Since we do not have walls to measure the pressure in the periodic boundary condition, we use a similar manner of pressure measurement in molecular dynamics simulations. We also introduce the method of pressure measurement below.

2.1. Interaction model

We calculate the direct interaction of each connection of particles, taking all mechanical interactions modeled by Dominik & Tielens (1997) and Wada et al. (2007) into account. The material parameters are the monomer radius r0, surface energy γ, Young’s modulus E, Poisson’s ratio ν, and the material density ρ0. Table 1 lists the values of the material parameters for ice and silicate.

We perform N-body simulations with ice particles except for one case with silicate particles. In protoplanetary disks, ice particles are the most dominant dust material beyond the snowline. Moreover, the computational time required for calculating ice particles is less than for silicate. Thus, we adopt ice particles in most simulations. We also treat a silicate case to compare with a previous study (Seizinger et al. 2012).

The critical displacement still shows a discrepancy between theoretical (ξcrit = 2 Å) and experimental (ξcrit = 32 Å) studies (Dominik & Tielens 1997; Heim et al. 1999). We adopt the same parameter as in Wada et al. (2011), ξcrit = 8 Å as a typical length for ice particles, and ξcrit = 20 Å for silicate particles to compare with Seizinger et al. (2012).

The parameter ξcrit is related to strength of rolling motion. The rolling motion between monomers is crucial in compression. The rolling energy Eroll is the energy required to rotate a particle around a connecting point by 90°. The rolling energy can be written as (1)(see Wada et al. 2007, for details). In the case of ice monomers, for example, Eroll = 4.37 × 10-9 erg for ξcrit = 8 Å.

We use a normalized unit of time in our simulations. For ice particles, the normalized unit of time is (2)which is a characteristic time, and it approximately represents the oscillation time of particles in contact at the critical collision velocity (see Wada et al. 2007, for details).

2.2. Damping force in normal direction

The normal force between two monomers is repulsive when the monomers are close or attractive when they are stretched out. Thus, normal oscillations occur at each connection. For realistic particles, these oscillations would dissipate because of viscoelasticity or hysteresis in the normal force (e.g., Greenwood & Johnson 2006; Tanaka et al. 2012). For such damping of normal oscillation, we add an artificial normal damping force to the particle interaction model, following the previous studies (Suyama et al. 2008; Paszun & Dominik 2008; Seizinger et al. 2012).

Assuming that two particles in contact have their position vectors x1 and x2, respectively, the contact unit vector nc is defined as (3)(see Fig. 2 in Wada et al. 2007). We introduce a damping force between contact particles in normal direction, defined as (4)where kn is the damping coefficient in normal direction and m0 the monomer mass. The adopted value of kn is on the order of 0.01. To show that the result is independent of the normal oscillation damping, we perform N-body simulations with the damping factor kn as a parameter.

The timescale of damping is (5)for kn = 0.01, it is much shorter than the simulation timescale, which is typically  ~107t0. We show that the obtained compressive strength is independent of the artificial normal damping force in our simulations (see Sect. 3.4).

2.3. Uniform compression by moving boundaries

We adopt the periodic boundary condition in our simulations. The aggregate in the computational region is surrounded by its copies, as shown in Fig. 1.

thumbnail Fig. 1

Schematic drawing of the periodic boundary condition. Each box illustrates a boundary box with a side length L for all directions. When the boundary starts to get closer, the aggregate sticks to the neighboring aggregates over the boundary and is compressed by them. It should be noted that this picture is illustrated in the 2D direction, but our simulations are performed in 3D.

Open with DEXTER

Initially, we set a cubic box whose sides are periodic boundaries with a size of L larger than the aggregate. Thus, the initial BCCA cluster is detached from its neighboring copies over the periodic boundaries. In our simulations, we gradually move the boundaries to the center of the aggregate to get closer to one another. As a result, the aggregate sticks to the neighboring copies and is compressed by them in a natural way. Therefore, the aggregate in the computational region corresponds to a small part of a whole large aggregate. In other words, although the number of particles in numerical simulations are limited because of computational cost, the periodic boundary condition enables us to investigate a large aggregate, such as a  ~ cm-sized dust aggregate in protoplanetary disks.

Another advantage of the periodic boundary condition is that we do not need to introduce the wall for compression. In the previous N-body simulations of static compression, dust aggregates are compressed by using the wall against the dust aggregate (Paszun & Dominik 2008; Seizinger et al. 2012). The wall itself may have some artificial effects on such experiments. For example, the wall moves in one direction and thus this may be different from isotropic compression. Besides, wall-particle interaction is different from particle-particle interaction, so it must be treated carefully. In contrast, the periodic boundary condition does not need walls for compression because a dust aggregate is compressed by the neighboring aggregate over the periodic boundary. In addition, the periodic boundaries in three directions make it possible to compress the aggregate isotropically. We calculate not only the interactions of particles in contact inside the computational region but also the interactions of the particles in contact across the periodic boundaries. Thus, no special treatment of interactions, which is wall-particle interactions in the case of simulations with walls in previous studies, is required when a particle crosses the periodic boundaries.

The computational cubic region has length L, and the coordinates in x,y, and z directions are set to be  − L/2 < x < L/2,  − L/2 < y < L/2, and  − L/2 < z < L/2, respectively. We adopt periodic boundary conditions for all directions to reproduce a part of a large aggregate, where L decreases with time t, L = L(t). The initial size of the box L0 is adopted as the maximum size of the dust aggregate in x, y, and z directions.

With the settings above, we move the boundaries of the computational region toward the center of the region. The velocity at the boundary is given by (6)where Cv is a dimensionless parameter (we call Cv the strain rate parameter hereafter). Owing to this definition of the boundary speed, the aggregate is compressed at a constant strain rate independent of the region scale L.

The box size decreases with the constant rate Cv in three directions. This corresponds to isotropic compression. Since , the box size is written as (7)Therefore, the whole time of compression is t0/Cv. Typically we chose Cv = 3 × 10-7 so the compression time is  ~3 × 106   t0 ~ 0.6 ms.

When a particle crosses a periodic boundary, the velocity should be treated carefully to reproduce the quasi-static compression with periodic boundary condition. Figure 2 illustrates how to calculate the velocity of particles across the periodic boundary.

thumbnail Fig. 2

Schematic drawing to illustrate how the particle velocity is calculated when a particle crosses a periodic boundary. For simplicity, we consider this situation in a 2D field, but we actually calculate this in a 3D situation. We consider that a dust particle is close to the boundary in the left figure. In the next time step, the particle crosses the boundary (dashed circle in the right figure). We put the particle on the other side of the boundary as expressed in Eqs. (8) and (10). The velocity component is converted as expressed in Eqs. (9) and (11). This treatment reproduces the isotropic compression in the velocity field well.

Open with DEXTER

When a particle goes out of the computational region across the boundary at x = L/2, we relocate the particle to the opposite side (i.e., from the boundary at x =  −L/2). In that case, the position of the particle in x direction is converted as (8)Since the two boundaries at x =  − L/2 and x = L/2 have a relative velocity of 2vb, the x-component of the velocity vx of the particle is also converted as (9)Owing to the conversion of vx, the velocity of particle against the boundary that the particle crosses does not change before and after the crossing. For a particle across the boundary at x =  − L/2, the position and the velocity are converted as \arraycolsep1.75ptWe also have the same treatments for particles across the boundaries at y =  ± L/2 and z =  ± L/2.

We introduce the constant strain rate at the boundaries for scaleless discussion. However, the initial aggregate is not moving. As the simulation starts, if all the particles in the aggregate are not moving, only the particles close to the boundaries have initial velocity. This is not a constant strain rate. To reproduce the scaleless constant strain rate initially, therefore, we first give all monomers the velocity smoothly connected to the boundary speed. The initial velocity is expressed as (12)where r is the position vector of the monomers.

thumbnail Fig. 3

Snapshots of the evolution of an aggregate under compression in the case of N = 16   384. The top three figures are 3D visualizations. They have the same scale with different time epochs. The white particles are inside a box enclosed by the periodic boundaries. The yellow particles are in neighboring boxes to the box of white particles. For visualization, we do not draw the copies on the back and front sides of the boundaries but only 8 copies of the white particles across the boundaries. Each bottom figure represents projected positions onto 2D plane of all particles in each corresponding top figure. The gray points in the bottom figures correspond to the positions of the white particles in the top figures, and the yellow points correspond to those of the yellow particles in the top figures. Scales are in μm.

Open with DEXTER

2.4. Pressure measurement

In previous studies, a dust aggregate is enclosed by walls, and the pressure is calculated by measuring the force exerted on the walls by the dust aggregate. In this work, a dust aggregate is compressed by themselves because of the periodic boundary condition. Therefore, we introduce another method of measuring the pressure on the aggregate. We calculate the pressure of the dust aggregate in the standard way in molecular dynamics simulations using the virial theorem as follows (e.g., Haile 1992).

thumbnail Fig. 4

a) Pressure P in [Pa] against filling factor φ. The ten thin solid lines show the results for the initial BCCA clusters with different initial random numbers and thick solid line shows the arithmetic average of the ten runs. b) Pressure P in [Pa] against filling factor φ. Same as the thick solid line in a) plotted with a dotted line of Eq. (25). The parameters are N = 16   384, Cv = 3 × 10-7, kn = 0.01, and ξcrit = 8   Å.

Open with DEXTER

We consider a virtual box that encloses the aggregate under consideration. We define the force acting from the walls of the virtual box on the particle i as Wi, and the sum of the forces from other particles on the particle i as Fi. The equation of motion of the particle i is given by (13)We take a scalar product of both sides of the equation with ri and take a long time average of both sides with time interval τ. The lefthand side becomes (14)The first term on the righthand side vanishes in the limit of τ → ∞. We define the taking-a-long-time average in t as  ⟨  ⟩ t. Taking a summation of all particles and a long time average of Eq. (13), we obtain (15)The first term on the righthand side is related to pressure P. The pressure is an average of all forces acting on the wall from all particles. Using the normal vector n of the wall surface directed outward, the force received by the wall that has an area dS is PndS. Therefore, (16)This equation is obtained by taking surface integral as (17)The translational kinetic energy K, averaged over a long time, is given by (18)Using K and P, Eq. (15) gives an expression of P as (19)We define the force from particle j on particle i as fi,j. Force Fi can be written as a summation of the force from another particle as (20)Using fi,j =  − fj,i, we finally obtain the pressure measuring formula as (21)The first term on the righthand side of the equation represents the translational kinetic energy per unit volume, and the second term represents the summation of the force acting at all connections per unit volume. This expression is useful for measuring the pressure of a dust aggregate under compression. We do not need to put any artificial object, such as walls, in simulations because Eq. (21) is totally expressed in terms of the summation of the physical quantities of each particle, which are the mass, the position, the velocity, and the force acting on the particle. In our calculations, we take an average of pressure for every 10 000 time steps, corresponding to 1000 t0 because we set 0.1 t0 as one time step in our simulation.

As mentioned in Sect. 2.2, the adopted damping force corresponds to rapid damping of normal oscillations. Thus, the kinetic energy of random motion rapidly dissipates. This corresponds to the static compression, and thus the compressive strength is determined by the second term of Eq. (21).

3. Results

The top three panels of Fig. 3 show snapshots of the evolution of an aggregate under compression in the case where N = 16   384, Cv = 3 × 10-7, kn = 0, and ξcrit = 8 Å. The top three panels have the same scale but different time epochs, which are t = 0, 1 × 106t0, and 2 × 106t0. The white particles are inside the computational region enclosed by the periodic boundaries, while the yellow particles are in the neighboring copy regions (for visualization, we do not draw particles on the front and backsides copy regions). The bottom three panels represent the projected positions onto the two-dimensional plane for the correspondent top three figures. We confirm that the dust aggregate is compressed by their copies from all directions. As the compression proceeds, the aggregate of white particles is compressed by the neighboring aggregate of yellow particles. We focus on how high pressure is generated by quasi-static compression in numerical simulations. Our numerical simulations have several parameters: the size of the initial BCCA cluster, the compression rate, the normal damping force, and the critical displacement (corresponds to the rolling energy). We investigate the dependence of the pressure on these parameters, by performing several runs with different parameter sets. Although we assume ice aggregates in most runs, we also investigate cases of silicate aggregates to compare them with previous studies.

3.1. Fiducial run: obtaining the compressive strength

We put a BCCA cluster as the initial aggregate. For each set of parameters, we randomly create ten BCCA clusters following Okuzumi et al. (2009) and take arithmetic averages of the ten simulations of the different initial clusters. The pressure is measured using Eq. (21) at each run. We define the filling factor of an aggregate as (22)where V0 is the monomer volume, N the number of monomers of the aggregate, and V the volume enclosed by the boundaries, which has a length of L. The filling factor also can be written as φ = ρ/ρ0. Figure 4 shows that the measured pressure as a function of the filling factor φ(t). The parameters of the simulations are N = 16   384, Cv = 3 × 10-7, kn = 0.01, and ξcrit = 8 Å. The corresponding Eroll is 4.74 × 10-9   erg for ξcrit = 8 Å. Each colored line in Fig. 4a shows each simulation with the different initial shape of the aggregate. Figure 4b shows the arithmetic average of the pressure measured in ten different runs. Each line shows in different ranges of φ. The lowest φ is determined with the largest size of the initial boundary boxes of the ten runs. We find that the compressive strength is reproduced well by (23)where P0 = 4.74 × 105 Pa. We analytically discuss why the compressive strength is proportional to φ3 in Sect. 4. In the high-density region (φ ≳ 10-1), the measured strength deviates from the line of P = P0φ3. This is because the dissipation mechanism changes in the high-density region (see Sect. 3.4). The deviation in the low-density region (φ ≲ 3 × 10-3) is partly caused by a finite boundary speed (or compression rate) as discussed in the next section. Another reason for the deviation in the low-density region is related to the density of the initial BCCA cluster. The filling factor of BCCA φBCCA is estimated as (24)where we use the radius and the volume of a BCCA cluster, and , respectively (e.g., Suyama et al. 2008). For N = 16   384, we obtain φBCCA ~ 3 × 10-3. In the early stage of compression, φ is lower than φBCCA because the initial BCCA clusters are apart from each other. This space between BCCA clusters would also cause the deviation from the line of P = P0φ3.

We now discuss the coefficient P0 of the compressive strength. Wada et al. (2008) show that Eroll is important in the collisional compressive strength. Thus, Eroll is expected to also be important in the static compressive strength. Considering that the characteristic volume is the monomer’s volume  ~ , we suppose , based on dimension analysis. Therefore, the compressive strength can be written as (25)We analytically discuss and confirm this equation in Sect. 4. We also plot this equation in Fig. 4b. This figure clearly shows that the result is well fit by Eq. (25).

We show that compressive strength is proportional to ξcrit, which is proportional to the rolling energy Eroll in Sect. 3.5. We also confirm that Eq. (25) is applicable to the case of different r0 in the silicate case.

3.2. Dependence on the boundary speed

To statically compress the aggregate, we should move the boundary at a low enough velocity not to create inhomogeneous structure. Figure 5 shows the dependency on the strain rate parameter.

thumbnail Fig. 5

Pressure P in [Pa] against filling factor φ with different strain rate parameter Cv. Each line shows the average of ten runs of the fixed strain rate: Cv = 1 × 10-7, 3 × 10-7, 1 × 10-6, 3 × 10-6, 1 × 10-5. The other parameters are the same for every ten runs: N = 16   384, kn = 0.01, and ξcrit = 8 Å. The dashed line is Eq. (25).

Open with DEXTER

Each line shows the average of ten runs. The fixed parameters are N = 16 384, kn = 0.01, and ξcrit = 8 Å. The strain rate parameter Cv is equal to 1 × 10-7, 3 × 10-7, 1 × 10-6, 3 × 10-6, and 1 × 10-5. The higher Cv, the higher pressure in the low-density region is required for compression. This is mainly caused by the ram pressure from the boundaries with high speed.

When the compression proceeds and the density becomes higher to reach the line of Eq. (25), the pressure follows the equation. From Fig. 5, Cv = 3 × 10-7 creates a sufficiently low boundary speed. The boundary speed can be calculated as a function of φ. Using Eq. (6) and , the velocity difference between a boundary and the next boundary, vd, can be written as (26)In the case of Cv = 3 × 10-7, vd = 12.7, 5.9, and 2.7 cm/s for φ = 10-3, 10-2, and 10-1, respectively.

Here, we discuss the velocity difference of boundaries, comparing with the effective sound speed of the aggregates. The effective sound speed can be estimated as (27)where we use Eq. (25). Using the rolling energy of ice particles, cs,eff is given by (28)Therefore, in the case of Cv = 3 × 10-7, vd is not low enough in the beginning of the simulation, where the aggregate has a low filling factor. However, the boundary velocity difference reaches lower than the effective sound speed when φ ≳ 10-2.

3.3. Dependence on the size of the initial BCCA cluster

To confirm that Eq. (25) is valid in the lower density region, we perform the simulations with the different number of particles, which is equivalent to the different sizes of the initial dust aggregates. Figure 6 shows dependence on the number of particles of the initial BCCA cluster.

thumbnail Fig. 6

Pressure P in [Pa] against filling factor φ with a different number of particles N. Each line shows the average of ten runs of the fixed number of particles: N = 1024,4096, and 16   384. The other parameters are Cv = 3 × 10-7, kn = 0.01, and ξcrit = 8 Å in the case of N = 1024, 4096, and Cv = 1 × 10-7, kn = 0.01, and ξcrit = 8   Å in the case of N = 16   384. The dashed line is Eq. (25).

Open with DEXTER

The initial numbers of particles are 1024, 4096, and 16 384. The other parameters are Cv = 3 × 10-7, kn = 0.01, and ξcrit = 8 Å in the case of N = 1024 and N = 4096, and Cv = 1 × 10-7, kn = 0.01, and ξcrit = 8 Å in the case of N = 16   384. We chose lower Cv in the case of N = 16   384 to investigate the strength in lower φ region. Each line represents the average of ten runs for each simulation as in Figs. 4b and 5. We draw the averaged line from lower φ than in Fig. 5. In such a low φ region, we consider that the pressure is zero for some runs because the aggregate is isolated from the copies of the aggregate over the periodic boundaries. Except for the initial deviation in low φ, all lines show good agreement with Eq. (25) where φ ≲ 0.1. The result agrees in lower φ for runs with larger N. Therefore, we conclude that the formula Eq. (25) is valid for φ ≲ 0.1.

3.4. Dependence on the normal damping force

As described in Sect. 2.2, we adopt the normal damping force to reduce the normal oscillations in addition to Wada et al. (2007). To confirm that this damping factor does not affect the simulation results, we set the damping factor kn as a parameter. Figure 7 shows dependence of pressure on the normal damping factor kn.

thumbnail Fig. 7

Pressure P in [Pa] against filling factor φ with different normal damping force. We put the same ten initial conditions varying the normal damping force with kn = 0, kn = 10-2, and kn = 101. Each line shows the result of one run. The other parameters are N = 16   384, Cv = 3 × 10-7, and ξcrit = 8 Å.

Open with DEXTER

The fixed parameters are N = 16   384 Cv = 3 × 10-7, and ξcrit = 8 Å. Each line represents the result of one run for kn = 0, 10-2, and 101, respectively. This figure clearly shows that the normal damping force does not affect the simulation results.

As mentioned in Sect. 3.1, the compressive strength in the low-density region (φ ≲ 0.1) is expected to be determined by the rolling motion. To confirm this, we calculate the total energy dissipations of all motions, which are normal damping, rolling, sliding, and twisting. Figure 8 shows the dissipated energy for each mechanism, the dissipated energies in the case without the normal damping and those in the case of kn = 0.01.

thumbnail Fig. 8

Energy dissipation of each dissipation mechanism in [erg] against filling factor φ. The solid lines show the result in the case without the normal damping and the dashed lines in the case of kn = 0.01 and The results in the case of kn = 10 are not plotted because they are the same as those in the case of kn = 0.01 and indistinguishable. The dissipation mechanisms are normal damping, rolling, sliding, and twisting. The dissipation energy by sliding motion is less than 10-9 erg.

Open with DEXTER

The dissipated energy in the case of kn = 10 is indistinguishable from those in the case of kn = 0.01, and thus we do not plot them. The dissipation energy of the sliding force is less than 10-9 erg, so it is not depicted in this figure. The dissipation by the rolling and twisting is almost the same in the cases with and without the normal damping. Thus, we confirm that the normal damping does not affect the compressive strength, although it dissipates the energy of the normal oscillations. Aside from the normal dissipation, the dominant dissipation mechanism is the rolling motion. This clearly shows that the static compression is determined by the rolling motion of each connection, as mentioned in Sect. 3.1. Where φ ≳ 0.1, the energy dissipation by twisting motion occurs. This is why Eq. (25) is valid until the filling factor reaches 0.1 as mentioned in Sect. 3.1. In the high-density region, where φ ≳ 0.1, another formulation is required but that is beyond the scope of this paper.

3.5. Dependence on the rolling energy

We also investigate the dependence of the compressive strength on ξcrit. Since Eroll is proportional to ξcrit, we investigate the dependence on the rolling energy in this section. Figure 9 shows the dependency on ξcrit.

thumbnail Fig. 9

Pressure P in [Pa] against filling factor φ with different critical displacement, ξcrit. We put the same ten initial conditions varying ξcrit with ξcrit = 32,16,8,4, and 2 Å, respectively. Each line shows the average of ten runs. The other parameters are N = 16   384 Cv = 3 × 10-7, and kn = 10-2.

Open with DEXTER

We vary ξcrit with 32, 16, 8, 4, and 2 Å. The fixed parameters are N = 16   384, Cv = 3 × 10-7, and kn = 10-2. This result shows that the compressive strength is almost the same in the low-density region. This is because the periodic boundary creates the additional voids as discussed in Sect. 3.1, so we should not focus on the low-density region. The lines in the case of ξcrit = 2,4, and 8 Å are on the corresponding lines of Eq. (25) where φ ≲ 0.1. The line in the case of ξcrit = 16 Å has a little deviation, and in the case of ξcrit = 32 Å it has a deviation from their corresponding lines of Eq. (25). The reason the lines in the case of ξcrit = 16,32 Å deviate from the corresponding lines of Eq. (25) is that the dissipation energy is dominated not by rolling motion but by a twisting motion as indicated in Fig. 10. This figure shows that dissipated energy of each dissipation mechanism.

thumbnail Fig. 10

Energy dissipation of each dissipation mechanism in [erg] against filling factor φ. Each panel represents the case of different ξcrit, which are 32, 16, and 8 Å, corresponding to Fig. 9.

Open with DEXTER

We show the results of the cases with ξcrit = 8, 16, and 32 Å. The normal damping does not contribute to the compressive strength as discussed in Sect. 3.4, and thus we focus on the rolling and twisting motions.

When ξcrit ≤ 8 Å, the dissipation energy is dominated by rolling motion. For ξcrit = 32 Å, on the other hand, the dissipation energy is dominated by a twisting motion. For ξcrit = 16 Å, the dissipation energy of rolling and twisting motion is comparable, and making it the marginal case. Thus, the reason Eq. (25) is not valid when ξcrit ≥ 16 Å is that the twisting motion is the dominant mechanism for determining the compressive strength. Therefore, we conclude that Eq. (25) is valid when ξcrit ≤ 8 Å.

3.6. Fractal structure

We also investigate how the fractal structure of the dust aggregate changes. Figure 11 shows how many particles are inside the distance rin for four snapshots.

thumbnail Fig. 11

Number of particles inside the radius rin against normalized radius rin/r0. This figure represents the fractal structure of the compressed aggregates in our simulation for various φ. We set a particle as an origin and count the number of particles inside r < rin, where r is the distance from the origin to each particle’s center. Then we count the same correlation of all particles as an origin and take their average (similar figure of Fig. 7 in the paper of Wada et al. 2008). Each line shows the result at the different time steps. The solid thick lines represents the structure of fractal dimension D = 2, and dashed lines represent D = 3 for each corresponding φ. The dotted line shows the number of particles in calculation. The region below this line corresponds to inside the periodic boundaries.

Open with DEXTER

We select one run from the case with N = 16   384, Cv = 3 × 10-7, kn = 10-2, and ξcrit = 8 Å. Each snapshot is when φ = 0.003,0.01,0.03 and 0.1, respectively. We take a particle as an origin and count the number of particles inside r < rin, where r is the length from the origin. Then we set the computational region as an origin for all the other particles inside and take an average of them. We obtain the same trend in several runs in the cases of different shapes of initial aggregates.

We also count particles beyond the periodic boundaries. In high rin, because of copies over the periodic boundary distributed as a fractal dimension of three. Therefore, where N(r < rin) ≳ 16   384, N must be . However, it is almost out of the range of Fig. 11. Figure 11 shows the number of particles in calculation, which is N = 16   384. The results over this line are affected by the periodic boundary condition and those below this line are in computational region. Thus, the results below the line represent the fractal structure inside the computational region and are not the artificial effect of the periodic boundary condition.

Since the initial aggregate is a BCCA cluster, N is proportional to . In the case of φ = 0.003, which is equivalent to φ of the initial BCCA cluster, as shown in Fig. 11. When the fractal dimension is three, N can be written as (29)where . We also plot this equation for each φ in Fig. 11, with good agreement on a large scale, while maintaining on a small scale.

Therefore, the structure evolution in the static compression is as follows. Initially, because the aggregate is a BCCA cluster. As compression proceeds, the fractal dimension D becomes three on a large scale, while it is two on a small scale. The transit scale from D = 2 to D = 3 becomes smaller as compression proceeds until D = 3 on any scale. This structure evolution means that the static compression reconstructs the aggregate first on a large scale when keeping the small-scale BCCA structure. This is the reason the rolling motion determines the compressive strength, as discussed in Sect. 4.

3.7. Silicate case: comparison with previous studies

The compressive strength has been investigated in the previous study (Seizinger et al. 2012). To investigate the connection of compressive strength from the low-density to the high-density regions, we perform simulations in the case of silicate with the same parameters of Seizinger et al. (2012). Figure 12 shows compression in the case of silicate whose monomer size is 0.6   μm.

thumbnail Fig. 12

Pressure P in [Pa] against filling factor φ. This figure is same as Fig. 4 but for the case of silicate particles (r0 = 0.6   μm).

Open with DEXTER

The parameters are N = 16   384, Cv = 3 × 10-7, and kn = 0.01. Figure 12a shows the results of ten runs with different initial aggregates and Fig. 12b shows their average. Using the rolling energy of silicate, which is Eroll = 1.42 × 10-8 erg, we also plot the line of Eq. (25) in Fig. 12b. Since t0 is given by 1.71 × 10-9 s for silicate aggregates, vd becomes 4.01 cm/s for φ = 10-2 with Cv = 3 × 10-7. This vd is larger than cs,eff( = 0.77 cm/s when φ = 0.01) for silicate aggregates, allowing the numerical results shown in Fig. 12 to deviate from the line of Eq. (25) in the low φ region. When vd = cs,eff, φ = 3.4 × 10-2, the compressive strength should obey Eq. (25) when φ ≳ 3.4 × 10-2. In the case of silicate, computational time is huge compared with ice particles. We take a relatively high value of the boundary speed to save on computational time. Therefore, the result deviates from Eq. (25) in the low-density region because of the high velocity. In other words, the compression is not static in the low density region. In the high-density region, on the other hand, the result is in good agreement with Eq. (25), suggesting that Eq. (25) is applicable to aggregates consisting of silicate particles with different r0.

Figure 13 compares our simulation results and Eq. (25) in the low-density region (φ < 0.1) with the results of Seizinger et al. (2012) and the fitting formula to experiments (Güttler et al. 2009). This figure corresponds to Fig. 4 in Seizinger et al. (2012). They performed similar N-body simulations to ours but using a BPCA aggregate composed of silicate particles as an initial condition. The compressive strength of our simulations shows good agreement with the same interaction model in Seizinger et al. (2012) with a little discrepancy: φ = 0.24 at P = 300 Pa in our simulations and φ = 0.21 at P = 300 Pa in Seizinger et al. (2012). The discrepancy, 13% in φ, may be caused by the difference in the initial aggregate or the pressure measurement method. The fitting formula of Güttler et al. (2009) suggests φ = 0.17 at P = 300 in the experiments. The discrepancy from our simulations is 29% in φ. In applicable uses of the static compression formula, we focus on obtaining φ with a given P.

4. Understanding the compressive strength formula

In this section, we analytically derive the compressive strength and confirm Eq. (25). First, we consider the structure of a fluffy aggregate in static compression in our simulations. As described in Sect. 2.3, we adopt the periodic boundary condition and put a BCCA cluster as the initial condition. This corresponds to a large aggregate that filled up with BCCA clusters three dimensionally. As compression proceeds, the initial BCCA cluster is compressed but the aggregate keeps smaller BCCA structure as confirmed in Sect. 3.6. Therefore, the aggregate in static compression always consists of BCCA clusters on some scale and filled up with them. Figure 14 illustrates the aggregate in static compression. The enclosed lines depict BCCA clusters on a small scale.

thumbnail Fig. 13

The filling factor φ against pressure P in [Pa]. This figure is same as Fig. 12, but plotted with a linear scale of φ and reversal of xy axis to compare with previous studies (see Fig. 4 in Seizinger et al. 2012). The dotted line is the result of numerical simulations in the high-density region (φ ≳ 0.1) in Seizinger et al. (2012) and the thin solid line is the fitting formula proposed by Güttler et al. (2009). Our results consistently connect to the previous simulations in the high-density region.

Open with DEXTER

thumbnail Fig. 14

Schematic drawing of compression of a dust aggregate consisting of a number of BCCA clusters. The left figure shows a dust aggregate consisting of many BCCA clusters and the BCCA clusters are distributed three dimensionally. Each enclosed line represents each region dominated by the BCCA clusters. The central figure is a BCCA cluster, receiving pressure from the next clusters. The BCCA cluster has a large void depicted in the central figure, and thus the void would be compressed, as expressed in the right figure. The required energy to compress the void is the energy to rotate the connection of monomers in contact. Therefore, the compression can be determined by the rolling motion of monomer connection on the connecting point of the subclusters.

Open with DEXTER

Next, we consider why the compressive strength can be determined by the rolling energy. The internal mass density and the volume filling factor of the aggregate are equal to those of the BCCA clusters. Compression of the whole aggregate proceeds by compression of each cluster. Therefore, the compressive strength of the whole aggregate would be determined by BCCA clusters. The righthand panel of Fig. 14 illustrates compression of one of BCCA clusters. The pressure on the BCCA cluster is exerted by neighbor clusters, which causes the compression of the BCCA cluster. The BCCA cluster can be divided further into two smaller subclusters because BCCA clusters are created by cluster-cluster aggregation. A large void exists between the two smaller clusters, and they are connected with one connection of monomers in contact. The compression of the BCCA cluster occurs by crashing the large void, which only requires rolling the monomers at the connection.

We now estimate the compressive strength. In static compression, the aggregate is compressed by external pressure. Each BCCA cluster feels a similar pressure, P. Using the pressure, the force on the BCCA cluster is approximately given by (30)Since the crashing of the large void is accompanied by the rolling of a pair of monomers in contact, the work required for the crashing is given by so-called the rolling energy of monomers, Eroll (Dominik & Tielens (1997) or see Eq. (1) for its definition). Therefore, the required force to compress the aggregate satisfies (31)Substituting Eq. (30), we further obtain the required pressure to compress the aggregate as (32)The radius of the BCCA clusters can be written by using the physical values of the whole aggregate. The internal density of the BCCA cluster is dependent on its radius. The BCCA cluster has the fractal dimension of two, and its radius is approximately given by rBCCA = N1/2r0, where N is the number of constituent monomers in the BCCA subcluster. The internal density of the BCCA cluster is evaluated as (33)Using Eqs. (32) and (33), we finally obtain the required pressure (or the compressive strength) as (34)This is the same as Eq. (25) obtained from our numerical simulations.

5. Summary

We investigated the static compressive strength of highly porous dust aggregates, whose filling factor φ is lower than 0.1. We performed numerical N-body simulations of static compression of highly porous dust aggregates. The initial dust aggregate is assumed to be a BCCA cluster. The particle-particle interaction model is based on Dominik & Tielens (1997) and Wada et al. (2007). We introduced a new method for compression and adopted the periodic boundary condition in order to compress the dust aggregate uniformly and naturally. Because of the periodic boundary condition, the dust aggregate in the computational region represents one part of a large aggregate, and thus we could investigate the compression of a large aggregate. The periodic boundaries move toward the center, and the distance between the boundaries becomes small. To measure the pressure of the aggregate, we adopted a similar manner to the one used in molecular-dynamics simulations. As a result of the numerical simulations, our main findings are as follows.

  • The relation between the compressionpressure P and the filling factor φ can be written as (35)where Eroll is the rolling energy of monomer particles and r0 the monomer radius. We defined the filling factor as φ = ρ/ρ0, where ρ is the mass density of the whole aggregate, and ρ0 is the material mass density. Equation (35) is independent of the numerical parameters: the number of particles, the size of the initial BCCA cluster, the boundary speed, and the normal damping force. We confirmed that Eq. (35) is applicable in different Eroll and r0. We also analytically confirmed Eq. (35).

  • Equation (35) is valid where φ ≲ 0.1 in the high-density region. In the low-density region, we confirmed that Eq. (35) is valid for φ ≳ 10-3 in the case of N = 16   384. From the results of different initial sizes of the aggregates, Eq. (35) is valid in the lower-density region in the case of the larger aggregates.

  • The initial BCCA cluster has a fractal dimension of two in the radius of the cluster, although the whole aggregate has a fractal dimension of three because of the periodic boundary. As compression proceeds, the fractal dimension inside the radius of the initial BCCA cluster becomes three, while the fractal dimension on a smaller scale keeps being two. This means that the initial set up, which is that the fractal dimension on a large scale is three and that on a small scale is two, reproduce the structure of a dust aggregate well in static compression as a consequence. This also supports the compressive strength being determined by BCCA structure on a small scale as shown in the analytical approach.

  • The static compression in the high-density region (φ ≳ 0.1) has been investigated in the silicate case in previous studies (e.g., Seizinger et al. 2012). We also performed the numerical simulations in the silicate case and confirmed that our results are consistent with those of previous studies in the high-density region.

The compressive strength formula allowed us to study how static compression affects the porosity evolution of dust aggregates in protoplanetary disks. In applications to dust compression in protoplanetary disks, we use the compressive strength formula to obtain φ with a given P. Moreover, the obtained compressive strength would be applicable to SPH simulations of dust collisions (Geretshauser et al. 2010, 2011). This application of the static compression process is important future work. In this work, we did not study shear or tensile strengths, but they are also worth investigating in future work.

Acknowledgments

We thank Kohji Tomisaka and Hiroki Senshu for fruitful discussions. We appreciate the careful reading and fruitful comments by the anonymous referee and the editor, Tristan Guillot. A.K. is supported by the Research Fellowship from the Japan Society for the Promotion of Science (JSPS) for Young Scientists (24·2120). S.O. acknowledges support by Grants-in-Aid for

JSPS Fellows (22·7006) from MEXT of Japan. K.W. acknowledges support by Grants-in- Aid for Scientific Research (24540459) from MEXT of Japan.

References

All Tables

Table 1

Material parameters in our simulation.

All Figures

thumbnail Fig. 1

Schematic drawing of the periodic boundary condition. Each box illustrates a boundary box with a side length L for all directions. When the boundary starts to get closer, the aggregate sticks to the neighboring aggregates over the boundary and is compressed by them. It should be noted that this picture is illustrated in the 2D direction, but our simulations are performed in 3D.

Open with DEXTER
In the text
thumbnail Fig. 2

Schematic drawing to illustrate how the particle velocity is calculated when a particle crosses a periodic boundary. For simplicity, we consider this situation in a 2D field, but we actually calculate this in a 3D situation. We consider that a dust particle is close to the boundary in the left figure. In the next time step, the particle crosses the boundary (dashed circle in the right figure). We put the particle on the other side of the boundary as expressed in Eqs. (8) and (10). The velocity component is converted as expressed in Eqs. (9) and (11). This treatment reproduces the isotropic compression in the velocity field well.

Open with DEXTER
In the text
thumbnail Fig. 3

Snapshots of the evolution of an aggregate under compression in the case of N = 16   384. The top three figures are 3D visualizations. They have the same scale with different time epochs. The white particles are inside a box enclosed by the periodic boundaries. The yellow particles are in neighboring boxes to the box of white particles. For visualization, we do not draw the copies on the back and front sides of the boundaries but only 8 copies of the white particles across the boundaries. Each bottom figure represents projected positions onto 2D plane of all particles in each corresponding top figure. The gray points in the bottom figures correspond to the positions of the white particles in the top figures, and the yellow points correspond to those of the yellow particles in the top figures. Scales are in μm.

Open with DEXTER
In the text
thumbnail Fig. 4

a) Pressure P in [Pa] against filling factor φ. The ten thin solid lines show the results for the initial BCCA clusters with different initial random numbers and thick solid line shows the arithmetic average of the ten runs. b) Pressure P in [Pa] against filling factor φ. Same as the thick solid line in a) plotted with a dotted line of Eq. (25). The parameters are N = 16   384, Cv = 3 × 10-7, kn = 0.01, and ξcrit = 8   Å.

Open with DEXTER
In the text
thumbnail Fig. 5

Pressure P in [Pa] against filling factor φ with different strain rate parameter Cv. Each line shows the average of ten runs of the fixed strain rate: Cv = 1 × 10-7, 3 × 10-7, 1 × 10-6, 3 × 10-6, 1 × 10-5. The other parameters are the same for every ten runs: N = 16   384, kn = 0.01, and ξcrit = 8 Å. The dashed line is Eq. (25).

Open with DEXTER
In the text
thumbnail Fig. 6

Pressure P in [Pa] against filling factor φ with a different number of particles N. Each line shows the average of ten runs of the fixed number of particles: N = 1024,4096, and 16   384. The other parameters are Cv = 3 × 10-7, kn = 0.01, and ξcrit = 8 Å in the case of N = 1024, 4096, and Cv = 1 × 10-7, kn = 0.01, and ξcrit = 8   Å in the case of N = 16   384. The dashed line is Eq. (25).

Open with DEXTER
In the text
thumbnail Fig. 7

Pressure P in [Pa] against filling factor φ with different normal damping force. We put the same ten initial conditions varying the normal damping force with kn = 0, kn = 10-2, and kn = 101. Each line shows the result of one run. The other parameters are N = 16   384, Cv = 3 × 10-7, and ξcrit = 8 Å.

Open with DEXTER
In the text
thumbnail Fig. 8

Energy dissipation of each dissipation mechanism in [erg] against filling factor φ. The solid lines show the result in the case without the normal damping and the dashed lines in the case of kn = 0.01 and The results in the case of kn = 10 are not plotted because they are the same as those in the case of kn = 0.01 and indistinguishable. The dissipation mechanisms are normal damping, rolling, sliding, and twisting. The dissipation energy by sliding motion is less than 10-9 erg.

Open with DEXTER
In the text
thumbnail Fig. 9

Pressure P in [Pa] against filling factor φ with different critical displacement, ξcrit. We put the same ten initial conditions varying ξcrit with ξcrit = 32,16,8,4, and 2 Å, respectively. Each line shows the average of ten runs. The other parameters are N = 16   384 Cv = 3 × 10-7, and kn = 10-2.

Open with DEXTER
In the text
thumbnail Fig. 10

Energy dissipation of each dissipation mechanism in [erg] against filling factor φ. Each panel represents the case of different ξcrit, which are 32, 16, and 8 Å, corresponding to Fig. 9.

Open with DEXTER
In the text
thumbnail Fig. 11

Number of particles inside the radius rin against normalized radius rin/r0. This figure represents the fractal structure of the compressed aggregates in our simulation for various φ. We set a particle as an origin and count the number of particles inside r < rin, where r is the distance from the origin to each particle’s center. Then we count the same correlation of all particles as an origin and take their average (similar figure of Fig. 7 in the paper of Wada et al. 2008). Each line shows the result at the different time steps. The solid thick lines represents the structure of fractal dimension D = 2, and dashed lines represent D = 3 for each corresponding φ. The dotted line shows the number of particles in calculation. The region below this line corresponds to inside the periodic boundaries.

Open with DEXTER
In the text
thumbnail Fig. 12

Pressure P in [Pa] against filling factor φ. This figure is same as Fig. 4 but for the case of silicate particles (r0 = 0.6   μm).

Open with DEXTER
In the text
thumbnail Fig. 13

The filling factor φ against pressure P in [Pa]. This figure is same as Fig. 12, but plotted with a linear scale of φ and reversal of xy axis to compare with previous studies (see Fig. 4 in Seizinger et al. 2012). The dotted line is the result of numerical simulations in the high-density region (φ ≳ 0.1) in Seizinger et al. (2012) and the thin solid line is the fitting formula proposed by Güttler et al. (2009). Our results consistently connect to the previous simulations in the high-density region.

Open with DEXTER
In the text
thumbnail Fig. 14

Schematic drawing of compression of a dust aggregate consisting of a number of BCCA clusters. The left figure shows a dust aggregate consisting of many BCCA clusters and the BCCA clusters are distributed three dimensionally. Each enclosed line represents each region dominated by the BCCA clusters. The central figure is a BCCA cluster, receiving pressure from the next clusters. The BCCA cluster has a large void depicted in the central figure, and thus the void would be compressed, as expressed in the right figure. The required energy to compress the void is the energy to rotate the connection of monomers in contact. Therefore, the compression can be determined by the rolling motion of monomer connection on the connecting point of the subclusters.

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.