Example: Planetary Gaps
In steady state the radial mass flux is constant
To create a custom gap profile in the gas, while still allowing for gas accretion, it is sufficient to simply force the inverse gap profile onto the viscosity
Gap Profiles
Since we are not self-consistently creating gaps, we have to impose a known gap profile. In this demonstration we are using the numerical fits to hydrodynamical simulations found by Kanagawa et al. (2017):
where
We therefore need to write a function that calculates the pertubation
[1]:
import numpy as np
def Kanagawa2017_gap_profile(r, a, q, h, alpha0):
"""Function calculates the gap profile according Kanagawa et al. (2017).
Parameters
----------
r : array
Radial grid
a : float
Semi-majo axis of planet
q : float
Planet-star mass ratio
h : float
Aspect ratio at planet location
alpha0 : float
Unperturbed alpha viscosity parameter
Returns
-------
f : array
Pertubation of surface density due to planet"""
# Unperturbed return value
ret = np.ones_like(r)
# Distance to planet (normalized)
dist = np.abs(r-a)/a
K = q**2 / (h**5 * alpha0) # K
Kp = q**2 / (h**3 * alpha0) # K prime
Kp4 = Kp**(0.25) # Fourth root of K prime
SigMin = 1. / (1 + 0.04*K) # Sigma minimum
SigGap = 4 / Kp4 * dist - 0.32 # Sigma gap
dr1 = (0.25*SigMin + 0.08) * Kp**0.25 # Delta r1
dr2 = 0.33 * Kp**0.25 # Delta r2
# Gap edges
mask1 = np.logical_and(dr1<dist, dist<dr2)
ret = np.where(mask1, SigGap, ret)
# Gap center
mask2 = dist < dr1
ret = np.where(mask2, SigMin, ret)
return ret
A planet of
[2]:
import dustpy.constants as c
rp = 10. * c.au
h = 0.05
q = 30. * c.M_earth / c.M_sun
alpha0 = 1.e-3
[3]:
r = np.geomspace(1.e0, 1.e3, 1000) * c.au
[4]:
import matplotlib.pyplot as plt
[5]:
fig = plt.figure(dpi=150)
ax = fig.add_subplot(111)
ax.semilogx(r/c.au, Kanagawa2017_gap_profile(r, rp, q, h, alpha0))
ax.axvline(rp/c.au, ls="--", c="C7", lw=1)
ax.set_xlim(1., 1000.)
ax.set_xlabel("Distance from star [AU]")
ax.set_ylabel(r"$\Sigma_\mathrm{g} \left( r \right)\ /\ \Sigma_0$")
fig.tight_layout()
plt.show()

The inverse of this pertubation needs to be imposed on the turbulent
Adding planets
For this example we want to add two planets: Jupiter and Saturn. We therefore add new groups and subgroups to organize them.
In addition to that we want to refine the grid close to the planet locations. For this we use the function similar to the one discussed in an earlier chapter.
[6]:
def refinegrid(ri, r0, num=3):
ind = np.argmin(r0 > ri) - 1
indl = ind-num
indr = ind+num+1
ril = ri[:indl]
rir = ri[indr:]
N = (2*num+1)*2
rim = np.empty(N)
for i in range(0, N, 2):
j = ind-num+int(i/2)
rim[i] = ri[j]
rim[i+1] = 0.5*(ri[j]+ri[j+1])
return np.concatenate((ril, rim, rir))
[7]:
from dustpy import Simulation
[8]:
sim = Simulation()
[9]:
ri = np.geomspace(1.e0, 1.e3, 100) * c.au
ri = refinegrid(ri, 5.2*c.au)
ri = refinegrid(ri, 9.6*c.au)
sim.grid.ri = ri
[10]:
sim.initialize()
[11]:
sim.addgroup("planets")
sim.planets.addgroup("jupiter")
sim.planets.addgroup("saturn")
[12]:
sim.planets.jupiter.addfield("a", 5.2*c.au, description="Semi-major axis [cm]")
sim.planets.jupiter.addfield("M", 1.*c.M_jup, description="Mass [g]")
[13]:
sim.planets.saturn.addfield("a", 9.6*c.au, description="Semi-major axis [cm]")
sim.planets.saturn.addfield("M", 95*c.M_earth, description="Mass [g]")
The planets are now set up.
[14]:
sim.planets.toc
Group
- jupiter: Group
- a: Field (Semi-major axis [cm])
- M: Field (Mass [g])
- saturn: Group
- a: Field (Semi-major axis [cm])
- M: Field (Mass [g])
Viscosity pertubation
We now have to write a funtion that returns the
The unperturbed
[15]:
alpha0 = 1.e-3
[16]:
from scipy.interpolate import interp1d
def alpha(sim):
"""Function returns the alpha viscosity profile to create planetary gaps in the gas surface density.
Parameters
----------
sim : Frame
Simulation frame
Returns
-------
alpha : array
alpha viscosity profile"""
# Unperturbed profile
alpha = np.ones_like(sim.gas.alpha)
for name, p in sim.planets.__dict__.items():
# Skip hidden fields
if name[0] == "_": continue
q = p.M / sim.star.M
# Aspect ratio
h = sim.gas.Hp / sim.grid.r
# Interpolate aspect ratio
f_h = interp1d(sim.grid.r, h)
hp = f_h(p.a)
# Inverse alpha profile
alpha /= Kanagawa2017_gap_profile(sim.grid.r, p.a, q, hp, alpha0)
return alpha0 * alpha
This function now produces the
[17]:
fig = plt.figure(dpi=150)
ax = fig.add_subplot(111)
ax.semilogx(sim.grid.r/c.au, alpha(sim))
ax.set_xlim(sim.grid.r[0]/c.au, sim.grid.r[-1]/c.au)
ax.set_ylim(8.e-4, 2.e-1)
ax.vlines(sim.grid.r/c.au, 1.e-4, 1.e0, lw=0.5, color="C7")
ax.set_xlabel("Distance from star [AU]")
ax.set_ylabel(r"$\alpha$ viscosity parameter")
fig.tight_layout()
plt.show()

We can now assign this profile to the
[18]:
sim.gas.alpha = alpha(sim)
Growing planets
Right now we could already start the simulation. But in this demonstration we want to let the planets grow linearly to their final mass. we therefore write a function that returns the planetary mass as a function of time.
[19]:
def Mplan(t, tini, tfin, Mini, Mfin):
"""Function returns the planetary mass.
Parameters
----------
t : float
Current time
tini : float
Time of start of growth phase
tfin : float
Time of end of growth phase
Mini : float
Initial planetary mass
Mfin : float
Final planetary mass"""
if t<=tini:
return Mini
elif t>=tfin:
return Mfin
else:
return (Mfin-Mini)/(tfin-tini)*(t-tini) + Mini
We want to start with planet masses of
[20]:
def Mjup(sim):
return Mplan(sim.t, 0., 1.e5*c.year, 1.*c.M_earth, 1.*c.M_jup)
def Msat(sim):
return Mplan(sim.t, 0., 1.e5*c.year, 1.*c.M_earth, 95.*c.M_earth)
These functions have to be assigned to the updaters of their respective fields.
[21]:
sim.planets.jupiter.M.updater.updater = Mjup
sim.planets.saturn.M.updater.updater = Msat
We also have to oragnize the update structure of the planets.
[22]:
sim.planets.jupiter.updater = ["M"]
sim.planets.saturn.updater = ["M"]
sim.planets.updater = ["jupiter", "saturn"]
And we have to tell the main simulation frame to update the planets group. This has to be done before the gas quantities, since Simulation.gas.alpha
needs to know the planetary masses. So we simply update it in the beginning.
[23]:
sim.updater = ["planets", "star", "grid", "gas", "dust"]
Since the
[24]:
sim.gas.alpha.updater.updater = alpha
To bring the simulation frame into consistency we have to update it.
[25]:
sim.update()
The planetary masses are now set to their initial values.
[26]:
msg = "Masses\nJupiter: {:4.2f} M_earth\nSaturn: {:4.2f} M_earth".format(sim.planets.jupiter.M/c.M_earth, sim.planets.saturn.M/c.M_earth)
print(msg)
Masses
Jupiter: 1.00 M_earth
Saturn: 1.00 M_earth
We can now start the simulation.
[27]:
sim.writer.datadir = "example_planetary_gaps"
[28]:
sim.run()
DustPy v1.0.6
Documentation: https://stammler.github.io/dustpy/
PyPI: https://pypi.org/project/dustpy/
GitHub: https://github.com/stammler/dustpy/
Please cite Stammler & Birnstiel (2022).
Checking for mass conservation...
- Sticking:
max. rel. error: 2.75e-14
for particle collision
m[114] = 1.93e+04 g with
m[116] = 3.73e+04 g
- Full fragmentation:
max. rel. error: 5.55e-16
for particle collision
m[55] = 7.20e-05 g with
m[55] = 7.20e-05 g
- Erosion:
max. rel. error: 1.78e-15
for particle collision
m[110] = 5.18e+03 g with
m[118] = 7.20e+04 g
Creating data directory example_planetary_gaps.
Writing file example_planetary_gaps/data0000.hdf5
Writing dump file example_planetary_gaps/frame.dmp
Writing file example_planetary_gaps/data0001.hdf5
Writing dump file example_planetary_gaps/frame.dmp
Writing file example_planetary_gaps/data0002.hdf5
Writing dump file example_planetary_gaps/frame.dmp
Writing file example_planetary_gaps/data0003.hdf5
Writing dump file example_planetary_gaps/frame.dmp
Writing file example_planetary_gaps/data0004.hdf5
Writing dump file example_planetary_gaps/frame.dmp
Writing file example_planetary_gaps/data0005.hdf5
Writing dump file example_planetary_gaps/frame.dmp
Writing file example_planetary_gaps/data0006.hdf5
Writing dump file example_planetary_gaps/frame.dmp
Writing file example_planetary_gaps/data0007.hdf5
Writing dump file example_planetary_gaps/frame.dmp
Writing file example_planetary_gaps/data0008.hdf5
Writing dump file example_planetary_gaps/frame.dmp
Writing file example_planetary_gaps/data0009.hdf5
Writing dump file example_planetary_gaps/frame.dmp
Writing file example_planetary_gaps/data0010.hdf5
Writing dump file example_planetary_gaps/frame.dmp
Writing file example_planetary_gaps/data0011.hdf5
Writing dump file example_planetary_gaps/frame.dmp
Writing file example_planetary_gaps/data0012.hdf5
Writing dump file example_planetary_gaps/frame.dmp
Writing file example_planetary_gaps/data0013.hdf5
Writing dump file example_planetary_gaps/frame.dmp
Writing file example_planetary_gaps/data0014.hdf5
Writing dump file example_planetary_gaps/frame.dmp
Writing file example_planetary_gaps/data0015.hdf5
Writing dump file example_planetary_gaps/frame.dmp
Writing file example_planetary_gaps/data0016.hdf5
Writing dump file example_planetary_gaps/frame.dmp
Writing file example_planetary_gaps/data0017.hdf5
Writing dump file example_planetary_gaps/frame.dmp
Writing file example_planetary_gaps/data0018.hdf5
Writing dump file example_planetary_gaps/frame.dmp
Writing file example_planetary_gaps/data0019.hdf5
Writing dump file example_planetary_gaps/frame.dmp
Writing file example_planetary_gaps/data0020.hdf5
Writing dump file example_planetary_gaps/frame.dmp
Writing file example_planetary_gaps/data0021.hdf5
Writing dump file example_planetary_gaps/frame.dmp
Execution time: 1:04:02
This is the result of our simulation.
[29]:
from dustpy import plot
[30]:
plot.panel(sim)

As can be seen, the core of Jupiter and Saturn carved gaps in the gas disk, which in turn affected the dust evolution.
To see that the planetary masses actually changed we can plot them over time.
[31]:
t = sim.writer.read.sequence("t")
Mjup = sim.writer.read.sequence("planets.jupiter.M")
Msat = sim.writer.read.sequence("planets.saturn.M")
[32]:
fig = plt.figure(dpi=150)
ax = fig.add_subplot(111)
ax.plot(t/c.year, Mjup/c.M_earth, label="Jupiter")
ax.plot(t/c.year, Msat/c.M_earth, label="Saturn")
ax.legend()
ax.set_xlim(t[0]/c.year, t[-1]/c.year)
ax.grid(visible=True)
ax.set_xlabel(r"Time [years]")
ax.set_ylabel(r"Mass [$M_\oplus$]")
fig.tight_layout()
plt.show()
