© 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 SiO₂ 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 SiO₂ 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 SiO₂ 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 $\begin{array}{ccc}\frac{\mathrm{d}\mathit{\rho }}{\mathrm{d}\mathit{t}}& \mathrm{=}& \mathrm{-}\mathit{\rho }\frac{\mathit{\partial }{\mathit{v}}_{\mathit{\alpha }}}{\mathit{\partial }{\mathit{x}}_{\mathit{\alpha }}}\mathit{,}\\ \frac{\mathrm{d}{\mathit{v}}_{\mathit{\alpha }}}{\mathrm{d}\mathit{t}}& \mathrm{=}& \frac{\mathrm{1}}{\mathit{\rho }}\frac{\mathit{\partial }{\mathit{\sigma }}_{\mathit{\alpha \beta }}}{\mathit{\partial }{\mathit{x}}_{\mathit{\beta }}}\mathit{,}\end{array}$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 |x^a − x^b| and is weighted by the kernel function W(|x^a − x^b|,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, ${\mathit{\sigma }}_{\mathit{\alpha \beta }}\mathrm{=}\mathrm{-}\mathit{p}{\mathit{\delta }}_{\mathit{\alpha \beta }}\mathrm{+}{\mathit{S}}_{\mathit{\alpha \beta }}\mathit{.}$(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 ${\mathit{ϵ}}^{\mathit{\alpha \beta }}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}}\left(\frac{\mathit{\partial }{{\mathit{x}}^{\mathrm{\prime }}}^{\mathit{\alpha }}}{\mathit{\partial }{\mathit{x}}^{\mathit{\beta }}}\mathrm{+}\frac{\mathit{\partial }{{\mathit{x}}^{\mathrm{\prime }}}^{\mathit{\beta }}}{\mathit{\partial }{\mathit{x}}^{\mathit{\alpha }}}\right)\mathit{,}$(4)where x′ are the coordinates of the deformed body. In particular, the relations that hold are $\begin{array}{ccc}\mathit{p}& \mathrm{=}& \\ {\mathit{S}}^{\mathit{\alpha \beta }}& \mathrm{=}& \end{array}$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 $\mathit{\phi }\mathrm{=}\mathrm{1}\mathrm{-}\mathrm{\Phi }\mathrm{=}\frac{\mathit{\rho }}{{\mathit{\rho }}_{\mathrm{s}}}\mathit{,}$(7)where ρ_s is the density of the matrix material. We use a constant value of ρ_s = 2000   kg   m^-3 for our SiO₂ material (e.g. Blum & Schräpler 2004). In particular, the elastic quantities are related to the filling factor by a power law $\mathit{K}\mathrm{\left(}\mathit{\phi }\mathrm{\right)}\mathrm{=}\mathrm{2}\mathit{\mu }\mathrm{\left(}\mathit{\phi }\mathrm{\right)}\mathrm{=}{\mathit{K}}_{\mathrm{0}}{\left(\frac{\mathit{\phi }}{{\mathit{\phi }}_{\mathrm{RBD}}}\right)}^{\mathit{\gamma }}\mathit{,}$(8)where γ = 4, K₀ 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 K₀ = 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, $\mathrm{\Sigma }\mathrm{\left(}\mathit{\phi }\mathrm{\right)}\mathrm{=}{\mathit{p}}_{\mathrm{m}}{\left(\frac{{\mathit{\phi }}_{\max }\mathrm{-}{\mathit{\phi }}_{\min }}{{\mathit{\phi }}_{\max }\mathrm{-}\mathit{\phi }}\mathrm{-}\mathrm{1}\right)}^{\mathrm{\Delta }\hspace{0.17em}\ln \mathrm{10}}\mathit{,}$(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 p_m = 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 p_m 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 $\mathit{T}\mathrm{\left(}\mathit{\phi }\mathrm{\right)}\mathrm{=}\mathrm{-}{\mathrm{10}}^{\mathrm{2.8}\mathrm{+}\mathrm{1.48}\mathit{\phi }} \mathrm{Pa}$(10)and the shear strength is the geometric mean of the compressive and tensile strength $\mathit{Y}\mathrm{\left(}\mathit{\phi }\mathrm{\right)}\mathrm{=}\sqrt{\mathrm{\Sigma }\mathrm{\left(}\mathit{\phi }\mathrm{\right)}\mathrm{\left|}\mathit{T}\mathrm{\right(}\mathit{\phi }\mathrm{\left)}\mathrm{\right|}}\mathit{.}$(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 $\mathit{p}\mathrm{\left(}\mathit{\phi }\mathrm{\right)}\mathrm{=}\{\begin{array}{c}\\ \mathrm{\Sigma }\mathrm{\left(}\mathit{\phi }\mathrm{\right)}& \mathrm{for} {\mathit{\phi }}_{\mathrm{c}}^{\mathrm{+}}\mathit{<}\mathit{\phi }\\ \mathit{K}\mathrm{\left(}{\mathit{\phi }}_{\mathrm{0}}^{\mathrm{\prime }}\mathrm{\right)}\mathrm{(}\mathit{\phi }\mathit{/}{\mathit{\phi }}_{\mathrm{0}}^{\mathrm{\prime }}\mathrm{-}\mathrm{1}\mathrm{)}& \mathrm{for} {\mathit{\phi }}_{\mathrm{c}}^{\mathrm{-}}\mathrm{\le }\mathit{\phi }\mathrm{\le }{\mathit{\phi }}_{\mathrm{c}}^{\mathrm{+}}\\ \mathit{T}\mathrm{\left(}\mathit{\phi }\mathrm{\right)}& \mathrm{for} \mathit{\phi }\mathit{<}{\mathit{\phi }}_{\mathrm{c}}^{\mathrm{-}}\mathit{,}\end{array}$(12)where ${\mathit{\phi }}_{\mathrm{c}}^{\mathrm{+}}\mathit{>}{\mathit{\phi }}_{\mathrm{c}}^{\mathrm{-}}$, and ${\mathit{\phi }}_{\mathrm{c}}^{\mathrm{+}}$ and ${\mathit{\phi }}_{\mathrm{c}}^{\mathrm{-}}$ are critical filling factors. The value of ${\mathit{\phi }}_{\mathrm{c}}^{\mathrm{+}}$ marks the transition between elastic and plastic compression and ${\mathit{\phi }}_{\mathrm{c}}^{\mathrm{-}}$ defines the transition between elastic and plastic tension. The filling factors in-between the critical values represent the elastic regime. The quantity ${\mathit{\phi }}_{\mathrm{0}}^{\mathrm{\prime }}$ 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 ${\mathit{\phi }}_{\mathrm{c}}^{\mathrm{+}}$ is exceeded, the pressure follows the much shallower compressive strength path Σ(φ) in the positive φ-direction, whereas if ${\mathit{\phi }}_{\mathrm{c}}^{\mathrm{-}}$ 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) ${\mathit{S}}_{\mathit{\alpha \beta }}\mathrm{\to }\mathit{f}{\mathit{S}}_{\mathit{\alpha \beta }}\mathit{.}$(13)The function f is defined by $\mathit{f}\mathrm{=}\min {}^{\mathrm{[}}\mathit{Y}^{\mathrm{2}}\mathit{/}\mathrm{3}{\mathit{J}}_{\mathrm{2}}\mathit{,}{\mathrm{1}}^{\mathrm{]}}$. An alternative approach (e.g. Wilkins 1964; Libersky et al. 1993) uses the expression $\mathit{f}\mathrm{=}\min \left[\sqrt{{\mathit{Y}}^{\mathrm{2}}\mathit{/}\mathrm{3}{\mathit{J}}_{\mathrm{2}}}\mathit{,}\mathrm{1}\right]$ but we expect this to lead to differences from our adopted approach that are marginal. In this relation, J₂ = 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 ${\mathit{\sigma }}_{\mathrm{dam}}\mathrm{=}\mathrm{(}\mathrm{1}\mathrm{-}\mathit{D}\mathrm{)}\mathit{\sigma }\mathit{.}$(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 $\frac{\mathrm{d}{\mathit{D}}^{\mathrm{1}\mathit{/}\mathrm{3}}}{\mathrm{d}\mathit{t}}\mathrm{=}{\mathit{n}}_{\mathrm{act}}\frac{{\mathit{c}}_{\mathrm{g}}}{{\mathit{R}}_{\mathrm{s}}}\mathit{,}$(15)where n_act accounts for crack accumulation and denotes the number of activated flaws, and R_s is the radius of a circumscribing sphere in which the crack evolves. It accounts for the maximum size of a damaged area. The quantity c_g 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 $\mathit{n}\mathrm{\left(}{\mathit{ϵ}}_{\mathrm{act}}\mathrm{\right)}\mathrm{=}{\mathit{k}}_{\mathrm{wb}}{\mathit{ϵ}}_{\mathrm{act}}^{{\mathit{m}}_{\mathrm{wb}}}\mathit{.}$(16)This distribution of flaws in a brittle material was proposed by Weibull (1939). It is based on two material parameters k_wb and m_wb, where k_wb 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 k_wb and m_wb are rarely available because they are difficult to measure. Unfortunately, small variations in k_wb and m_wb 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 SiO₂ 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 $\mathit{n}\mathrm{\left(}\mathit{\phi }\mathrm{\right)}\mathrm{=}\frac{\mathrm{1}}{\sqrt{\mathrm{2}\mathit{\pi }}{\mathit{\phi }}_{\mathit{\sigma }}}\exp \left(\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}\frac{\mathit{\phi }\mathrm{-}{\mathit{\phi }}_{\mathit{\mu }}}{{\mathit{\phi }}_{\mathit{\sigma }}}\right)\mathit{,}$(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 R_c 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 R_c 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 R_c. 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 ${\begin{array}{}{\mathrm{˜}}_{}\\ {\mathit{\phi }}_{\mathrm{i}}\end{array}}^{\mathit{a}}$ following the Gaussian. The same ${\begin{array}{}{\mathrm{˜}}_{}\\ {\mathit{\phi }}_{\mathrm{i}}\end{array}}^{\mathit{a}}$ is assigned to all particles within a radius of typical clump, which is fixed to R_c = 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 φ_σ.

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.

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.

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.

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 r_t = 10   cm and r_p = 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 v₀ = 10   m   s^-1 and v₀ = 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.

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.

For v₀ = 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 m_SPH). This class is mainly described by global quantities. Fragments with masses m_SPH < m_frag < m₂ 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 v₀ = 10   m   s^-1 and v₀ = 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 m₁. For φ_σ > 0.01, a second largest fragment and a power-law population appears. With increasing inhomogeneity m₁ rapidly decreases, while the mass of the second largest fragment m₂ only slightly decreases. The masses m₁ and m₂ become comparable in the end. For high inhomogeneity values, most of the mass is stored in the power-law population m_pw, as we discuss below. However, for high inhomogeneities m_pw seems to remain fairly constant, while the mass of the sub-resolution population m_sr increases. For the higher collision value (bottom plot of Fig. 4), the evolution is similar but starts with m_pw ≠ 0. The mass m₁ constantly decreases, but not as rapidly as in the low velocity case. Moreover, m₂ 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 m_sr 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 v₀ = 10   m   s^-1 is shown in the upper plot and the high velocity case with v₀ = 12.5   m   s^-1 in the lower plot. For the former, a large amount of energy E₁ is stored in the largest fragment, but a considerable amount E_sr is also stored in the sub-resolution population. As E₁ 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 E_pw 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 E_sr 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: E₁, E₂, and E_pw behave similar to the low velocity case but with E₂ ≠ 0 and E_pw ≠ 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.

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.

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.

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 v₀ = 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 v₀ = 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 v₀ = 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 (v₀ = 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 v₀ = 10   m   s^-1 (top) and v₀ = 12.5   m   s^-1 (bottom). The figure shows the ratio of the cumulative mass m_cum to the fragment mass m_frag, which are both normalised by the total mass of the power-law population m_pw given in Table 2. The cumulative distribution can be accurately (e.g. Güttler et al. 2010) described by a power law ${\mathit{m}}_{\mathrm{cum}}\mathrm{\left(}{\mathit{m}}_{\mathrm{frag}}\mathrm{\right)}\mathrm{=}{\mathrm{\int }}_{\mathrm{0}}^{{\mathit{m}}_{\mathrm{frag}}}\mathit{n}\mathrm{\left(}\mathit{m}\mathrm{\right)}\hspace{0.17em}\mathit{m}\hspace{0.17em}\mathrm{d}\mathit{m}\mathrm{=}{\left(\frac{{\mathit{m}}_{\mathrm{frag}}}{{\mathit{\mu }}_{\mathrm{pw}}}\right)}^{\mathit{\kappa }}\mathit{,}$(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 m_pw. 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 (v₀ = 12.5   m   s^-1) are systematically below the κ values of the lower velocity (v₀ = 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.

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.

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.

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 v₀ = 10   m   s^-1 and at φ_σ ~ 0.01 for v₀ = 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 φ_σ.

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 κ.

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.

5. Discussion and outlook

We have presented the first simulations of inhomogeneous porous pre-planetesimals, where the inhomogeneity of SiO₂ 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 R_c 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

All Tables

All Figures

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.

In the text

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.

In the text

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.

In the text

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.

In the text

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.

In the text

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.

In the text

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.

In the text

	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 κ.
In the text

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.

In the text