Collisional Disruption of Planetesimals in the Gravity Regime with iSALE Code: Comparison with SPH code for Purely Hydrodynamic Bodies
Ryo Suetsugu, Hidekazu Tanaka, Hiroshi Kobayashi, Hidenori Genda
aa r X i v : . [ a s t r o - ph . E P ] M a y Collisional Disruption of Planetesimals in the Gravity Regimewith iSALE Code: Comparison with SPH code for PurelyHydrodynamic Bodies
Ryo Suetsugu , Hidekazu Tanaka , Hiroshi Kobayashi , Hidenori Genda
1. School of Medicine, University of Occupational and Environmental Health, Iseigaoka,Yahatanishi-ku, KitaKyushu, 807-8555, Japan2. Astronomical Institute, Tohoku University, Aramaki, Aoba-ku, Sendai, 980-8578, Japan3. Department of Physics, Graduate School of Science, Nagoya University, Furo-cho,Chikusa-ku, Nagoya, 464-8602, Japan4. Earth-Life Science Institute, Tokyo Institute of Technology, Ookayama, Meguro-ku,Tokyo, 152-8550, Japan [email protected]
Received ; acceptedNot to appear in Nonlearned J., 45. 2 –
ABSTRACT
In most of the previous studies related to collisional disruption of planetesi-mals in the gravity regime, Smoothed Particle Hydrodynamics (SPH) simulationshave been used. On the other hand, impact simulations using grid-based hydro-dynamic code have not been sufficiently performed. In the present study, weexecute impact simulations in the gravity regime using the shock-physics codeiSALE, which is a grid-based Eulerian hydrocode. We examine the dependenceof the critical specific impact energy Q ∗ RD on impact conditions for a wide rangeof specific impact energy ( Q R ) from disruptive collisions to erosive collisions, andcompare our results with previous studies. We find collision outcomes of theiSALE simulation agree well with those of the SPH simulation. Detailed analysismainly gives three results. (1) The value of Q ∗ RD depends on numerical resolution,and is close to convergence with increasing numerical resolution. The differencein converged value of Q ∗ RD between the iSALE code and the SPH code is within30%. (2) Ejected mass normalized by total mass ( M ej /M tot ) generally dependson various impact conditions. However, when Q R is normalized by Q ∗ RD that iscalculated for each impact simulation, M ej /M tot can be scaled by Q R /Q ∗ RD , andis independent of numerical resolution, impact velocity and target size. (3) Thissimilarity law for Q R /Q ∗ RD is confirmed for a wide range of specific impact energy.We also derive a semi-analytic formula for Q ∗ RD based on the similarity law andthe crater scaling law. We find that the semi-analytic formula for the case witha non-porous object is consistent with numerical results. Subject headings:
Impact processes, Cratering, Planetary formation
1. INTRODUCTION
Collisions are one of the most important processes in planet formation becauseplanetary bodies in the Solar System are thought to have experienced a lot of collisionsduring the accretion process (e.g., Lissauer 1993). Thus, collisional processes have beenexamined extensively. Roughly speaking, collisional outcomes can be classified intodisruptive collisions and erosive collisions by the specific impact energy Q R , given by Q R = (cid:18) M tar V + 12 M imp V (cid:19) /M tot = (cid:18) M R v (cid:19) /M tot , (1)where M tar and M imp are the mass of the target and the impactor ( M tar > M imp , M tot = M tar + M imp ), respectively, and V tar and V imp are the velocities of the target andthe impactor in the frame of the center of mass when the two objects contact each other,respectively, M R is the reduced mass, given by M imp M tar /M tot , and v imp is the impactvelocity ( v imp = V imp − V tar for negative V tar ). In particular, the specific impact energyrequired to disperse the largest body such that it has exactly half its total mass after thecollision is called the critical specific impact energy Q ∗ RD . In the case of Q R > Q ∗ RD , collisionsbetween planetesimals are regarded as disruptive collisions, while they are non-disruptivecollisions for Q R ≪ Q ∗ RD , whose mass ejection is small (hereafter called erosive collisions).The values of Q ∗ RD have been investigated by laboratory experiments and numericalsimulations (e.g., Benz & Asphaug 1999; Nakamura et al. 2009). When the target is smallenough to neglect the effect of the target’s gravity, the critical specific impact energy ismainly estimated by laboratory experiments (Housen & Holsapple 1999; Nakamura et al.2009). As target size increases, collision outcomes gradually become dominated bythe gravity of the target. However, direct experimental measurements of a large scalecollision are difficult to carry out in the laboratory. Thus, the values of Q ∗ RD for largetargets ( & & km/s): Lagrangianhydrocode such as Smoothed Particle Hydrodynamics (SPH) methods (Love & Ahrens1996; Melosh & Ryan 1997; Benz & Asphaug 1999; Jutzi et al. 2010; Genda et al. 2015;Jutzi 2015; Movshovitz et al. 2016; Genda et al. 2017), or a hybrid code of Eulerianhydrocode and N-body (Leinhardt & Stewart 2009). These numerical simulations showedthe dependence of the value of Q ∗ RD on various impact conditions such as target size, impactvelocity, material properties, and impact angle. For example, the value of Q ∗ RD in thegravity regime increases nearly monotonically with the size of the target because collisionalfragments are more easily bound by the gravitational force of the target. The critical specificimpact energy also depends on the material property (e.g. material strength, porosity, andfriction) of the impactor and the target (Leinhardt & Stewart 2009; Jutzi et al. 2010; Jutzi2015). Notably, the friction significantly dissipates impact energy (Kurosawa & Genda2018), which tends to hinder the disruption of the target. The value of Q ∗ RD then reachesabout 10 times the value of Q ∗ RD without the friction (Jutzi 2015). Moreover, recent impact 4 –simulations show that Q ∗ RD depends not only on impact conditions, but also on numericalresolution (Genda et al. 2015, 2017). Genda et al. (2015) performed SPH simulation atvarious numerical resolutions, and showed that Q ∗ RD at high numerical resolution is ratherlow compared to the case of low resolution.In addition to the critical specific impact energy, the understanding of erosivecollisions is also important in relation to the formation of planetary bodies. In most of theprevious studies, the contribution of erosive collision to growth of the planets has beenunderestimated because the amount of mass ejected by erosive collision is much smallerthan the total mass. However, some previous studies showed that erosive collision also playsan important role in planetary accretion (Kobayashi & Tanaka 2010; Kobayashi et al.2010, 2011). Kobayashi & Tanaka (2010) assumed a simple fragmentation model describingboth disruptive collisions and erosive collisions, and investigated mass depletion timein a collision cascade based on analytic consideration and numerical simulation. Theyshowed that erosive collisions occur much more frequently than disruptive collisions andthe mass depletion time is mainly determined by erosive collisions. Recently, the validityof the simple fragmentation model was examined by Genda et al. (2017), who performedimpact simulations for a wide range of specific impact energy using the SPH method withself-gravity and without material strength (i.e. a purely hydrodynamic case), and showedthat the fragmentation model is consistent with collisional outcomes of simulations withina factor of two. They also showed that the ejected mass normalized by the total mass canbe scaled by Q R /Q ∗ RD for their parameter range.However, almost all high velocity collisions have been examined by the SPH method.Another common hydrodynamic simulation, whose computational domain is discretized bygrids, has also been carried out (e.g., Leinhardt & Stewart 2009). However, the grid-basedcode is only used for the shock deformation immediately after collision, and a large partof the disruption is calculated by N-body simulation. Thus, impact simulation using thegrid-based code has not been sufficiently performed, though it is important to examine theproblem with a different numerical approach.In this study, we perform impact simulations in the gravity regime by using shock-physicscode iSALE (Amsden et al. 1980; Collins et al. 2004; W¨unnemann et al. 2006; Collins et al.2016), which is a grid-based Eulerian hydrocode, and has been widely distributed toacademic users in the impact community. This code has been used to understand variousimpact phenomena: crater formation (Collins et al. 2008; Cremonese et al. 2012), impactjetting (Johnson et al. 2015; Wakita et al. 2017; Kurosawa et al. 2018), pairwise collisionsof planetesimals with/without self-gravity (Davison et al. 2010, 2012) and comparison withexperimental data (Nagaki et al. 2016; Kadono et al. 2018). We examine the dependenceof Q ∗ RD on numerical resolution and impact conditions for a wide range of specific impactenergy from disruptive collisions to erosive collisions, and compare our results with previousstudies. Furthermore, using numerical results obtained by the iSALE code and the craterscaling law, we derive a semi-analytic formula for Q ∗ RD . 5 –In Section 2, we present methods for impact simulations and analysis. We show ournumerical outcomes of simulations in the case of disruptive collisions in Section 3. InSection 4, we establish a similarity law for Q R /Q ∗ RD for a wide range of impact energy, andderive a semi-analytic formula for Q ∗ RD . We discuss effects of oblique collisions and materialproperties in our results in Section 5. Section 6 summarizes our results.
2. NUMERICAL METHODS
In this study, we examine collisions between planetesimals using shock-physics codeiSALE-2D, the version of which is iSALE-Chicxulub. The iSALE-2D is an extension ofthe SALE hydrocode (Amsden et al. 1980). To simulate hypervelocity impact processesin solid materials, SALE was modified to include an elasto-plastic constitutive model,fragmentation models, and multiple materials (Melosh et al. 1992; Ivanov et al. 1997).More recent improvements include a modified strength model (Collins et al. 2004), and aporosity compaction model (W¨unnemann et al. 2006; Collins et al. 2011).The iSALE-2D supports two types of equation of state: ANEOS (Thompson & Lauson1972; Melosh 2007) and Tillotson equation of state (Tillotson 1962). These equations ofstate have been widely applied in previous studies including planet- and planetesimal-sizecollisional simulations (e.g. Canup & Asphaug 2001; Canup 2004; Fukuzaki et al. 2010;´Cuk & Stewart 2012; Sekine & Genda 2012; Hosono et al. 2016; Wakita et al. 2017). In oursimulation, we use the Tillotson equation of state for basalt because almost all previousstudies related to collisional disruption have used the Tillotson equation of state, whichallows us to directly compare our results with theirs. The Tillotson equation of statecontains ten material parameters, and the pressure is expressed as a function of the densityand the specific internal energy; all of which are convenient when used in works regardingfluid dynamics. Although the Tillotson parameters for basalt of the iSALE-2D are set toexperimental values, we used the parameter sets of basalt referenced in previous works(Benz & Asphaug 1999; Genda et al. 2015, 2017).We employ the two-dimensional cylindrical coordinate system and perform head-onimpact simulations between two planetesimals (Figure 1). We assumed that planetesimalsare not differentiated. Planetesimals are also assumed to be composed of basalt. Fornominal cases, the radius of the target R tar and the impact velocity of the impactor v imp arefixed at 100 km and 3 km/s, respectively. We also examine the dependence of collisionaloutcome on target size and impact velocity in Section 3.2. To carry out impact simulationswith various impact energy Q R , we changed the radius of the impactor R imp . For example, R imp = 14 - 21 km (i.e., Q R ≃
12 - 41 kJ/kg). In this study, we consider four cases with thenumber of cells per target radius ( n tar = 100 , , n v × n h , see Figure 1) is changed dependingon n tar . For example, ( n v × n h ) = (450 × , (900 × , (1800 × × n tar = 100, 200, 400, and 800, respectively. In the case of R tar = 100 km, the values ofthe spatial cell size for each numerical resolution ∆ x (= R tar /n tar ) are ∆ x = 1000, 500, 250,and 125 m, and the size of the computational domain is fixed at ( n v ∆ x, n h ∆ x ) =(450 km,450 km).The aim of this study is to make a direct comparison of collisional outcomes betweendifferent numerical codes (SPH code and iSALE code). Therefore, although the iSALE-2Dcan deal with the effects of material strength, damage, and porosity of the target and theimpactor, these effects are not taken into account in the present work; that is, the fluidmotion is purely hydrodynamic. The self-gravity is calculated by the algorithm in theiSALE-2D based on a Barnes-Hut type algorithm, which can reduce the computational costof updating the gravity field. In most of our calculations, the opening angle θ , which isthe ratio of mass length-scale to separation distance, is set to 1.0. Although the value ofthe opening angle is rather large, we confirmed that the difference in ejected mass betweencases for θ = 1 . θ = 0 . ∼ R tar = 300 km because the gravity of such a largeplanetesimal is relatively strong (Section 3.2). Then, the calculation time of a single impactsimulation mainly depends on the numerical resolution. For example, in our calculationwith n tar = 800, the calculation time is a few weeks for R tar = 100 km, while it is twomonths for R tar = 30 km, using a computer with an Intel R Core TM i7-4770K Processor (3.50GHz). Other input parameters for iSALE simulation are summarized in Table 1.In order to obtain the value of Q ∗ RD , we need to estimate the mass of ejecta caused bya single disruptive collision, M ej , defined as M ej = M tot − M lrg , (2)where M lrg is the mass of the largest body resulting from the impact. The mass of thelargest body M lrg is determined by the following procedure. First, we find the groups of cellsin contact with cells of non-zero densities and compare the total masses of the cells in therespective groups. The largest total mass is regarded as the preliminary mass of the largestbody M lrg . This procedure roughly corresponds to a friends-of-friends algorithm to identifycollisional fragments used in previous studies with a SPH method (e.g., Benz & Asphaug1999; Genda et al. 2015, 2017). Next, if the constituent cells in the tentative largest bodyare gravitationally bound, the largest body is determined. Otherwise, we find the cellswhere the specific kinetic energies are larger than the specific gravitational potential energyof the tentative largest body and remove them from the largest body. This procedure isiteratively performed until the number of the cells in the largest mass converges. We regardthe converged M lrg as the mass of the largest body. 7 –
3. DISRUPTIVE COLLISIONS AND THE CRITICAL SPECIFIC IMPACTENERGY Q ∗ RD Q ∗ RD at a nominal case We examine the critical specific impact energy Q ∗ RD for a nominal case with the targetradius R tar = 100 km and the impact velocity v imp = 3 km/s. Figure 2 shows a timeseries of a simulation of an impactor’s head-on collision ( R imp = 16 km) with the target( R tar = 100 km) at an impact velocity of v imp = 3 km/s. We adopt the number of cells pertarget radius n tar = 800. The color contour represents the specific kinetic energy. When t = 0 s, the impactor starts colliding with the target. The shock wave generated by theimpact propagates through the target and arrives at the antipode of the impact point at t ≃
50 s. Most of the ejecta formed by the collision have relatively large amounts of kineticenergy, and continue escaping from the gravity of the target even when t &
300 s.In Figure 3(a), the black curve represents the time evolution of the ejected mass ( M ej )normalized by the total mass of colliding bodies ( M tot ) in the impact simulation shown inFigure 2. Immediately after the collision, the value of M ej /M tot increases rapidly. However,it converges within a short time ( t .
170 s) and eventually becomes nearly constant( M ej /M tot ≃ .
49) at t = 350 s. Figure 3(a) also shows the resolution dependence of the timeevolution of M ej /M tot . The general feature of the ejected mass for each numerical resolutionis similar to that with n tar = 800. However, the converged values of M ej /M tot increase asnumerical resolution increases, which was also observed in the results obtained by using theSPH code (Genda et al. 2015). Figure 3(b) shows the case with more destructive collision( R imp = 19 km). Although the differences in the values of M ej /M tot between numericalresolutions become smaller, the basic feature is similar to the cases with R imp = 16 km.Figure 4(a) shows M ej /M tot for various numerical resolutions as a function of Q R .The values of M ej /M tot are listed in Table 2. In all the cases, we find that the ejectedmass increases as impact energy increases. We also find that M ej /M tot for each impactenergy increases with n tar , and the differences in M ej /M tot between numerical resolutionsbecome smaller in the case of higher resolution. The vertical dashed lines in Figure 4(a)represent the critical specific impact energy Q ∗ RD for each numerical resolution, which can becalculated by the linear interpolation of the two data sets of Q R across M ej /M tot = 0 . Q ∗ RD = 24 . , . , .
73 and 18.99 kJ/kg for n tar = 100,200, 400, and 800, respectively. Figure 4(b) shows Q ∗ RD as a function of the spatial cell size∆ x divided by the diameter of the target D tar . Here, ∆ x/D tar corresponds to (2 n tar ) − .The values of Q ∗ RD decrease monotonically with decreasing ∆ x/D tar and would be closeto convergence. Even for the highest resolution simulation ( n tar = 800), the value of Q ∗ RD does not fully converge. However, Genda et al. (2015, 2017) showed the dependence of Q ∗ RD can be approximated well by a linear function of the inverse of the numerical resolution. 8 –Therefore, we also fit our result with the following relation, Q ∗ RD = Q ∗ RD , + b (∆ x/D tar ) , (3)where Q ∗ RD , and b are fitting parameters. As a result, we find that the fitting parametersare determined to be Q ∗ RD , = 18 .
25 kJ/kg and b = 1241 kJ/kg in the case of R tar = 100km and v imp = 3 km/s. The value of fitting parameter Q ∗ RD , in Equation (3) correspondsto the converged Q ∗ RD in the limit of ∆ x/D tar → n tar → ∞ ). From these results, weconclude that ejecta masses and Q ∗ RD obtained by iSALE code also depend on numericalresolution as well as those in SPH simulations by Genda et al. (2015, 2017).We compare these results with Genda et al. (2015) in detail. In the case of R tar = 100km and v imp = 3 km/s, Genda et al. (2015) carried out SPH simulations for four differentnumerical resolutions ( N SPH = 5 × , 1 × , 5 × , and 5 × , where N SPH is thenumber of SPH particles in the target) and estimated the converged values of the head-oncritical specific impact energy. Although their hydrocode is not based on two-dimensionalEulerian code such as iSALE-2D, but three-dimensional Lagrangian code, we simply assumethat the number of cells on target diameter 2 n tar is equivalent to N / . Under theassumption, the case of lowest numerical resolution ( n tar = 100) in this study roughlycorresponds to the case of highest resolution ( N SPH = 5 × ) in Genda et al. (2015) .Figure 5 is the same as Figure 4(b), but head-on Q ∗ RD obtained from SPH simulations(Genda et al. 2015) are also plotted. The converged values in both Genda et al. (2015) andin our study are similar values within a 30% range of each other. The slope from numericaldata obtained by the iSALE-2D is somewhat steeper than in the case of SPH simulations.Thus, for the same numerical resolution, the value of Q ∗ RD estimated by SPH simulations iscloser to each of their converged values.Recently, Genda et al. (2017) found that the dependence of M ej /M tot on Q R /Q ∗ RD isindependent of numerical resolutions: once the converged value (or reasonable value of Q ∗ RD )is obtained from very high-resolution simulations, the general behavior of M ej /M tot is givenby low resolution simulations. Thus, we examine whether the ejected mass due to disruptivecollision can be scaled by Q R /Q ∗ RD because they mainly focus on erosive collisions. When the target consists of N SPH
SPH particles, the relationship between N SPH and thenumber of SPH particles for the target radius n SPH is given as N / = (4 π/ / n SPH ≃ . n SPH . However, it is not clear whether n SPH is equal to n tar because of the difference inthe numerical scheme. Thus, we assume N / ≃ . n SPH ≃ n tar . In comparison with SPH code, this correspondence does not change even if we useimpactor resolution. In Genda et al. (2015), the impactor is composed of ∼ × SPHparticles for the highest resolution case ( N SPH = 5 × ). Under the assumption that2 n imp ≃ N / , imp (where N SPH , imp is the number of SPH particles for the impactor), thenumber of cells per impactor radius n imp is estimated to be 13.6, which is consistent withthe value of n imp used by the lowest resolution case (see Table 2). 9 –Figure 6 shows the same results as shown in Figure 4(a), but Q R is normalized by eachcalculated value of Q ∗ RD . We find that our numerical data can be clearly scaled by Q R /Q ∗ RD .Asterisks in Figure 6 represent collision outcomes of SPH simulations for N SPH = 5 × (Genda et al. 2015). Our results agree well with their results despite the different numericalschemes. Therefore, this would suggest that even numerical scheme-dependence of ejectedmass can be scaled by Q R /Q ∗ RD . Stewart & Leinhardt (2009) derived an empirical universallaw given by M lrg /M tot = − . Q R /Q ∗ RD −
1) + 0 . M ej M tot = 0 . Q R Q ∗ RD . (4)Equation (4) is also drawn in Figure 6. Our results are in a good agreement withEquation (4) for Q R /Q ∗ RD = 0 . In Section 3.1, we obtained Q ∗ RD with high accuracy for the case with R tar = 100 kmand v imp = 3 km/s. Here, we examine the dependences of Q ∗ RD on the target size and theimpact velocity.Figure 7 shows the dependence of Q ∗ RD on ∆ x/D tar for different target sizes. The valuesof Q ∗ RD for R tar = 30 and 300 km are estimated by the linear interpolation (Equation (3))in the same way as the case with R tar = 100 km. In order to compare the convergence of Q ∗ RD for each target size, Q ∗ RD is normalized by the potential energy at the target’s surface U tar ( ≡ GM tar /R tar , where G is the gravitational constant) (e.g., Movshovitz et al. 2016),whose values are given by U tar = 0 .
68, 7.55, and 67.9 kJ/kg for R tar = 30, 100, and 300km, respectively. Using the least-squares fit to Q ∗ RD /U tar for each target size, we estimatethe values of Q ∗ RD , /U tar for R tar = 30 and 300 km are 4.16 and 1.90, which correspondto Q ∗ RD , = 2 .
83 and 129.01 kJ/kg, respectively. Therefore, Q ∗ RD , /U tar is proportional to R − . , which leads to Q ∗ RD , ∝ R . for constant density and impact velocity. We alsofind that the slopes of the lines become steeper in the case of smaller target size (see alsoGenda et al. 2017). This is explained by low resolutions of impactors for small targetsbecause small mass ratios between impactors and targets are required to find Q ∗ RD for smalltargets. This result also suggests that the value of Q ∗ RD for a large target is closer to theconverged value in the case where the resolution is the same.Figure 8 shows the dependence of the critical specific impact energy for head-oncollision on target sizes. In Figure 8, the diamond symbols represent Q ∗ RD , obtained fromthe present work and the other symbols are the head-on critical specific impact energyfrom previous works. The critical specific impact energy in the gravity regime is knownto increase with R tar . In fact, Q ∗ RD , estimated by the present work also increases as 10 –target radius increases. We find that the values of Q ∗ RD , become rather low comparedto the critical specific impact energy obtained by some previous works that considermaterial strength and/or damage. On the other hand, our results are roughly consistentwith the values of Q ∗ RD obtained by a purely hydrodynamic target (Genda et al. 2015;Movshovitz et al. 2016), although their numerical scheme is different from ours. Thus, itseems that the dependence of the values of the critical specific impact energy on numericalmethods becomes small in the case of high numerical resolution.Figure 9 shows the dependence of Q ∗ RD , on impact velocity v imp ( v imp = 3, 5, 7, and9 km/s). In order to examine impact velocity dependence, the target size is fixed at 100km. The values of Q ∗ RD , somewhat depend on impact velocities in the range of v imp = 3- 9 km/s. When the impact velocities are rather high, impact energies would be easilyconverted into internal energy (Housen & Holsapple 1990). As a result, with increasingimpact velocity, the kinetic energy of ejecta decreases, which leads to the disruption of thetarget being hindered. From our numerical results, Q ∗ RD , is proportional to v . .Here, we compare our results with the scaling law. In the case of the gravity regime, M lrg /M tot is described by a function of the scaling parameter (Housen & Holsapple 1990),Π G = ( ρG ) − µ/ R − µ tar v µ − Q R , (5)where ρ is the density of the target, and µ is a constant value dependent on material(Table 3). Therefore, M lrg M tot = F (Π G ) , (6)where F (Π G ) is a monotonically increasing function of Π G . For M lrg /M tot = 0 .
5, Π G isequal to a constant value Π ∗ G and also Q R = Q ∗ RD . Hence, from Equation (5), Q ∗ RD can begiven as (Leinhardt & Stewart 2012; Movshovitz et al. 2016), Q ∗ RD = Π ∗ G ( ρG ) µ/ R µ tar v − µ imp . (7)Using Equations (5) and (7), we haveΠ G = Π ∗ G Q R Q ∗ RD . (8)From Equation (7), the dependence of Q ∗ RD on R tar and v imp can be approximated by apower-law given by Q ∗ RD ∝ R p tar v q imp (9)where p and q are fitting parameters. Based on our numerical results shown in Figures 8and 9, the dependence becomes p = 1 .
66 and q = 0 .
48. According to the scaling law, the 11 –values of p and q depend on the value of µ (see Equation (7)). Then, from Figure 8, we canestimate the value of µ to be 0.55, which is in excellent agreement with the value of µ if thetarget is composed of non-porous material (e.g. water and rock) (see Table 3). However, thevalue of µ obtained by Figure 9 becomes smaller than the case with non-porous material.On the other hand, Melosh & Ryan (1997) analytically examined the dependence of Q ∗ RD on R tar and v imp , and they obtained p = 1 . q = 0 .
5. Thus, the target size dependenceagrees well with the scaling law, while the velocity dependence is consistent with analyticconsideration.
4. CONNECTION BETWEEN EROSIVE AND DISRUPTIVE COLLISIONS
In this section, we examine Q R -dependence of the ejecta mass for a wide range ofimpactor radii R imp , including erosive collisions. The numerical methods are basically thesame as those described in Section 2, but slightly modified. In the previous sections, thenumber of cells for the target n tar was fixed because we only consider disruptive collisions inwhich the target is entirely and largely deformed. However, in the case of erosive collision,large deformation due to the impact appears near the impact point and its area depends onimpactor size. Thus, in this section, we fixed the number of cells for the impactor n imp (seeFigure 1), while the value of n tar is changed depending on the impactor size. Typically, weset n imp = 20.First, we examine the dependence of erosive collisions on numerical resolution in thesame way as we did in Section 3. Figure 10(a) shows that mass ejected by erosive collisionsin the case of R tar = 100 km and v imp = 3 km/s as a function of Q R . Blue circles and greensquares represent collision outcomes for n imp = 10, and 20, respectively. We find that themasses ejected by erosive collisions also depend on numerical resolution and become largerdue to the higher numerical resolution. Figure 10(b) shows ejected mass as a function of Q R normalized by each calculated value of Q ∗ RD . We make new estimates for the values of Q ∗ RD based on numerical data in Figure 10(a); Q ∗ RD = 30 . n imp = 10 and20, respectively. Lower resolution simulations result in larger Q ∗ RD . However, we find that M ej /M tot is nicely scaled by Q R /Q ∗ RD (Figure 10(b)). Asterisks in Figure 10(b) representnumerical results obtained by SPH impact simulations with the same impact conditions(Genda et al. 2017). We find that their features agree very well with our results. However,we also note that the material properties of planetesimals are neglected in both simulations.Since the material properties affect collisional outcomes (see Section 5), it is also necessaryto compare outcomes between different impact simulations including material properties,which is out of the scope of this study.Next, we examine the dependence of erosive collisions on target size and impactvelocity. Figure 11 shows the dependence of the ejected mass for five different cases ofimpact conditions as a function of impact energy normalized by each calculated value of 12 – Q ∗ RD . Although the slope of numerical data for the case with low impact velocity ( v imp = 2km/s) is slightly different from the others, M ej /M tot is clearly scaled by Q R /Q ∗ RD . We alsoconfirm that our results are consistent with the dependence on R tar and v imp obtained byGenda et al. (2017). From Equation (8), the ratio Q R /Q ∗ RD is equal to Π G / Π ∗ G . Thus wecan say that the scaling in Figure 11 is equivalent to the scaling law of M lrg (= M tot − M ej )for the gravity regime proposed by Housen & Holsapple (1990). For erosive collisions with Q R /Q ∗ RD ≪
1, the obtained ejected mass is fitted by a linear relation, M ej M tot = ψ (cid:18) Q R Q ∗ RD (cid:19) . (10)From Figures 10(b) and 11, the non-dimensional parameter ψ is obtained as 0.4, independentof the numerical resolution, the target size, and the impact velocity. The value of ψ depends on the impact angle and becomes larger with oblique collisions. According to SPHsimulations by Genda et al. (2017), ψ is 1.2 for the typical oblique (45 ◦ ) collision.The ejected masses from erosive collisions are also described by the crater scaling law(Holsapple 1993; Housen & Holsapple 2011). When the densities of the target and theimpactor are the same, the total mass of fragments with velocity greater than v can begiven by (Housen & Holsapple 2011) M ( v ) M imp = 3 k π C µ (cid:18) vv imp (cid:19) − µ , (11)where k and C are constants whose values are dependent on target material (Table 3). Afragment with a velocity higher than the escape velocity of the target v esc ( ≡ p GM/R tar )is not bound by the target’s gravity. Therefore, M ( v esc ) would correspond to the ejectedmass M ej obtained in this study. In the case of erosive collisions, since the mass of theimpactor is significantly lower than that of the target, the specific impact energy is given as(i.e. classical definition of specific impact energy) Q R ≃ . M imp v /M tot . (12)Using Equations (12) and (11) with M ej = M ( v esc ), M ej can be written as M ej M tot = 3 k π C µ (cid:18) v imp v esc (cid:19) µ − (cid:18) Q R v (cid:19) . (13)Equation (13) is essentially the crater scaling law. Thus we assume that the value of Q R is considerably smaller than Q ∗ RD .Then, since Equation (13) should be equal to Equation (10), we obtain the followingsemi-analytic formula for Q ∗ RD as Q ∗ RD = ψπ k C − µ v (cid:18) v imp v esc (cid:19) − µ , (14) 13 –or Q ∗ RD = ψ k (cid:18) π (cid:19) (3 µ +2) / C − µ ( ρG ) µ/ R µ tar v − µ imp . (15)From Equation (7), Π ∗ G is written byΠ ∗ G = ψ k (cid:18) π (cid:19) (3 µ +2) / C − µ . (16)Although Equation (15) has the same functional form derived by previous studies(e.g. Housen & Holsapple 1990; Leinhardt & Stewart 2012; Movshovitz et al. 2016), thedifference is that, except for ψ , the value of Π ∗ G is determined not by numerical data, butby experimental data. For example, in the case of rock (see, Table 3), Q ∗ RD is written as, Q ∗ RD = 3 . × (cid:18) R tar
100 km (cid:19) . (cid:18) v imp / s (cid:19) . J / kg , (17)where we assume ψ = 0 . Q ∗ RD on target size obtained by Equation (15) forseveral materials. In the case of non-porous material (water and rock), the values of Q ∗ RD for rock are larger than they are in the case with water. This is because the gravitationalpotential increases as the density increases. In comparison with the critical specific valuesobtained by iSALE code and SPH code, the values of Q ∗ RD for the case with rock are roughlyconsistent with the numerical results. In the case of oblique collision we set ψ = 1 .
2, butwe assume that material parameters are unchanged because it is generally thought that thecrater size hardly depends on impact angle, except for in the case of a very low impactangle (Melosh 2011). The value of Q ∗ RD at R tar = 100 km agrees with the numerical dataobtained by Genda et al. (2015), while deviation between semi-analytic and numericalresults for R tar = 10 km becomes large.For sand, the slope becomes smaller compared with non-porous material because thevalue of µ depends on the porosity of the target, and it decreases with increasing degreesof porosity (Table 3). Also, energy dissipation by compaction takes place within poroustargets such as sand. Thus, the value of Q ∗ RD for sand is much higher than it is for the casewith non-porous material. The circles in Figure 12 represent numerical results obtainedby Jutzi et al. (2010), who performed oblique impact simulation using the SPH method,including the effect of porosity. Indeed, their numerical data also reflects the effect ofthe porosity: a small slope and large critical specific impact energy. However, we findthat semi-analytic results are significantly different from Jutzi et al. (2010). Since thedependence of material parameters on the impact angle is small, the deviation betweensemi-analytic and numerical result would be caused by the value of ψ . Therefore, furtherstudies are needed to clarify the effect of material properties on non-dimensional parameter ψ . 14 –
5. DISCUSSION5.1. Effect of oblique impacts
So far, we have focused on the case of head-on collisions with the impact angle θ = 0 ◦ .Here, we discuss the effect of oblique impacts on our results.In the previous section, we have obtained the linear relation between ejecta mass andimpact energy (Equation (10)). Even if ejected mass is averaged over impact angles, thelinear relation holds: ¯ M ej /M tot ∝ Q R / ¯ Q ∗ RD where overlines denote angle-averaged quantities(Genda et al. 2017). Thus, the linear relation for the oblique impact case is written by¯ M ej M tot = ¯ ψ Q R ¯ Q ∗ RD , (18)where ¯ ψ is angle-averaged ψ whose value is 0.88 for a purely hydrodynamic case (seeGenda et al. 2017).As we have shown, in the framework of the crater scaling law, the total mass offragments with velocity greater than v can be scaled by v/v imp (Equation (11)). On theother hand, although oblique impact experiments have not been sufficiently performed, itis thought that the total mass of fragments formed by an oblique impact M ( v, θ ) can bescaled by the normal component of the impact velocity, v imp cos θ (Housen & Holsapple2011). Under this approximation, the crater scaling law for the oblique impact is given by M ( v, θ ) M imp = 3 k π C µ (cid:18) vv imp cos θ (cid:19) − µ . (19)Since the total ejecta mass at an oblique impact with the angle θ , M ej ( θ ), is equal to M ( v esc , θ ), the total ejecta mass M ej ( θ ) is obtained from Equation (19) as M ej ( θ ) M imp = M ej (0 ◦ ) M imp (cos θ ) µ . (20)Using the probability distribution for impact angle (Shoemaker 1962), angle-averagedejected mass ¯ M ej /M imp is given by (see also Genda et al. 2017)¯ M ej M imp = Z ◦ ◦ M ej (0 ◦ ) M imp sin θ (cos θ ) µ +1 dθ = 22 + 3 µ M ej (0 ◦ ) M imp . (21)The amount of mass ejected by head-on erosive collision is reduced by the factor 2 / (2 + 3 µ )due to the effect of impact angle. In the case of a purely hydrodynamic body, the factorbecomes 0.548, which is consistent with numerical results obtained by Genda et al. (2017). 15 –Using Equations (13), (18), and (21), ¯ Q ∗ RD can be given by,¯ Q ∗ RD = 2 + 3 µ ψ k (cid:18) π (cid:19) (3 µ +2) / C − µ ( ρG ) µ/ R µ tar v − µ imp . (22)¯ Q ∗ RD for a rocky planetesimal without material strength is written as,¯ Q ∗ RD = 1 . × (cid:18) R tar
100 km (cid:19) . (cid:18) v imp / s (cid:19) . J / kg . (23)Equation (23) shows that the value of ¯ Q ∗ RD is four times as large as that of head-on Q ∗ RD (see Equation (17)). Indeed, Genda et al. (2017) examined the dependence of Q ∗ RD for thecase of a purely hydrodynamic body on impact angle, and estimated the value of ¯ Q ∗ RD tobe 1 . × J/kg, which is close to the value of ¯ Q ∗ RD obtained by the above semi-analyticformula (Equation (23)). In this study, we had assumed that target and impactor planetesimals are purelyhydrodynamic bodies including self-gravity, but they would have material strength, frictionand compaction. Recently, Jutzi (2015) examined the dependence of collision outcomeson material properties of the target, and showed that the mass of the largest body afterdisruptive collisions becomes considerably large compared to the case with a purelyhydrodynamic body. Although the effect of the material properties has not been takeninto account in this study so far, the iSALE-2D can deal with several energy dissipationmechanisms inside planetary bodies.As a demonstration of this effect, we further perform a simulation of collisionsbetween planetesimals including the effects of material strength and damage; the parametervalues are listed in Table 4. Figure 13 shows the difference in M ej /M tot between thepurely hydrodynamic case and the case with material strength and damage. We usethe same impact conditions as in Figure 2, except for the setting for material propertiesof planetesimals. Although the impact energy is unchanged, the ejected mass becomessignificantly smaller due to energy dissipation. The value of M ej /M tot eventually convergesto M ej ≃ . M tot at t = 350 s. This result means that the value of Q ∗ RD substantiallyincreases due to material strength and damage (Jutzi 2015). In addition to such effects,the effect of porosity also plays an important role for the determination of Q ∗ RD . It isgenerally thought that the value of Q ∗ RD becomes large via impact energy dissipation due tocompaction (Jutzi 2015). However, since the effect of porosity depends on the density andthe strength of the target, Q ∗ RD would become small due to inefficient reaccumulation of aporous target after collision (Jutzi et al. 2010). Therefore, we will investigate the effects ofmaterial properties on Q ∗ RD in the future. 16 –
6. SUMMARY
In this study, we have performed head-on impact simulations of purely hydrodynamicplanetesimals in the gravity regime by using shock-physics code iSALE-2D, and have madea comparison of collisional outcomes between the SPH code and iSALE-2D code. We foundour numerical simulation results agree well with those obtained by the SPH simulation.Detailed analysis gives the following three results. • The value of Q ∗ RD depends on numerical resolution. With decreasing the spatial cellsize ∆ x , the differences in the values of Q ∗ RD between numerical resolutions linearlydecrease and would be close to convergence. Thus, the converged value of Q ∗ RD at∆ x → Q ∗ RD for each numerical resolution. • The converged value Q ∗ RD , obtained by the iSALE code is similar to the case ofthe SPH code, and the difference in Q ∗ RD , between them is within a 30% range ofvariation. • The relationship between ejected mass normalized by total mass ( M ej /M tot ) andimpact energy Q R generally depends on various impact conditions. However, when Q R is scaled by Q ∗ RD that is calculated for each impact simulation, the relationship isindependent of numerical resolution, impact velocity and target size. This similaritylaw for Q R /Q ∗ RD is confirmed for a wide range of specific impact energy from disruptivecollisions to erosive collisions.Using the above similarity law and the crater scaling law, we obtained a semi-analyticformula for the critical specific impact energy Q ∗ RD (Equation (15)). In the case of anon-porous object, the values of Q ∗ RD estimated by the semi-analytic formula agree withnumerical results obtained by the iSALE code and SPH code. However, the values of Q ∗ RD for porous objects are inconsistent with numerical results from SPH simulation taking intoaccount the effect of porosity (Jutzi et al. 2010). Thus, the value of ψ would depend onmaterial properties, as the value of µ depends on the porosity of the target.As mentioned above, most of our results can reproduce the results obtained fromSPH simulations by Genda et al. (2015, 2017) despite different numerical methods.Especially, the correspondence of Q ∗ RD would help us better understand planet formation.Kobayashi et al. (2010) assumed a simple fragmentation model describing both disruptivecollision and erosive collision, and analytically derived the final mass of protoplanetsformed in the protoplanetary disk. According to the analytic formula, the mass of formedprotoplanets is proportional to Q ∗ RD0 . : a factor of two difference in Q ∗ RD directly affects themass of protoplanets by a factor of 1.8. On the other hand, the final mass of protoplanetscan be determined by a balance between the growth time of embryos and the mass depletiontime in collision cascades using the simple fragmentation model. The mass depletion timeis dominated by erosive collisions where the ejected mass is nicely scaled by Q R /Q ∗ RD . 17 –Thus, the depletion time also depends on the value of Q ∗ RD (In fact, the depletion time isproportional to Q ∗ RD0 . (Kobayashi & Tanaka 2010)). As a result, the determination ofthe values of Q ∗ RD can provide constraints on the formation of planetary bodies.In this study, we showed the correspondence of Q ∗ RD for the case of purely hydrodynamicbodies. To determine a more realistic value of Q ∗ RD for various types of planetesimals, how Q ∗ RD depends on material properties needs to be clarified. Therefore, we will investigate theeffects of material properties in the future.We thank Thomas Davison and an anonymous reviewer for their useful comments onour manuscript. We gratefully acknowledge the developers of iSALE, including GarethCollins, Kai W¨unnemann, Boris Ivanov, Jay Melosh, Dirk Elbeshausen, and ThomasDavison. This work was supported by JSPS KAKENHI Grant Numbers JP26287101 andJP15H03716. H.K. was supported by Grant-in-Aid for Scientific Research (JP17K05632,JP17H01105, JP17H01103) and JSPS Core-to-Core Program “International Networkof Planetary Sciences”. H.G. was also supported by JSPS KAKENHI Grant NumberJP17H02990. 18 – REFERENCES
Amsden, A. A., Ruppel, H. M., Hirt, C. W., 1980. SALE: A simplified ALE computerprogram for fluid flow at all speeds. Los Alamos National Laboratories Report,LA-8095, Los Alamos, New Mexico, 101p.Benz, W., Asphaug, E., 1999. Catastrophic disruptions revisited. Icarus 142, 5-20Canup, R., 2004. Simulations of a late lunar-forming impact. Icarus 168, 433-456Canup, R., Asphaug, E., 2001. Origin of the Moon in a giant impact near the end of theEarth’s formation. Nature 412, 708-712Cintala, M. J., Berthoud, L., H¨orz, F., 1999. Ejection-velocity distributions from impactsinto coarse-grained sand. Meteo. Planet Sci. 34, 605-623Collins, G. S., Elbeshausen, D., W¨unnemann, K., Davison, T. M., Ivanov, B. A, Melosh, H.J., 2016. iSALE: A multi-material, multi-rheology shock physics code for simulatingimpact phenomena in two and three dimensions. iSALE-Dellen manual.Collins, G. S., Kenkmann, T., Osinski, G. R., W¨unnemann, K., 2008. Mid-sized complexcrater formation in mixed crystalline-sedimentary targets: Insight from modelingand observation. Meteo. Planet Sci. 43, 1955-1977Collins, G. S., Melosh, H. J., Ivanov, B. A., 2004. Modeling damage and deformation inimpact simulations. Meteo. Planet. Sci. 39, 217-231Collins, G. S., Melosh, H. J., W¨unnemann, K., 2011. Improvements to the epsilon-alphaporous compaction model for simulating impacts into high-porosity solar systemobjects. Int. J. Impact Eng. 38, 434-439Cremonese, G., Martellato, E., Marzari, Kuhrt, E., Scholten, F., Preusker, F., W¨unnemann,K., Borin, P., Massironi, M., Simioni, E., Ip, W., The Osiris Team., 2012. Hydrocodesimulations of the largest crater on asteroid Lutetia. Planet. Space Sci. 66, 147-154´Cuk, M., Stewart, S. T., 2012. Making the Moon from a fast-spinning Earth: a giant impactfollowed by resonant despinning. Science 338, 1047-1052Davison, T. M., Collins, G. S., Ciesla, F. J., 2010 Numerical modelling of heating in porousplanetesimal collisions. Icarus 208, 468-481Davison, T. M., Ciesla, F. J., Collins, G. S., 2012 Post-Impact Thermal Evolution of PorousPlanetesimals. Geochim. Cosmochim. Acta 95, 252-269Fukuzaki, S., Sekine, Y., Genda, H. Sugita, S., Kadono, T., Matsui, T., 2010. Impact-induced N production from ammonium sulfate: Implications for the origin andevolution of N in Titan’s atmosphere. Icarus 209, 715-722 19 –Gault, D. F., Quaide, W. L., Oberbeck, V. R., 1963. Spray Ejected from the Lunar Surfaceby Meteoroid Impact. NASA Tech. Note D-1767Genda, H., Fujita, T., Kobayashi, H., Tanaka, H., Abe, Y., 2015. Resolution dependence ofdisruptive collisions between planetesimals in the gravity regime. Icarus 262, 58-66Genda, H., Fujita, T., Kobayashi, H., Tanaka, H., Suetsugu, R., Abe, Y., 2017. Impacterosion model for gravity-dominated planetesimals. Icarus 294, 234-246Holsapple, K., 1993. The scaling of impact processes in planetary sciences. Annu. Rev.Earth Planet. Sci. 21, 333-373Hosono, N., Saitoh, T.R., Makino, J., Genda, H., Ida, S., 2016. The giant impact simulationswith density independent smoothed particle hydrodynamics. Icarus 271, 131-157Housen, K., Holsapple, K., 1990. On the fragmentation of asteroids and planetary satellites.Icarus 84, 226-253Housen, K., Holsapple, K., 1999. Scale effects in strength-dominated collisions of rockyasteroids. Icarus 142, 21-33Housen, K., Holsapple, K., 2011. Ejecta from impact craters. Icarus 211, 856-875Ivanov, B. A., Deniem, D., Neukum, G., 1997. Implementation of dynamic strength modelsinto 2D hydrocodes: Applications for atmospheric breakup and impact cratering.Int. J. Impact Eng., 20, 411-430Johnson, B. C., Minton, D. A., Melosh, H. J., Zuber, M. T., 2015. Impact jetting as theorigin of chondrules. Nature 517, 339-341Jutzi, M., Michel, P., Benz, W., Richardson, D. C., 2010. Fragment properties at thecatastrophic disruption threshold: The effect of the parent bodyfs internal structure.Icarus 207, 54-65Jutzi, M., 2015. SPH calculations of asteroid disruptions: The role of pressure dependentfailure models. Planet. Space Sci. 107, 3-9Kadono, T., Suzuki, A., Araki, S., Asada, T., Suetsugu R., Hasegawa, S., 2018 Investigationof impact craters on flat surface of cylindrical targets based on experiments andnumerical simulations, submitted to Planet. Space Sci.Kobayashi, H., Tanaka, H., 2010. Fragmentation model dependence of collision cascades.Icarus 206, 735-746Kobayashi, H., Tanaka, H., Krivov, A. V., Inaba, S., 2010. Planetary growth with collisionalfragmentation and gas drag. Icarus 209, 836-847 20 –Kobayashi, H., Tanaka, H., Krivov, A. V., 2011. Planetary core formation with collisionalfragmentation and atmosphere to form gas giant planets. Astrophys. J. 738, 35-45Kurosawa, K., Okamoto, T., Genda, H., 2018. Hydrocode modeling of the spallation processduring hypervelocity impacts: Implications for the ejection of Martian meteorites.Icarus 301, 219-234.Kurosawa, K., Genda, H., 2018. Effects of friction and plastic deformation inshock-comminuted damaged rocks on impact heating. Geophys. Res. Lett. 45,10.1002/2017GL076285Leinhardt, Z. M., Stewart, S. T., 2009. Full numerical simulations of catastrophic smallbody collisions. Icarus 199, 542-559Leinhardt, Z. M., Stewart, S. T., 2012. Collisions between gravity-dominated bodies. I.Outcome regimes and scaling laws. Astrophys. J. 745, 79-105Lissauer J. J., 1993, Planet formation. Annu. Rev. Astron. Astrophys. 31, 129-174Love, S. G., Ahrens, T. J., 1996. Catastrophic impacts on gravity dominated asteroids.Icarus 124, 141-155Melosh, H. J., 2007. A hydrocode equation of state for SiO . Meteorit. Planet. Sci. 42,2079-2098Melosh, H. J., 2011. Plnanetary surface Processes. Cambridge Univ. Press, CambridgeMelosh, H.J., Ryan, E.V., 1997. Asteroids: Shattered but Not Dispersed. Icarus 129,562-564Melosh, H. J., Ryan, E. V., Asphaug, E., 1992. Dynamic fragmentation in impacts -Hydrocode simulation of laboratory impacts. J. Geophys. Res., 97, 14735-14759Movshovitz, N., Nimmo, F., Korycansky, D. G., Asphaug, E., Owen, J. M., 2016. Impactdisruption of gravity-dominated bodies: New simulation data and scaling. Icarus275, 85-96Nagaki, K., Kadono, T., Sakaiya, T., Kurosawa, K., Hironaka, Y., Shigemori, K., Arakawa,M., 2016. Recovery of entire shocked samples in a range of pressure from ∼
100 GPato Hugoniot elastic limit. Meteo. Planet. Sci. 51, 1153-1162Nakamura, A., Hiraoka, K., Yamashita, Y., Machii, N., 2009. Collisional disruptionexperiments of porous targets. Planet. Space Sci. 57, 111-118Sekine, Y., Genda, H., 2012. Giant impacts in the Saturnian system: A possible origin ofdiversity in the inner mid-sized satellites. Planet. Space Sci. 63, 133-138 21 –Shoemaker, E. M., 1962. Interpretation of lunar craters. In: Kopal, Z. (Ed.), Physics andAstronomy of the Moon. Academic Press, New YorkStewart, S. T., Leinhardt,Z. M., 2009. Velocity-Dependent Catastrophic Disruption Criteriafor Planetesimals. Astrophys. J. 691, L133-L137Schmidt, R. M., Hausen, K. T., 1987. Some recent advances in the scaling of impact andexplosion cratering. Int. J. Impact Eng. 5, 543-560Thompson, S. L., Lauson, H. S., 1972. Improvements in the CHART D radiation-hydrodynamic code III: revised analytic equations of state. Sandia NationalLaboratory Report, SC-RR-710714Tillotson, J.H., 1962. Metallic equations of state for hypervelocity impact. General AtomicRept. GA-3216, 1-142Wakita, S., Matsumoto, Y., Oshino, S., Hasegawa, Y., 2017. Planetesimal collisions as achondrule forming event. Astrophys. J. 834, 125, (8 pp).W¨unnemann, K., Collins, G. S., Melosh, H. J., 2006. A strain-based porosity model for usein hydrocode simulations of impacts and implications for transient crater growth inporous targets. Icarus 180, 514-527This manuscript was prepared with the AAS L A TEX macros v5.2. 22 – ∆ x ImpactorTarget R tar D tar R imp n tar n imp n v n h v imp Fig. 1.— Schematic illustration of our numerical setting. In a two-dimensional cylindricalcoordinate system, an impactor with velocity v imp has a head-on collision with a target. Theradii of the target and impactor are R tar and R imp , respectively. 23 –Fig. 2.— Time series of a simulation of a head-on impact between a target with R tar = 100km and an impactor with R imp = 16 km. The impact velocity and the number of cells pertarget radius are set to be 3 km/s, and 800, respectively. The color contour represents thespecific kinetic energy [J/kg]. 24 – M e j / M t o t Time [s]aR imp =16km n tar =100n tar =200n tar =400n tar =800 0 0.2 0.4 0.6 0.8 1 0 50 100 150 200 250 300 350 M e j / M t o t Time [s]bR imp =19km n tar =100n tar =200n tar =400n tar =800
Fig. 3.— The time evolution of ejected mass normalized by the total mass ( M ej /M tot )for R imp = 16 km (a) and 19 km (b). Blue, green, red and black represent the cases with n tar = 100, 200, 400, and 800, respectively. Black curve in Panel (a) corresponds to thecollision shown in Figure 2. 25 – M e j / M t o t Q R [kJ/kg]a n tar =100n tar =200n tar =400n tar =800
20 25 0 0.002 0.004 0.006 Q R D * [ k J / k g ] ∆ x/D tar b Fig. 4.— Determination of Q ∗ RD value and its dependence on the numerical resolution. Panel(a) shows normalized mass of ejecta due to disruptive collisions as a function of the impactenergy. Diamonds, squares, triangles, and circles represent the cases with n tar = 100, 200,400, and 800, respectively. The data listed in Table 2 are used. The vertical dashed linesrepresent the critical specific impact energies Q ∗ RD for each resolution. Panel (b) shows thecritical specific impact energy Q ∗ RD as a function of the spatial cell size divided by the target’sdiameter (∆ x/D tar ). The dashed line is the least-squares fit to the results (Equation (3)). 26 – Q R D * [ k J / k g ] ∆ x/D tar This studyGenda et al. (2015)
Fig. 5.— Comparison of resolution dependence of Q ∗ RD between this study and Genda et al.(2015). The lines are the least-squares fit to the results (Equation (3)). 27 –
0 0.5 1 1.5 2 2.5 M e j / M t o t Q R /Q RD * n tar =100Genda et al. (2015)n tar =200n tar =400n tar =800 Fig. 6.— Same as Figure 4(a), but Q R is normalized by each Q ∗ RD . Asterisks representnumerical data obtained by Genda et al. (2015). Dashed line represents M ej /M tot obtainedby the empirical law (Equation (4)). 28 – Q R D * / U t a r ∆ x/D tar R tar = 30kmR tar =100kmR tar =300km Fig. 7.— Critical specific impact energy Q ∗ RD normalized by the potential energy of thetarget U tar as a function of ∆ x/D tar in the case of v imp = 3 km/s. Circles, squares andtriangles represent R tar = 30, 100 and 300 km, respectively. The lines are the least-squaresfit to the results based on Equation (3). 29 – Target radius (R tar ) [m]
Melosh & Ryan1997Benz & Asphaug1999Genda et al. 2015Movshovitz et al. 2016
Leinhardt & Stewart 2009This study C r i t i c a l s p e c i fi c i m p a c t e n e r g y ( Q R D * o r Q D * ) [ J / k g ] Fig. 8.— Dependence of Q ∗ RD obtained by head-on impact simulations as a function of targetradius. The diamond symbols represent our results, which are converged values ( Q ∗ RD , ) forthe case with 3 km/s. Other symbols are the values of critical specific impact energy ( Q ∗ RD orclassical definition of Q ∗ D (= 0 . M imp v /M tot )) obtained by various previous studies. Graydata symbols represent cases for planetesimals with material strength and/or damage. Blacksymbols represent purely hydrodynamic cases. The dashed line is the least-squares fit to ourresults. 30 – Q R D , * [ J / k g ] v imp [m/s] Fig. 9.— Dependence of converged values Q ∗ RD , on impact velocity ( v imp = 3, 5, 7, and 9km/s) in the case of R tar = 100 km. The dashed line is the least-squares fit to the results. 31 – -3 -2 -1 -1 M e j / M t o t Q R [kJ/kg]a =10n imp =20n imp Genda et al. (2017) -3 -2 -1 -2 -1 M e j / M t o t Q R /Q RD * =10n imp =20n imp Genda et al. (2017) b Fig. 10.— (a) Dependence of normalized mass of ejecta on numerical resolution in the caseof R tar = 100 km and v imp = 3 km/s. Circles and squares represent outcomes of n imp = 10and 20, respectively. Asterisks represent numerical data obtained by Genda et al. (2017).Panel (b) is the same as panel (a), but Q R is normalized by each Q ∗ RD . Dashed and dottedlines represent M ej /M tot obtained by Equation (10) with ψ = 0 . ψ = 1), respectively. 32 – -3 -2 -1 -2 -1 M e j / M t o t Q R /Q RD * =3km/s tar =100kmR imp v =3km/s tar =30kmR imp v =3km/s tar =300kmR imp v =2km/s tar =100kmR imp v =5km/s tar =100kmR imp v Fig. 11.— Dependence of M ej /M tot on Q R /Q ∗ RD for five different impact conditions in thecase of n imp = 20. Squares, circles, triangles, inverted triangles, and diamonds represent theresults for ( v imp , R tar )=(3 km/s, 30 km), (3 km/s, 100 km), (3 km/s, 300 km), (2 km/s,100 km), and (5 km/s, 100 km), respectively. Dashed and dotted lines represent M ej /M tot obtained by Equation (10) with ψ = 0 . ψ = 1), respectively. 33 – C r i t i c a l s p e c i fi c i m p a c t e n e r g y [ J / k g ] Target radius [m]Jutzi et al. 2010This studyGenda et al. 2015SandRock Water
Fig. 12.— Comparison of Q ∗ RD between semi-analytic formula and numerical results. Blue,red, and green lines correspond to Q ∗ RD for the case with water, rock and sand, respectively.Solid and dotted lines represent head-on collisions ( ψ = 0 .
4) and oblique (45 ◦ ) collisions( ψ = 1 . Q ∗ RD for three different numerical data. Filled and opensymbols represent head-on collisions and oblique (45 ◦ ) collisions, respectively. 34 – M e j / M t o t Time [s]
Fig. 13.— Comparison of M ej /M tot between the purely hydrodynamic case and the casewith material properties. Except for the setting for material properties, impact conditionsin Figure 2 are used. Solid and dotted curves represent the case with material strength anddamage (see also Table 4), and the case of a purely hydrodynamic body shown in Figure 2,respectively. 35 –Table 1: iSALE input parametersDescription ValuesCell per target radius ( n tar ) 100 200 400 800Horizontal cells ( n h ) 450 900 1800 3600Vertical cells ( n v ) 450 900 1800 3600Setup type PLANET a Surface temperature [K] 293 b Gradient type SELF a Gradient dimension 3Self-gravity accuracy parameter ( θ ) 1.0 or 0.5 ( R tar = 300 km)Self-gravity update frequency 10Volume fraction cutoff 10 − Density cutoff [kg/m ] 5 b Velocity cutoff 1.697 × v impb Linear term of artificial viscosity 0.24 b Quadratic term of artificial viscosity 1.2 ba See Collins et al. (2016) for detailed description. b Default parameter values of the iSALE code (See Collins et al.(2016)).
Table 2: Collision outcomes R tar = 100km, M tar = 1 . × kg, v imp = 3 km/s M imp [kg] Q R [J/kg] n tar = 100 n tar = 200 n tar = 400 n tar = 800 n imp M ej /M tot n imp M ej /M tot n imp M ej /M tot n imp M ej /M tot . × . ×
14 1 . × −
28 2 . × −
56 2 . × −
112 2 . × − . × . ×
15 2 . × −
30 2 . × −
60 3 . × −
120 3 . × − . × . ×
16 2 . × −
32 4 . × −
64 4 . × −
128 4 . × − . × . ×
17 4 . × −
34 5 . × −
68 5 . × −
136 5 . × − . × . ×
18 5 . × −
36 6 . × −
72 6 . × −
144 6 . × − . × . ×
19 6 . × −
38 7 . × −
76 7 . × −
152 7 . × − . × . ×
20 8 . × −
40 8 . × −
80 8 . × −
160 8 . × − . × . ×
21 8 . × −
42 9 . × −
84 9 . × −
168 9 . × −
37 –Table 3: Constants in erosive anddisruptive collisionsConstants Water Rock Sand µ a b c C a b c k a b c ρ [kg/m ] 1000 a b c Π ∗ G /ψ a Schmidt & Hausen (1987);Housen & Holsapple (2011) b Gault et al. (1963) c Cintala et al. (1999)Table 4: iSALE material parametersDescription Values a Poisson’s ratio 0.25Specific heat capacity [J/(kg K)] 10 Strength model Rock b Cohesion (undamaged) [Pa] 2 . × Coefficient of internal friction (undamaged) 2 . × Cohesion (damaged) [Pa] 10 Coefficient of internal friction (damaged) 0.4Limiting strength at high pressure (damaged) [Pa] 2 . × Damage model Ivanov c Minimum failure strain 10 − Damage model constant 10 − Threshold pressure for damage model [Pa] 3 × Parameter values generated by iSALEMat. b See Ivanov et al. (1997) for detailed description. cc