© ESO, 2013

1. Introduction

The formation of planetesimals, which are kilometer-sized objects that are massive enough for gravity to come into play, constitutes a key step in the core accretion scenario for planet formation proposed by Pollack et al. (1996). However, the earlier growth from mm to km sized bodies has not yet been fully understood. Understanding the interplay between porosity, impact velocity, and the size and structure of colliding aggregates in the meter-size regime is crucial to unraveling the process of planetesimal formation. Obviously, this size regime renders laboratory experiments impossible. Thus, for years to come astrophysicists will have to rely on computer simulations to obtain the necessary insight into this complex process.

Molecular dynamics simulations featuring detailed micro-mechanical interactions have been employed to study collisions of sub-mm sized dust and ice aggregates (e.g., Dominik & Tielens 1997; Wada et al. 2007; Ringl et al. 2012). Owing to the high computational demand, a different approach is necessary for the mm to meter size regime. In this regime, smoothed particle hydrodynamics (SPH) simulations have often been utilized to model preplanetesimal collisions (e.g., Sirono 2004; Schäfer et al. 2007; Geretshauser et al. 2011). SPH constitutes a continuum approach that is able to simulate the collisional behavior of macroscopic aggregates including physical processes such as compaction or fragmentation. As a continuum approach, SPH requires various material parameters such as the compressive, tensile, and shear strength, so that a proper calibration is necessary (Geretshauser et al. 2010). Typically, the calibration process is based on comparison with results from laboratory experiments (Güttler et al. 2009). However, not all material parameters have been obtained in this way. For instance, the shear strength has only been based on theoretical models and estimations so far.

Only a few laboratory experiments have been performed to investigate the mechanical properties of porous dust aggregates. Blum & Schräpler (2004) measured the tensile strength of highly porous dust aggregates generated by random ballistic deposition (RBD). Since they used monodisperse, spherical silica grains, their experiments are comparable to our simulations. It has been shown both theoretically (Bertini et al. 2009) and experimentally (Blum et al. 2006) that the mechanical properties depend on the shape and size distribution of the grains. Recently, Meisner et al. (2012) presented results from various experiments on the mechanical properties of irregularly shaped quartz aggregates.

Determining material parameters directly from molecular dynamics simulations of porous dust aggregates constitutes a tempting alternative. Paszun & Dominik (2008) presented the first attempt to obtain the compressive strength from ab initio simulations. A few years later, Seizinger et al. (2012) studied the compressive strength in greater detail and especially revealed the differences between static and dynamic compaction processes.

The aim of the present work is to extend this approach to determine the tensile and shear strength of porous aggregates. Together with the compressive strength we then can provide continuum simulations (such as SPH) with a complete parameter set describing the transition from elastic to plastic deformation.

2. Interaction model

We simulated the behavior of porous dust aggregates with a molecular dynamics approach. In this work we use cuboidal shaped aggregates with an edge length of 30 to 60   μm. Depending on its volume filling factor, such an aggregate consists of up to 6 × 10⁴ micron-sized spherical grains (monomers).

The interaction between individual monomers is based on the work of Dominik & Tielens (1997). Based on earlier theoretical work by Johnson et al. (1971); Dominik & Tielens (1995, 1996), they developed a detailed micromechanical model of the interaction of two microscopic spherical grains. When two monomers touch each other, surface forces allow for the creation of an adhesive contact. Upon deformation of these contacts caused by the relative motion of the monomers, kinetic energy is dissipated. This approach is favorable to our purpose because the process of internal restructuring of the aggregate is modeled far more realistically than in simpler hard-sphere models.

Later, Wada et al. (2007) derived almost the same interaction laws from corresponding potentials. Seizinger et al. (2012) calibrated the interaction model by comparison with laboratory experiments on the compression of porous dust aggregates (Güttler et al. 2009). They observed that the original model of Dominik & Tielens (1997) was too soft. To increase the strength of the aggregates, Seizinger et al. (2012) introduced the rolling and sliding modifiers m_r and m_s that modify the strength of the corresponding type of interaction. By increasing the rolling interaction by a factor of 8 and the sliding interaction by a factor of 2.5, they achieved excellent agreement between numerical simulations and laboratory results.

In this work we use the modified interaction model presented by Seizinger et al. (2012) with m_r = 8 and m_s = 2.5 unless stated otherwise. The material parameters are listed in Table 1.

3. Sample generation

In this section we briefly summarize our sample-generation methods. Here we use cuboidal-shaped dust aggregates (also referred to as dust cakes) of different porosities. In principle, we may employ different methods to generate these samples. As shown in in Fig. 1 the relation between the volume filling factor and the average coordination number n_c depends on the selected generation method (Seizinger & Kley 2013).

The volume filling factor φ is given by $\begin{array}{ccc}\mathit{\phi }\mathrm{=}\frac{\mathit{N}{\mathit{V}}_{\mathrm{p}}}{\mathit{V}}\mathit{,}& & \end{array}$(1)where N denotes the number of monomers, V_p is the volume of an individual monomer, and V is the total volume occupied by the sample. The coordination number, n_c, denotes the number of contacts a monomer has established with its neighbors. The mean n_c of the sample is calculated by averaging the number of contacts of each monomer.

Fig. 1 Relation between the volume filling factor φ and the average coordination number nc. Figure taken from Seizinger & Kley (2013).

3.1. Static compaction

Static compaction constitutes a method that very closely resembles the generation process of dust cakes in laboratory experiments. Initially, a dust cake is built by RBD where single monomers are successively dropped onto the existing sample. They come to rest right at the spot where they hit the existing sample. This growth process results in fractal, highly porous aggregates with a volume filling factor of φ = 0.15 (Watson et al. 1997).

In the second step, the RBD cake is put into a box and compacted until the desired filling factor is reached. However, the compaction must occur slowly enough to avoid inhomogeneities (Seizinger et al. 2012). For the material/monomer size used in this work, a typical speed of the compacting wall is 1   cm   s^-1. The simulation time needed to compact an aggregate is proportional to the desired size and filling factor. Since the number of monomers increases with the size and compactness of the final dust cake, the computational effort per integration step rises as well. Thus, generating samples by static compaction can become a computationally expensive, time-consuming procedure.

Elastic loading of dust cakes compressed to filling factors above ≈0.45 constitutes another setback of this method (Seizinger & Kley 2013). Because of the elastic energy stored in the contacts between monomers the aggregate will start to expand once the confining walls of the compaction box are removed. This effect will alter the results of a measurement of the tensile or shear strength significantly.

Up to a filling factor of ≈0.58 dust cakes may be stabilized in the following way. After slightly disturbing the positions of monomers, the aggregate is kept in a box until the energy induced by the disturbance is damped away (Seizinger & Kley 2013).

3.2. Ballistic aggregation and migration

The generation procedure of ballistic aggregation and migration (BAM) has been proposed by Shen et al. (2008). A larger aggregate is generated by successively shooting in single monomers from random directions onto the existing aggregate.

In Seizinger & Kley (2013) we employed three different methods to select the final position of a monomer hitting the aggregate:

• 1.

  Select the position closest to the spot, where the monomerimpacts the aggregate (referred to as “shortest migration”).

• 2.

  Select the position randomly from all available possibilities (referred to as “random migration”).

• 3.

  Select the position that is closest to the center of mass (referred to as “center migration”).

The volume filling factor of the generated aggregate depends on the selection mechanism. For a given coordination number, aggregates generated with the “shortest migration” feature the lowest filling factor, whereas the “center migration” method leads to the most compact aggregates.

The relation between the filling factor and the coordination number is displayed in Fig. 1, for the different preparation methods. For comparison we show the results for the hexagonal close packing as well, see Seizinger & Kley (2013) for more details. After generating a sufficiently large aggregate, a cuboidal shaped dust cake of the desired size is sliced out to be used for the subsequent numerical experiments.

4. Tensile strength

4.1. Setup

Fig. 2 Example of the tension required to pull the plates apart from each other for two cubical samples of different porosity. As the pull distance increases the required force decreases because cracks form in the sample.

Fig. 3 Outcome of a typical pull experiment on a cubical sample agglomerate. Left: since the adhesion between particle-wall contacts is stronger than between particles, the uppermost layer of particles is ripped off when pulling the plates away from each other. Right: adhesion between particles that are close to one of the plates has been artificially increased. The red dotted line indicates where the additional gluing effect sets in.

In principle, the procedure for determining the tensile strength of a given sample is simple. After attaching two plates at the top and bottom of the cubic aggregate, the plates are pulled apart at constant speed. During this process the pulling force exerted on these plates is measured. At first, the force will increase with the distance between the plates. If a certain distance is exceeded, cracks will form. Thus, the strength of the sample is reduced and the force required to pull the two plates apart decreases. An example showing a typical relation between the pull distance and tension is shown in Fig. 2. It is very similar to experimental data (see Blum & Schräpler 2004, Fig. 4).

The force is determined by summing up the individual interaction forces of all monomers that are in contact with the wall. In the case of the tensile strength, only the component in the direction of the motion of the wall (which is equivalent to the normal vector of the wall) is taken into account. To allow easier comparison between different sized samples, we use their base area to normalize the values and plot the corresponding tension instead of the force.

In accordance with Blum et al. (2006), we define the maximum tension that is measured during a run as the tensile strength of the sample. The displacement at which force/pressure peaks depends on the porosity of the sample. In a sample with a high filling factor, monomers are fixed tightly, which hampers internal restructuring. In contrast, in a fluffy sample individual chains of monomers can be unfolded, and thus the material can be stretched out significantly before the formation of cracks sets in.

To model this setup in our simulations, the sample is put into a box of flat walls. Before slowly moving the top and bottom walls away from each other, we must ensure that enough monomers are in contact with these walls. For this purpose the top and bottom walls are slowly pushed inwards a short distance. For samples with filling factors below φ = 0.2, we use a value of one monomer radius, whereas for more compact samples, we decrease the distance to 0.5−0.1 monomer radii.

Another problem arises when pulling the two plates away from each other. The critical force F_c required to break a contact between two monomers is given by F_c = 3πγR, where γ denotes the surface energy and R the reduced radius of the two particles (Johnson et al. 1971). The wall is modeled as a particle of infinite radius, which means that the reduced radius of a particle-wall contact equals twice the reduced radius of a particle-particle contact (Seizinger et al. 2012). Thus, contacts between two particles can be broken more easily than particle-wall contacts. As a result, the monomers that are in contact with one of the plates tend to get ripped off the remaining sample (see left panel of Fig. 3).

To counter this effect we artificially increase the strength of the adhesion between two monomers depending on the distance to the plates. To achieve this “gluing effect” the force/potential of the normal interaction that is responsible for the adhesion is multiplied by a gluing factor κ. To avoid discontinuities in the particle interaction, κ is interpolated linearly depending on the distance to the closest plate. Above a threshold of eight particle radii κ is set to 1 and thus the default JKR interaction is used. As shown in the righthand panel of Fig. 3, this mechanism leads to the rupture somewhere in the middle of the sample rather than just tearing off the upper- or lowermost layer of monomers. Because the first cracks will form where the aggregate is weakest, the exact location is random owing to the inhomogeneous structure of the aggregate.

We tested different maximum values of κ and found that a value of 2 is sufficient for our purpose. For κ < 1.5, samples do not break in the center anymore. On the other hand, higher values do not alter the measured tensile strength significantly (see Fig. 4).

Fig. 4 Tensile strength for different wall glue factors κ and normal interaction modifiers mn using the model from Seizinger et al. (2012) (mr = 8, ms = 2.5). The dotted lines represent linear fits for filling factors below 0.5. For comparison, we also performed simulations using the model employed by Wada et al. (2007) (mr = ms = 1). All samples are BAM cakes with an edge length of 50 × 50 × 30   μm.

4.2. Results

Apart from the wall-gluing factor κ, there are several other parameters whose influence has to be studied. To determine the influence of the rolling and sliding interaction, we performed a series of simulations using m_r = m_s = 1, which is equivalent to the model of Wada et al. (2007). In this and all of the following simulations, we used a wall-gluing factor of κ = 2. Apparently, internal restructuring that is governed by rolling and sliding does not play a major role when determining the tensile strength (see the purple squares in Fig. 4 showing results for the model from Wada et al. 2007).

Contacts between the monomers have to be broken while the sample aggregate is torn apart. Therefore, we expect that the measured tensile strength depends on the number of contacts that have to be broken, as well as on the critical force F_c that is necessary to break a contact between individual monomers. To check this hypothesis, we alter the strength of the normal force by multiplying it by the normal force modifier m_n. Indeed, when doubling the strength of the normal interaction (and thus F_c) by setting m_n = 2, we observe a steeper increase in the tensile strength with the filling factor (blue asterisks in Fig. 4). When determining a linear fit for filling factors below 0.5 to the κ = 3,m_n = 1 and κ = 3,m_n = 2 simulations we get slopes of 7.9   kPa and 15.3   kPa, respectively. Their ratio of 15.3/7.9 = 1.94 is very close to the value of the normal interaction modifier m_n = 2. This strongly suggests that the pull-off force F_c is critical for the measured value of the tensile strength.

Independent of κ, m_n, or the rolling and sliding interaction, there is a striking drop in the tensile strength for filling factors above 0.5. To unravel its cause we first used a different type of sample. Much to our surprise, the critical filling factor at which the tensile strength drops is different for each type of aggregate and close to the maximum filling factor that may be achieved by a given generation method (see Fig. 5). Apparently, the micro-mechanical behavior is not governed by the volume filling factor alone.

Fig. 5 Comparison of the relation between the filling factor φ and the tensile strength of different sample types. Most noticeably, the measured tensile strength always drops when a certain type-specific filling factor is exceeded.

Therefore we take the average coordination number n_c into account. From the relation between the filling factor and the coordination number of the different sample types (see Fig. 1), we find that the drop in the measured tensile strength of BAM aggregates coincides with a value of n_c around ≈4.5 (see Fig. 6). This points to the influence of the coordination number on the micro-mechanical properties of the sample aggregates. If the average coordination number is low, the majority of monomers may react to an external stress by rearranging themselves. Thus, a large number of monomers participate in absorbing the external stress, and the aggregate exhibits rather ductile behavior. As the coordination number increases, monomers are fixated in their positions more tightly. For compact aggregates the monomers cannot rearrange themselves freely any longer, which means that fewer monomers have to absorb the applied strain. Therefore, the aggregates become brittle.

The effect of brittleness can clearly be seen in the case of the hexagonal lattice aggregates. Because of their regular, crystal-like structure (see Seizinger & Kley 2013, Fig. 1a), their capability of internal restructuring is very limited. Thus, contacts break very easily when external strain is applied. As a result, the measured tensile strength is considerably lower than for BAM or static compaction aggregates.

Fig. 6 Comparison of the relation between tensile strength and initial coordination number nc of different sample types. The tensile strength of the different BAM aggregates drops for values of nc around 4.5. For the hexagonal close-packing aggregates we see a clear drop for nc → 6.

Additionally, the pressure exerted on the sample when slowly pressing on the top and bottom walls suffices to disrupt very compact aggregates. Moving the top and bottom walls inwards by a distance of 0.5 to 0.15 particle radii is necessary to establish a firm contact between the walls and the sample. For a BAM (center migration) aggregate with φ = 0.59, for example, the average coordination number decreased from 5.94 to 4.95 after moving the top wall inwards by a distance of only 0.1 monomer radii. This means that the strength of the sample is lowered during the preparation process. This raises the question of whether the transition from ductile to brittle behavior or the disturbance when affixing the plates is the dominant effect.

Recently, Kataoka et al. (2013) presented a different approach to determining the compressive strength of highly porous (φ < 0.1) dust aggregates by using periodic boundary conditions. A similar approach might work for the tensile strength as well and would avoid the problem of attaching the plates to a highly compact sample without lowering its strength. Luckily, measuring the tensile strength of hexagonal lattice aggregates also allows us to circumvent this problem. Because of their regular grid structure, the contact between the wall and all particles of the top or bottom layer is established without compacting the sample. Nevertheless, we observe a significant drop in the tensile strength for n_c → 6 (see Fig. 6). This observation allows us to conclude that the disruption caused by affixing the plates only plays a secondary role.

Fig. 7 Dissipated energy (solid line) and total number of broken contacts (dotted line) of BAM and static compaction (SC) aggregates while pulling the plates apart. The size of both samples is 50 × 50 × 30   μm, and their initial filling factor is φ = 0.28. The critical force Fc and critical distance δc at which a contact breaks are used to normalize the energies. For comparison, the corresponding tension (dashed-dotted line) is also plotted.

To explain the discrepancy between BAM and static compaction aggregates as shown in Fig. 5, the process of pulling the plates apart is analyzed more closely. Comparing the dissipated energy of the different types of aggregates, we clearly see that tearing the static compaction aggregate apart requires more energy (Fig. 7). Thus, a higher value of the tensile strength is measured. Tracking the number of contacts that have been broken since the start of the simulation provides us with a hint at an explanation of this observation: Compared to the BAM aggregate, fewer than half as many contacts break while pulling the static compaction aggregate apart, although more kinetic energy is dissipated.

This observation might be surprising at first sight, yet it illustrates the importance of the internal structure for the outcome of the measurements. In reaction to the applied strain, the internal structure of the aggregates changes, where inelastic rolling accounts for ≈80% and the breaking of contacts for only ≈3% of the total dissipated energy for both aggregate types. In the case of the static compaction aggregate, the final number of broken contacts is reached much more quickly after pulling the sample roughly 7.5   μm apart. Likewise, the peak of the measured tension is reached earlier. At this point, a lot more energy has been dissipated by internal restructuring (mainly by inelastic rolling) than for the BAM aggregate.

However, for the BAM aggregates, the total number of broken contacts increases much faster. This means that the internal structure of this type of aggregate allows for less restructuring before contacts start to break. Since roughly 97% of the dissipated energy is required for restructuring rather than breaking the contacts, we measure a lower tensile strength for BAM aggregates.

We think that the reason behind this observation lies in the fact that the tensile strength test is the reversal of the compression process by which the static compaction aggregates have been generated. While being slowly compacted to the desired filling factor, enough time had been given to the fractal chains of the initial RBD aggregates to rearrange themselves. As a result, the monomers of the static compaction aggregates adopt a structure that is favorable to withstand an external load.

While pulling the plates apart, the total number of broken contacts (Fig. 7) sometimes decreases. This happens when the connection between two already stretched out parts of the sample breaks and both parts “snap back”.

We also checked whether there are any preferred directions resulting from the generation process. For this reason we rotated cubic static compaction samples by 90° before determining the tensile strength. Reassuringly, we measured the same values (not shown in this work) and may rule out that the direction of the compaction induces any preferred direction in the structure of the sample.

Fig. 8 Tensile strength measurements for BAM (center migration) aggregates of different size. The error bars have been determined by averaging the results of six runs with different samples.

Fig. 9 Comparison of the results of simulations with static compaction and BAM (center migration) aggregates to laboratory experiments by Blum & Schräpler (2004); Blum et al. (2006). The fitting curve has been obtained by combining both data sets for filling factors φ < 0.5.

As a last step we varied the geometry of the samples to check that the size influences the results. As seen in Fig. 8, the results do not vary significantly if we alter the size of the samples. In Fig. 9 we compare our results with laboratory experiments performed by Blum & Schräpler (2004) and Blum et al. (2006). Considering that their samples have been produced by static compaction, our results show good agreement with their data for filling factors below φ = 0.5. Because of the drop in the tensile strength for n_c → 6 explained in the previous paragraphs, we cannot compare our simulations to laboratory experiments for higher filling factors. Earlier compression simulations have already indicated that our physical model may no longer be valid for highly compact aggregates (Seizinger et al. 2012). Luckily, the filling factor regime relevant to the growth processes of planetesimals is below 0.5 (e.g., Teiser et al. 2011).

Fig. 10 Tensile strength for aggregates composed of differently sizes of monomers. The default monomer radius is rp = 0.6   μm, which is also used in the corresponding laboratory experiments.

To determine the fitting curve depicted in Fig. 9, only data points for filling factors below 0.5 have been taken into account. Because of the significant difference between the static compaction and BAM aggregates, we determined two fit curves T_SC(φ) and T_BAM(φ), respectively. We obtain $\begin{array}{ccc}{\mathit{T}}_{\mathrm{SC}}\mathrm{\left(}\mathit{\phi }\mathrm{\right)}\mathrm{=}\mathrm{13.4}\hspace{0.17em}{\mathit{\phi }}^{\mathrm{1.62}}\hspace{0.17em}\mathrm{kPa}\mathit{,}& & \end{array}$(2)and $\begin{array}{ccc}{\mathit{T}}_{\mathrm{BAM}}\mathrm{\left(}\mathit{\phi }\mathrm{\right)}\mathrm{=}\mathrm{11.2}\hspace{0.17em}{\mathit{\phi }}^{\mathrm{1.88}}\hspace{0.17em}\mathrm{kPa}\mathit{.}& & \end{array}$(3)Based on their generation process, the laboratory samples should resemble the static compaction cakes. Indeed, for φ ≈ 0.2, our simulations agree very well with laboratory experiments (Blum & Schräpler 2004; Blum et al. 2006). However, for higher filling factors, the laboratory results lie somewhere between the static compaction and BAM results. In private conversation Jürgen Blum (Braunschweig) pointed out that creating more compact samples in the lab sometimes proved to be a difficult task. Thus, we also determined a fit T(φ) to the combined results of the static compaction and BAM aggregates. Using the combined data points from both aggregate types shown in Fig. 9 for a single fit T(φ), we find, for values of φ < 0.5,$\begin{array}{ccc}\mathit{T}\mathrm{\left(}\mathit{\phi }\mathrm{\right)}\mathrm{=}\mathrm{12.6}\hspace{0.17em}{\mathit{\phi }}^{\mathrm{1.77}}\hspace{0.17em}\mathrm{kPa}\mathit{.}& & \end{array}$(4)

4.3. Influence of the monomer size

Most other numerical simulations dealing with silicates have been performed with monomers with a diameter of 1.2 and 1.5   μm because these sizes have been used in laboratory experiments with spherical silicate grains. Out of curiosity we varied the size of the monomers. In general, it can be said that according to JKR-theory, the adhesion forces increase as grains get smaller. Indeed, we find that the tensile strength depends closely on the size of the monomers (see Fig. 10). Our interaction model has not been calibrated for monomer radii other than r_p = 0.6   μm. Therefore the rolling and sliding modifiers m_r and m_s may not have the correct values to properly describe restructuring processes. However, in the case of the tensile strength, this problem does not arise, as we have already seen that it is mainly governed by the normal interaction (see Fig. 4).

From Fig. 4 we see that the tensile strength scales linearly with the normal interaction. Altering the monomer size by a factor of 2 for the transition from 0.6   μm to 1.2   μm sized monomers, the change in the measured tensile strength differs from 2. At first glance this may seem odd since the critical pull-off force F_c depends linearly on the monomer radius. However, the dependence of the normal force acting upon the monomers before they are separated on the monomer radius is nonlinear (see Seizinger et al. 2012, Eqs. (2) and (3)). Nevertheless, these simulations confirm the importance of the pull-off force F_c that has already been shown in Fig. 4.

The results clearly demonstrate the effect of the stickiness of the single monomers on the tensile strength. For future work it would be interesting to perform simulations with aggregates composed of monomers of different sizes.

5. Shear strength

5.1. Setup

Fig. 11 A snapshot taken during a shear strength test using an aggregate with an edge length of 30 ×   30   ×   40   μm. The upper plate is slowly moving to the right. Adhesion between particles that are close to one of the plates has been artificially increased.

The shear strength of a porous aggregate is determined in a similar way to the tensile strength. As before, two plates are attached to the top and bottom of the sample. During the shearing motion of the plates, the force acting on them is tracked. However, for the shear strength, the direction of motion is perpendicular to the normal of the wall, i.e. tangential to a cuboid surface. During the simulation the vertical positions of the top and bottom, walls remain fixed to keep the filling factor constant. This means, similar to the “fixed walls” used in the work of Seizinger et al. (2012), that the normal component of the force acting upon the walls is ignored. As before, the initial base area of the sample is used to normalize the force.

The setup closely resembles the tensile strength setup. A firm contact between the sample and the plates is achieved in the same way as described in Sect. 4.1. To prevent the monomers that are in contact with the moving wall from being torn away from the sample, an additional “gluing effect” is applied to particles that are close to one of the plates. A snapshot taken during a typical simulation is depicted in Fig. 11.

While the top plate is slowly moving, shearing sets in. With increasing pulling distance, cracks will form and reduce the strength of the sample. Thus, we expect a similar shape to the one for the tensile strength if we plot the tension acting on the moving plate with respect to the displacement. Indeed, the example shear strength curve shown in Fig. 12 resembles the curves shown in Fig. 2. Similar to the tensile strength case, we define the shear strength as the maximum tension that is measured during the simulation. Again, the higher the porosity of a sample, the larger the necessary displacement at which the force peaks.

Fig. 12 Example of the tension (force per area) acting on the upper plate during the shearing motion.

5.2. Results

Fig. 13 Comparison of the relation between the filling factor φ and shear strength of different sample types and sizes. The error bars have been determined by averaging the results from six different samples. The black dotted line has been obtained by fitting a power law to the results from BAM and static compaction aggregates with an edge length of 40 × 40 × 50   μm.

Owing to the computational demand of the simulations, the size of our samples is limited to values below 0.1   mm. To study the dependency of our results on the sample size, we prepared both BAM and static compaction aggregates with different edge lengths. For each data point, six different samples with equal statistical properties have been generated.

Some of the results are shown in Fig. 13. As we can see, the results of the different sample sizes do not alter significantly. In order to check whether the length of the sample in the direction of the shearing motion is sufficient we also performed simulations for sample sizes of 80 × 40 × 50   μm and 120 × 40 × 50   μm. Owing to the huge number of particles, these simulations took several weeks. Therefore we restricted the values of the filling factor to φ = 0.33 and φ = 0.49. The deviation to the values obtained from the smaller 40 × 40 × 50   μm aggregates was ≈8–10% for φ = 0.33 and ≈2–3% for φ = 0.49. Thus, we may draw the conclusion that the samples are in fact sufficiently large. With the exception of the most compact samples (φ = 0.59), the error bars obtained by averaging the results from the six samples are very small.

Interestingly, we do not observe a significant difference between the static compaction and BAM aggregates as in the case of the tensile strength. As explained in Sect. 4.2, owing to their generation process, the internal structure of the static compaction aggregates is more favorable to counteracting external loading/tension. However, this does not apply to shearing motion that results in different kinds of deformation compared to the tensile strength test. Therefore, the two types of samples exhibit similar values for the shear strength.

To provide SPH simulations with a model for the shear strength that is easy to implement, we describe the dependency of shear strength S on the filling factor φ with a power law S(φ) = aφ^b. Using the results from BAM and static compaction aggregates of 40 × 40 × 50   μm edge length, we obtained (see Fig. 13) $\begin{array}{ccc}\mathit{S}\mathrm{\left(}\mathit{\phi }\mathrm{\right)}\mathrm{=}\mathrm{21.7}{\mathit{\phi }}^{\mathrm{2.65}}\hspace{0.17em}\mathrm{kPa}\mathit{.}& & \end{array}$(5)

5.3. Comparison with the SPH continuum model

One objective of the present investigations has been to compare the resulting strengths with those adopted in the SPH simulations by Geretshauser et al. (2010). To model shear failure, Sirono (2004) introduced a von Mises yielding criterion in his SPH simulations. The required shear strength is in principle equivalent to the shear strength obtained in our calculations. Güttler et al. (2009) calibrated their SPH model in an extensive process by comparing simulation results with laboratory experiments. They find a new representation for the dynamic compressive strength, $\begin{array}{ccc}{\mathrm{\Sigma }}_{\mathrm{SPH}}\mathrm{\left(}\mathit{\phi }\mathrm{\right)}\mathrm{=}\mathrm{13}{\left(\frac{{\mathit{\phi }}_{\mathrm{2}}\mathrm{-}{\mathit{\phi }}_{\mathrm{1}}}{{\mathit{\phi }}_{\mathrm{2}}\mathrm{-}\mathit{\phi }}\mathrm{-}\mathrm{1}\right)}^{\mathrm{0.58}\mathrm{·}\ln \mathrm{10}}\hspace{0.17em}\mathrm{kPa}& & \end{array}$(6)with φ₁ = 0.12 and φ₂ = 0.58, and they chose the tensile strength according to Blum & Schräpler (2004), $\begin{array}{ccc}{\mathit{T}}_{\mathrm{SPH}}\mathrm{\left(}\mathit{\phi }\mathrm{\right)}\mathrm{=}{\mathrm{10}}^{\mathrm{2.8}\mathrm{+}\mathrm{1.48}\hspace{0.17em}\mathit{\phi }}\hspace{0.17em}\mathrm{Pa}\mathit{.}& & \end{array}$(7)For the shear strength, no experimental data have been available. Therefore the shear strength was adopted following Sirono (2004) according to $\begin{array}{ccc}{\mathit{S}}_{\mathrm{SPH}}\mathrm{=}\sqrt{{\mathrm{\Sigma }}_{\mathrm{SPH}}{\mathit{T}}_{\mathrm{SPH}}}\mathit{.}& & \end{array}$(8)In Fig. 14 the fit curves of tensile strength T(φ), Eq. (4), and shear strength S(φ), Eq. (5), are compared to the corresponding values of the SPH model, T_SPH(φ), Eq. (7), and S_SPH(φ), Eq. (8). As can be seen, the tensile strength curves match rather well. This emphasizes that the present molecular dynamics method is well suited to modeling highly porous aggregates. The shear strength curves, however, differ by nearly one order of magnitude. This indicates that the approach of Sirono (Eq. (8)) for the SPH shear strength, which is not based directly on laboratory experiments, might be inappropriate. But during the calibration process it was found already that the SPH simulation results for the chosen reference problems only depend weakly on the exact values of the shear strength (Güttler et al. 2009).

Fig. 14 Comparison of the fit curves for the compressive Σ, tensile T, and shear strength S derived in this work, and the corresponding functions ΣSPH, TSPH, and SSPH adopted in the SPH code by Geretshauser et al. (2010). The compressive strength Σ(φ) has already been determined in earlier work (Seizinger et al. 2012).

6. Conclusions

This work supports the observation of Seizinger & Kley (2013) that the sample generation method influences its mechanical behavior significantly. Whereas the bouncing behavior of microscopic dust aggregates differs little for BAM and static compaction aggregates, they do behave differently when external strain is applied (see Fig. 5). It is important to keep this in mind when comparing numerical simulations to laboratory results.

Observing the transition from ductile to brittle behavior for coordination numbers of ≈6 is very interesting. It certainly influences the outcome of collisions as well. For brittle aggregates, fragmentation will play a significantly larger role.

In this work we determined simple power laws to describe the relation between the tensile strength (see Eq. (4)) or shear strength (see Eq. (5)) and the porosity. In combination with earlier work on the the compressive strength (Seizinger et al. 2012), it provides a complete description when the inelastic regime is entered upon deformation of porous dust aggregates. Since the dissipation of the kinetic impact energy is critical, this knowledge is crucial for continuum simulations of collisions of macroscopic porous aggregates.

Comparing with a special SPH model, we find that our tensile strength agrees well with the tensile strength adopted in the SPH code. The same holds for the compressive strength as found in earlier work (Seizinger et al. 2012). However, the shear strength differs significantly. Future analysis has to show whether our improved relation for the shear strength will have a fundamental impact on the SPH simulation results, or whether the shear strength only alters details in the simulations, as might be indicated by previous work.

Acknowledgments

A. Seizinger acknowledges support through the German Research Foundation (DFG) grant KL 650/16. The authors acknowledge support through DFG grant KL 650/7. Additional support through the German Research Foundation (DFG) through grant KL 650/11 within the Collaborative Research Group FOR 759: The formation of Planets: The Critical First Growth Phase is acknowledged. We thank the anonymous referee for pointing out possible misunderstandings and helping to improve the quality of the paper.

References

All Tables

All Figures

	Fig. 1 Relation between the volume filling factor φ and the average coordination number nc. Figure taken from Seizinger & Kley (2013).
In the text

	Fig. 2 Example of the tension required to pull the plates apart from each other for two cubical samples of different porosity. As the pull distance increases the required force decreases because cracks form in the sample.
In the text

Fig. 3 Outcome of a typical pull experiment on a cubical sample agglomerate. Left: since the adhesion between particle-wall contacts is stronger than between particles, the uppermost layer of particles is ripped off when pulling the plates away from each other. Right: adhesion between particles that are close to one of the plates has been artificially increased. The red dotted line indicates where the additional gluing effect sets in.

In the text

Fig. 4 Tensile strength for different wall glue factors κ and normal interaction modifiers mn using the model from Seizinger et al. (2012) (mr = 8, ms = 2.5). The dotted lines represent linear fits for filling factors below 0.5. For comparison, we also performed simulations using the model employed by Wada et al. (2007) (mr = ms = 1). All samples are BAM cakes with an edge length of 50 × 50 × 30   μm.

In the text

	Fig. 5 Comparison of the relation between the filling factor φ and the tensile strength of different sample types. Most noticeably, the measured tensile strength always drops when a certain type-specific filling factor is exceeded.
In the text

	Fig. 6 Comparison of the relation between tensile strength and initial coordination number nc of different sample types. The tensile strength of the different BAM aggregates drops for values of nc around 4.5. For the hexagonal close-packing aggregates we see a clear drop for nc → 6.
In the text

Fig. 7 Dissipated energy (solid line) and total number of broken contacts (dotted line) of BAM and static compaction (SC) aggregates while pulling the plates apart. The size of both samples is 50 × 50 × 30   μm, and their initial filling factor is φ = 0.28. The critical force Fc and critical distance δc at which a contact breaks are used to normalize the energies. For comparison, the corresponding tension (dashed-dotted line) is also plotted.

In the text

	Fig. 8 Tensile strength measurements for BAM (center migration) aggregates of different size. The error bars have been determined by averaging the results of six runs with different samples.
In the text

	Fig. 9 Comparison of the results of simulations with static compaction and BAM (center migration) aggregates to laboratory experiments by Blum & Schräpler (2004); Blum et al. (2006). The fitting curve has been obtained by combining both data sets for filling factors φ < 0.5.
In the text

	Fig. 10 Tensile strength for aggregates composed of differently sizes of monomers. The default monomer radius is rp = 0.6   μm, which is also used in the corresponding laboratory experiments.
In the text

	Fig. 11 A snapshot taken during a shear strength test using an aggregate with an edge length of 30 ×   30   ×   40   μm. The upper plate is slowly moving to the right. Adhesion between particles that are close to one of the plates has been artificially increased.
In the text

	Fig. 12 Example of the tension (force per area) acting on the upper plate during the shearing motion.
In the text

	Fig. 13 Comparison of the relation between the filling factor φ and shear strength of different sample types and sizes. The error bars have been determined by averaging the results from six different samples. The black dotted line has been obtained by fitting a power law to the results from BAM and static compaction aggregates with an edge length of 40 × 40 × 50   μm.
In the text

	Fig. 14 Comparison of the fit curves for the compressive Σ, tensile T, and shear strength S derived in this work, and the corresponding functions ΣSPH, TSPH, and SSPH adopted in the SPH code by Geretshauser et al. (2010). The compressive strength Σ(φ) has already been determined in earlier work (Seizinger et al. 2012).
In the text