5. Dust Coagulation
DustPy calculates the collisional evolution of dust by solving the Smolukowski equation
\(\begin{split} \frac{\partial}{\partial t} f\left( m \right) = &\int\limits_0^\infty \int\limits_0^{m'} K\left( m, m', m'' \right) f\left( m' \right) f\left( m'' \right) R\left( m', m'' \right) \mathrm{d}m'' \mathrm{d}m'\\ &-f\left( m \right) \int\limits_0^\infty f\left( m' \right) R\left( m, m' \right) \mathrm{d}m'. \end{split}\)
Where \(f\left( m \right)\) is the dust mass distribution, \(R\) the collision rates, and \(K\) defines the collision outcomes.
The first integral sums over all possible particle collisions and \(K\left( m, m', m'' \right)\) defines the amount that gets added to \(f\left( m \right)\) if particles \(m'\) and \(m''\) collide. The second integral removes that particles from the distribution that collided with other particles.
This equation can be generalized to
\(\frac{\partial}{\partial t} f\left( m \right) = \int\limits_0^\infty \int\limits_0^{m'} M\left( m, m', m'' \right) f\left( m' \right) f\left( m'' \right) R\left( m', m'' \right) \mathrm{d}m'' \mathrm{d}m'\)
where the matrix \(M\) has both the positive and the negative terms. Also note that the integral over \(m''\) only goes up to mass \(m'\). Since the collisions are symmetric this saves some computation.
For details on the numerical implementation, we refer to Brauer et al. (2008) and Birnstiel et al. (2010).
The standard model of DustPy knows three different collision outcomes: sticking, erosion, and full fragmentation.
[1]:
from dustpy import Simulation
sim = Simulation()
sim.initialize()
The parameters that controll the collisional behavior are pre-calculated when calling Simulation.initialze() and stored in Simulation.dust.coagulation.
[2]:
sim.dust.coagulation
[2]:
Group (Coagulation quantities)
------------------------------
A : Field (Fragment normalization factors), constant
eps : Field (Remnant mass distribution), constant
lf_ind : Field (Index of largest fragment), constant
phi : Field (Fragment distribution), constant
rm_ind : Field (Smaller index of remnant), constant
stick : Field (Sticking matrix), constant
stick_ind : Field (Non-zero elements of sticking matrix), constant
-----
These fields are strictly tied to the mass grid. Neither those fields nor the mass grid should be changed after calling Simulation.initialize().
Sticking
Sticking means that two particles collide and form a new larger particles conserving the total mass. There are two complications when implementing this numerically:
First, since we need to cover a high dynamic range from particles masses around \(10^{-15}\,\mathrm{g}\) to more than \(10^{15}\,\mathrm{g}\) while still having a reasonable number of mass bins, the mass grid has to be logarithmic. This has the disadvantage that the new particle mass resulting from a collision will not directly lie on the mass grid, but will fall in between two grid masses. The resulting mass has to be therefore distributed between the two adjacent grid cells. This means there will be density added to a mass bin that is larger than the actual mass of the resulting particle. This can have unwanted side effects if the mass grid is too coarse. Please see Drążkowska et al. (2014) and the notebook with the analytical coagulation kernels for details on the minimum mass resulution that should be used.
Second, if the masses of the two colliding particles differ by more than \(15\) orders of magnitude (double-precision floating-point format) the resulting mass is (for a computer) identical to the mass of the larger particle.
[3]:
m1 = 1.e-16
m2 = 1.
m_tot = m1 + m2
m_tot == m2
[3]:
True
The mass error if this happens is small for a single collision, but it effectively halts the growth of larger particles that primarily grow by sweeping up many smaller particles. There are techniques to prevent this from happening by re-arranging sums. Please have a look at apendices A and B of Brauer et al. (2008) for a detailed description of this mechanism.
The matrix \(M\) has for sticking at maximum four non-zero values for any pair of particle collision: two negative terms for the colliding particles that get removed from the distribution and two positive terms for the two mass bins that border the mass of the resulting particle. If two particles of identical mass collide \(M\) has only three non-zero values: one negative and two positives.
Since \(M\) is in the case of sticking a sparse matrix that contains mostly zeroes, we only store the non-zero elements in Simulation.dust.coagulation.stick and Simulation.dust.coagulation.stick_ind.
The first field is the four non-zero elements of \(M\), while the latter field stores the positions of those elements within \(M\).
Let’s look at a collision of two particles i and j.
[4]:
i = 5
j = 3
The non-zero elemets of \(M\) are
[5]:
sim.dust.coagulation.stick[:, j, i]
[5]:
[-1. -1. 0.76265439 0.23734561]
at the positions
[6]:
sim.dust.coagulation.stick_ind[:, j, i]
[6]:
[3 5 6 7]
As you can see, we remove one particle each from bins 3 and 5 and add about \(75\,\%\) of a particle to mass bin 6 and about \(25\,\%\) of a particle to mass bin 7.
To check for mass conservation we have to multiply the non-zero elements in \(M\) with their respective masses and sum up all elements.
[7]:
m_tot = sim.grid.m[i] + sim.grid.m[j]
err = (sim.dust.coagulation.stick[:, j, i] * sim.grid.m[sim.dust.coagulation.stick_ind[:, j, i]]).sum()
rel_err = err/m_tot
print("Mass error: {:7.1e} g\nRel. error: {:7.1e}".format(err, rel_err))
Mass error: 6.5e-27 g
Rel. error: 8.2e-16
The relative mass error is of the order of machine precision. This as good as it can get.
Note: Since the collsions are symmetric (collision of i with j is identical to collision of j with i), we do not use collisions with i>j.
[8]:
sim.dust.coagulation.stick[:, i, j]
[8]:
[0. 0. 0. 0.]
If two identical particles collide there are only three non-zero elements.
[9]:
sim.dust.coagulation.stick[:, j, j]
[9]:
[-2. 0.90784249 0.09215751 0. ]
[10]:
sim.dust.coagulation.stick_ind[:, j, j]
[10]:
[ 3 5 6 -1]
In this case we remove two particles from mass bin \(3\), add \(90\,\%\) of a particle to mass bin \(5\) and \(10\,\%\) of a particle to mass bin 6. And index of \(-1\) means that it should be ignored.
In any case the sum over Simulation.dust.coagulation.stick for any possible collision should be \(-1\). This means that for every collision the total number of particles is reduced by \(1\). Two particles collide and form one new particle.
[11]:
sim.dust.coagulation.stick[:, j, i].sum()
[11]:
-1.0000000000000004
Fragmentation & Erosion
If the collision velocity exceeds the fragmentation velocity particles fragment. Here we distinguish two cases: Either the particles are roughly equal sized, then both particles fully fragment. Or one particle is significantly larger than the other, then only the smaller particle fragments, but it is chipping off some mass of the larger particle. The latter case is called erosion.
The threshold of these two cases can be set with Simulation.ini.dust.erosionMassRatio. By default this is \(10\), meaning both particles fully fragment if their mass ratio is smaller than \(10\). Otherwise we have erosion. This value has to be set, before calling Simulation.initialize().
For fragmentation we follow the method described in Rafikov et al. (2020), who developed an \(\mathcal{O} \left( N_m^2 \right)\) algorithm for fragmentation. We modified it, however, to make it strictly mass conserving.
In either case, full fragmentation and erosion, we produce a distribution of fragments of kind
\(n\left( m \right) \mathrm{d}m \propto m^\gamma \mathrm{d}m\).
\(\gamma\) has to be set in Simulation.ini.dust.fragmentDistribution before calling Simulation.initialize(). By default it is \(-\frac{11}{6}\) taken from Dohnanyi (1969). From this we can calculate a normalized distribution of fragments with a given maximum fragment size.
This, for example, is the fragment distribution up to mass bin \(15\). It is similar to the \(\varphi\) in Rafikov et al. (2020)
[12]:
sim.dust.coagulation.phi[15, :]
[12]:
[0.04013542 0.04239721 0.04478647 0.04731037 0.0499765 0.05279288
0.05576798 0.05891073 0.06223059 0.06573754 0.06944212 0.07335547
0.07748936 0.0818562 0.08646913 0.09134202 0. 0.
0. 0. 0. 0. 0. 0.
0. 0. 0. 0. 0. 0.
0. 0. 0. 0. 0. 0.
0. 0. 0. 0. 0. 0.
0. 0. 0. 0. 0. 0.
0. 0. 0. 0. 0. 0.
0. 0. 0. 0. 0. 0.
0. 0. 0. 0. 0. 0.
0. 0. 0. 0. 0. 0.
0. 0. 0. 0. 0. 0.
0. 0. 0. 0. 0. 0.
0. 0. 0. 0. 0. 0.
0. 0. 0. 0. 0. 0.
0. 0. 0. 0. 0. 0.
0. 0. 0. 0. 0. 0.
0. 0. 0. 0. 0. 0.
0. 0. 0. 0. 0. 0. ]
that is normalized to \(1\).
[13]:
sim.dust.coagulation.phi[15, :].sum()
[13]:
0.9999999999999999
That means it can be simply multiplied by any desired total fragment mass.
You may have noticed that the exponent of the fragment distribution is not identical to \(\gamma\) from above. Our \(\varphi\) here is already tranformed to code units, meaning the mass density integrated over the mass bin. Since the mass bin widths increase with increasing masses, this changes the slope.
In a collision resulting in full fragmentation, the largest fragment size is equal to the larger collision partner. In an erosive collision the largest fragment has the size of the smaller collision partner. The index of the largest fragment is stored in Simulation.dust.coagulation.lf_ind.
[14]:
sim.dust.coagulation.lf_ind[j, i]
[14]:
5
In this example, this is equal to the index of the larger particle
[15]:
sim.dust.coagulation.lf_ind[j, i] == i
[15]:
True
suggesting that this is a collision that would result in full fragmentation (if the collision velocity is high enough).
Again, note that this is symmetric and the transpose would return \(-1\) indicating that it is not used.
[16]:
sim.dust.coagulation.lf_ind[i, j]
[16]:
-1
Simulation.dust.coagulation.A stores the total fragment mass resulting from a collision of particles i and j. In a collision resulting in full fragmentation the total fragment mass is identical to the total mass of both collision partners.
[17]:
m_frag = sim.grid.m[i] + sim.grid.m[j]
sim.dust.coagulation.A[j, i] == m_frag
[17]:
True
In an erosive collision the total mass of fragments is given by
\(m_\mathrm{fragments} = \left( 1 + \chi \right)m_\mathrm{small}\),
where \(\chi\) is the amount of mass that gets excavated from the larger collision partner in units of the smaller collision partner. This value can be set in Simulation.ini.dust.excavatedMass and is by default \(1\). This value has to be set before calling Simulation.initialize().
[18]:
i = 50
j = 3
m_frag = 2. * sim.grid.m[j]
sim.dust.coagulation.A[j, i] == m_frag
[18]:
True
In an erosive collision event the larger particle loses mass, which has to be taken into account. Since the new mass of the larger particle after an erosive collision will most likely not directly fall onto a mass bin, it has to be distributed between the two adjacent mass bins. Simulation.dust.coagulation.rm_ind is the lower mass index into which the remnant mass is distributed, while Simulation.dust.coagulation.eps is the fraction of mass that is distributed into that bin. That means
that the fraction of \(1\) - Simulation.dust.coagulation.eps is added into mass bin Simulation.dust.coagulation.rm_ind + \(1\).
In that way it is possible to check for mass conservation in erosive events.
[19]:
# Total mass of collision partners
m_tot = sim.grid.m[i] + sim.grid.m[j]
# Index of larges fragment
k_lf = sim.dust.coagulation.lf_ind[j, i]
# Smaller index of remnant mass
k_rm = sim.dust.coagulation.rm_ind[j, i]
# Fraction of mass in k_rm
eps = sim.dust.coagulation.eps[j, i]
# Mass of fragments
m_frag = (sim.dust.coagulation.A[j, i] * sim.dust.coagulation.phi[k_lf, :]).sum()
# Mass of remnant particle
m_remnant = eps*sim.grid.m[k_rm] + (1.-eps)*sim.grid.m[k_rm+1]
err = m_frag + m_remnant - m_tot
rel_err = err/m_tot
print("Mass error: {:7.1e} g\nRel. error: {:7.1e}".format(err, rel_err))
Mass error: 1.7e-21 g
Rel. error: 1.2e-16
As you can see this is accurate up to machine precision.
When starting a simulation with Simulation.run(), DustPy will perform a quick check for mass conservation by calculating the relative error for every possible particle collision and every possible collision outcome. If you implemented a custom collision model, you also have to modify the mass conservation check for accurate results.
Probabilities
To decide wether a collision event results in sticking or fragmentation, the fragmentation probability is calculated with a smooth transition between sticking and fragmentation around the fragmentation velocity. This probability is stored in Simulation.dust.p.frag.
The sticking probability is simply given by
\(p_\mathrm{stick} = 1 - p_\mathrm{frag}\).
[20]:
sim.dust.p
[20]:
Group (Probabilities)
---------------------
frag : Field (Fragmentation probability)
stick : Field (Sticking probability)
-----
The sum of both is always \(1\) in the default setup. The probabilities at the boundaries are always \(0\).
[21]:
p_tot = sim.dust.p.stick + sim.dust.p.frag
(p_tot[1:-1, ...]==1).all()
[21]:
True
Bouncing
Bouncing is by default not implemented. If you want to include bouncing collisions into your simulations you have to set functions for calculating \(p_\mathrm{stick}\) and \(p_\mathrm{frag}\), such that the sum of both is smaller than \(1\). If neither sticking, nor fragmentation occurs, the mass distribution does not change, i.e., we have bouncing.
Collision Rates
The collision rates \(R\) as defined above are given by
\(\begin{split} R_\mathrm{stick} &= k \cdot p_\mathrm{stick}\\ R_\mathrm{frag} &= k \cdot p_\mathrm{frag} \end{split}\)
where \(k\) is given by
\(k = \frac{1}{1+\delta_{ij}}\frac{\pi\left( a_i + a_j \right)^2}{\sqrt{2\pi\left( H_i^2 + H_j^2 \right)}} v_\mathrm{rel}\),
which is the product of geometrical crosssection and relative velocity devided by a normalization factor that accounts for the two-dimensional structure of the disk. See (A.14) in appendix 2 of Birnstiel et al. (2010) for details on this.
\(k\) is stored for every possible particle collision at any location in the disk in Simulation.dust.kernel.
Coagulation Sources
The parameters stored in Simulation.dust.coagulation only have to be calculated once in the beginning of the simulation. But to calculate the coagulation sources the current densities and relative velocities are needed. Since these change over time, they have to be calculated on the fly.
This is done by summing over all possible collisions, evaluating the collision outcome and by multiplying with their respective collision rates. The result is stored in Simulation.dust.S.coag.
[22]:
sim.dust.S.coag
[22]:
[[ 0.00000000e+00 0.00000000e+00 0.00000000e+00 ... 0.00000000e+00
0.00000000e+00 0.00000000e+00]
[-6.23483044e-08 -6.59511337e-08 -4.67031877e-08 ... 0.00000000e+00
0.00000000e+00 0.00000000e+00]
[-4.88452102e-08 -5.16638811e-08 -3.65876621e-08 ... 0.00000000e+00
0.00000000e+00 0.00000000e+00]
...
[ 0.00000000e+00 0.00000000e+00 0.00000000e+00 ... 0.00000000e+00
0.00000000e+00 0.00000000e+00]
[ 0.00000000e+00 0.00000000e+00 0.00000000e+00 ... 0.00000000e+00
0.00000000e+00 0.00000000e+00]
[ 0.00000000e+00 0.00000000e+00 0.00000000e+00 ... 0.00000000e+00
0.00000000e+00 0.00000000e+00]]