A Dust Aggregate Model
Based onNumerical Simulations of Aggregate Collisions

ダストアグリゲイト同士のDEM(?)衝突シミュレーション
和田 浩二
千葉工業大学 惑星探査研究センター
田中秀和1，陶山徹2，木村宏3，山本哲生1
1北大低温研，2新潟県立自然科学館，3CPS

A collision of BCCAs
8192+8192 ice particles (r=0.1µm, ξc = 8Å)Collision velocity = 22 m/s

Background
Collisional growth of dust(< µm)?

Planetesimal formation(> km)

Structure evolution of dust aggregates in protoplanetary disks

When and how are aggregates compressed and/or disrupted ?
Numericalsimulation of dust aggregate collisions!
ダストアグリゲイト衝突の数値シミュレーション！

©NASABallistic Cluster-Cluster Aggregation (BCCA)
In the early growth stage, undeformed BCCAs are formed
because of their low collision velocity (< mm/s)
 A series of hit-and sticks ofcomparable aggregates Fluffy structure (fractal dimension < ~ 2)

How are the BCCA structures compressed?
Dominik & Tielens 1997;
Wada et al. 2007, 2008; Suyama et al. 2008

Background
Collision velocity of dust
in protoplanetary disks< several 10 m/se.g., ＜～50 m/s (Hayashi model, without turbulence)

Is it possible for dust to grow through collisions?
To what extent is dust compressed?
Experimental: Blum & Wurm 2000, Wurm et al. 2005
Numerical: Dominik & Tielens 1997, Wada et al. 2008; 2009

Ballistic Particle-Cluster Aggregation (BPCA)
 Formed by one-by-one sticking of monomers

 Compact structure (fractal dimension ~ 3)Dust should be compact
in high velocity collisions causing their disruption
Collisions of BPCA clusters
implication for growth and disruption of compact dust

Objective
To construct a structural evolution model of dust aggregatesby numerical simulations of aggregate collisionsCollisions of BCCA & BPCA clusters
 Compression process（BCCAs）

Gyration radius → Degree of compression

 High-velocity collisions（BPCAs）

Number of particles in the largest remnant→ Growth efficiency
Coordination numbers in the largest remnant→ Degree of compression

Simulation Method

Grain interaction model

Johnson, Kendall and Roberts (1971)
Johnson (1987), Chokshi et al. (1993)
Dominik and Tielens (1995,96)
Wada et al. (2007)

Elastic spheres having surface energy
NormalSliding

rollingtwistingξδsϕ

JKR theoryUnδ

Critical sticking velocity:exp.~10×theo.!?

(Dominik & Tielens 1996)(Dominik & Tielens 1995)

JKR and rolling resistance have been tested with experiments using ~1µm
SiO2 particles. (Heim et al. 1999; Poppe et al. 2000; Blum & Wurm 2000)

Grain interaction model

Johnson, Kendall and Roberts (1971)
Johnson (1987), Chokshi et al. (1993)
Dominik and Tielens (1995,96)
Wada et al. (2007)

Elastic spheres having surface energy
NormalSliding

twistingrollingξδsContact & Separations, ξ, ϕ > critical displacements
Critical slide
Critical roll
Critical twist

ϕEnergy dissipation

scrit ~ 1.5 Å (for 0.2 µm quartz)
ξcrit ~ 2 Å (or ~30 Å (Heim et al.,1999))
ϕcrit ~ 1°

Ebreak: Energy to break a contact
Eroll : Energy to roll a pair of gains by 90˚A classical study

Dominik and Tielens (1997)
Each grain motion is directly calculated,taking into account particle interactions
modeling grain interactions seriously

Limitations:
D&T “recipe”
 2-D, Head-on collision~ nk Eroll

Max. compression> 10 nk Ebreak

Catastrophic disruption

Eimpact =  Small size （40＋40 grains）
Confirmed by experiments(Blum & Wurm 2000)

 Initial structure: only 1 type

Eroll : Energy to roll a grain by 90°
Ebreak : Energy to break a contact
nk : Number of contacts in initial aggregatesCollisions between BCCA clusters: Compression process

Initial Conditions and Parameters
Collisions of BCCA clusters

Results areaveraged

 BCCA clusters are• composed of 512, 2048, or 8192 particles (10 types randomly produced)• impacted by head-on collision

512+5122048+2048

particle : radius = 0.1 µm,Ice(E = 7 GPa, ν = 0.25, γ = 100 mJ/m2)
SiO2 (E = 54 GPa, ν = 0.17, γ = 25 mJ/m2)
Critical rolling displacement : ξcrit = 2, 8, 30 Å8192+8192

Example of simulations
Ice, 8192 + 8192, ξcrit = 8 Å

Eimpact ~ 0.7 Eroll

Eimpact ~ 0.3 nk Eroll

Eimpact ~ 13 nk EbreakVimpact = 0.2 m/s

Vimpact = 17 m/s

Vimpact = 39 m/s

Numerical Results on
Gyration Radius

Gyration radius rg : compression process
Eimpact ~ 0.01 Eroll
Impact velocity: 0.024 m/s

Eimpact ~ 0.19 N Eroll
Impact velocity: 13 m/s

Ice, 8192 + 8192, ξcrit = 8 Årg =12xi − x g∑
N ixg: center of massGyration radius rg : compression processrg / (r1 N1/3)

10SiO2, N = 1024, ξc = 2 Åξc = 8 Åξc = 16 Å
Ice, N = 1024, ξc = 2 Åξc = 8 Åξc = 16 Å
Ice, N = 4096, ξc = 8 Å
N = 16384, ξc = 8 Å1
1e-6 1e-5 1e-4 0.001 0.01 0.1

Eimpact / (N Eroll)1Eimpact ~ (0.1 – 1) N Eroll
Max. compression
Consistent withDominik & Tielens (1997)

10100

Gyration radius rg : compression processrg / (r1 N1/2.5)

5

SiO2, N = 1024, ξc = 2 Åξc = 8 Åξc = 16 ÅIce, N = 1024, ξc = 2 Åξc = 8 Åξc = 16 ÅIce, N = 4096, ξc = 8 Å
N = 16384, ξc = 8 Å

Scaled by
Eimpact / (N Eroll)rg is normalized by

r 1N

11e-6 1e-5 1e-4 0.001 0.01 0.1

Eimpact / (N Eroll)

1

10100

12.5

 11
Scaled by− No restructuringslope =− Eimpact / (N Eroll) d f dc rg,BCCA
df :fractal dimension of BCCA
dc: fractal dimension
compression
r ofismax.
normalized
bygrg,compEimp ~ NEroll

Max. compression& Disruption

Eimp ~ Erolllog rgGyration radius rg : compression process

log Eimp

r 1N

12.5

Gyration radius rg : compression processrg / (r1 N1/2.5)

5

SiO2, N = 1024, ξc = 2 Åξc = 8 Åξc = 16 ÅIce, N = 1024, ξc = 2 Åξc = 8 Åξc = 16 ÅIce, N = 4096, ξc = 8 ÅN = 16384, ξc = 8 ÅScaled by
Eimpact / (N Eroll)rg is normalized byr 1N

11e-6 1e-5 1e-4 0.001 0.01 0.1

Eimpact / (N Eroll)

1

10100

12.5

Gyration radius rg : compression processrg / (r1 N1/2.5)

5

SiO2, N = 1024, ξc = 2 Åξc = 8 Åξc = 16 Å
Ice, N = 1024, ξc = 2 Åξc = 8 Åξc = 16 Å
Ice, N = 4096, ξc = 8 Å
N = 16384, ξc = 8 Å

Scaled by
Eimpact / (N Eroll)rg is normalized by

r 1N

12.5

Not fully compressed Eimpact ≈ 0.81
 NEroll r N 2.5
rg1− 0.1

11e-6 1e-5 1e-4 0.001 0.01 0.1

Eimpact / (N Eroll)

110100(df = 2, dc =2.5)

df : fractal dimension of BCCA
dc : fractal dimension of max. compression

The number of particles N(<r) within r in an aggregate

rIce, 8192 + 8192, ξcrit = 8 Å

1d = 2.5

N / (r/r1 )2.5

Increasing V
0.1

0.010.001
110100r/r

1000

Successive collisions in a BCCA mode
Suyama et al. 2008
Fractal dimension ~2.5
Decrease in density

CG by Dr. T. Takeda, 4D2Uproject, NAOJ

Summary of Compression Process
3D BCCA clusters (df ~ 2) are not fully compressed•Fractal dimension for max. compression : dc ~ 2.5•Gyration radius：

1- 0.1rg ~ r1N 2.5 [ Eimpact / (nkEroll) ]rg ~ r1N 1/dc[ Eimpact / (NEroll) ] - (1/df – 1/dc)
Successive collisions also lead to dc~2.5
The results for single collisions are applicable.

Collisions between BPCA clusters: High-velocity collisions

Initial Conditions and Parameters
Collisions of BPCA clusters
 BPCA clusters are:Results areaveraged• composed of 500, 2000, or 8000 particles (3 types randomly produced)• Impact parameter： b(defined by using characteristic radius rc）

b=00.19
0.39
0.49
0.58
0.69
0.77
0.89
1.00
l

b=ｒｃ,bｒｃ,arc =lrc ,a + rc ,b

5 1( xi − xG ) 2∑
3N

8000+8000

Ice (E = 7.0×1010 Pa, ν = 0.25, γ = 100 mJ/m2, R = 0.1µm) , critical rolling displace. ξcrit = 8ÅImpact velocity vimp= 6 - 260 m/s

A collision of BPCAs
8000+8000 ice particles (r=0.1µm, ξc = 8Å)Collision velocity = 57 m/s

CG by Dr. T. Takeda,
4D2Uproject, NAOJ

Collisions of BPCA clusters

N=8000+8000, ice, ξc = 8Å, vimp = 70 m/s ( Eimp = 42 NEbreak )b=0

b = 0.39

b = 0.69

b = 1.00

Degree of Disruption:
Growth Efficiency

Largest fragment mass Nlarge: growth efficiency
vimp [m/s] (ice)110100f ≡ Nlarge / Ntotal

f = Nlarge / Ntotal

b = 0.00: growth efficiency
f > 0.5 + growth

8000+80000.8f < 0.50.62000+2000

dependent on N

0.4500+500
0.2

00.1- growth

110Eimpact / (N Ebreak)

100

1000

Largest fragment mass Nlarge: growth efficiency
vimp [m/s] (ice)110100f = Nlarge / Ntotal

b = 0.190.8: growth efficiency
f > 0.5 + growthf < 0.5

0.6

0.4

0.2

00.1

f ≡ Nlarge / Ntotal

110Eimpact / (N Ebreak)

100

1000

- growth

Largest fragment mass Nlarge: growth efficiency
vimp [m/s] (ice)110100f = Nlarge / Ntotal

b = 0.39
0.8: growth efficiencyf > 0.5 + growthf < 0.50.60.40.200.1

f ≡ Nlarge / Ntotal

110Eimpact / (N Ebreak)

100

1000

- growth

Largest fragment mass Nlarge: growth efficiency
vimp [m/s] (ice)110100f = Nlarge / Ntotal

b = 0.49
0.8: growth efficiency
f > 0.5 + growthf < 0.50.60.40.200.1f ≡ Nlarge / Ntotal

110Eimpact / (N Ebreak)

100

1000

- growth

Largest fragment mass Nlarge: growth efficiency
vimp [m/s] (ice)110100f = Nlarge / Ntotal

b = 0.580.8: growth efficiency
f > 0.5 + growthf < 0.50.6

0.40.200.1

f ≡ Nlarge / Ntotal

110Eimpact / (N Ebreak)

100

1000

- growth

Largest fragment mass Nlarge: growth efficiency
vimp [m/s] (ice)110100f = Nlarge / Ntotal

b = 0.690.8: growth efficiency
f > 0.5 + growthf < 0.5

0.6

0.4

0.2

00.1

f ≡ Nlarge / Ntotal

110Eimpact / (N Ebreak)100

1000- growthLargest fragment mass Nlarge: growth efficiency
vimp [m/s] (ice)110100f = Nlarge / Ntotal

b = 0.770.8: growth efficiency
f > 0.5 + growthf < 0.50.60.40.200.1

f ≡ Nlarge / Ntotal

110Eimpact / (N Ebreak)100

1000- growthLargest fragment mass Nlarge: growth efficiency
vimp [m/s] (ice)110100f = Nlarge / Ntotal

b = 0.89
0.8: growth efficiency
f > 0.5 + growthf < 0.50.60.40.200.1

f ≡ Nlarge / Ntotal

110Eimpact / (N Ebreak)

100

1000

- growth

Largest fragment mass Nlarge: growth efficiency
vimp [m/s] (ice)110100f = Nlarge / Ntotal

b = 1.000.8f ≡ Nlarge / Ntotal: growth efficiency
f > 0.5 + growthf < 0.5

0.6

0.4- growth

Offset collisions

0.2

independent of N
00.1

110

Eimpact / (N Ebreak)

100

1000

Largest fragment mass Nlarge: growth efficiency
vimp [m/s] (ice)110100f = Nlarge / Ntotal

8000+8000
0.8: growth efficiency
f > 0.5 + growthf < 0.50.60.4- growth

Averageweighted by b20.200.1

f ≡ Nlarge / Ntotal

110Eimpact / (N Ebreak)

100

1000

Growth efficiency averaged
vimp [m/s] (ice)110

Averaged for b2100

f = Nlarge / Ntotal

8000+8000
0.8: growth efficiency
f > 0.5 + growthf < 0.50.60.40.200.1

f ≡ Nlarge / Ntotal

110

Eimpact / (N Ebreak)

100

1000

- growth

Growth efficiency averaged
vimp [m/s] (ice)f = Nlarge / Ntotal

110Averaged for b2100f ≡ Nlarge / Ntotal: growth efficiency
f > 0.5 + growth~ 60 m/s

0.8

f < 0.50.6- growth

0.4

small dependence on N
0.200.1

1

10Eimpact / (N Ebreak)1001000

Degree of compression:
Coordination number

Coordination number Nc
Number of particles in contact with a particlee.g.,
Nc = 4
Nc = 2 for BCCA and BPCA
Max. Nc = 12 for close-packing
An index of compression:
The more compact are aggregates, the larger Nc is.

What value of Nc is achieved at BPCA collisions?

Mean Coordination Number

Coordination number Nc @ BPCA collisions
410ucol [m/s] 100 Ice, 8Å, 8000+8000

3.6

3.2

2.8

2.4

20.1

1

10

100

Eimpact / (N Ebreak)1000

Why Nc = 4 ?
Particles are stable enough with Nc = 4 in 3D:

1D

2D

3D

SummaryDust aggregates remain fluffly only through collisions.
Fractal dimension ~ 2.5Coordination number < 4Very fluffy planetesimals could be formed !?
~10-4 g/cc (Suyama et al. 2008)
Other compression processes are required.

Icy aggregates can grow at collision velocity ~ 50 m/s.
Planetesimals can be formed through collisions of dust.