Tensile Strength of Porous Dust Aggregates

TL;DR

This study numerically investigates the tensile strength of porous dust aggregates, relating it to initial volume filling factors, and applies the findings to estimate monomer sizes in comet 67P, enhancing understanding of cometary material properties.

Contribution

The paper introduces a numerical and semi-analytical model for tensile strength of dust aggregates, linking it to initial volume filling factors and applying it to cometary data.

Findings

Maximum tensile stress occurs when volume filling factor halves.

Derived relation between tensile strength and volume filling factor.

Estimated monomer size range for comet 67P to match observed tensile strength.

Abstract

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 $10^{-2}$ 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 $\mathrm{\mu m}$ to reproduce its tensile strength using our model.

Tables

Table 1. Table 1: Material parameters of ice (Israelachvili, 1992 ; Dominik & Tielens, 1997 ) . The parameters of silicate are selected according to Seizinger et al. ( 2013 ) .

Material	Ice	Silicate
Monomer radius $r_0$ [$\mu m$]	0.1	0.6
Material density ${\rho }_0$ [g cm-3]	1.0	2.65
Surface energy $\gamma$ [mJ m-2]	100	20
Poisson’s ratio $\nu$	0.25	0.17
Young’s modulus $E$ [GPa]	7	54
Critical rolling displacement ${\xi }_{\mathrm{crit}}$ [Å]	8	20

Equations

$$E_{\mathrm{roll}}$$

$$\frac{1}{R}=\frac{1}{r_1}+\frac{1}{r_2}.$$

$$F_{\mathrm{c}}$$

$${\mathbf{n}}_{\mathrm{c}}=\frac{{\mathbf{x}}_1-{\mathbf{x}}_2}{|{\mathbf{x}}_1-{\mathbf{x}}_2|}$$

$${\mathbf{F}}_{\mathrm{damp}}=-k_{\mathrm{n}}\frac{m_0}{t_{\mathrm{c}}}({\mathbf{n}}_{\mathrm{c}}\cdot {\mathbf{v}}_{\mathrm{r}}){\mathbf{n}}_{\mathrm{c}},$$

$$t_{\mathrm{c}}=0.95\left(\frac{r_0^{7/6}{\rho }_0^{1/2}}{{\gamma }^{1/6}{E^{\ast }}^{1/3}}\right),$$

$$\frac{1}{E^{\ast }}=\frac{1-{\nu }_1^2}{E_1}+\frac{1-{\nu }_2^2}{E_2}.$$

$$c_{\mathrm{s},\mathrm{eff}}\sim \sqrt{\frac{P}{\rho }},$$

$$P=\frac{E_{\mathrm{roll}}}{r_0^3}{\left(\frac{\rho }{{\rho }_0}\right)}^3$$

$$c_{\mathrm{s},\mathrm{eff}}$$

$$L_x=L_0+2v_{\mathrm{b}}t.$$

$$x$$

$$v_x$$

$$x$$

$$v_x$$

$$m_0\frac{{\mathrm{d}}^2{x}_i}{\mathrm{d}{t}^2}=W_{x,i}+F_{x,i},$$

$$\frac{m_0}{\tau }{\int }_0^{\tau }{x}_i\frac{{\mathrm{d}}^2{x}_i}{\mathrm{d}{t}^2}dt=\frac{m_0}{\tau }{\left[x_i\frac{\mathrm{d}{x}_i}{\mathrm{d}t}\right]}_0^{\tau }-\frac{m_0}{\tau }{\int }_0^{\tau }\frac{\mathrm{d}{x}_i}{\mathrm{d}t}\frac{\mathrm{d}{x}_i}{\mathrm{d}t}dt.$$

$${⟨\sum _{i=1}^N\frac{m_0}{2}{\left(\frac{\mathrm{d}{x}_i}{\mathrm{d}t}\right)}^2⟩}_t=-\frac{1}{2}{⟨\sum _{i=1}^N{x}_i{W}_{x,i}⟩}_t-\frac{1}{2}{⟨\sum _{i=1}^N{x}_i{F}_{x,i}⟩}_t.$$

$$K_x={⟨\sum _{i=1}^N\frac{m_0}{2}{\left(\frac{\mathrm{d}{x}_i}{\mathrm{d}t}\right)}^2⟩}_t,$$

$${⟨\sum _{i=1}^N{x}_i{W}_{x,i}⟩}_t=L_x{P}_x{L}_0^2=P_{x}V,$$

$$P_x=-\frac{2K_x}{V}-\frac{1}{V}{⟨\sum _{i=1}^N{x}_i{F}_{x,i}⟩}_t.$$

$$F_{x,i}=\sum _{j\ne i}{f}_{x,i,j},$$

$$P_x=-\frac{2K_x}{V}-\frac{1}{V}{⟨\sum _{i<j}(x_i-x_j)f_{x,i,j}⟩}_t$$

$$\varphi =\frac{(4/3)\pi r_0^{3}N}{V}$$

$$P_{x,\max }=P_0{\varphi }_{\mathrm{init}}^{1.8},$$

$$P_{x,\max }\propto \frac{F_{\mathrm{c}}}{r_{\mathrm{agg}}^2}.$$

$$r_{\mathrm{agg}}\propto N_{\mathrm{agg}}^{1/D}{r}_0,$$

$${\varphi }_{\mathrm{init}}=N_{\mathrm{agg}}{\left(\frac{r_0}{r_{\mathrm{agg}}}\right)}^3,$$

$$r_{\mathrm{agg}}\propto r_0{\varphi }_{\mathrm{init}}^{-1/(3-D)}.$$

$$P_{x,\max }\sim C\frac{F_{\mathrm{c}}}{r_0^2}{\varphi }_{\mathrm{init}}^{2/(3-D)},$$

$$P_{x,\max }$$

Full text

planets and satellites: formation — protoplanetary disks — methods: numerical — methods: analytical

††journal: ApJ

1 Introduction

Planetesimal formation is one of the most important and unsolved issues of planet formation theory. In protoplanetary disks, sub-$\mathrm{\mu}$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., Goldreich & Ward, 1973), streaming instability (e.g., Youdin & Goodman, 2005; Johansen et al., 2007, 2011), 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., Kozasa et al., 1992; Ossenkopf, 1993; Dominik & Tielens, 1997; Blum & Wurm, 2000; Wada et al., 2007, 2008; Suyama et al., 2008). The sub-$\mathrm{\mu}$m-sized constituent grains are called monomers. Finally, it is found that planetesimals form via direct coagulation (e.g., Okuzumi et al., 2012; Kataoka et al., 2013a).

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., Fulle et al., 2016), and especially it is remarkable that its tensile strength was estimated. The tensile strength of 67P for its surface is 3–150 Pa (Groussin et al., 2015; Basilevsky et al., 2016), while for bulk comet 10–200 Pa (Hirabayashi et al., 2016). 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. Blum & Schräpler (2004) 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 (SiO2) spheres with $0.76\mathrm{\ \mu m}$ 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. Blum et al. (2006) 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. Meisner et al. (2012) used dust aggregates consisted of quartz (crystallized SiO2) monomers with a size range from 0.1 $\mathrm{\mu m}$ to 10 $\mathrm{\mu m}$ and measured the tensile strength using the Brazilian disc test (e.g., Li & Wong, 2013). Gundlach et al. (2018) also performed the Brazilian disc test to measure the tensile strength of dust aggregates composed of polydisperse spherical ice (H2O) monomers and monodisperse spherical silica monomers. They used silica monomers whose radii are $0.15\mathrm{\ \mu m}$, $0.50\mathrm{\ \mu m}$, and $0.75\mathrm{\ \mu m}$ to investigate the monomer radius dependence. Moreover, they succeeded in measuring the tensile strength of ice dust aggregates whose monomer radius is 2.4 $\mathrm{\mu m}$ on average.

On the other hand, there is only one numerical study about the tensile strength of dust aggregates. Seizinger et al. (2013) performed three-dimensional simulations to reproduce the experimental results by Blum & Schräpler (2004) and Blum et al. (2006). They used dust aggregates whose volume filling factor ranges from 0.15 to 0.6 and monomers are silicate spheres with 0.6 $\mathrm{\mu m}$ radius. In their simulations, a $\mathrm{\mu m}$-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 Dominik & Tielens (1997). In addition, they introduced the rolling and sliding modifiers to make numerical simulations correspond with experimental results (Seizinger et al., 2012). 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 (Blum & Schräpler, 2004; Blum et al., 2006; Seizinger et al., 2013). The size of the used dust aggregates ranges from $\mathrm{\mu m}$ 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 $10^{-2}$ 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 (Blum & Schräpler, 2004; Blum et al., 2006; Gundlach et al., 2018) and numerical simulations (Seizinger et al., 2013), 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 Dominik & Tielens (1997) and Wada et al. (2007) 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 Kataoka et al. (2013b), which is described in Section 2.3. At the boundaries of the calculation box, we set moving periodic boundaries in the $x$-axis direction and fixed periodic boundaries in the $y$- and $z$-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 Dominik & Tielens (1997) and Wada et al. (2007). Based on the JKR theory (Johnson et al., 1971) and the following studies by Dominik & Tielens (1995, 1996), Dominik & Tielens (1997) carried out two-dimensional simulations of monomer interactions. To expand into three-dimensional simulations, Wada et al. (2007) tested their recipe, and then Wada et al. (2008) 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 $r_{0}$, material density $\rho_{0}$, surface energy $\gamma$, Poisson’s ratio $\nu$, Young’s modulus $E$, and the critical rolling displacement $\xi_{\mathrm{crit}}$. These parameters of ice and silicate are listed in Table 1. To compare our results with those by Seizinger et al. (2013), we set the same values for silicate.

If a rolling displacement exceeds the critical one $\xi_{\mathrm{crit}}$, a monomer begins to roll inelastically. The critical rolling displacement has different values between the theoretical one ($\xi_{\mathrm{crit}}=2$ Å, Dominik & Tielens, 1997) and the experimental one ($\xi_{\mathrm{crit}}=32$ Å, Heim et al., 1999). We adopt $\xi_{\mathrm{crit}}=8$ Å as a fiducial value and investigate the dependence of our results on $\xi_{\mathrm{crit}}$ in Section 3.3.

The rolling energy $E_{\mathrm{roll}}$ needed to rotate a monomer around its connection point by 90∘ is described as

[TABLE]

where $R$ is the reduced monomer radius (Wada et al., 2007). The reduced radius $R$ of monomer radii $r_{1}$ and $r_{2}$ is defined as

[TABLE]

Here, the reduced monomer radius is $R=r_{0}/2$ because we assume no size distribution of monomers.

In our simulations, the maximum force needed to separate two sticking monomers (breaking) is

[TABLE]

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 (Greenwood & Johnson, 2006; Tanaka et al., 2012). Therefore, we add an artificial damping force in the normal direction (Suyama et al., 2008; Paszun & Dominik, 2008; Seizinger et al., 2012; Kataoka et al., 2013b). 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 $\bm{x}_{1}$ and $\bm{x}_{2}$, and velocities $\bm{v}_{1}$ and $\bm{v}_{2}$, respectively, the contact unit vector $\bm{n}_{\mathrm{c}}$ is defined as

[TABLE]

(Wada et al., 2007). The damping force applied to each monomer is introduced as

[TABLE]

where $k_{\mathrm{n}}$ is the damping coefficient, $m_{0}$ is the monomer mass, $t_{\mathrm{c}}$ is the characteristic time, and $\bm{v}_{\mathrm{r}}$ is the relative velocity (Kataoka et al., 2013b). When we calculate the damping force experienced by the monomer ($\bm{x}_{1},\bm{v}_{1}$), $\bm{v}_{\mathrm{r}}=\bm{v}_{2}-\bm{v}_{1}$ is the relative velocity of the other monomer ($\bm{x}_{2},\bm{v}_{2}$). We adopt $k_{\mathrm{n}}=0.01$ as a fiducial value.

The characteristic time is given by Wada et al. (2007) as

[TABLE]

where $E^{\ast}$ is the reduced Young’s modulus of monomers 1 and 2 defined as

[TABLE]

In our simulation, the Young’s modulus $E_{1}=E_{2}=E$ and the Poisson’s ratio $\nu_{1}=\nu_{2}=\nu$ are uniform.

2.3 Initial Dust Aggregates

The initial dust aggregates are statically and isotropically compressed ballistic cluster-cluster aggregations (BCCAs) investigated by Kataoka et al. (2013b). 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 $x$-axis direction and fixed boundaries in the $y$- and $z$-axis directions to measure the tensile strength of dust aggregates. The initial calculation box is a cube whose length on each side is $L_{0}$. The length in the $y$- and $z$-axis directions does not change, while the length in the $x$-axis direction $L_{x}$ increases. Therefore, the coordinates in the $x$-, $y$-, and $z$-axis directions are $-L_{x}/2<x<L_{x}/2$, $-L_{0}/2<y<L_{0}/2$, and $-L_{0}/2<z<L_{0}/2$, respectively.

The velocity at the boundary in the $x$-axis direction $v_{\mathrm{b}}>0$ 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 $c_{\mathrm{s,eff}}$ is described as

[TABLE]

where $P$ and $\rho$ 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

[TABLE]

(Kataoka et al., 2013b). From Equation (8) and (9), we can obtain the effective sound speed of dust aggregates as

[TABLE]

where $\phi=\rho/\rho_{0}$ is the volume filling factor of dust aggregates. Since $v_{\mathrm{b}}$ is independent of time $t$, the length in the $x$-axis direction $L_{x}$ can be written as

[TABLE]

We treat the coordinates and velocity of a monomer across a periodic boundary as follows. When a monomer passes the moving periodic boundary at $x=L_{x}/2$, its position $x$ and velocity $v_{x}$ are converted as

[TABLE]

On the other hand, in the case of the moving periodic boundary at $x=-L_{x}/2$, its position and velocity are converted as

[TABLE]

At $y=\pm L_{0}/2$ and $z=\pm L_{0}/2$, 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 Kataoka et al. (2013b) because we have no walls. This is different from Seizinger et al. (2013), who measured the tensile stress considering the force exerted on walls.

The tensile stress is calculated only in the $x$-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 $i$ in the $x$-axis direction is described as

[TABLE]

where $W_{x,i}$ is the force exerted from the walls of the virtual box on the monomer $i$ and $F_{x,i}$ is the total force from other monomers on the monomer $i$. We multiply Equation (16) by $x_{i}$ and take a long-time average with the time interval $\tau$. The left-hand side of Equation (16) becomes

[TABLE]

The first term on the right-hand side of Equation (17) becomes zero when $\tau\to\infty$. Here, we define the long-time average as $\langle\rangle_{t}$ and take a summation of all monomers of Equation (16). Then, Equation (16) can be written as

[TABLE]

The left-hand side of Equation (18) can be defined as

[TABLE]

which is the time-averaged kinematic energy in the $x$-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 $x$-axis direction $P_{x}$. Since the virtual box is the same as the calculation box, we obtain

[TABLE]

where $V=L_{x}L_{0}^{2}$ is the volume of the calculation box. Therefore, Equation (18) gives an expression of the tensile stress $P_{x}$ as

[TABLE]

The total force from other monomers on the monomer $i$ can be described as

[TABLE]

where $f_{x,i,j}$ is the force from the monomer $j$ on the monomer $i$ in the $x$-axis direction. Thus, Equation (21) can be written as

[TABLE]

because of the relation that $f_{x,i,j}=-f_{x,j,i}$. Equation (23) is different from that of Kataoka et al. (2013b) because we consider the tensile stress only in the $x$-axis direction.

We take an average of the tensile stress $P_{x}$ 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 $0.7t_{\mathrm{c}}=1.9\times 10^{-11}$ s, and therefore 10,000 time-steps correspond to $1.9\times 10^{-7}$ 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 (Kataoka et al., 2013b). 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 $\phi_{\mathrm{init}}$. Finally, we stretch it statically and one-dimensionally. Figure 1 shows the overview of our simulations.

3 Results

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 $N$, boundary velocity $v_{\mathrm{b}}$, and the damping coefficient $k_{\mathrm{n}}$ in Section 3.2. Finally, in Section 3.3, we investigate the dependence on physical parameters: the initial volume filling factor $\phi_{\mathrm{init}}$, monomer radius $r_{0}$, surface energy $\gamma$, and the critical rolling displacement $\xi_{\mathrm{crit}}$.

3.1 Fiducial Run

We measure the tensile stress of 10 runs for the fiducial parameter set. The fiducial values are $N=16384$, $v_{\mathrm{b}}=10\mathrm{\ cm\ s^{-1}}$, $k_{\mathrm{n}}=0.01$, $\phi_{\mathrm{init}}=0.1$, $r_{0}=0.1\mathrm{\ \mu m}$, $\gamma=100\mathrm{\ mJ\ m^{-2}}$, and $\xi_{\mathrm{crit}}=8$ Å. Figure 2 shows three snapshots of a fiducial run. Each particle represents a 0.1 $\mathrm{\mu m}$-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 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

[TABLE]

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 $\phi$ decreases and the tensile stress $P_{x}$ 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 $N$, we plot tensile stress when $N=2^{10}=1024$, $N=2^{12}=4096$, $N=2^{14}=16384$, and $N=2^{16}=65536$ in Figure 4(a). Changing $N$ corresponds to changing the size of the calculation box. Obviously, the tensile strength has no dependence on $N$. Because of the smoothness of the tensile stress plot and calculation costs, we set $N=16384$ as the fiducial value.

Figure 4(b) shows tensile stress when boundary velocity $v_{\mathrm{b}}=1\mathrm{\ cm\ s^{-1}}$, $10\mathrm{\ cm\ s^{-1}}$, and $100\mathrm{\ cm\ s^{-1}}$. All boundary velocities are less than the effective sound speed of dust aggregates $c_{\mathrm{s,eff}}$ (Equation (10)). There is no difference among the three boundary velocities. Therefore, we can conclude that dust aggregates are stretched statically. We set $v_{\mathrm{b}}=10\mathrm{\ cm\ s^{-1}}$ as the fiducial value considering sampling rates of tensile stress and calculation costs.

Tensile stress with various damping coefficients is plotted in Figure 4(c). No damping force corresponds to $k_{\mathrm{n}}=0$. We change the strength of damping force from weak damping ($k_{\mathrm{n}}=0.01$) to strong damping ($k_{\mathrm{n}}=1$). Undoubtedly, there is no dependence on the damping force in this range. We use $k_{\mathrm{n}}=0.01$ 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 5. The calculated tensile strength is proportional to $\phi_{\mathrm{init}}^{1.8}$ from the fitting. The tensile strength $P_{x,\mathrm{max}}$ can be described with the initial volume filling factor $\phi_{\mathrm{init}}$ as

[TABLE]

where $P_{0}\sim 6\times 10^{5}$ Pa in this case. The analytical interpretation of Equation (25) is discussed in Section 4.1.

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 $\mathrm{\mu m}$ and 0.9 $\mathrm{\mu m}$. Figure 6 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 $\gamma=50\mathrm{\ mJ\ m^{-2}}$ and $25\mathrm{\ mJ\ m^{-2}}$ and plot it in Figure 7. 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 $\xi_{\mathrm{crit}}$ in Figure 8. The critical rolling displacement is changed from $\xi_{\mathrm{crit}}=2\textrm{\ \AA}$

citep[e.g.,][]Dominik1997 to $\xi_{\mathrm{crit}}=32\textrm{\ \AA}$ (e.g., Heim et al., 1999). Tensile stress has a marginal dependence on $\xi_{\mathrm{crit}}$ 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 $\xi_{\mathrm{crit}}=8\textrm{\ \AA}$ in our simulations.

4 Discussions

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 (Blum & Schräpler, 2004; Blum et al., 2006; Seizinger et al., 2013; Gundlach et al., 2018) 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 $P_{x,\mathrm{max}}$ and $\phi_{\mathrm{init}}$ can be derived by considering the maximum force needed to separate two sticking monomers $F_{\mathrm{c}}$ and the radius of a dust aggregate $r_{\mathrm{agg}}$. When the tensile stress has a maximum value, the force $F_{\mathrm{c}}$ is applied on a connection between two monomers of a dust aggregate. This means that

[TABLE]

The radius of a dust aggregate is given as

[TABLE]

where $D$ and $N_{\mathrm{agg}}$ are the fractal dimension and the number of monomers of a dust aggregate, respectively. The initial volume filling factor is described as

[TABLE]

and then, the radius of a dust aggregate is obtained as

[TABLE]

From Equation (26), the tensile strength can be written as

[TABLE]

where $C$ is a constant. The fractal dimension $D$ of BCCAs is $\sim 1.9$ (Mukai et al., 1992; Okuzumi et al., 2009).

Using the fitting result of Equations (3) and (25), we obtain

[TABLE]

In the case of ice monomer whose radius is 0.1 $\mathrm{\mu m}$, we find $C=0.12\pm 0.01$.

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 9. The curves in Figure 9 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 $\xi_{\mathrm{crit}}$-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. $F_{\mathrm{c}}$.

To confirm that $D\sim 1.9$ on a small scale of a dust aggregate in our simulations, we calculate the number of monomers inside the radius $r_{\mathrm{in}}$ for five snapshots of a fiducial run and plot it in Figure 10. We take the snapshots during the continuous strain of the dust aggregate. The parameters of the fiducial run are $N=16384$, $v_{\mathrm{b}}=10\mathrm{\ cm\ s^{-1}}$, $k_{\mathrm{n}}=0.01$, $\phi_{\mathrm{init}}=0.1$, $r_{0}=0.1\mathrm{\ \mu m}$, $\gamma=100\mathrm{\ mJ\ m^{-2}}$, and $\xi_{\mathrm{crit}}=8$ Å. The method to count the number of monomers $N(r<r_{\mathrm{in}})$ is as follows. At first, we set a monomer in the calculation box as the center and count $N(r<r_{\mathrm{in}})$ including monomers outside the periodic boundaries. Next, we take an average of $N(r<r_{\mathrm{in}})$ for all monomers in the calculation box.

Also, we plot $N(r<r_{\mathrm{in}})$ as a function of $r_{\mathrm{in}}/r_{0}$ when $D=2$ and $D=3$ in Figure 10. This relationship is derived by Equation (27) as

[TABLE]

On a small scale that $N(r<r_{\mathrm{in}})\lesssim 20$, all snapshot results mostly correspond to the relationship when $D=2$. When the scale becomes large enough, dust aggregates have a fractal dimension of three.

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 1 and summarize the results in Figure 11. Both results by experiments (Blum & Schräpler, 2004; Blum et al., 2006; Gundlach et al., 2018) and simulations (this work; Seizinger et al., 2013) correspond very well. This means that there is little influence of artificial adhesion force introduced by Seizinger et al. (2013) and small aggregates whose sizes are from $\mathrm{\mu m}$ to mm. The right panel of Figure 11 is derived by Equation (31), i.e., $P_{x,\mathrm{max}}\propto r_{0}^{-1}$.

Gundlach et al. (2018) measured the tensile strength of ice aggregates with the monomer radius of $2.38\pm 1.11\mathrm{\ \mu m}$, which is shown in Figure 12. 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 $r_{0}=1.27\mathrm{\ \mu m}$, $r_{0}=2.38\mathrm{\ \mu m}$, and $r_{0}=3.49\mathrm{\ \mu m}$. The experimental value is lower than the theoretical lines. From their low value of the tensile strength, Gundlach et al. (2018) inferred that the specific surface energy of ice, $\gamma_{\mathrm{ice}}$, has a very value of $=0.02\mathrm{\ J\ m^{-2}}$ at low temperatures ($\lesssim 150$ K). However, Gundlach & Blum (2015) also estimated as $\gamma_{\mathrm{ice}}=0.19\mathrm{\ J\ m^{-2}}$ from the constant sticking threshold velocity of $\sim 10\mathrm{\ m\ s^{-1}}$ for $T<210$ K in their ice impact experiments. The reason for low tensile strength in Gundlach et al. (2018) is thus unknown. We also confirm that the experimental value is higher than the compressive strength theoretically expected by Kataoka et al. (2013b) (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 (Kofman et al., 2015), while the surface porosity is 0.87 (Fornasier et al., 2015). 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 $\mathrm{g\ cm^{-3}}$ (Pätzold et al., 2016). In consideration of the high dust-to-water mass ratio of 67P (e.g., Fulle et al., 2016), its main material is not H2O ice, but silicate. Assuming that the bulk density is 0.533 $\mathrm{g\ cm^{-3}}$ and the material density is 2.65 $\mathrm{g\ cm^{-3}}$ (Table 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 $P_{x,\mathrm{max}}=3$–200 Pa (see Section 1), $\gamma=20\mathrm{\ mJ\ m^{-2}}$, and $\phi_{\mathrm{init}}=0.20$ in Equation (33), we find that the monomer radius of 67P has to be 3.3–220 $\mathrm{\mu m}$. 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., Blum et al., 2014, 2017) 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.

5 Conclusions

We investigated the tensile strength of porous dust aggregates whose initial volume filling factors $\phi_{\mathrm{init}}$ are from $10^{-2}$ 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 Dominik & Tielens (1997) and Wada et al. (2007). The initial dust aggregates are statically and isotropically compressed BCCAs investigated by Kataoka et al. (2013b). At boundaries of the calculation box, we set moving periodic boundaries in the $x$-axis direction and fixed periodic boundaries in the $y$- and $z$-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 $P_{x,\mathrm{max}}$ can be written as

[TABLE]

where $F_{\mathrm{c}}$ is the maximum force needed to separate two sticking monomers, $\phi_{\mathrm{init}}$ is the initial volume filling factor, $r_{0}$ is the monomer radius, and $\gamma$ is the surface energy.

• •

We analytically confirmed the dependence in Equation (33). It is found that the tensile strength is determined by monomer-connection breaking. This is consistent with the fact that $P_{x,\mathrm{max}}$ is proportional to $F_{\mathrm{c}}$ as shown in Equation (33).

• •

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 (Mukai et al., 1992; Okuzumi et al., 2009).

• •

Equation (33) is consistent with the previous experimental (Blum & Schräpler, 2004; Blum et al., 2006; Gundlach et al., 2018) and numerical (Seizinger et al., 2013) studies of silicate dust aggregates, while it is inconsistent with the previous experimental study of ice dust aggregates (Gundlach et al., 2018).

• •

We estimated that the monomer radius of comet 67P has to be 3.3–220 $\mathrm{\mu m}$ 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 $\mathrm{\mu m}$ is inconsistent with the typical radius of dust monomers in the interstellar medium: sub-$\mathrm{\mu m}$. To reduce the monomer radius of 67P, other mechanisms, such as sintering (e.g., Sirono & Ueno, 2017), to decrease the tensile strength of dust aggregates are needed.

We thank Satoshi Okuzumi for fruitful discussions.