© The Authors 2023

1 Introduction

Agglomeration of dust particles in accretion disks around young stars is the first step toward the formation of planets in a planetary system (Weidenschilling 1977; Weidenschilling & Cuzzi 1993; Wurm & Blum 1998; Blum 2006; Blum et al. 2006; Birnstiel et al. 2010a). To better understand the growth processes of dust aggregates, observational studies such as Birnstiel et al. (2010b) and laboratory experiments such as those reported by Blum & Wurm (2000), Blum et al. (2006), and Gundlach & Blum (2015) are important. In addition, numerical simulations based on granular-mechanics codes such as Dominik & Tielens (1997), Wada et al. (2007, 2008, 2009), Ringl et al. (2012), and Umstätter & Urbassek (2021) were performed to understand the collision behavior of granular aggregates. Much of this simulational work by Umstätter & Urbassek (2020, 2021), Wada et al. (2007, 2008), for example, focuses on head-on collisions, and only relatively fewer studies (Wada et al. 2009; Ringl et al. 2012; Hasegawa et al. 2021) consider noncentral collisions. The recent work by Arakawa et al. (2022b) and Arakawa et al. (2022a) focuses on identifying the role of viscous dissipation and rolling friction on growth and fragmentation in aggregate collisions. Wada et al. (2013) and Hasegawa et al. (2021) extend the simulations to include the case of nonequal masses of the colliding aggregates.

Previous simulational studies obtained relevant results on the outcomes of collisions between granular aggregates. Paszun & Dominik (2009) showed how aggregate collisions change the form and filling factor of aggregates and found that the size distribution of ejecta follows a power law. At the same time, Wada et al. (2009) studied fractal aggregates containing up to around 8000 grains and demonstrated that these can survive and even grow by collisions at relatively high velocities. Later, Wada et al. (2013) extended these studies to collisions between aggregates of unequal masses and showed that for highly unequal mass ratios, the tendency for growth is even increased. Hasegawa et al. (2021) used aggregates containing up to 260,000 grains and examined the role of mass transfer between the colliding aggregates.

Collisions may lead to aggregate growth by the merging of colliding aggregates or to aggregate destruction by fragmentation of the aggregates. Such processes change the size distribution of aggregate ensembles, such as protoplanetary disks, with time and the information obtained from simulational studies can hence be used to investigate the evolution of protoplanetary disks. Monte-Carlo-based research such as by Ormel et al. (2009) and Drazkowska et al. (2013) has been devoted to study the dust evolution and requires information on the outcome of aggregate collisions over a large portion of the pertinent parameter space: collision velocity, impact parameter, aggregate porosity, etc.

This present study aims to identify the collisional outcome of silica dust grains as a function of these parameters to provide insight into the collision dynamics. Our study supplements the results of previous work – such as by Wada et al. (2009, 2013), Paszun & Dominik (2009) and Hasegawa et al. (2021) – by focusing on developing a systematics of collision outcomes as a function of the aggregate porosity. As a result, we derived a map of collision outcomes as a function of the impact parameter and collision velocity and calculated the impact-parameter-averaged velocity dependence of fusion and fragmentation probabilities. In addition, the dependence of these quantities on the aggregate porosity is evaluated.

2 Methods

2.1 Granular-mechanics model

The model employed in this work continues earlier work (Ringl & Urbassek 2012) and is based on the prolific work by Chokshi et al. (1993) and Dominik & Tielens (1995, 1996, 1997). We repeat some details here for completeness and specialize the formulae to the case of grains of equal radius r and identical material properties, which is relevant for the present study.

Two grains interact when their distance is smaller than the sum of their radii, that is to say their overlap, $$\delta =2r-\left|x_i-x_j\right|,$$(1)

is positive. Grains then are subject to a normal force along the line connecting their respective centers (i.e., along the normal direction n). This normal force consists of an elastic repulsion that follows Hertz’ law, a viscous drag force proportional to the normal component of their relative velocity and an adhesive force which we assume to be constant. Using Young’s modulus Y and the Poisson number ν of the material as well as the surface energy γ, the normal force is $$\begin{array}{l}{F}_n=-2\pi \gamma rn\hfill \\ +\max \left(0,\frac{2}{3}\frac{Y}{1-v^2}\sqrt{\frac{r}{2}}{\delta }^{3/2}-A\sqrt{\frac{\delta r}{2}}{\upsilon }_n\cdot n\right)n,\hfill \end{array}$$(2)

where the maximum function was used to avoid the frictional force leading to an attractive contribution (cf. Pöschel & Schwager 2005).

We furthermore defined the relative contact velocity $${\upsilon }_c={\upsilon }_i-{\upsilon }_j-\left(r-\frac{\delta }{2}\right)\left({\omega }_i+{\omega }_j\right)\times n,$$(3)

its tangential contribution $${\upsilon }_t={\upsilon }_c-\left({\upsilon }_c\times n\right)n,$$(4)

where the tangential direction is given by $$t=\frac{{\upsilon }_t}{\left|{\upsilon }_t\right|}.$$(5)

The relative angular velocity is $${\omega }_{\text{rel}}={\omega }_i-{\omega }_j,$$(6)

its twisting component $${\omega }_{\text{twist}}=\left({\omega }_{\text{rel}}\cdot n\right)n,$$(7)

and its rolling component $${\omega }_{\text{roll}}={\omega }_{\text{rel}}-{\omega }_{\text{twist}}.$$(8)

These definitions enabled us to define a sliding force $$F_{\text{slid}}=-t\min \left({\eta }_{\mathrm{tang}}\left|{\upsilon }_t\right|,\frac{Y}{2\left(1+v\right)}\frac{\delta r}{8\pi }\right),$$(9)

where the minimum function selects either a constant force according to Dominik & Tielens (1996, 1997) or a velocity-proportional force to avoid unphysical oscillations at small relative velocities (cf. Ringl & Urbassek 2012). A similar minimum function was already implemented by Haff & Werner (1986; Pöschel & Schwager 2005). The constant η_tang was chosen depending on the time step ∆t of our simulations, $${\eta }_{\mathrm{tang}}=0.05\frac{m}{\text{Δ}t},$$(10)

where the mass m of the grains was used. The constant branch in the equation of the sliding force, Eq. (9), is proportional to the square of the contact radius, $a=\sqrt{\delta r/2}$. The tangential sliding force led to a sliding torque $$d_{\text{slid}}=-rn\times F_{\text{slid}},$$(11)

where the minus sign had to be used, as n points from the center of grain i away from the point of contact. Rolling motion led to an asymmetry in the adhesive neck at the point of contact (Dominik & Tielens 1995, 1997). We assumed the rolling torque to be constant with the average torque for rolling over a critical inelastic distance given in Dominik & Tielens (1995), unless the relative angular velocity was too small (to avoid unphysical oscillations, as above). We used $$d_{\text{roll}}=-\min \left({\eta }_{\text{roll}}\left|{\omega }_{\text{roll}}\right|,4\pi \gamma r{\xi }_{\text{crit}}\right)\frac{{\omega }_{\text{roll}}}{\left|{\omega }_{\text{roll}}\right|},$$(12)

where the constant branch in the minimum function is given by the adhesive force and the critical rolling distance ξ_crit. The constant η_roll was chosen as $${\eta }_{\text{roll}}=0.4r^2{\eta }_{\mathrm{tang}}$$(13)

(Ringl & Urbassek 2012). Twisting motion also led to a torque, which again selects between the constant value given in Dominik & Tielens (1996) and Dominik & Tielens (1997) and a velocity-proportional term in the range of small velocities, $$d_{\text{twist}}=-\min \left({\eta }_{\text{twist}}\left|{\omega }_{\text{twist}}\right|,\frac{Y{\left(\delta r/2\right)}^{3/2}}{\left(1+v\right)12\pi }\right)\frac{{\omega }_{\text{twist}}}{\left|{\omega }_{\text{twist}}\right|},$$(14)

where η_twist = η_roll. The constant branch of the twisting torque is proportional to the third power of the contact radius a.

We note that our approach ignores tangential elastic forces, which might apply for small tangential displacements of contacting grains and which are included in codes such as those by Dominik & Nübold (2002), Wada et al. (2007), and Seizinger et al. (2012). These small displacements are believed to be relevant for aggregate restructuring in low-velocity collisions, but hardly at the high collision velocities that are needed for aggregate fragmentation. At high velocities, dissipative sliding friction – apart from normal viscoelastic energy dissipation – governs the processes occurring under fragmentation (Umstätter & Urbassek 2021). We note that Arakawa et al. (2022a) recently used the code developed by Wada et al. (2007) to study noncentral aggregate collisions. They find that – using our collision parameters – the processes for velocities >0.25 m s⁻¹ are governed by tangential friction, rather than elastic forces. Snapshots 1a–d in that paper do indeed correspond nicely with our results in Fig. 2 below for a velocity of 5 m s⁻¹, where we observe a breakup into three major fragments. The difference in our velocities and those used by Arakawa et al. (2022a) almost perfectly corresponds to the differences in the choice of parameters leading to a different critical sliding velocity (a factor 12.6). Furthermore, above a velocity of roughly 5 m s⁻¹, the relevance of fusion processes vanishes (Fig. 6); taking the factor of 12.6 difference in velocities into account, this corresponds again with Arakawa et al. (2022a) who find a transition from positive growth efficiency to negative growth efficiency between 50 and 60 m s⁻¹; readers can refer to Sect. 3.2 and Appendix B for more details.

2.2 Grains

For this work we chose grains of radius r = 0.76 μm with material parameters according to silica: γ = 0.025 J m⁻², Y = 54 GPa, ν = 0.17, and ρ = 2000 kg m⁻³ (Ringl & Urbassek 2012). We chose the critical rolling displacement according to the experiments of Heim et al. (1999) as ξ_crit = 32 Å and the damping constant in the viscous friction of the normal motion as A = 0.5 ns (Ringl & Urbassek 2012). Such grains relax to an equilibrium overlap of δ_eq = 3.01 Å. To separate two grains from equilibrium overlap, a relative velocity in the normal direction of υ = 0.1717 m s⁻¹ – the so-called fragmentation velocity, υ_fr – is necessary. We note that without viscous friction (defined by the parameter A), a lower velocity of υ = 0.1531 m s⁻¹ is sufficient to separate the grains.

2.3 Aggregate construction

Spherical aggregates are characterized by their radius R and the number N of grains they contain. The filling factor Φ of an aggregate is defined as the ratio of the volume filled by the grains, N(4π/3)r³, and the aggregate volume, (4π/3)R³; the aggregate porosity is 1 − Φ.

We constructed spherical aggregates as follows. We requested that a number N of grains be in the aggregate and that it has a filling factor Φ; this determined the radius R of the aggregate as $$R={\left(\frac{N}{\text{Φ}}\right)}^{1/3}r.$$(15)

We placed the first grain at the origin. Then we iterated the following procedure. We randomly picked a grain that had already been placed and attached the next grain to it in a random direction. If the new grain was wholly inside the aggregate radius, that is |x| ≤ R − r, and if it did not overlap with any other grain than the one it had been attached to, it was kept, or else deleted. This procedure was repeated until the total number of grains reached N. Such an aggregate has a coordination number of 2. The aggregates were then relaxed such that all the grains that are in contact with each other reached equilibrium overlap. We constructed aggregates with N = 20 000 grains each, with four different filling factors. The four aggregates are summarized in Table 1.

2.4 Collision setup

Figure 1 shows our basic setup of a collision. To start such a simulation, we took two identical relaxed aggregates centered around the origin and translated them into two dimensions to separate them along the direction of their collision and perpendicular to that direction to introduce a nonzero impact parameter B. For this work, we used the reduced impact parameter $$b=\frac{B}{R}.$$(16)

We defined the collision angle θ = sin⁻¹[B/(2R)] = sin⁻¹(b/2) via the impact parameter and the angle ψ of the collision plane to be the angle of the plane perpendicular to the line connecting the aggregate centers when they just touch. From geometry, we have $$\psi ={90}^{°}-\theta ={90}^{°}-{\sin }^{-1}\left(b/2\right).$$(17)

For each filling factor, we performed collisions for eight impact parameters b = 0−1.75 in steps of ∆b = 0.25 and nine relative velocities υ = 2.5, 3.75, 5, 6.25, 7.5, 8.75, 10, 12.5, and 15 m s⁻¹, leading to a total number of 288 collisions. We ran the simulations for a total duration of 200 μs with a time step of ∆t = 50 ps.

Fig. 1 Schematic diagram showing the moment when the two granular aggregates just touch. The aggregates are equal in size with radius R and collide with a relative velocity υ. The point of contact is denoted by “C.” The solid green line connects the centers of the aggregates. The dashed gray lines are drawn from the centers of the aggregates along the direction of motion. The impact parameter of the collision, B, is the perpendicular distance between the gray lines, shown by black arrows. The collision angle is denoted by θ. The dashed black line represents the edge-on view of the common tangent plane of the two aggregates on the point of contact; it is denoted as the “collision plane.” The angle between the collision plane and the direction of motion is denoted as ψ.

2.5 Detection of fragments

The two colliding aggregates generate a number of resulting fragments. We used the algorithm by Stoddard (1978); two grains with centers at x_i and x_j belong to the same cluster if |x_i − x_j| ≤ 2r. This allowed us to determine the post-collisional aggregates and the number of grains they contain. With N₁ we denote the number of grains in the largest aggregate and with N₂ the number of grains in the second-largest aggregate.

3 Results and discussion

3.1 Collision outcomes

Figure 2 gives a synopsis of the final collision outcomes for the simulations performed for the filling factor 0.28. Some velocities were left out in order to make the presentation more readable. The corresponding synopses for the other filling factors are included in Appendix A.

An inspection of the synopsis reveals that a variety of outcomes are possible.

1. The two aggregates can merge and form a united cluster (example: υ = 2.5 m s⁻¹, b = 0); this outcome is denoted as fusion, but it has also been given other names in the past such as agglomeration, coalescence, or merging.

2. Another outcome is the scattering of the two aggregates off each other (example: υ = 2.5 m s⁻¹, b = 1.75); this outcome happens at large impact parameters (typically b > 1) and is denoted as sliding. It is important to note that the two aggregates typically do not leave the collision intact but are surrounded by a fragment cloud which increases in mass for higher velocities. Furthermore, the two surviving aggregates keep more or less their initial direction.

3. We classify the remaining cases as fragmentation; other names such as erosion or shattering might be chosen here. In the fragmentation outcomes, there might be several large fragments, such as in υ = 7.5 m s⁻¹ and b = 0.5, or a multitude of monomer grains might be visible, such as in υ = 15 m s⁻¹ and b = 0.25; we do not subdivide these outcomes any further (cf. the discussion in Sect. 3.6).

The origin of the term sliding for the second collision class is shown in Fig. 3, which provides snapshots of the time evolution of the collision with b = 1.25 and υ = 5 m s⁻¹. In this collision, the two aggregates do indeed slide along the collision plane; sliding friction leads to a small amount of mass transfer between the aggregates as well as to the ejection of monomers and small fragments, but the major part of the noncollided parts of the aggregates move onward, on only slightly disturbed trajectories, giving rise to two large fragments.

We note that in laboratory and microgravity experiments, the more or less intact scattering repulsion of aggregates is observed even at small impact parameters (Weidling et al. 2009, 2012; Heißelmann et al. 2010; Jankowski et al. 2012; Weidling & Blum 2015; Schräpler et al. 2022), corresponding to the “backscattering” or bouncing of the two aggregates. Such a behavior has not been observed in the present study nor, to our knowledge, in other studies of aggregate collisions, except for unrealistically large filling factors (Wada et al. 2011) or in aggregates composed of grains of very different masses such as dust-covered chondrules (Umstätter et al. 2019). We note, however, that such a bouncing behavior is well-known for solid or liquid clusters (Kalweit & Drikakis 2006; Sommerfeld & Kuschel 2016; Nietiadi & Urbassek 2022).

We thus base our classification on only three collision classes: fusion, sliding, and fragmentation. This appears sufficient for the analysis of dust aggregate collisions since here mostly aggregate growth (incorporated in fusion) or destruction (synonymous with fragmentation) is relevant; the sliding class basically keeps the aggregate sizes unchanged. We note that collisions of liquid droplets – of interest in technological applications such as spraying – have been classified similarly; however, the fragmentation regime is denoted as “stretching separation,” where a multitude of small droplets are produced after the colliding droplets separate again (Ashgriz & Poo 1990; Estrade et al. 1999; Sommerfeld & Kuschel 2016). Also collisions of solid clusters, analyzed with the help of molecular dynamics simulations, lead to a similar classification (Kalweit & Drikakis 2006; Nietiadi & Urbassek 2022).

Fig. 2 Tableau of the final state after the collision of two aggregates of filling factor 0.28 colliding with velocity υ (in units of m s−1) and reduced impact parameter b. Grains were projected into the plotted plane in order to allow for a better perception of the monomer distribution. Grains are colored according to their initial aggregate affiliation. The size of the subpanels changes with the collision velocity.

Fig. 3 Series of snapshots showing the time evolution of a sliding collision: b = 1.25 and υ = 5 m s−1. Grains are colored according to their initial aggregate affiliation.

3.2 Collision classes

In this section, we set up an algorithm to classify collisions into the three classes qualitatively described in Sect. 3.1. We base the classification on the sizes of the two largest clusters found at the end of the collision, N₁ and N₂, where N₁ ≥ N₂. The clusters were identified with Stoddard’s cluster algorithm as described in Sect. 2.5.

1. The analysis of fusion outcomes is based on the size of the largest cluster, N₁. In Fig. 4a, a contour plot of N_fu = N₁/N in the b-υ plane is shown; it has been obtained by interpolation of the simulation data calculated at the discrete b and υ values where simulations have been performed. A comparison with the visual collision outcomes of Fig. 2 shows that a threshold value of 1.5 can indeed be used to identify the fusion regime, such that fusion is characterized by $$N_{\text{fu}}=N_1/N\ge \mathrm{1.5.}$$(18)

   A threshold value of 1.5 means that 75% of all the grains initially present in the two aggregates have combined into the largest post-collision aggregate. We note that for large b, the contour lines are rather narrowly spaced such that in this region any changes in the threshold value would only negligibly influence the boundary of the fusion regime. For small b, however, changes in the threshold value have a higher influence; this is due to the large fragment cloud developing for higher velocities.

2. We base the analysis of sliding on the second-largest fragment, N₂. A plot of the contours of N_sl = N₂/N is shown in Fig. 4b. Keeping in mind that the largest fragment is larger than N₂, it necessarily is N₂/N ≤ 1. Again comparing this with the collision outcomes of Fig. 2, a threshold value of 0.75 for N₂/N appears to describe sliding events well and we characterize the sliding regime by $$N_{\text{sl}}=N_2/N\ge \mathrm{0.75.}$$(19)

   It is important to note that with this definition, sliding events cannot simultaneously belong to the fusion class since N₂/N ≥ 0.75 implies N₁/N ≤ 1.25 by mass conservation.

   Asymmetric collision outcomes – such as N₁ = 1.25N and N₂ = 0.75N – are formally included in this definition, but Fig. 2 shows that this hardly occurs. Rather, collisions with N₂ < 0.75N are characterized by a large mass loss (erosion) of both colliding aggregates and they are hence better classified as fragmentation rather than sliding.

3. In this study, we shall treat all events that are neither sliding nor fusing as fragmentation. This allows us to classify all collision outcomes uniquely into one of the three classes.

We note that an assessment of the fragmentation class could be based on the quantity N_fr = (N₁ + N₂)/N. Small values of N_fr indicate a small mass in the two largest clusters and hence a high degree of fragmentation. Figure 4c shows that a threshold value of N_fr = (N₁ + N₂)/N ≤ 1.5 characterizes a large part of the fragmentation regime; however, the “central part” of the classification diagram – around υ = 7.5 m s⁻¹ and b = 0.75 – is not included. One might argue that this central part could be classified by a regime of its own (“several large fragments,” see Sect. 3.6), but we failed to find a quantifier for this regime.

Earlier attempts to classify collision outcomes of granular or atomic aggregates were also based on the sizes of the largest clusters, N₁ and N₂ (Svanberg et al. 1998), or linear combinations of these were used. Also in the study of collisions between granular aggregates, the sizes of the largest surviving clusters (or fragments) have been used early on to analyze collisions (Wada et al. 2009, 2013; Paszun & Dominik 2009; Hasegawa et al. 2021). The study of collisions between partners of unequal masses made the analysis of the two largest clusters particularly relevant (Wada et al. 2013; Hasegawa et al. 2021). In a prominent work on atomic-cluster collisions, Kalweit & Drikakis (2006) proposed the quantifier X = (N₁ − N₂)/(2N) and used a threshold value of X = 0.15 to identify fusion (X > 0.15); this appears to be less restrictive than our criterion. Several studies of granular collisions (Ringl et al. 2012; Umstätter & Urbassek 2021) used the fragmentation parameter Ns = 1 − (N₁ + N₂)/(2N) which is closely connected to N_fr. We note that the study of the outcomes of mass-asymmetric collisions is more complex since its discussion is usually based on the initial sizes of both colliding aggregates (Wada et al. 2013; Hasegawa et al. 2021).

Our results agree qualitatively well with the previous study of Ringl et al. (2012), which used slightly smaller aggregates (N = 7200) but investigated impact parameters only up to b ≤ 0.9. In that study, however, the sliding class was virtually nonexistent as it becomes only noticeable for larger impact parameters.

Contour plots of the three quantifiers for the other filling factors investigated are found in Appendix A.

For the sake of completeness, we mention that discussions of the collision-induced changes in aggregates also use the notion of a “growth efficiency” (Seizinger et al. 2013; Wada et al. 2013; Planes et al. 2021). In the case of a massive target impacted by a light projectile, it indicates whether the target grew by assembling mass from the projectile. For the case of two equal-mass aggregates, its meaning is unclear since there is no distinction between “target” and “projectile”: if one of them grows, the other one shrinks, and there is no net growth induced by the collision. Quantitatively, the growth efficiency is defined in our case as $$f=\frac{N_1}{N}-1,$$(20)

and f > 0 indicates growth. If these data are plotted (see Appendix A), the growth regime f > 0 looks similar to the fusion regime described above. The growth velocity, υ_gr, is defined as the velocity where f = 0; it is the maximum velocity at which growth occurs. As the data in Appendix A show, for central collisions, it is f > 0 for all velocities studied and υ_gr > 15 m s⁻¹; the only exception is for the most compact aggregates, Φ = 0.32, where υ_gr = 13.9 m s⁻¹.

We can compare these data with previously published simulations. Ringl et al. (2012) studied collisions of aggregates (filling factor of 0.20) containing 7200 grains, which are identical to those considered here, and they obtained υ_gr = 17 m s⁻¹ for b = 0 and 5.1 m s⁻¹ for b = 1; this is in fair agreement with our present data of 20 m s⁻¹ (extrapolated) and 5 m s⁻¹ (see Fig. B.1a). In order to compare this with other grain material, we scaled our velocity with the so-called fragmentation velocity, υ_fr, which is defined as the smallest velocity at which a bound pair of grains dissociates. For the silica grains considered here, it is υ_fr = 0.17 m s⁻¹ (cf. Sect. 2.2 and Ringl et al. 2012). Wada et al. (2009) studied collisions of fractal aggregates composed of 8000 ice grains. For central collisions, they found a growth velocity of υ_gr = 31 υ_fr, which is somewhat below our value of around (100−120)υ_fr (cf. Fig. B.1a and Ringl et al. 2012). However, Wada et al. (2009) noted that υ_gr increases with aggregate size. Also, the structure of these fractal aggregates – constructed by ballistic particle-cluster aggregation – differs from the more compact clusters studied here.

Fig. 4 Classification of the collisions occurring with impact parameter b and velocity υ into the collision classes’ fusion, sliding, and fragmentation. Data are for aggregates of Φ = 0.28. Color shading was added according to the class quantifiers (a) Nfu, (b) Nsl, and (c) Nfr. The contour lines delineating the respective collision classes are marked in each plot.

3.3 Influence of filling factor Φ

Figure 5 summarizes the classification of all collision outcomes for the four filling factors studied here. The contour lines of N_fu = 1.5 and N_sl = 0.75 were used for the classification of the outcomes, as introduced in Sect. 3.2; the contour lines of N_fr = 1.5 were additionally plotted.

The boundaries of the collision classes given here are subject to some uncertainty. (i) The threshold values defined in Sect. 3.2 have some degree of arbitrariness. (ii) To enhance statistical reliability, more collisions could be simulated with identical b and υ, in which the aggregate structure was varied, for instance by rotating the aggregates with respect to each other.

We see in Fig. 5 that with increasing Φ, the fragmentation regime shrinks and its boundary moves toward higher υ. The sliding regime increases for larger Φ, extending both toward higher υ and lower b. The fusion regime shows the least changes with increasing Φ. In detail, however, for larger Φ, the fusion regime extends toward higher b values in the lower υ domain, but the fusion boundary tends to decrease toward lower υ for near-central collisions.

The strongest influence of Φ is observed in the central part of the classification diagram –at υ ~ 7.5 m s⁻¹ and b ~ 0.75. In this region, the three regimes meet. With an increasing filling factor, the sliding regime grows at the expense of fragmentation.

Previous studies of the influence of the aggregate filling factor on aggregate collisions were mostly restricted to central collisions. Gunkelmann et al. (2016) show that for filling factors below 0.11, the fragmentation regime extends to increasingly smaller collision velocities. Wada et al. (2009) used aggregates built with various ballistic growth processes; including also non-central collisions in their study, they demonstrate this same trend – that with increasing porosity, aggregates fragment at lower velocity. Their findings are thus in line with our results.

Fig. 5 Classification of the collision outcomes in the b–υ parameter space for all the Φ values (0.20, 0.24, 0.28, and 0.32) investigated. The defining contour lines of Nfu, Nsl, and Nfr, as introduced in Sect. 3.2 have been added.

3.4 Average over impact parameter

For most applications, the impact parameter is not known, and thus averages of collision properties over the impact parameter are of interest. We present here the probabilities of fusion, sliding, and fragmentation at fixed velocity υ. The fusion probability is given by $$\langle p_{\text{fu}}\rangle =\frac{{\displaystyle {\int }_0^{b_{\max }}2\pi b{p}_{\text{fu}}\left(b,\upsilon \right)\text{d}b}}{\pi b_{\max }{}^2},$$(21)

where p_fu(b, υ) assumes the value 1 or 0, depending on whether the collision (b, υ) was considered to belong to the fusion class or not.

The average sliding probability is defined analogously, with p_sl(b, υ) = 1 if the collision (b, υ) is in the sliding class, and 0 otherwise. The average fragmentation probability is given by p_fr(b, υ) = 1 − p_fu(b, υ) − p_si(b, υ). The b integration in Eq. (21) was performed using the trapezoidal rule and assuming that for the maximum impact parameter, b_max = 2, it is p_sl(b, υ) = 1.

Figure 6 displays the resulting probabilities for aggregates with the various filling factors simulated. We see that at the lowest velocities fusion dominates, in agreement with the results of Fig. 5 that show fusion prevailing for all but the largest impact parameters. With increasing velocity, sliding becomes the dominant collision class with a maximum at around 7.5–10 m s⁻¹. At even higher velocities, fragmentation gains importance. In summary, the averaged probabilities show the pattern expected from the detailed impact-parameter-resolved collisions discussed above.

Similarly, the growth efficiency, Eq. (20), can be averaged over the impact parameters (see Fig. 7). The impact-parameter-averaged growth efficiency, 〈f〉, shows a steep decline in the velocity range of 2.5–5 m s⁻¹ where the fusion probability strongly drops. For higher velocities, 〈f〉 becomes negative and continues declining, albeit at a smaller rate, as fragmentation gains in importance, 〈f〉. This behavior qualitatively coincides with that found in previous studies (Hasegawa et al. 2021). The impact-parameter-averaged growth velocity, 〈υ_gr〉, amounts to 6.3 (7.5, 8.8, 8.5) m s⁻¹ for aggregates with filling factor Φ = 0.20 (0.24, 0.28, 0.32). The trend of increasing growth velocity for an increasing filling factor corroborates the fact that aggregates of higher porosity show an increased trend to fragmentation.

These data can be compared to the study by Wada et al. (2009) using fractal aggregates of 8000 ice grains, which was recently extended to 260 000 grains (Hasegawa et al. 2021). These studies find that the critical collision energy for growth amounts to (10–15) N_totE_break, where N_tot = 2N and $E_{\text{break}}=\left(m/4\right){\upsilon }_{\text{fr}}^2$. Relating the critical impact energy to the average growth velocity, 〈υ_gr〉, via E = Nm〈υ_gr〉²/4, the abovementioned studies thus find 〈υ_gr〉 = (4.5 − 5.5)υ_fr. This value is considerably smaller than our result of 6.3 m s⁻¹ (corresponding to 37υ_fr) for the smallest filling factor simulated. We assume that this difference is due to the fluffy nature of the aggregates used in the studies of ice grains as opposed to the compact aggregates used in the present study.

Fig. 6 Velocity dependence of the averaged fusion, sliding, and fragmentation probabilities for the filling factors simulated in this study.

Fig. 7 Impact-parameter-averaged growth efficiency, 〈f〉, for the filling factors simulated in this study. The line 〈f〉 = 0 separates the collisions where aggregates grow on average from those where they erode.

3.5 Emission angle of fragment cloud

At higher velocity, the collision outcome is accompanied by a cloud of small ejecta, which are mostly monomers (Ringl et al. 2012). This cloud is clearly visible in Fig. 2 at velocities υ ≳ 5 m s⁻¹, even if the outcome is not in the fragment class. The fragment cloud shows, at its center, a constant slope with respect to the collision velocity; we denote this emission angle by α. For central collisions (b = 0), it is α = 90°, and for sliding collisions (b = 1.75), it is α = 0°.

One may therefore suspect that α is connected to the angle ψ that the collision plane makes to the collision velocity (Fig. 1). Such a relationship would originate if grains were ejected during the collision from the collision plane – since here shear forces are the strongest – and in a direction parallel to the plane. Figure 8 displays the correlation between these two angles α and ψ. The angle α was determined from the inertia tensor of the monomer distribution in the central region where the fragment cloud shows a constant slope. We see that this expectation is not fulfilled well; apart from central collisions (α = ψ = 90°), it is always α < ψ, that is to say the fragment cloud is more inclined toward the “forward” direction – the direction of the velocity vector – than the collision plane. This fact demonstrates that fragment ejection is not entirely in directions parallel to the collision plane, but the aggregate momentum also influences the ejection direction. We note that for more compact aggregates with a large value of Φ, α is better aligned in the direction of ψ; this appears plausible since aggregate cohesion reduces the tendency to eject fragments in the forward direction.

At higher impact parameters, b ≳ 0.75, the fragment cloud shows a deviation from its linear shape at its outermost ends; this gives it an S-like morphology. These outermost parts of the fragment cloud are aligned in the direction of the collision velocity. They correspond to ejecta that were emitted early in the collision and proceed in a forward direction; in other words, these are fragments that were stripped off in these high-impact parameter collisions – their inertia lets them fly on trajectories along their initial momenta. A comparison of the morphologies for various filling factors (see Appendix C) shows that the relevance of this forward component of the ejection cloud decreases for high filling factors.

An evaluation of the velocity dependence of α (see Appendix C) shows that α does not depend on υ. This is in line with the idea that the ejection angle is connected to the (velocity-independent) angle of the collision plane.

The emission angle does not change with time as we can see when continuing our simulation for another 200 μs (see Appendix C). This appears plausible since fragment collisions, which could change the cloud shape, become less frequent at late times since the density in the fragment clouds decreases with time.

Fig. 8 Variation of the emission angle α of the fragment cloud with the angle of the collision plane, ψ, for various Φ. Impact parameters b corresponding to ψ, Eq. (17), are shown. The dashed black line shows α = ψ, which characterizes fragment ejection in directions parallel to the collision plane. Data are for υ = 15 m s−1.

Fig. 9 Series of snapshots showing the time evolution of a collision generating three large fragments: b = 1 and υ = 5 m s−1. Grains are colored according to their initial aggregate affiliation.

3.6 Several large fragments

Close to the line separating fusion from fragmentation – and in particular in the “center” of the classification diagram, where the three regimes meet – typically several large fragments appear. The occurrence of three fragments is particularly intriguing such as for υ = 5 m s⁻¹ and with b = 1 in Fig. 2. These three fragments correspond to a fused product aggregate (containing, however, only a small fraction of the total number of grains), plus two sheared off remains of (partly) sliding aggregates. Figure 9 displays a series of snapshots demonstrating the formation of these three fragments. We note that the occurrence of such three-fragment outcomes is typical, not only for the present code (see Fig. 1 of Ringl et al. 2012), but also for other granular-mechanics codes, such as Fig. 17 of Hasegawa et al. (2021). It would appear tempting to introduce a collision class of its own for such few-fragment cases, but we did not succeed in finding an appropriate quantifier.

4 Summary

In this work, we have studied collisions of equal-mass silica-grain aggregates as a function of the impact parameter (b), collision velocity (υ), and filling factor (Φ). The main conclusions of this work can be stated as follows:

• The collision outcomes can be classified as fusion, sliding, and fragmentation events, in agreement to earlier studies of atomic-cluster (Kalweit & Drikakis 2006) and granular-aggregate collisions (Ringl et al. 2012). Fusion mainly occurs for low velocities, sliding for large impact parameters, and fragmentation for high velocities, but the details depend on the exact values of impact parameter b, velocity υ, and filling factor Φ;

• We base a quantitative description of these classes on the numbers of grains contained in the largest (N₁) and second-largest (N₂) aggregates present after the collision. We define descriptors of fusion: N_fu = N₁/N ≥ 1.5; sliding, N_sl = N₂/N ≥ 0.75; and fragmentation otherwise. These descriptors allow us to characterize the collision events obtained for aggregates with all filling factors studied;

• The boundaries of the three collision classes change with Φ. With an increasing filling factor, the sliding probability – and to a lesser degree also the fusion probability at small velocities – increases at the expense of the fragmentation probability;

• Sliding only occurs for noncentral collisions as these occur most frequently in experiment, simulation data for central collisions and their analyses should be treated with caution;

• When collisions are averaged over the impact parameter, the probability for fusion is dominant at low velocities, but quickly vanishes for higher velocity, above around 7.5 m s⁻¹. The probability for sliding then dominates but decreases with increasing velocity, such that at the highest velocities investigated here, fragmentation becomes most prominent.

Our results extend previous simulation results on the influence of aggregate porosity (Gunkelmann et al. 2016) to noncentral collisions. As we demonstrate in the present study, only noncentral collisions allow for aggregate sliding and are therefore essential for an understanding of the collision behavior of aggregates.

Acknowledgements

Simulations were performed at the High Performance Cluster Elwetritsch (RHRK, TU Kaiserslautern, Germany).

Appendix A Simulation results for other filling factors

A.1 Simulation results for filling factor 0.20

Fig. A.1 Simulation results for the quantifiers (a) fusion, (b) sliding, and (c) fragmentation for collisions of two aggregates of filling factor 0.20 colliding with velocity υ and impact parameter b.

Fig. A.2 Tableau of the final state of two aggregates of filling factor 0.20 colliding with velocity υ and impact parameter b.

A.2 Simulation results for filling factor 0.24

Fig. A.3 Simulation results for the quantifiers (a) fusion, (b) sliding, and (c) fragmentation for collisions of two aggregates of filling factor 0.24 colliding with velocity υ and impact parameter b.

Fig. A.4 Tableau of the final state of two aggregates of filling factor 0.24 colliding with velocity υ and impact parameter b.

A.3 Simulation results for filling factor 0.32

Fig. A.5 Simulation results for the quantifiers (a) fusion, (b) sliding, and (c) fragmentation for collisions of two aggregates of filling factor 0.32 colliding with velocity υ and impact parameter b.

Fig. A.6 Tableau of the final state of two aggregates of filling factor 0.32 colliding with velocity υ and impact parameter b.

Appendix B Growth efficiency

Fig. B.1 Dependence of the growth efficiency, f = (N1 − N)/N, on impact parameter b and collision velocity υ for filling factors of (a) 0.20, (b) 0.24, (c) 0.28, and (d) 0.32.

Appendix C Fragment cloud

Fig. C.1 Variation of the emission angle α of the fragment cloud with velocity υ for all impact parameters studied. Data are for Φ = 0.28.

Fig. C.2 Change of the emission angle α of the fragment cloud with time. Data are for Φ = 0.28. (a) shows the changes between times t = 100 μs and t = 200 μs for several angles of the collision plane ψ for a collision velocity of 15 m/s; (b) shows the changes up to a time of t = 400 μs for the case b = 1 and υ = 7.5 m/s.

References

All Tables

All Figures

Fig. 1 Schematic diagram showing the moment when the two granular aggregates just touch. The aggregates are equal in size with radius R and collide with a relative velocity υ. The point of contact is denoted by “C.” The solid green line connects the centers of the aggregates. The dashed gray lines are drawn from the centers of the aggregates along the direction of motion. The impact parameter of the collision, B, is the perpendicular distance between the gray lines, shown by black arrows. The collision angle is denoted by θ. The dashed black line represents the edge-on view of the common tangent plane of the two aggregates on the point of contact; it is denoted as the “collision plane.” The angle between the collision plane and the direction of motion is denoted as ψ.

In the text

Fig. 2 Tableau of the final state after the collision of two aggregates of filling factor 0.28 colliding with velocity υ (in units of m s−1) and reduced impact parameter b. Grains were projected into the plotted plane in order to allow for a better perception of the monomer distribution. Grains are colored according to their initial aggregate affiliation. The size of the subpanels changes with the collision velocity.

In the text

	Fig. 3 Series of snapshots showing the time evolution of a sliding collision: b = 1.25 and υ = 5 m s−1. Grains are colored according to their initial aggregate affiliation.
In the text

	Fig. 4 Classification of the collisions occurring with impact parameter b and velocity υ into the collision classes’ fusion, sliding, and fragmentation. Data are for aggregates of Φ = 0.28. Color shading was added according to the class quantifiers (a) Nfu, (b) Nsl, and (c) Nfr. The contour lines delineating the respective collision classes are marked in each plot.
In the text

	Fig. 5 Classification of the collision outcomes in the b–υ parameter space for all the Φ values (0.20, 0.24, 0.28, and 0.32) investigated. The defining contour lines of Nfu, Nsl, and Nfr, as introduced in Sect. 3.2 have been added.
In the text

	Fig. 6 Velocity dependence of the averaged fusion, sliding, and fragmentation probabilities for the filling factors simulated in this study.
In the text

	Fig. 7 Impact-parameter-averaged growth efficiency, 〈f〉, for the filling factors simulated in this study. The line 〈f〉 = 0 separates the collisions where aggregates grow on average from those where they erode.
In the text

	Fig. 8 Variation of the emission angle α of the fragment cloud with the angle of the collision plane, ψ, for various Φ. Impact parameters b corresponding to ψ, Eq. (17), are shown. The dashed black line shows α = ψ, which characterizes fragment ejection in directions parallel to the collision plane. Data are for υ = 15 m s−1.
In the text

	Fig. 9 Series of snapshots showing the time evolution of a collision generating three large fragments: b = 1 and υ = 5 m s−1. Grains are colored according to their initial aggregate affiliation.
In the text

	Fig. A.1 Simulation results for the quantifiers (a) fusion, (b) sliding, and (c) fragmentation for collisions of two aggregates of filling factor 0.20 colliding with velocity υ and impact parameter b.
In the text

	Fig. A.2 Tableau of the final state of two aggregates of filling factor 0.20 colliding with velocity υ and impact parameter b.
In the text

	Fig. A.3 Simulation results for the quantifiers (a) fusion, (b) sliding, and (c) fragmentation for collisions of two aggregates of filling factor 0.24 colliding with velocity υ and impact parameter b.
In the text

	Fig. A.4 Tableau of the final state of two aggregates of filling factor 0.24 colliding with velocity υ and impact parameter b.
In the text

	Fig. A.5 Simulation results for the quantifiers (a) fusion, (b) sliding, and (c) fragmentation for collisions of two aggregates of filling factor 0.32 colliding with velocity υ and impact parameter b.
In the text

	Fig. A.6 Tableau of the final state of two aggregates of filling factor 0.32 colliding with velocity υ and impact parameter b.
In the text

	Fig. B.1 Dependence of the growth efficiency, f = (N1 − N)/N, on impact parameter b and collision velocity υ for filling factors of (a) 0.20, (b) 0.24, (c) 0.28, and (d) 0.32.
In the text

	Fig. C.1 Variation of the emission angle α of the fragment cloud with velocity υ for all impact parameters studied. Data are for Φ = 0.28.
In the text

	Fig. C.2 Change of the emission angle α of the fragment cloud with time. Data are for Φ = 0.28. (a) shows the changes between times t = 100 μs and t = 200 μs for several angles of the collision plane ψ for a collision velocity of 15 m/s; (b) shows the changes up to a time of t = 400 μs for the case b = 1 and υ = 7.5 m/s.
In the text