EDP Sciences
Highlight
Free Access
Issue
A&A
Volume 536, December 2011
Article Number A104
Number of page(s) 11
Section Planets and planetary systems
DOI https://doi.org/10.1051/0004-6361/201117645
Published online 19 December 2011

© ESO, 2011

1. Introduction

Terrestrial planets are thought to form by core accretion in protoplanetary accretion discs (e.g. Lissauer & Stewart 1993). The discs consist of both hydrogen and helium gas and a small fraction of dust, which in the early stages exists as micron sized dust grains (e.g. Dullemond et al. 2007). Interaction between the gas and dust induces relative motions between the dust grains and causes collisions among them (e.g. Weidenschilling & Cuzzi 1993). Mutual sticking by van der Waals forces occurs and by this process millimetre- to centimetre-sized, highly porous pre-planetesimals form (e.g. Blum & Wurm 2008). Additional growth by means of collision is endangered by several obstacles: rebound, material loss by accretion or photoevaporation, and fragmentation. Zsom et al. (2010) found that the growth of pre-planetesimals to kilometre-sized planetesimals may be halted by mutual compaction and bouncing, which leads finally to a population of dust and pebble-sized objects in the protoplanetary disc. Weidenschilling (1977) found that metre-sized bodies have the highest drift speed towards the host star, where the objects get accreted or photoevaporated. For these objects, the drift timescale is much shorter than the time needed for a planet to form. One of the most serious obstacles to planetesimal formation is the fragmentation barrier: with increasing aggregate size the mutual collision velocities also increase and catastrophic disruption might become more common than mass gain (e.g. Brauer et al. 2008).

If by some mechanism a sufficient population of kilometre-sized objects is formed, their self-gravity acts as an accretion mechanism and ensures further growth to planets. In this gravity dominated regime, kilometre-sized objects and yet larger ones are referred to as planetesimals.

To investigate planetesimal formation, global dust coagulation models (e.g. Dullemond & Dominik 2005; Zsom et al. 2010; Birnstiel et al. 2010) numerically compute the evolution of dust aggregates in a protoplanetary disc starting from micron-sized dust grains. These models require detailed collision statistics evaluated for the outcome of two-body pre-planetesimal collisions as functions of important parameters such as mass ratio, porosity, impact parameter, and relative velocity of the collision partners. In the past decade, many experimental studies have been carried out to collect these data (e.g. Blum & Wurm 2008; Güttler et al. 2010, for summaries). However, laboratory setups become difficult and infeasible for centimetre-sized pre-planetesimals and larger ones since they have to be carried out under protoplanetary disc conditions, i.e. in a vacuum and under micro-gravity. In this regime, collision data has to be acquired by numerical simulations.

In their pioneering work, Benz & Asphaug (1994, 1995) developed a solid body smoothed particle hydrodynamics (SPH) implementation for brittle material based on the dynamical fracture model by Grady & Kipp (1980), which evolves the propagation of cracks starting from an initial Weibull (1939) distribution of flaws. On these grounds, Benz (2000) investigated collisions between rocky non-porous pre-planetesimals in their typical collision velocity regime from 5   m   s-1 to 40   m   s-1 and found only disruption. However, experiments (e.g. Blum & Wurm 2008) indicate that pre-planetesimals are fluffy, porous objects rather than solid rocks. For this reason, Sirono (2004) developed an SPH implementation of a porosity model based on porosity-dependent material quantities such as bulk and shear moduli as well as compressive, shear, and tensile strengths. In his simulations of porous ice aggregates, he found that owing to their porous structure the kinetic energy can be dissipated effectively, permitting collisional growth. Sirono’s approach included a damage model based again on the dynamical fracture model of Grady & Kipp (1980). In addition, a damage restoration model tries to capture damage healing by compression. However, Schäfer et al. (2007) found that Sirono’s damage model is not applicable to their SPH simulations of porous ice because it includes damage by compression. In contrast, materials such as porous ice and porous dust form new molecular bondings when they are compressed. Thus, their strengths increase but do not decrease as predicted by the Sirono damage model. On the basis of the implementation of Benz & Asphaug (1994, 1995), Jutzi et al. (2008) combined the damage model of Grady & Kipp (1980) and the p − α model of Herrmann (1969) to simulate porous brittle media. The model was calibrated with laboratory high velocity impact experiments (Jutzi et al. 2009). However, the model by Jutzi et al. (2008) is also not applicable to porous protoplanetary dust because it includes damage by compression and has no suitable model of damage healing.

To simulate SiO2 dust, which is used as protoplanetary dust analogue in many experimental studies, we developed a suitable porosity model based on the ideas of Sirono (2004), that was implemented into a solid body SPH code to numerically investigate pre-planetesimal collisions (Geretshauser et al. 2010). With the aid of laboratory benchmark experiments (see also Güttler et al. 2009), the porosity model was calibrated and tested for the simulation of SiO2 dust. In particular, the simulation reproduced the compaction, bouncing, and fragmentation behaviour of homogeneous dust aggregates of various porosities in a quantitatively correct way. However, this approach does not include an explicit damage model.

In this paper, we present an approach to explicitly resolve inhomogeneous structures of dust aggregates. Since an inhomogeneous density distribution is one of the main sources of fracture in porous aggregates, in this work we approximate the damage evolution by the evolution of density. The inhomogeneity is described by a Gaussian density distribution and the typical size of a clump. The aggregate inhomogeneity of SiO2 dust can easily be measured in the laboratory by X-ray tomography (see Güttler et al. 2009, Fig. 10). In contrast, the Grady & Kipp (1980) damage model relies on the correct parameters for the initial Weibull distribution of flaws, which are difficult to determine in laboratory experiments. Our approach of resolving inhomogeneous structures is closely connected to the porosity model, where the strength quantities are porosity dependent. Therefore, any fluctuation in porosity induces a fluctuation in strength over the aggregate, which is utilised to model the evolution of the material under stress.

To classify the collisional outcome of our simulations involving inhomogeneous aggregates, we utilise the four-population model. This classification scheme was introduced by Geretshauseret al. (2011) to enable a precise mapping of all collisional outcomes involving any combination of sticking, bouncing, and fragmentation processes. Four fragment populations are distinguished according to their masses: both the largest and the second largest fragment as well as a set of fragments whose mass distribution can be approximated by a power law. To take into account the resolution limit, the sub-resolution population contains all fragments that consist of a single SPH particle.

The outline of this paper is the following. In Sect. 2, we describe the SPH method with extensions for the simulation of solid bodies and the applied porosity model. Section 3 firstly describes a previous study of a damage model, then we present our approach based on the inhomogeneity of dust aggregates and its implementation. The results of collisions with inhomogeneous aggregates are discussed in Sect. 4. We summarise our findings in Sect. 5, where we also discuss the results and provide an outlook on future work.

2. Solid body SPH and porosity model

For our simulations, we apply the smoothed particle hydrodynamics (SPH) numerical method originally developed by Lucy (1977) and Gingold & Monaghan (1977) for astrophysical flows. In subsequent studies, SPH was extended to model the elastic and plastic deformations of solid materials (e.g. Libersky & Petschek 1991; Benz & Asphaug 1994; Randles & Libersky 1996). Our parallel SPH code parasph (Hipp & Rosenstiel 2004; Schäfer 2005; Schäfer et al. 2007; Geretshauser et al. 2010) is based on these concepts. In the SPH scheme, the time evolution of a solid body is simulated by means of the Lagrangian equations of continuum dynamics which represent the continuum and momentum equations, respectively. The quantities ρ, vα, and σαβ have their usual meanings, i.e. density, velocity, and the stress tensor. The Greek indices denote the spatial components and the Einstein sum convention applies to them. Since our equation of state (EOS) is energy-independent, we do not solve the energy equation, hence omit to state it here. Within the framework of SPH, the continuous quantities are discretised into mass packages, which we refer to as “particles”. These represent the sampling points of the method and interact with each other within a limited spatial range, called the smoothing length h. The influence of a particle b on a particle a depends on the particle distance |xa − xb| and is weighted by the kernel function W(|xa − xb|,h). This is usually a spherically symmetric function with compact support, which is differentiable at least to first order. We utilise the cubic spline kernel of Monaghan & Lattanzio (1985).

The particles of the SPH scheme move according to the Lagrangian form of the equations of continuum mechanics. Thus, they represent a natural frame of reference for any deformation of the simulated solid body. In the context of fragmentation, two fragments are considered to be separate once their constituent subsets of particles no longer interact, i.e., when the fragments are separated by more than 2   h. The fragments may be spread out over an unlimited computational domain. This is why computations are carried out at the particle positions.

The stress tensor in the momentum equation can be divided into the term accounting for pure hydrostatic pressure p and the term representing pure shear, which is given by the traceless deviatoric stress tensor Sαβ. Hence, (3)In contrast to viscous fluids, p and Sαβ do not depend on the velocity gradient but on the deformation of a solid body given by the strain tensor (4)where x′ are the coordinates of the deformed body. In particular, the relations that hold are where the quantity d denotes the dimension and K and μ are the bulk modulus and shear modulus, respectively, for elastic deformation. To obey the principle of frame invariance, rotation terms have to be added to the deviatoric stress tensor. For this, we apply the Jaumann rate form, which was previously used by Benz & Asphaug (1994) and Schäfer et al. (2007), and described by Geretshauser et al. (2010) in more detail.

The above relations (Eqs. (5) and (6)) describe the time evolution of a solid body in the elastic regime, i.e., for reversible deformations. With respect to pre-planetesimal collisions, we are concerned with highly porous material that undergoes elastic and plastic deformations. Both aspects are described by means of the porosity model of Geretshauser et al. (2010), which is based on the ideas of Sirono (2004) and can be categorised as a plasticity-based porosity model.

In this approach, the quantities for the elastic regime (bulk and shear moduli), as well as the strength quantities (compressive, shear, and tensile strength), depend on the filling factor, which accounts for the sub-resolution porosity. The filling factor φ is related to the density ρ and the porosity Φ by (7)where ρs is the density of the matrix material. We use a constant value of ρs = 2000   kg   m-3 for our SiO2 material (e.g. Blum & Schräpler 2004). In particular, the elastic quantities are related to the filling factor by a power law (8)where γ = 4, K0 is the bulk modulus of an uncompressed random ballistic deposition (RBD) dust sample with φRBD = 0.15 (e.g. Blum & Schräpler 2004), such that K0 = K(φRBD) = 4.5   kPa.

For the compressive, tensile, and shear strength we adopt some empirical relations (see Güttler et al. 2009; Geretshauser et al. 2010; and Geretshauser et al. 2011, for details). These relations allow us to reproduce the compaction, bouncing, and fragmentation behaviour to high accuracy in a velocity range of the order of 0.1 to 10   m   s-1 (Geretshauser et al. 2010). Geretshauser et al. (2011) showed that all collision types found in laboratory experiments can be reproduced. In the same reference, the laboratory velocity thresholds for the transitions between sticking, bouncing, and fragmentation behaviour resembled the simulation results for collisions of medium porosity aggregates with up to 27.5   m   s-1. Therefore, we expect these relations to be valid for collisions below the sound speed of the dust material, which is  ~ 30   m   s-1 (Güttler et al. 2009). For the compressive strength, (9)where φmin + ε < φ < φmax and ε = 0.005. The quantities φmin = 0.12 and φmax = 0.58 are the minimum and maximum filling factors, respectively, in the compressive strength relation. Both values can be exceeded by the material. The power of the compressive strength relation is ln(10) times the parameter Δ with Δ = 0.58. The quantity pm = 260   Pa accounts for the mean pressure of the compressive strength relation. According to Geretshauser et al. (2010, 2011), these parameters are valid for impacts below  ~30   m   s-1. The influence of Δ and pm was studied by Geretshauser et al. (2010). For φ ≤ φmin + ε, the compressive strength relation is continuously extended by the constant function Σ(φ) = Σ(φmin + ε) and for φmax ≤ φ, we set Σ(φ) = ∞. The tensile strength is given by (10)and the shear strength is the geometric mean of the compressive and tensile strength (11)We note that the three strength quantities are filling-factor dependent. This is because as the filling factor increases the monomers of the dust material are packed more closely. Thus, each monomer establishes more bonds, which increases the strength. The filling-factor dependence of the strength quantities is exploited by the approach presented in this paper. Together the bulk modulus, the compressive strength, and the shear strength yield the full equation of state for pure hydrostatic compression or tension (12)where , and and are critical filling factors. The value of marks the transition between elastic and plastic compression and defines the transition between elastic and plastic tension. The filling factors in-between the critical values represent the elastic regime. The quantity is the reference filling factor that represents the filling factor of the porous material at vanishing external pressure. Once the material undergoes plastic deformation, the reference filling factor is reset using the middle relation of Eq. (12). Details of this pressure reduction process can be found in Geretshauser et al. (2010). In the elastic regime (middle line of Eq. (12)), the pressure varies with density according to the elastic path given by Eq. (8). If is exceeded, the pressure follows the much shallower compressive strength path Σ(φ) in the positive φ-direction, whereas if is exceeded, the pressure follows the shallower tensile strength path T(φ) in the negative φ-direction (cf. Fig. 1 in Geretshauser et al. 2010). In both cases, the absolute value of the pressure, which is computed according to the elastic equation of state, is reduced to the compressive or tensile strength, respectively. This accounts for the stress relieved by plastic deformation. Details of the implementation are described by Geretshauser et al. (2010).

For the plastic deformation by pure shear, the deviatoric stress tensor Sαβ is reduced according to the von Mises plasticity given by (e.g. Benz & Asphaug 1994) (13)The function f is defined by . An alternative approach (e.g. Wilkins 1964; Libersky et al. 1993) uses the expression but we expect this to lead to differences from our adopted approach that are marginal. In this relation, J2 = 0.5   SαβSαβ and Y = Y(φ) are the second irreducible invariant of the deviatoric stress tensor and the shear strength, respectively.

3. Resolving inhomogeneous structures

3.1. Previous approach

In contrast to ductile media, brittle materials, such as basalt, granite, or porous pumice, do not rupture by plastic flow. This is because the material is not completely homogeneous but contains little flaws. These are little defects in the medium. With increasing strain, cracks develop originating from these flaws and start to pervade the solid body. In brittle media, stress is relieved by developing cracks.

To model damage, Grady & Kipp (1980) introduce a scalar parameter 0 ≤ D ≤ 1, where D = 0 represents undamaged and D = 1 fully disintegrated material. The stress tensor (Eq. (3)) is modified according to (14)Thus, damaged material may feel less stress than undamaged material. The local damage D is defined as the ratio of the volume defined by the growing crack to the volume in which the crack is growing (see Grady & Kipp 1980). The local time evolution of the damage D is given by (15)where nact accounts for crack accumulation and denotes the number of activated flaws, and Rs is the radius of a circumscribing sphere in which the crack evolves. It accounts for the maximum size of a damaged area. The quantity cg is the crack growth velocity. The concept of flaw activation originates from the idea that cracks do not start to grow for any strain applied but are activated for some strain threshold ϵact. The number of flaws n per unit volume with ϵact is given mostly by a power law (16)This distribution of flaws in a brittle material was proposed by Weibull (1939). It is based on two material parameters kwb and mwb, where kwb is the number of flaws per unit volume. The flaw distribution according to Eq. (16) is set as an initial condition for the material before the simulation.

A disadvantage is that the material parameters kwb and mwb are rarely available because they are difficult to measure. Unfortunately, small variations in kwb and mwb lead to large differences in the activation thresholds and simulation results (e.g. Schäfer 2005; Jutzi et al. 2009).

Because the damage model by Grady & Kipp (1980) is only suitable for brittle material, it is not applicable to pre-planetesimal material represented by SiO2 dust, which has properties that are between those of a ductile and a brittle material. Therefore, it is desirable to develop a simple damage approximation consistent with the applied porosity model. To summarise, three ingredients are essential for a damage model: (1) a distribution of defects, (2) the stress reduction due to damage, and (3) the typical size of a damaged area.

3.2. The approach of resolving inhomogeneous structures

In contrast to the damage model presented in Sect. 3.1, our approach does not explicitly consider the evolution of defects, hence is not strictly a damage model. However, it is designed to behave as an approximation to a damage model in some aspects, which are described in this section. As a significant advantage, our approach relies on quantities that can be more easily measured by laboratory experiments than the parameters of the model by Grady & Kipp (1980).

The starting point for our approach is the inhomogeneous nature of the material. Macroscopic dust aggregates created by the RBD method are not completely homogeneous (Güttler et al. 2009). The filling factor instead was found to follow a Gaussian distribution around a median of φμ ~ 0.15 with φσ = 0.013.

According to the porosity model of Sect. 2, regions of lower filling factor also represent regions of weaker compressive, tensile, and shear strengths, whereas regions of higher filling factor are stronger. Therefore, as an analogue to the Weibull distribution we propose an initial distribution of the filling factor given by the Gaussian function (17)where n(φ) is the number density for a filling factor φ, φσ is the width of the distribution function, and φμ is the median filling factor. The pressure limits of compressive, tensile, and shear strengths Σ(φ), T(φ), and Y(φ), respectively, mark the transition from the elastic to the plastic hydrostatic regime. Therefore, they can be regarded as analogues to the activation threshold ϵact. Hence, we can associate with each φ analogues to the activation threshold, which accounts for the first criterion of a damage model (see end of Sect. 3.1) that is based on the distribution of defects.

Damaged regions in this simple inhomogeneity scheme are therefore represented by areas of low filling factor. Under constant tension, φ becomes increasingly smaller with time in these regions and, owing to the φ-dependence of either the tensile strength T(φ) or shear strength Y(φ), the strength decreases in these regions as in Eq. (14) with increasing D. This fulfils the second criterion: the stress reduction.

Instead of artificially introducing a damage parameter D and an activation threshold ϵact, the filling factor φ naturally takes over this twofold role. On the one hand, it represents a quantity which determines the activation thresholds for the plastic regime, and on the other, it also represents a damage parameter that reduces the strength quantities. Therefore, in the inhomogeneity approach, the time evolution of the approximated damage is given by the time evolution of φ or equivalently the density. The propagation speed of the defects is, intrinsically, the sound speed.

We introduce a typical size Rc of an inhomogeneous clump to comply with the third criterion. In terms of this parameter, an absolute length scale is introduced into the inhomogeneity approach.

A great advantage of our approach is that the material parameters φσ, φμ, and Rc can be determined in laboratory measurements, e.g. by X-ray tomography (Güttler et al. 2009). However, no systematic study has been performed yet. A single measurement for a random ballistic deposition dust sample, which has the lowest filling factor (φ = 0.15) and the highest degree of homogeneity, has found a standard deviation of φσ ~ 0.013 (Fig. 10 in Güttler et al. 2009). There exist no data on either the typical clump size or clump size distribution. Therefore, as a first approximation, in Sect. 4 we study the effect of φσ on the outcome of pre-planetesimal collisions for a fixed clump size Rc. The influence of the latter will be studied in future work.

3.3. Implementation issues

The inhomogeneity is imposed on the simulated aggregate as an initial condition. Initially, we generate a homogeneous aggregate by assigning a constant initial filling factor φi. As a second step, φi is modified according to a Gaussian distribution with a standard deviation φσ, which we will simply call either the “inhomogeneity” or “degree of inhomogeneity” in this paper. For the distribution, we randomly select a particle a and set its density to a new following the Gaussian. The same is assigned to all particles within a radius of typical clump, which is fixed to Rc = 9.75   mm, the equivalent of one smoothing length h. Together with a spatial distance of 2.6 mm between the SPH particles, one clump is resolved by  ~290 particles, which is well above the required resolution (Geretshauser et al. 2010). Since there are no empirical data for the typical clump size yet, as a first approximation we assume that the clump size is constant. With the chosen size, the simulated spherical decimetre bodies can be regarded as being assembled from roughly spherical centimetre aggregates in sticking processes with different velocities similar to those described by Güttler et al. (2010) and Geretshauser et al. (2011). However, the effect of the typical clump size and resolution on the resulting fragment distribution will be studied in future work.

It is very likely that dust aggregates are compacted during their evolution in the protoplanetary disc (Zsom et al. 2010), in particular if the aggregates are assembled by smaller units in sticking processes. For this reason, we do not consider very fluffy aggregates with φi = 0.15 as produced by random ballistic deposition. We instead simulate spheres whose filling factor lies between the minimum (φmin = 0.12) and maximum (φmax = 0.58) filling factor of the compressive strength relation (Eq. (9)), i.e. φi = 0.35.

The mass of all individual SPH particles, which is a numerical parameter in the SPH equations, takes a constant value in the simulation. This value was computed to be consistent with the initial homogeneous aggregate (here φi = 0.35). Therefore, the modification of the density introduces a slight inconsistency in the SPH particle distribution. This means that the density value of an arbitrary particle differs slightly from the value computed from the particle number density using the standard SPH sum approach. Potentially, this could introduce pressure fluctuations causing spurious particle motion. However, in our algorithm the density is determined by solving the continuity equation and not by evaluating the number density of the particles (e.g. Monaghan 2005). As a consequence this inconsistency is marginal. Moreover, the time evolution of a inhomogeneous dust aggregate at rest was simulated and no spurious particle motion has been detected. In addition, the aggregate displayed no signs of instabilities within simulation times much longer than the collision timescale. Consequently, the presented approach is assumed to yield stable aggregates. Evaluating the generated SPH particle distribution yields the resulting filling factor distribution: the number density n(φ) of the volume fraction with filling factor φ is shown in Fig. 1 as a binned and cumulative plot for various φσ.

thumbnail Fig. 1

Initial filling factor distributions of the target and projectile as binned (top) and cumulative diagram (bottom). The number density n(φ) of volumes with filling factor φ is plotted over φ. A larger degree of inhomogeneity is characterised by the increasing standard deviation φσ of the Gaussian (Eq. (17)) around the initial filling factor φi of the homogeneous aggregate. Here, we choose φi = 0.35. The values of φσi can be found in Table 1.

Open with DEXTER

thumbnail Fig. 2

Interior and exterior view of the targets for different inhomogeneities with the filling factor colour coded. The inhomogeneity increases from the homogeneous target a) to the most inhomogeneous target f). In particular, the standard deviations are φσ = 0 a), φσ = 0.01 b), φσ = 0.02 c), φσ = 0.03 d), φσ = 0.04 e), and φσ = 0.05 f). All targets have a similar pattern that originates from considering interacting particle packages in the implementation. With increasing φσ the maximum and minimum filling factor of the different spots increase and decrease, respectively. The projectiles have a similar appearance. The result of collisions among these aggregates with 10   m   s-1 are depicted in Fig. 3.

Open with DEXTER

The inhomogeneous clumps with different filling factors are clearly visible in Fig. 2, where the density structure of the target is colour coded for different φσ. The packages are picked randomly but the random pattern is the same for all aggregates to ensure that the results are comparable. Owing to the increasing φσ, the regions with φ ≠ φi gain higher maximum and lower minimum filling factors. Consequently, the inhomogeneity pattern for φσ1 and φσ2 (see Table 1) is barely visible in Fig. 2 (first row, second and third). In contrast, for φσ5 the difference between the extreme values of the filling factor create a clearly visible pattern with high contrast.

The standard deviation of the Gaussian filling factor curve φσ is supplied to the distribution routine as an input parameter. The range of φσ begins slightly below the empirical value for a homogeneous random ballistic deposition aggregate (φσ = 0.013). There is no empirical data on the inhomogeneity of aggregates resulting from subsequent collisions. However, owing to the collisional history of an aggregate with macroscopic pre-planetesimals it is very likely that realistic aggregates have a higher degree of inhomogeneity. Therefore, we increase φσ in steps of 0.01 up to φσ = 0.050. The real standard deviation computed from the full width at half maximum (FWHM) is smaller than the φσ supplied to the SPH particle initialisation algorithm (algorithmic value). The FWHM values are listed in Table 1. For the sake of simplicity, we use the algorithmic value this paper.

Table 1

Selection of the standard deviations φσ as they are used in the inhomogeneity algorithm.

4. Simulation results

We now perform test simulations using our inhomogeneity approach. The initial setup is the same as in most simulations of Geretshauser et al. (2011). It consists of two colliding spheres with φi = 0.35. The radii of the target and projectile are rt = 10   cm and rp = 6   cm, and they consist of 238 238 and 51 477 SPH particles, respectively, with masses of 1.23 × 10-5   kg each. The SPH particles are set on a cubic grid separated by a distance of 2.6 mm. We consider two impact velocities of v0 = 10   m   s-1 and v0 = 12.5   m   s-1. For homogeneous aggregates, the smaller velocity leads to the sticking of the projectile to the target (see Fig. 3a). In contrast, the higher velocity results in fragmentation (see also Geretshauser et al. 2011). For these two collision velocities, the inhomogeneity of both the target and the projectile is defined by the standard deviation of the Gaussian φσ, which is varied from 0 to 0.05 yielding the distributions depicted in Fig. 1. The resulting targets are shown in Fig. 2. The projectiles display a similar density pattern.

thumbnail Fig. 3

Outcome of a collision between a target and projectile with the same φσ for different inhomogeneities. In all cases, the target and projectile radii are rt = 10   cm and rp = 6   cm, respectively, and the collision velocity is v0 = 10   m   s-1. Analogous to the initial targets in Fig. 2, the standard deviations for the Gaussian of the initial targets are φσ = 0 a), φσ = 0.01 b), φσ = 0.02 c), φσ = 0.03 d), φσ = 0.04 e), and φσ = 0.05 f). The collision outcome is shown in the impact direction. In the homogeneous case a) and for small inhomogeneities b), the target stays intact and forms one massive object with the projectile. For φσ ≥ 0.02, the target fragments c)–f). The fragment sizes decrease with increasing φσ and at the same time the number of fragments increases.

Open with DEXTER

For v0 = 10   m   s-1, the final fragment distribution is shown in Fig. 3 from the impact direction. For homogeneous aggregates (a) and the small standard deviation (φσ = 0.01, b), the target remains intact and the target and projectile form one large aggregate. For φσ = 0.01, some small fragments are visible. For φσ ≥ 0.02, the target fragments and with greater inhomogeneity the largest fragments of the distribution decrease in size.

The first impression from our visual inspection of the collision products is confirmed by the more detailed analysis. For this purpose, we utilise the four-population model (Geretshauser et al. 2011). According to this classification scheme of collision outcome, four classes of fragments are distinguished according to their mass: the largest fragment (index “1”) and the second largest fragment (index “2”) are characterised by single values of mass m, kinetic energy E, and other parameters. To account for the resolution limit, the sub-resolution population (index “sr”) is introduced, which contains all fragments consisting of a single SPH particle (mass mSPH). This class is mainly described by global quantities. Fragments with masses mSPH < mfrag < m2 form a class whose cumulative mass distribution is described by a power-law (index “pw”). This class is characterised by global quantities and distribution functions.

The evolution of the masses of each population with φσ is shown in Fig. 4 as a fraction of the total mass. The upper and the lower plots display the results for v0 = 10   m   s-1 and v0 = 12.5   m   s-1, respectively. The exact values can be found in Table 2. For the lower collision velocity and φσ ≤ 0.01, nearly all mass is stored in the largest fragment m1. For φσ > 0.01, a second largest fragment and a power-law population appears. With increasing inhomogeneity m1 rapidly decreases, while the mass of the second largest fragment m2 only slightly decreases. The masses m1 and m2 become comparable in the end. For high inhomogeneity values, most of the mass is stored in the power-law population mpw, as we discuss below. However, for high inhomogeneities mpw seems to remain fairly constant, while the mass of the sub-resolution population msr increases. For the higher collision value (bottom plot of Fig. 4), the evolution is similar but starts with mpw ≠ 0. The mass m1 constantly decreases, but not as rapidly as in the low velocity case. Moreover, m2 at first increases, then slightly decreases until the largest and second largest fragment become comparable. The mass of the power-law population rapidly increases for 0.01 ≤ φσ ≤ 0.04. For larger φσ, the power-law population only slightly increases, whereas msr increases at a higher rate. The mass evolution with increasing inhomogeneity can be summarised as follows: (1) inhomogeneity leads to fragmentation at collision velocities, at which homogeneous aggregates do not fragment, i.e. inhomogeneous aggregates are weaker; (2) for the two investigated velocities a larger φσ leads to a decrease in the mass of the largest and second largest fragments; (3) a higher degree of inhomogeneity causes an increase in the masses of both the power-law and sub-resolution population.

We turn to the analysis of the distribution of the residual total energy after the collision, which is shown in Fig. 5 as a fraction of the initial kinetic energy. The residual total energy is the sum of the translational energy and the energy stored in the internal degrees of freedom (e.g. rotation, vibration) of the respective populations. The tabulated values can be found in Table 2. In Fig. 5, the low velocity case with v0 = 10   m   s-1 is shown in the upper plot and the high velocity case with v0 = 12.5   m   s-1 in the lower plot. For the former, a large amount of energy E1 is stored in the largest fragment, but a considerable amount Esr is also stored in the sub-resolution population. As E1 decreases, the contribution of the second largest fragment at first increases, then decreases slightly with larger φσ. As for the mass, the energy contribution of the power-law population Epw increases with a higher degree of inhomogeneity. Remarkably, the total energy fraction of these three populations remains almost constant. In contrast, the energy contribution of the sub-resolution population Esr increases more strongly than its mass contribution in Fig. 4. The overall residual energy increases with a higher degree of inhomogeneity and this energy excess is stored mostly in the sub-resolution population. For the high velocity case (bottom plot of Fig. 5), this is even more evident: E1, E2, and Epw behave similar to the low velocity case but with E2 ≠ 0 and Epw ≠ 0 for φσ = 0. For a higher degree of inhomogeneity, less energy can again be dissipated and the total residual energy of the largest three populations remains nearly constant for increasing φσ. As a consequence, the residual energy excess is mainly stored in the sub-resolution population for large values of φσ. The conclusions of our energy analysis are: (1) with increasing impact velocity less energy can be dissipated by the system and the residual energy increases as previously found by Geretshauser et al. (2011); (2) with increasing inhomogeneity, the effectiveness of the energy dissipation decreases for the two investigated velocities; (3) the energy contribution of all populations show a similar dependence on φσ to their mass contributions: the fractions of energy stored in the largest and second largest fragment follow a decreasing trend with increasing inhomogeneity. The contribution of power-law and sub-resolution populations increase; (4) while the overall energy contribution of the largest fragment populations remains nearly constant with φσ, the fraction stored in the sub-resolution population drastically increases, such that the energy that cannot be dissipated owing to the increasing inhomogeneity is stored in rapidly moving single SPH particles.

thumbnail Fig. 4

The total masses of the four populations with inhomogeneities measured by φσ. The collision velocities are v0 = 10   m   s-1 (top) and v0 = 12.5   m   s-1 (bottom). The masses of the largest fragment m1, of the second largest fragment m2, of the power-law population mpw, and of the sub-resolution population msr are stacked and normalised by the total mass of the system mtot. In both cases, m1 and m2 decrease with increasing φσ. Conversely, mpw and msr increase. For the homogeneous case (φσ = 0), sticking is found for the low velocity collision, whereas for v0 = 12.5   m   s-1m2 ≠ 0 and mpw ≠ 0, which indicates fragmentation.

Open with DEXTER

thumbnail Fig. 5

The energy of each population for increasing inhomogeneity parameter φσ and collision velocity v0 = 10   m   s-1 (top) and v0 = 12.5   m   s-1 (bottom). The energy of each population is the sum of translational, rotational and vibrational energy and normalised by the initial kinetic energy Einit. The two latter energy contributions are negligible for the head-on collisions presented here. Owing to the strong correlation with mass, this diagram follows a similar trend to Fig. 4. The energy contributions are divided up into E1, E2, Epw, and Esr for the largest and second largest fragment, the power-law and sub-resolution population, respectively. In both velocity cases, both E1 and E2 decrease with increasing φσ. In contrast, Epw and Esr increase. In particular, the high fraction of kinetic energy stored in the sub-resolution population is remarkable.

Open with DEXTER

In comparison to the initial filling factor distribution of Fig. 1, the final filling factor distribution of the cumulated mass is shifted towards higher filling factors. Figure 6 reveals that in most cases more than 85% of the fragment mass is stored in filling factors  ≳ 42% (except for v0 = 12.5   m   s-1 and φσ = 0.050). With the same exception, the distributions look quite similar in both velocity cases. There is a slight tendency for there to be more mass at lower filling factors in the collision with v0 = 12.5m   s-1. From Fig. 6, it can also be seen that the broader the initial filling factor distribution, the broader the final filling factor distribution. For the v0 = 10.0   m   s-1 case, the two most homogeneous aggregates hardly produce any fragments. This causes a jump in the curves without any intermediate steps. The distribution for the highest collision velocity (v0 = 12.5   m   s-1) and the highest degree of inhomogeneity (φσ = 0.050) is particularly interesting since  ~30% of the fragment mass possess filling factors  <0.35, which is a larger fraction than for the initial aggregate. This implies that rarefaction occurs during the collision. On the microscopic scale, this is equivalent to monomer re-ordering and stretching of monomer chains. We note that a tiny fraction of particles acquire overdensities (φ > 1) in the collision. In the most violent case, this applies to single SPH particles, which in total represent 0.03% of the total mass of the system, which is negligible. These overdensities indicate that the monomers are breaking up on the microscopic scale.

We have so far considered all members of the four-population model. We now constrain our analysis to the power-law population. Its mass distribution is shown in the cumulative plot of Fig. 7 again for v0 = 10   m   s-1 (top) and v0 = 12.5   m   s-1 (bottom). The figure shows the ratio of the cumulative mass mcum to the fragment mass mfrag, which are both normalised by the total mass of the power-law population mpw given in Table 2. The cumulative distribution can be accurately (e.g. Güttler et al. 2010) described by a power law (18)where n(m)   dm is the number density of fragments of mass m, which is normalised by the total mass of the power-law population mpw. The quantity μpw is the normalised mass of the most massive member of the power-law population and κ is the power-law index or fragmentation parameter. We fitted the cumulative mass distributions of Fig. 7 for all fragment masses with Eq. (18). Finding an explanation for the kink in the cumulative mass curves is the subject of ongoing research. This deviation from a power law might be either a physical or a resolution effect (see also Geretshauser et al. 2011). The results are listed in Table 3. Figure 8 shows that the power-law index κ decreases slightly as the degree of inhomogeneity increases for both collision velocities, which leads to increasingly shallower slopes in Fig. 7 owing to the production of smaller fragments with increasing φσ. In addition, the κ values for the higher impact velocity (v0 = 12.5   m   s-1) are systematically below the κ values of the lower velocity (v0 = 10.0   m   s-1), except for small values of φσ for which there is only a small number of fragments (see Table 2) and hence the quality of the statistics is insufficient. This finding indicates that a larger fraction of less massive fragments are produced at higher impact velocities, which is expected.

thumbnail Fig. 6

Cumulative mass of the fragment masses over their average filling factor φ. The mass is normalised by the total fragment mass mtot. The collision velocities are v0 = 10   m   s-1 (top) and v0 = 12.4   m   s-1 (bottom). Compared to the initial distribution of filling factors (Fig. 1), the curves are shifted significantly towards higher filling factors. This reflects the compression taking place during the collision.

Open with DEXTER

thumbnail Fig. 7

Cumulative mass distribution of a power-law population for different inhomogeneity parameters φσ and collision velocities v0 = 10   m   s-1 (top) v0 = 12.5   m   s-1 (bottom). In both velocity cases both the mass of the largest member of the power-law population and the power-law slope decrease. This indicates the fragmentation to smaller aggregates for increasing φσ. For equal φσ, the slopes are shallower for the higher velocity. The power-law fit parameters can be found in Table 3.

Open with DEXTER

As shown in Fig. 9, the mass of the largest member of the power-law population, given by μpw, at first increases for larger values of φσ. It has a maximum at φσ ~ 0.02 for v0 = 10   m   s-1 and at φσ ~ 0.01 for v0 = 12.5   m   s-1, and then decreases in an exponential fashion. When comparing both velocity cases (see Table 2), we find that more fragments are produced in the high velocity case. This was also observed by Geretshauser et al. (2011). The influence of the inhomogeneity on the power-law population can be summarised as follows: (1) the overall mass of the power-law distribution is increased as discussed above; (2) a higher degree of inhomogeneity leads to the production of smaller fragments; (3) the largest member of the power-law population reaches its highest mass at small φσ, after which its mass decreases with increasing φσ.

thumbnail Fig. 8

Fragmentation parameter κ of the power-law population. For increasing inhomogeneity, indicated by the standard deviation φσ, the fragmentation parameter decreases. This indicates shallower slopes of the power-law mass distribution which is interpreted as an increasing fraction of smaller fragments. For φσ ≥ 0.02, the higher impact velocity results in smaller values of κ.

Open with DEXTER

thumbnail Fig. 9

Fragmentation strength of the power-law population. The quantity μpw represents the normalised mass of the most massive member of the power-law population and serves as an indicator of the fragmentation strength. For increasing inhomogeneity, represented by larger values of φσ, μpw at first increases and then decreases again. The increase may be caused by our low quality statistics owing to the small number of fragments.

Open with DEXTER

Table 2

Simulation results presented according to the four-population classification.

Table 3

Fit values of the power-law population.

5. Discussion and outlook

We have presented the first simulations of inhomogeneous porous pre-planetesimals, where the inhomogeneity of SiO2 dust aggregates is explicitly resolved. This inhomogeneity can in principal also be determined in laboratory experiments (Güttler et al. 2009). The approach is based on the concept that according to our porosity model (Geretshauser et al. 2010) inhomogeneities in filling factor cause fluctuations in the compressive, shear, and tensile strength in the aggregate. These fluctuations can be regarded as flaws in the material. In contrast to damage models designed for brittle material (Grady & Kipp 1980), we did not explicitly evolve the propagation of these flaws. Instead, the defects in the dust material, which behaves more like a fluid, were approximated by the time evolution of the filling factor or – equivalently – the density. The inhomogeneity of an aggregate was assumed in the initial SPH particle distribution by using a Gaussian distribution to model the filling factor. The inhomogeneity is characterised by the standard deviation φσ of the Gaussian and the typical radius Rc of a clump.

Using this inhomogeneity approach, we have performed test simulations of collisions between dust aggregates of intermediate porosity with fixed typical clump size. Two collision velocities have been chosen: one below and one above the fragmentation threshold for homogeneous aggregates. The results have been analysed using the four-population model of Geretshauser et al. (2011). For the lower collision velocity, we have shown that inhomogeneity leads to fragmentation. For both velocities, the masses of largest and second largest fragment are lower for a higher degree of inhomogeneity whereas the masses of the power-law and sub-resolution population are higher. Focussing on the power-law population, the number of fragments and the fraction of small fragments increase with increasing φσ. These findings show that the inhomogeneity approach has passed the first test in explicitly resolving the heterogeneity of porous objects. In a future study, we will address the influence of the typical clump size on the outcome of pre-planetesimal collisions.

Our findings indicate that inhomogeneous dust aggregates are weaker than their homogeneous equivalents. A slight inhomogeneity is sufficient to cause a catastrophic disruption instead of growth as the result of a dust aggregate collision. Therefore, inhomogeneity might explain the higher velocity thresholds for fragmentation in simulations than the a lower value found in laboratory experiments (Geretshauser et al. 2011).

Furthermore, macroscopic dust aggregates in protoplanetary discs are produced by the subsequent impacts of smaller aggregates at different impact velocities, i.e. pre-planetesimals have a collision history, thus are very likely to be inhomogeneous. This has a negative effect on planetesimal formation by subsequent collisions because pre-planetesimals might become weaker as the number of collisions increases. As a result, their threshold velocity for catastrophic disruption might be smaller. However, in a future study we will investigate the influence of the clump size of the stability of the pre-planetesimals. If the clumps are small relative to the pre-planetesimal size, this might increase the stability of the aggregates.

With the model presented in this paper, we are able to include these features. Further studies might be carried out that investigate the inhomogeneity created by subsequent multiple impacts. The result could be classified according to the inhomogeneity model that we have presented here.

To date it has been found that the simulations of collisions between homogeneous aggregates carried out by means of the porosity model by Geretshauser et al. (2010) depend only on the filling factor and produce the same fragment distribution for an identical set of input parameters. By randomly assigning an inhomogeneity pattern it is now possible to profoundly investigate the statistics of a fragment distribution. Statistical fluctuations in the quantities of the four-population model can now be estimated in particular for simulations with a small number of fragments.

In high velocity grazing collisions, the target begins to rotate. For highly porous aggregates, which have a low tensile strength, high spinning rates tend to lead to fragmentation of the target. With increasing inhomogeneity, this might also be true for aggregates with medium and low porosity. A quantitative investigation of this effect can be carried out by means of the presented approach.

The filling factor distribution of dust aggregates can be determined in the laboratory by X-ray tomography measurements (Güttler et al. 2009). These empirical data can be directly implemented into our inhomogeneity model. In addition the successful material calibration (Geretshauser et al. 2010), the obtained inhomogeneity parameters can be used to improve the realistic simulation of porous dust aggregates.

Acknowledgments

R.J.G. and R.S. wish to thank S. Arena, W. Benz, and F. Meru for many illuminating discussions. We are obliged to our referee Gareth Collins for his profound revision and constructive comments which helped to improve and clarify the manuscript. The SPH simulations were performed on the university cluster of the computing centre of Tübingen, the bwGriD clusters in Karlsruhe, Stuttgart, and Tübingen. Computing time was also provided by the High Performance Computing Centre Stuttgart (HLRS) on the national supercomputer NEC Nehalem Cluster under project grant SPH-PPC/12848. This project was funded by the Deutsche Forschungsgemeinschaft within the Forschergruppe “The Formation of Planets: The Critical First Growth Phase” under grant Kl 650/8-1.

References

Appendix A: List of symbols

Table A.1

List of symbols.

All Tables

Table 1

Selection of the standard deviations φσ as they are used in the inhomogeneity algorithm.

Table 2

Simulation results presented according to the four-population classification.

Table 3

Fit values of the power-law population.

Table A.1

List of symbols.

All Figures

thumbnail Fig. 1

Initial filling factor distributions of the target and projectile as binned (top) and cumulative diagram (bottom). The number density n(φ) of volumes with filling factor φ is plotted over φ. A larger degree of inhomogeneity is characterised by the increasing standard deviation φσ of the Gaussian (Eq. (17)) around the initial filling factor φi of the homogeneous aggregate. Here, we choose φi = 0.35. The values of φσi can be found in Table 1.

Open with DEXTER
In the text
thumbnail Fig. 2

Interior and exterior view of the targets for different inhomogeneities with the filling factor colour coded. The inhomogeneity increases from the homogeneous target a) to the most inhomogeneous target f). In particular, the standard deviations are φσ = 0 a), φσ = 0.01 b), φσ = 0.02 c), φσ = 0.03 d), φσ = 0.04 e), and φσ = 0.05 f). All targets have a similar pattern that originates from considering interacting particle packages in the implementation. With increasing φσ the maximum and minimum filling factor of the different spots increase and decrease, respectively. The projectiles have a similar appearance. The result of collisions among these aggregates with 10   m   s-1 are depicted in Fig. 3.

Open with DEXTER
In the text
thumbnail Fig. 3

Outcome of a collision between a target and projectile with the same φσ for different inhomogeneities. In all cases, the target and projectile radii are rt = 10   cm and rp = 6   cm, respectively, and the collision velocity is v0 = 10   m   s-1. Analogous to the initial targets in Fig. 2, the standard deviations for the Gaussian of the initial targets are φσ = 0 a), φσ = 0.01 b), φσ = 0.02 c), φσ = 0.03 d), φσ = 0.04 e), and φσ = 0.05 f). The collision outcome is shown in the impact direction. In the homogeneous case a) and for small inhomogeneities b), the target stays intact and forms one massive object with the projectile. For φσ ≥ 0.02, the target fragments c)–f). The fragment sizes decrease with increasing φσ and at the same time the number of fragments increases.

Open with DEXTER
In the text
thumbnail Fig. 4

The total masses of the four populations with inhomogeneities measured by φσ. The collision velocities are v0 = 10   m   s-1 (top) and v0 = 12.5   m   s-1 (bottom). The masses of the largest fragment m1, of the second largest fragment m2, of the power-law population mpw, and of the sub-resolution population msr are stacked and normalised by the total mass of the system mtot. In both cases, m1 and m2 decrease with increasing φσ. Conversely, mpw and msr increase. For the homogeneous case (φσ = 0), sticking is found for the low velocity collision, whereas for v0 = 12.5   m   s-1m2 ≠ 0 and mpw ≠ 0, which indicates fragmentation.

Open with DEXTER
In the text
thumbnail Fig. 5

The energy of each population for increasing inhomogeneity parameter φσ and collision velocity v0 = 10   m   s-1 (top) and v0 = 12.5   m   s-1 (bottom). The energy of each population is the sum of translational, rotational and vibrational energy and normalised by the initial kinetic energy Einit. The two latter energy contributions are negligible for the head-on collisions presented here. Owing to the strong correlation with mass, this diagram follows a similar trend to Fig. 4. The energy contributions are divided up into E1, E2, Epw, and Esr for the largest and second largest fragment, the power-law and sub-resolution population, respectively. In both velocity cases, both E1 and E2 decrease with increasing φσ. In contrast, Epw and Esr increase. In particular, the high fraction of kinetic energy stored in the sub-resolution population is remarkable.

Open with DEXTER
In the text
thumbnail Fig. 6

Cumulative mass of the fragment masses over their average filling factor φ. The mass is normalised by the total fragment mass mtot. The collision velocities are v0 = 10   m   s-1 (top) and v0 = 12.4   m   s-1 (bottom). Compared to the initial distribution of filling factors (Fig. 1), the curves are shifted significantly towards higher filling factors. This reflects the compression taking place during the collision.

Open with DEXTER
In the text
thumbnail Fig. 7

Cumulative mass distribution of a power-law population for different inhomogeneity parameters φσ and collision velocities v0 = 10   m   s-1 (top) v0 = 12.5   m   s-1 (bottom). In both velocity cases both the mass of the largest member of the power-law population and the power-law slope decrease. This indicates the fragmentation to smaller aggregates for increasing φσ. For equal φσ, the slopes are shallower for the higher velocity. The power-law fit parameters can be found in Table 3.

Open with DEXTER
In the text
thumbnail Fig. 8

Fragmentation parameter κ of the power-law population. For increasing inhomogeneity, indicated by the standard deviation φσ, the fragmentation parameter decreases. This indicates shallower slopes of the power-law mass distribution which is interpreted as an increasing fraction of smaller fragments. For φσ ≥ 0.02, the higher impact velocity results in smaller values of κ.

Open with DEXTER
In the text
thumbnail Fig. 9

Fragmentation strength of the power-law population. The quantity μpw represents the normalised mass of the most massive member of the power-law population and serves as an indicator of the fragmentation strength. For increasing inhomogeneity, represented by larger values of φσ, μpw at first increases and then decreases again. The increase may be caused by our low quality statistics owing to the small number of fragments.

Open with DEXTER
In the text

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

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

Initial download of the metrics may take a while.