Tensile Strength of Porous Dust Aggregates
Comets are thought to have information about the formation process of our solar system. Recently, detailed information about comet 67P/Churyumov-Gerasimenko has been found by a spacecraft mission Rosetta. It is remarkable that its tensile strength was estimated. In this paper, we measure and formulate the tensile strength of porous dust aggregates using numerical simulations, motivated by porous dust aggregation model of planetesimal formation. We perform three-dimensional numerical simulations using a monomer interaction model with periodic boundary condition. We stretch out a dust aggregate with a various initial volume filling factor between and 0.5. We find that the tensile stress takes the maximum value at the time when the volume filling factor decreases to about a half of the initial value. The maximum stress is defined to be the tensile strength. We take an average of the results with 10 different initial shapes to smooth out the effects of initial shapes of aggregates. Finally, we numerically obtain the relation between the tensile strength and the initial volume filling factor of dust aggregates. We also use a simple semi-analytical model and successfully reproduce the numerical results, which enables us to apply to a wide parameter range and different materials. The obtained relation is consistent with previous experiments and numerical simulations about silicate dust aggregates. We estimate that the monomer radius of comet 67P has to be about 3.3–220 to reproduce its tensile strength using our model.
0000-0003-1844-5107]Misako Tatsuuma ††email@example.com \move@AU\move@AF\@affiliationDepartment of Astronomy, Graduate School of Science, The University of Tokyo, 7-3-1 Hongo, Bunkyo-ku, Tokyo 113-0033, Japan \move@AU\move@AF\@affiliationDivision of Theoretical Astronomy, National Astronomical Observatory of Japan, 2-21-1 Osawa, Mitaka, Tokyo 181-8588, Japan
0000-0003-4562-4119]Akimasa Kataoka \move@AU\move@AF\@affiliationDivision of Theoretical Astronomy, National Astronomical Observatory of Japan, 2-21-1 Osawa, Mitaka, Tokyo 181-8588, Japan
Astronomical Institute, Graduate School of Science, Tohoku University, 6-3 Aramaki, Aoba-ku, Sendai 980-8578, Japan
Planetesimal formation is one of the most important and unsolved issues of planet formation theory. In protoplanetary disks, sub-m-sized dust grains are believed to coagulate, settle to the disk midplane as they grow, and form km-sized planetesimals. There are several scenarios about the planetesimal formation such as gravitational instability (e.g., Goldreich1973), streaming instability (e.g., Youdin2005; Johansen2007; Johansen2011), and direct coagulation. In the direct coagulation scenario, dust grains grow larger by pairwise collisions. Recently, it has been proposed that dust grains become not compact but porous by pairwise collisions, and properties of these fluffy dust aggregates have been investigated theoretically and experimentally (e.g., Kozasa1992; Ossenkopf1993; Dominik1997; Blum2000; Wada2007; Wada2008; Suyama2008). The sub-m-sized constituent grains are called monomers. Finally, it is found that planetesimals form via direct coagulation (e.g., Okuzumi2012; Kataoka2013L).
In recent years, physical properties of comets have been investigated by observation and exploration. Comets are the most primitive bodies in our solar system and are thought to be leftover planetesimals. In 2014, a spacecraft Rosetta reached comet 67P/Churyumov-Gerasimenko (hereinafter 67P). This mission was the first one to orbit and land onto a comet. There are many unexpected results about 67P (e.g., Fulle2016), and especially it is remarkable that its tensile strength was estimated. The tensile strength of 67P for its surface is 3–150 Pa (Groussin2015; Basilevsky2016), while for bulk comet 10–200 Pa (Hirabayashi2016). This tensile strength depends on composition and formation process of comets, i.e., planetesimals.
There are several experimental studies about the tensile strength of dust aggregates. Blum2004 directly measured the tensile strength of dust aggregates whose volume filling factors are 0.2 and 0.54. They used dust aggregates consisted of monodisperse silica (SiO) spheres with radius. In their experiments, a mm-sized dust aggregate was attached to two plates at its top and bottom, and then the two plates were pulled apart. Blum2006 conducted the same experiments using dust aggregates which have volume filling factors of 0.23, 0.41, and 0.66. In addition to the monodisperse spherical silica monomers, they used irregularly shaped diamond monomers with a narrow size distribution and irregular silica monomers with a wide size distribution. Meisner2012 used dust aggregates consisted of quartz (crystallized SiO) monomers with a size range from 0.1 to 10 and measured the tensile strength using the Brazilian disc test (e.g., Li2013). Gundlach2018 also performed the Brazilian disc test to measure the tensile strength of dust aggregates composed of polydisperse spherical ice (HO) monomers and monodisperse spherical silica monomers. They used silica monomers whose radii are , , and to investigate the monomer radius dependence. Moreover, they succeeded in measuring the tensile strength of ice dust aggregates whose monomer radius is 2.4 on average.
On the other hand, there is only one numerical study about the tensile strength of dust aggregates. Seizinger2013 performed three-dimensional simulations to reproduce the experimental results by Blum2004 and Blum2006. They used dust aggregates whose volume filling factor ranges from 0.15 to 0.6 and monomers are silicate spheres with 0.6 radius. In their simulations, a -sized cubic aggregate was attached to two plates, which is the same as previous experiments except for the size of dust aggregates; a mm-sized dust aggregate was used in the previous experiments. The interaction between two monomers is mainly based on Dominik1997. In addition, they introduced the rolling and sliding modifiers to make numerical simulations correspond with experimental results (Seizinger2012). To avoid for monomers being peeled off the plate, they also used artificial adhesion force as “gluing effect.” Although their results correspond well with the laboratory ones, the influences of their artificial adhesion force and small aggregates should also be checked. They also obtained a fitting formula of the tensile strength as a function of the filling factor of dust aggregates. However, their formula does not include the dependence on the monomer size and material.
In this work, we numerically investigate the tensile strength of dust aggregates composed of single-sized spherical monomers. In the previous works, a dust aggregate was attached to two plates, and then they were pulled apart (Blum2004; Blum2006; Seizinger2013). The size of the used dust aggregates ranges from to mm, while planetesimals are km-sized. To unravel the planetesimal formation mechanism, it is important to investigate the tensile strength of dust aggregates whose size is larger than km. Therefore, we use the periodic boundary condition to remove effects of plates. Moreover, we perform simulations using dust aggregates whose volume filling factors are lower than those of the previous works. Then, we find a power-law relation between the tensile strength and initial volume filling factor of dust aggregates whose filling factors range from to 0.5. We will also construct a theoretical model to explain the power-law dependence. This model reproduces the dependence on all material parameters in our simulations.
This paper is constructed as follows. In Section 2, we describe settings of our simulations, which include the model of monomer interactions, initial dust aggregates, periodic boundary condition, how to measure the tensile strength without plates, and an overview of our simulations. Then, we summarize our results of fiducial runs and the investigation of parameter dependence in Section 3. There are three numerical parameters: the number of monomers, boundary velocity, and the strength of the damping force, while four physical parameters: the initial volume filling factor, monomer radius, surface energy, and the critical rolling displacement. We also find an analytical expression of the tensile strength of dust aggregates, compare our results with previous experiments (Blum2004; Blum2006; Gundlach2018) and numerical simulations (Seizinger2013), and apply the analytical expression to comet 67P in Section 4. Finally, we conclude this work and discuss future works in Section 5.
2 Simulation Settings
We perform three-dimensional numerical simulations to measure the tensile strength of dust aggregates consisting of spherical monomers. In this section, we describe settings of our simulations. First, we introduce the monomer interaction model based on Dominik1997 and Wada2007 in Section 2.1. In Section 2.2, we explain the damping force in the normal direction. The initial conditions of our simulations are statically and isotropically compressed dust aggregates investigated by Kataoka2013, which is described in Section 2.3. At the boundaries of the calculation box, we set moving periodic boundaries in the -axis direction and fixed periodic boundaries in the - and -axis directions, which is explained in Section 2.4. Thus, we can simulate one-direction stretching of dust aggregates. The details of the calculation method of tensile stress, which is the same as molecular dynamics, is summarized in Section 2.5. In Section 2.6, we describe an overview of our simulations.
2.1 Monomer Interaction Model
We calculate interactions of each connection between two monomers using the theoretical model by Dominik1997 and Wada2007. Based on the JKR theory (Johnson1971) and the following studies by Dominik1995; Dominik1996, Dominik1997 carried out two-dimensional simulations of monomer interactions. To expand into three-dimensional simulations, Wada2007 tested their recipe, and then Wada2008 conducted three-dimensional simulations of dust aggregate collisions. In the model, there are four kinds of interactions named normal (sticking and breaking), sliding, rolling, and twisting. The material parameters needed to describe the interactions are the monomer radius , material density , surface energy , Poisson’s ratio , Young’s modulus , and the critical rolling displacement . These parameters of ice and silicate are listed in Table 2.1. To compare our results with those by Seizinger2013, we set the same values for silicate.
If a rolling displacement exceeds the critical one , a monomer begins to roll inelastically. The critical rolling displacement has different values between the theoretical one ( Å, Dominik1997) and the experimental one ( Å, Heim1999). We adopt Å as a fiducial value and investigate the dependence of our results on in Section 3.3.
The rolling energy needed to rotate a monomer around its connection point by 90 is described as
where is the reduced monomer radius (Wada2007). The reduced radius of monomer radii and is defined as
Here, the reduced monomer radius is because we assume no size distribution of monomers.
In our simulations, the maximum force needed to separate two sticking monomers (breaking) is
Table 0. \Hy@raisedlink\hyper@@anchor\@currentHrefMaterial parameters of ice (Israelachvili1992; Dominik1997). The parameters of silicate are selected according to Seizinger2013.
|Monomer radius ||0.1||0.6|
|Material density [g cm]||1.0||2.65|
|Surface energy [mJ m]||100||20|
|Young’s modulus [GPa]||7||54|
|Critical rolling displacement [Å]||8||20|
2.2 Damping Force in Normal Direction
The force in the normal direction induces oscillation at each connection between two monomers. In reality, the oscillation would attenuate because of viscoelasticity or hysteresis of monomers (Greenwood2006; Tanaka2012). Therefore, we add an artificial damping force in the normal direction (Suyama2008; Paszun2008; Seizinger2012; Kataoka2013). The dependence of our results on the damping force is investigated in Section 3.2.
We describe the damping force as follows. In the case that two contacting monomers have their position vectors and , and velocities and , respectively, the contact unit vector is defined as
(Wada2007). The damping force applied to each monomer is introduced as
where is the damping coefficient, is the monomer mass, is the characteristic time, and is the relative velocity (Kataoka2013). When we calculate the damping force experienced by the monomer (), is the relative velocity of the other monomer (). We adopt as a fiducial value.
The characteristic time is given by Wada2007 as
where is the reduced Young’s modulus of monomers 1 and 2 defined as
In our simulation, the Young’s modulus and the Poisson’s ratio are uniform.
2.3 Initial Dust Aggregates
The initial dust aggregates are statically and isotropically compressed ballistic cluster-cluster aggregations (BCCAs) investigated by Kataoka2013. We set these initial conditions to simulate the planetesimal formation mechanism. The calculation boundary is treated periodically, thus we do not have to consider the aggregate radius.
2.4 One-Direction Stretching by Moving Boundaries
We set moving boundaries in the -axis direction and fixed boundaries in the - and -axis directions to measure the tensile strength of dust aggregates. The initial calculation box is a cube whose length on each side is . The length in the - and -axis directions does not change, while the length in the -axis direction increases. Therefore, the coordinates in the -, -, and -axis directions are , , and , respectively.
The velocity at the boundary in the -axis direction has to be constant and less than the effective sound speed of dust aggregates for statical stretching. We investigate the dependence on the velocity in Section 3.2. The effective sound speed of dust aggregates is described as
where and are the pressure and mean internal density of dust aggregates, respectively. Because the initial dust aggregates are statically and isotropically compressed, their pressure is given as
where is the volume filling factor of dust aggregates. Since is independent of time , the length in the -axis direction can be written as
We treat the coordinates and velocity of a monomer across a periodic boundary as follows. When a monomer passes the moving periodic boundary at , its position and velocity are converted as
On the other hand, in the case of the moving periodic boundary at , its position and velocity are converted as
At and , the coordinates of a monomer are converted similarly, but its velocity is not changed.
2.5 Tensile Stress Measurement
We calculate tensile stress in the same way as Kataoka2013 because we have no walls. This is different from Seizinger2013, who measured the tensile stress considering the force exerted on walls.
The tensile stress is calculated only in the -axis direction as follows. At first, we assume a virtual box, which is the same as the calculation box. The equation of motion of the monomer in the -axis direction is described as
where is the force exerted from the walls of the virtual box on the monomer and is the total force from other monomers on the monomer . We multiply Equation (16) by and take a long-time average with the time interval . The left-hand side of Equation (16) becomes
The first term on the right-hand side of Equation (17) becomes zero when . Here, we define the long-time average as and take a summation of all monomers of Equation (16). Then, Equation (16) can be written as
The left-hand side of Equation (18) can be defined as
which is the time-averaged kinematic energy in the -axis direction of all monomers. The first energy term on the right-hand side of Equation (18) is related to the tensile stress in the -axis direction . Since the virtual box is the same as the calculation box, we obtain
where is the volume of the calculation box. Therefore, Equation (18) gives an expression of the tensile stress as
The total force from other monomers on the monomer can be described as
where is the force from the monomer on the monomer in the -axis direction. Thus, Equation (21) can be written as
because of the relation that . Equation (23) is different from that of Kataoka2013 because we consider the tensile stress only in the -axis direction.
We take an average of the tensile stress for every 10,000 time-steps, at least. In some simulations, the tensile stress fluctuates, and thus we take a longer time average to smooth it (see Section 3.1 for details). One time-step in our simulation is s, and therefore 10,000 time-steps correspond to s.
2.6 Overview of Our Simulations
The overview of our numerical simulations is as follows. First, we randomly create a BCCA to change the initial condition. Next, we compress it statically and isotropically (Kataoka2013). The compression of a BCCA corresponds to the formation of a planetesimal. It is necessary for the BCCA to be attached to all boundaries so that we can stretch it. Then, we stop compression and define this volume filling factor as the initial one . Finally, we stretch it statically and one-dimensionally. Figure 2.6 shows the overview of our simulations.
We perform 10 simulations with different initial dust aggregates for every parameter set. First, we perform fiducial runs to investigate what occurs in our stretching simulations in Section 3.1. Then, we show that the results do not depend on any numerical parameters, such as the number of particles , boundary velocity , and the damping coefficient in Section 3.2. Finally, in Section 3.3, we investigate the dependence on physical parameters: the initial volume filling factor , monomer radius , surface energy , and the critical rolling displacement .
3.1 Fiducial Run
We measure the tensile stress of 10 runs for the fiducial parameter set. The fiducial values are , , , , , , and Å. Figure 3.1 shows three snapshots of a fiducial run. Each particle represents a 0.1 -radius ice monomer. The light gray monomers are in the calculation box with periodic boundaries, while the dark gray monomers are in the neighbor boxes. The box with white lines shows the final state of the calculation box. By stretching the dust aggregate, the chain-like structure appears.
Figure 3.1 shows the time evolution of tensile stress of 10 fiducial runs averaged for every 10,000 time-steps (left) and 200,000 time-steps (right). The volume filling factor at each time-step is calculated as
and the tensile stress is calculated according to Equation (23). We choose the number of time-steps when the rate of change of the volume filling factor does not exceed 10%. As tensile displacement increases, the volume filling factor decreases and the tensile stress increases. The maximum value of tensile stress is called the tensile strength. To calculate the tensile strength for every parameter set, we find 10 maximum values of the obtained tensile stress and take an average of them.
3.2 Numerical Parameter Dependence
To investigate the dependence on the number of particles , we plot tensile stress when , , , and in Figure 3.2(a). Changing corresponds to changing the size of the calculation box. Obviously, the tensile strength has no dependence on . Because of the smoothness of the tensile stress plot and calculation costs, we set as the fiducial value.
Figure 3.2(b) shows tensile stress when boundary velocity , , and . All boundary velocities are less than the effective sound speed of dust aggregates (Equation (10)). There is no difference among the three boundary velocities. Therefore, we can conclude that dust aggregates are stretched statically. We set as the fiducial value considering sampling rates of tensile stress and calculation costs.
Tensile stress with various damping coefficients is plotted in Figure 3.2(c). No damping force corresponds to . We change the strength of damping force from weak damping () to strong damping (). Undoubtedly, there is no dependence on the damping force in this range. We use for all the other simulations.
3.3 Physical Parameter Dependence
We measure the tensile strength of dust aggregates which have various initial volume filling factors as shown in Figure 3.3. The calculated tensile strength is proportional to from the fitting. The tensile strength can be described with the initial volume filling factor as
To investigate the dependence of tensile strength on the monomer radius, we perform simulations in the case of ice monomers whose radii are 0.3 and 0.9 . Figure 3.3 shows the summary of the monomer radius dependence. The plotted dashed lines are based on Equation (31), which is an analytical expression of tensile strength (see Section 4.1). It is confirmed that the tensile strength is in inverse proportion to the monomer radius.
To clarify the surface energy dependence, we calculate the tensile strength when and and plot it in Figure 3.3. Other parameters are the same as the fiducial values. The dashed lines represent Equation (31) (see Section 4.1). Obviously, the tensile strength is in proportion to the surface energy.
Finally, we investigate the dependence of tensile stress on the critical rolling displacement in Figure 3.3.
The critical rolling displacement is changed from
citep[e.g.,]Dominik1997 to (e.g., Heim1999). Tensile stress has a marginal dependence on because the main mechanism of displacement is rolling (see Section 4.1). On the other hand, tensile strength, which is the maximum value of tensile stress, has no difference. We can conclude that the tensile strength is almost the same even if the critical rolling displacement has uncertainty. Therefore, we fix in our simulations.
Now, we discuss the obtained physical parameter dependence of the tensile strength of dust aggregates (Section 3.3) and apply our results to previous studies of experiments, numerical simulations, and comet 67P. In Section 4.1, we find an analytical expression of the tensile strength using material parameters: the initial volume filling factor, monomer radius, and the surface energy. Then, we compare our results with previous experiments and numerical simulations about silicate dust aggregates (Blum2004; Blum2006; Seizinger2013; Gundlach2018) in Section 4.2. Finally, we apply our interpretation to comet 67P in Section 4.3.
4.1 Semi-Analytical Model of Tensile Strength
The relationship between and can be derived by considering the maximum force needed to separate two sticking monomers and the radius of a dust aggregate . When the tensile stress has a maximum value, the force is applied on a connection between two monomers of a dust aggregate. This means that
The radius of a dust aggregate is given as
where and are the fractal dimension and the number of monomers of a dust aggregate, respectively. The initial volume filling factor is described as
and then, the radius of a dust aggregate is obtained as
From Equation (26), the tensile strength can be written as
where is a constant. The fractal dimension of BCCAs is (Mukai1992; Okuzumi2009dustagg).
In the case of ice monomer whose radius is 0.1 , we find .
We confirm Equation (31) from the perspective of energy dissipation. All energy dissipations, which are caused by the normal, sliding, rolling, twisting, and damping force in the normal direction, are plotted in Figure 4.1. The curves in Figure 4.1 run time-wise from right to left and arise during the stretching of a dust aggregate. The main energy dissipation mechanism is the rolling, which corresponds to the -dependence of tensile stress (Section 3.3). The energy dissipation by the normal arises when the tensile stress has a maximum value. This energy dissipation is caused by connection breaking between two contacting monomers. For this reason, tensile strength is determined by the connection breaking, i.e. .
To confirm that on a small scale of a dust aggregate in our simulations, we calculate the number of monomers inside the radius for five snapshots of a fiducial run and plot it in Figure 4.1. We take the snapshots during the continuous strain of the dust aggregate. The parameters of the fiducial run are , , , , , , and Å. The method to count the number of monomers is as follows. At first, we set a monomer in the calculation box as the center and count including monomers outside the periodic boundaries. Next, we take an average of for all monomers in the calculation box.
4.2 Comparison with Previous Studies
To compare our results with previous experiments and numerical simulations, we perform simulations with the parameters of silicate listed in Table 2.1 and summarize the results in Figure 4.2. Both results by experiments (Blum2004; Blum2006; Gundlach2018) and simulations (this work; Seizinger2013) correspond very well. This means that there is little influence of artificial adhesion force introduced by Seizinger2013 and small aggregates whose sizes are from to mm. The right panel of Figure 4.2 is derived by Equation (31), i.e., .
Gundlach2018 measured the tensile strength of ice aggregates with the monomer radius of , which is shown in Figure 4.2. We do not know which monomer radius determines the tensile strength of dust aggregates when the monomer radius has a size distribution. Therefore, we plot dashed lines derived by Equation (31) when , , and . The experimental value is lower than the theoretical lines. From their low value of the tensile strength, Gundlach2018 inferred that the specific surface energy of ice, , has a very value of at low temperatures ( K). However, Gundlach2015 also estimated as from the constant sticking threshold velocity of for K in their ice impact experiments. The reason for low tensile strength in Gundlach2018 is thus unknown. We also confirm that the experimental value is higher than the compressive strength theoretically expected by Kataoka2013 (Equation (9)).
4.3 Application to Comet 67P
We can estimate the monomer radius of comet 67P using Equation (31) as follows. From the exploration of 67P, it was found that the micro-porosity of 67P is 0.75–0.85 (Kofman2015), while the surface porosity is 0.87 (Fornasier2015). In other words, the volume filling factor of 67P is about 0.13–0.25. The averaged nucleus bulk density of 67P was estimated to be 0.533 (Patzold2016). In consideration of the high dust-to-water mass ratio of 67P (e.g., Fulle2016), its main material is not HO ice, but silicate. Assuming that the bulk density is 0.533 and the material density is 2.65 (Table 2.1), we can estimate that the volume filling factor of 67P is about 0.20, which is consistent with the values of 0.13–0.25. Substituting –200 Pa (see Section 1), , and in Equation (33), we find that the monomer radius of 67P has to be 3.3–220 . To explain the low tensile strength of 67P using only our model, we have to consider a larger radius than that of interstellar materials.
Another idea to decrease the tensile strength is assuming an aggregate of aggregates (e.g., Blum2014; Blum2017) instead of a simple aggregate of monomers as discussed in this paper. The aggregate-of-aggregate model may explain the low strength with small monomers, but this is beyond the scope of this paper.
We investigated the tensile strength of porous dust aggregates whose initial volume filling factors are from to 0.5. We performed three-dimensional numerical simulations with periodic boundary condition to measure the tensile stress of dust aggregates. The monomer interaction model is based on Dominik1997 and Wada2007. The initial dust aggregates are statically and isotropically compressed BCCAs investigated by Kataoka2013. At boundaries of the calculation box, we set moving periodic boundaries in the -axis direction and fixed periodic boundaries in the - and -axis directions. The calculation method of the tensile stress is the same way as molecular dynamics. In our simulations, we created a BCCA at first, compressed it three-dimensionally, and then stretched it one-dimensionally. For every parameter set, we conducted 10 stretching-simulations with different initial dust aggregates, found 10 maximum values of the obtained tensile stress, and took an average of them, which is called the tensile strength. Our main findings of the tensile strength of porous dust aggregates are as follows.
As a result of numerical simulations, we found that the tensile strength can be written as
where is the maximum force needed to separate two sticking monomers, is the initial volume filling factor, is the monomer radius, and is the surface energy.
It is confirmed that the energy dissipation during a stretching simulation support the dependence in Equation (33). The energy dissipation caused by monomer-connection breaking arises when the tensile stress has a maximum value. Also, the dependence on the initial volume filling factor corresponds to the fractal dimension of BCCAs, which is about 1.9 (Mukai1992; Okuzumi2009dustagg).
Equation (33) is consistent with the previous experimental (Blum2004; Blum2006; Gundlach2018) and numerical (Seizinger2013) studies of silicate dust aggregates, while it is inconsistent with the previous experimental study of ice dust aggregates (Gundlach2018).
We estimated that the monomer radius of comet 67P has to be 3.3–220 using Equation (33). Assuming that the main material of 67P is silicate and its volume filling factor is about 0.20, we obtained the monomer radius to reproduce its tensile strength of 3–200 Pa.
From the point of view of planet formation, the conclusion that the monomer radius of 67P has to be 3.3–220 is inconsistent with the typical radius of dust monomers in the interstellar medium: sub-. To reduce the monomer radius of 67P, other mechanisms, such as sintering (e.g., Sirono2017), to decrease the tensile strength of dust aggregates are needed.
We thank Satoshi Okuzumi for fruitful discussions.