Simulations of Dynamical Gas-Dust Circumstellar Disks: Going Beyond the Epstein Regime
Olga P. Stoyanovskaya, Fedor A. Okladnikov, Eduard I. Vorobyov, Yaroslav N. Pavlyuchenkov, Vitaliy V. Akimkin
SSimulations of Dynamical Gas–Dust Circumstellar Disks:Going Beyond the Epstein Regime
O. P. Stoyanovskaya, F. A. Okladnikov, E. I. Vorobyov,Ya. N. Pavlyuchenkov, V. V. AkimkinLavrentyev Institute of Hydrodynamics,Siberian Branch of the Russian Academy of Sciences,Novosibirsk, Russia,Institute of Computational Technology,Siberian Branch of the Russian Academy of Sciences,Novosibirsk, Russia,Novosibirks State University, Novosibirsk, RussiaSouthern Federal University, Institute of Physics, Rostov-on-Don, Russia,Department of Astrophysics, University of Vienna, Vienna, AustriaInstitute of Astronomy, Russian Academy of Sciences,Moscow, RussiaFebruary 19, 2021
Abstract
In circumstellar disks, the size of dust particles varies from submicron to several centimeters, whileplanetesimals have sizes of hundreds of kilometers. Therefore, various regimes for the aerodynamic dragbetween solid bodies and gas can be realized in these disks, depending on the grain sizes and velocities:Epstein, Stokes, and Newton, as well as transitional regimes between them. This means that simulations ofthe dynamics of gas–dust disks require the use of a drag coefficient that is applicable for a wide range for sizesand velocities for the bodies. Furthermore, the need to compute the dynamics of bodies of different sizes inthe same way imposes high demands on the numerical method used to find the solution. For example, inthe Epstein and Stokes regimes, the force of friction depends linearly on the relative velocity between thegas and bodies, while this dependence is non-linear in the transitional and Newton regimes. On the otherhand, for small bodies moving in the Epstein regime, the time required to establish the constant relativevelocity between the gas and bodies can be much less than the dynamical time scale for the problem—thetime for the rotation of the disk about the central body. In addition, the dust may be concentrated inindividual regions of the disk, making it necessary to take into account the transfer of momentum betweenthe dust and gas. It is shown that, for a system of equations for gas and monodisperse dust, a semi-implicitfirst-order approximation scheme in time in which the interphase interaction is calculated implicitly, whileother forces, such as the pressure gradient and gravity are calculated explicitly, is suitable for stiff problemswith intense interphase interactions and for computations of the drag in non-linear regimes. The piece-wisedrag coefficient widely used in astrophysical simulations has a discontinuity at some values of the Mach andKnudsen numbers that are realized in a circumstellar disk. A continuous drag coefficient is presented, whichcorresponds to experimental dependences obtained for various drag regimes.
Simulations of the dynamics of circumstellar gas–dust disks are necessary to improve our understanding ofmechanisms for the formation of planets, and have been carried out in many studies, such as [5, 49, 10, 33, 6,13, 43, 35, 14]. Thus far, models for the dynamics of gas–dust disks have been developed in which the gas istaken to be a carrier phase, and solid particles to be a dispersed phase. In this case, the dynamics of a dustcloud can be described like the dynamics of a continuous medium (the basis for applying this approximation ispresented in Appendix A). If we suppose that the solid phase is represented by particles of a single size at each1 a r X i v : . [ a s t r o - ph . E P ] F e b oint in space, the continuity equation and equation of motion of the gas and dust in the disk will have theform ∂ρ g ∂t + ∇ ( ρ g v ) = S g , ρ g (cid:20) ∂v∂t + ( v · ∇ ) v (cid:21) = −∇ p + ρ g g − f D + f g , (1) ∂ρ d ∂t + ∇ ( ρ d u ) = S d , ρ d (cid:20) ∂u∂t + ( u · ∇ ) u (cid:21) = ρ d g + f D + f d , (2)where ρ g and ρ d are the volume densities of the gas and dust, v and u the velocities of the gas and dust, p the gas pressure, g the acceleration due to the self-gravitation of the gas and dust, as well as the gravity ofthe central star, S g , S d the sources and sinks for the gas and dust, f g , f d the forces acting on the gas and dust,apart from the pressure force, gravitation, and drag, f D the drag force per unit volume of the medium, f D = n d F D , (3)where n d is the volume number density of dust particles, F D the force of friction per particle, F D = 12 C D sρ g (cid:107) v − u (cid:107) ( v − u ) , (4)where s is the area of the “frontal” cross section of a particle (for example, for a sphere of radius ), and C D isthe dimensionless drag coefficient, which is a function of two dimensionless quantities – the Mach number,Ma = (cid:107) v − u (cid:107) c s (5)and the Knudsen number, Kn = λa , (6)Here, λ is the mean free path of a gas molecule c s is a sound speed in pure gas.The dust and solid bodies in the circumstellar disk comprise objects with various sizes, from submicrondust particles to planetesimals up to hundreds of kilometers. Large and small bodies interact differently withthe gas. Small dust particles intensively exchange momentum with the gas, and their steady-state or terminalvelocities relative to the gas are small. Large bodies can have velocities that are appreciably different from thegas velocity. This means that Ma, Kn and || f D || in the disk vary by orders of magnitude. Consequently, fornumerical simulations of the dynamics of the two-phase medium in the disk, it is necessary to use a coefficient offriction that is applicable over a broad range of values of Ma and Kn. Therefore, 2 presents an analysis of variousforms of drag coefficients and shows important differences between them and their advantages when applied inmodels of a disk. Furthermore, || f D || in Eqs. (1) and (2) appreciably exceeds the sum of the other forces forsmall bodies, but is comparable to the other forces for large bodies. In general, f D depends non-linearly on therelative velocity between the gas and bodies. These factors impose high demands on the numerical method usedto solve the equations of motion (1) and (2). Section 3.1 presents a review of this problematic situation andpublished results of studies on methods for the simulation of the dynamics of a gas–dust medium with intenseinterphase interactions. Section 3.2 presents a semi-implicit scheme that is first order in time, and Section 3.3shows that this scheme satisfies the necessary requirements for its application in simulations of circumstellardisks. The conclusions are presented in Section 4. The applicability of the model to a continuous–mediumdescription of the dynamics of the dust in the disk is justified in Appendix A, where alternative approaches arealso presented. Appendix B describes the simple, one-dimensional model for a circumstellar disk used for oursimulations. In the case of a compressible gas, the drag coefficient C D in (4) is a function of the two independent variablesMa and Kn. The characteristic form of C D is presented in Fig. 1. In a regime where the gas flows aroundthe bodies like a continuous medium, the drag coefficient is determined by the Reynolds number Re (see, forexample, [4]), which is derived from Ma and Kn: Re = 4 MaKn . (7)2igure 1: Drag coefficient as a function of the Reynolds and Mach numbers for the free-molecular flow regime(Epstein regime) and the case when the matter flows around a body like a continuous medium (the Stokes,transition, and Newton regimes). The boundary between free-molecular flow and continuous-medium flow isdetermined by the Knudsen number, and the boundaries of the continuous-medium flow regime by the Reynoldsnumber. Stream lines around spherical particles are shown for the continuous-medium flow regimes. Thefunctional dependence of the drag force F D on the particle radius a and relative velocity between the particlesand gas ∆ V = (cid:107) v − u (cid:107) are presented for all regimes. 3igure 2: Regions of determination of the piecewise function C D (8), corresponding to the various flow regimes.I—Epstein, II—Stokes, III—transition (non-linear Stokes), IV—Newton. The green arrows indicate the bound-aries at which a function is continuous, and the red x’s boundaries where a function is discontinuous.Particles with radii less than 2 . λ , interact with the gas in a free-molecular-flow regime, or Epstein regime [15].The drag coefficients for such bodies do not depend on Re, but do depend on Ma. This dependence correspondsto the standard astrophysical coefficient proposed in [46]: C D =
83 Ma − , Kn − <
94 (Epstein regime; I),24Re − , Re < − . , < Re <
800 (transition regime or non-linear Stokes regime; III),0 . , Re >
800 (Newton regime; IV). (8)An approximation for the Stokes regime, the transition regime (also known as the non-linear Stokes regime),and the Newton regime in (8) was proposed in [32], and was applied in [47] to the aerodynamics of solid bodiesin a disk. This approximation was expanded to encompass the Epstein regime in [46]. The coefficient (8) isused in [33, 23] and other studies in the framework of modern models for circumstellar disks.Figure 2 presents the regions I–IV within which the drag coefficient C D (as a function of Kn − and Ma) iscontinuously differentiable. At the boundaries I–II, II–III, and III–IV, C D is continuous, whereas this coefficientis discontinuous at the boundaries I–III and I–IV. The presence of a discontinuity is a serious inadequacyfor simulations, since this can lead to the appearance of artificial singularities in the solution. We convincedourselves that the conditions for which a particle makes the transition from the Epstein regime I to the transitionregime (non-linear Stokes regime) III, which is accompanied by a discontinuity of the coefficient of friction (8).is indeed realized in simulations of the dynamics of a disk. This conclusion was based on the simulations of thedynamics of a viscous, self-gravitating, gaseous disk with a dust component using the FEOSAD model, which isdescribed in detail in [43]. The disk model used differs from the basic model presented in [43]in the additionalconsideration of the influence of the dust dynamics on the gas dynamics, and also in the absence of the restrictionon the size of the dust a < . λ , which artificially maintains the particles in the Epstein regime. We used dragcoefficient (8) and the numerical scheme from Section 3.2 for our simulations. We used a constant value for theShakura–Syunyaev viscosity over the disk, α = 10 − , and the fragmentation velocity v frag = 30 m s − . ThisShakura–Syunyaev viscosity corresponds to a disk with suppressed magnetorotational instability, in which thetransport of mass and angular momentum occurs mainly via gravitational torques in a non-axially symmetric4igure 3: Values of Ma and Re (left column) and of Ma and Kn − (central column) for dust grains in all cellsof circumstellar disks with ages of 141 000 years (upper row) and 707000 years (lower row). The various dragregimes from (8) are shown in the central panels by Roman numerals and in the right panels by different colors.I—Epstein (white), II—Stokes (blue), III—transition (non-linear Stokes) (green), IV—Newton (red). Blackindicates drag cells of the disk. The standard astrophysical drag coefficient (8) is discontinuous at the boundarybetween the white and green region.disk. A fragmentation velocity v frag determines the maximum particle size. The fragmentation velocity chosencorresponds to an enhanced probability of adhesion, as is realized when the particles are covered in a thick layerof ice.We computed the quantities Kn, Ma and Re in each computational cell for the evolution of the disk overone million years. Figure 3 presents the set of values of these quantities in the disk in the Ma, Re (left panels)and Kn − , Ma (central panels) axis for ages of 141000 and 707000 years. The central panels show that thedust and gas interact in the free-molecular flow Epstein regime in the most part of the disk, however all fourof the drag regimes indicated in (8) are realized in the disk. At a time of 707000 years (central panel in thesecond row), particles are present in the disk at the boundary of regions I–III, where the drag coefficient (8)undergoes a discontinuity. The right lower panel in Fig.3 where different colours indicate different regions in thedisk with different drag regimes, shows that this discontinuity occurs at a distance of 6 au from the star. Thedistribution of physical quantities determining the drag regime in the central region of the circumstellar diskat the same times are shown in Fig. 4. The black contour in the upper panels, which present the Knudsen andMach numbers, shows the boundary where a = 4 λ/
9. In a disk with an age of 141000 years, this boundary passesthrough the region Ma ≈ .
01, corresponding to a continuous coefficient of friction, while it passes through theregion Ma ≈ / < × and Ma < ≤
1) 5igure 4: Distributions of physical quantities determining the drag regimes in the central regions of circumstellardisks with ages of 141000 years (left six panels) and 707000 years (right six panels). The upper panels presentthe Knudsen and Mach numbers and the lower panels present the surface densities of gas and grown dust in gcm − , The central panels show the maximum grain size a in cm and the relative velocity between the gas anddust ∆ V in km s − . The black contours in the upper panels show the boundary at which a = 4 λ/
9. For a diskwith an age of 141000 years, this passes through the region Ma ≈ .
01, providing a continuous drag coefficient,while it passes through the region Ma ≈ /
9, leading to a discontinuous drag coefficient for a disk with an ageof 707000 years. 6igure 5: Henderson drag coefficient (10)-(13) and the standard astrophysical expression (8) for various Ma asa function of Re. C D = 24Re + S (cid:18) .
33 + 3 . − . T d /T g . T d /T g exp( − .
247 Re S ) (cid:19) + 0 . S (cid:18) − exp (cid:18) − MaRe (cid:19)(cid:19) (9)+ exp (cid:18) − . √ Re (cid:19) . . (cid:16) . . √ Re (cid:17) . . √ Re + 0 . + 0 . (10)while for supersonic regimes (Ma ≥ . C D = 0 . . + 1 . (cid:114) MaRe (cid:32) S + 1 . S (cid:115) T d T g − S (cid:33) . (cid:114) MaRe , (11)where S = Ma (cid:114) γ , (12) T d /T g is the ratio of the dust and gas temperatures. Intermediate values for the drag coefficient (when 1 < Ma < .
75) can be obtained using a linear interpolation: C D = C D (Ma = 1) + 43 (Ma −
1) ( C D (Ma = 1 . − C D (Ma = 1)) , (13)where C D (Ma = 1) can be found from (10), by substituting Ma = 1, and C D (Ma = 1 .
75) can be found from(11), by substituting Ma = 1 .
75. Here and below, we assume γ = 1 . T d /T g = 1.The values of the standard astrophysical coefficient (8) and the Henderson coefficient (10)–(13) for the rangeof parameters of a circumstellar disk 10 − ≤ Ma ≤ − ≤ Re ≤ are presented in Figs.5 and 6. It7igure 6: Henderson drag coefficient (10)–(13) and the standard astrophysical expression (8) for various Re asa function of Ma.follows from the upper row of panels in Fig. 5 that both coefficients are similar when Ma ≤ .
01 for any value ofRe. This is confirmed by the upper panels of Fig. 6. The values of the two drag coefficients differ quantitativelywhen 0 . ≤ Ma ≤ ≥
1. First, a discontinuity of the first type is visible in thestandard astrophysical coefficient in the central and right panels of the lower row of Fig. 5. Furthermore, thelower right panel shows that the Henderson coefficient falls with growth in Re, while the standard astrophysicalcoefficient grows. Second, the central and right panels of the lower row in Fig. 6 show that the Hendersoncoefficient has a feature (a bump near Ma = 1, well known in aeromechanics), which is absent in the simplerastrophysical approximation for the coefficient.Thus, when simulating the dynamics of a gas–dust medium with Ma > /
9, it is recommended to use onlythe Henderson coefficient. However, if the treatment is restricted to Ma ≤ /
9, the standard astrophysicalcoefficient is less computationally expensive and convenient for programming.
We define the time scale for the velocity of a spherically symmetrical particle to relax to the velocity of theambient gas as follows: t stop = m p (cid:107) v − u (cid:107)(cid:107) F D (cid:107) = πa ρ s (cid:107) v − u (cid:107) C D πa ρ g (cid:107) v − u (cid:107) = 83 aρ s C D ρ g (cid:107) v − u (cid:107) , (14)where m p is the mass of the particle and ρ s is the material density of the grains.For dust grains with sizes of about one micron in an orbit with a radius of 100 au, the estimated relaxation,or stopping, time t stop in the disk is of order 100 s (see, for example, [40, 22]), while the disk dynamics require8imulations over 10 yrs ≈ s or more. This means that the equations of motion of the gas and dust havestiff relaxation terms associated with intense interphase interactions. This makes simulations the dynamics ofgas–dust disks a challenging task in modern computational astrophysics (see, for example, [16]).Problems with stiff relaxation terms have been actively studied in applied mathematics [11]. In additionto the dynamics of multiphase media, such problems also arise in plasma physics, models for elasticity withmemory, problems in kinetic theory, and so forth. When explicit methods are used to obtain stable numericalsolutions, the time step must satisfy the condition τ ≤ t stop , which is unacceptably expensive for simulationsof multi-physics systems of equations in the two-dimensional and three-dimensional cases, even on modernsupercomputers.In particular, when modeling the dynamics of a gaseous, self-gravitating disk (e.g., [42, 39]), the time stepindicated by the Courant condition is several Earth days. When necessary, this step can be made up to a factorof 100 shorter, but the time for the solutions to be obtained then increases by a factor of 100–1000, making itinfeasible to carry out series of computational experiments.There is no restriction on the time step from the stability point of view when using implicit methods, butrestrictions can arise due to the required accuracy [31]. For example, in an early study, Jim and Livermore[20] showed that, in general, in the presence of stiff relaxation terms, even high-accuracy grid methods canreproduce asymptotic behavior of the solutions only if a very small spatial step is used. Examples of themanifestation of such effects in simulations of the dynamics of circumstellar disks can be found in [40, 43, 23].In particular, the diffusion overdissipation of solutions obtained in simulations of the dynamics of a gas–dustmedium using smoothed particle hydrodynamics (SPH) with computations of the interphase interactions usinga classical approach [29] were studied in [24, 23]. An empirical criterion was derived, according to which theoverdissipation of the solutions is at an acceptable level if the smoothing radius h < c s t stop is used. Accordingto this criterion, simulations of the dynamics of micron-sized particles in a circumstellar disk requires a spatialresolution h ≈
10 km, whereas the radius of the disk is of order 10 km. It is clear that such a spatial resolutionis beyond the capabilities of modern computational technology.In a number of cases, a search for a numerical solution of a stiff problem can be substituted with the solutionto a simpler problem obtained from the original one via expansion in a small parameter (this general idea ispresented in [11], and its application to simulations of a gas–dust medium in a circumstellar disk in [21, 1, 2],for example). The simplified problem becomes non-stiff, but this approach is justified only if the relaxationtime is much smaller than the time scale for the motion of the carrier phase, and it cannot be used to modelthe transition regime when the relaxation time is close to the dynamical time scale. Moreover, this approachmay involve a change in the type of equations considered.Therefore, an alternative to adopting a simplified problem is the development of implicit and semiimplicitnumerical methods that are able to precisely reproduce the equilibrium values, even with crude spatial resolution.Attempts are currently being made to find a general approach to the construction of highaccuracy methodsfor problems with relaxation terms (see, e.g., [3]). However, in most studies of the numerical solution of suchsystems, the methods are constructed taking into account the specifics of the problem to be solved and specificforms for the “stiff relaxation term”. In particular, when simulating the dynamics of a gas–dust medium, thesystem has not one, but two stiff relaxation terms—in the equations of motion for the gas and dust [27, 24, 19,48]). On average over the disk, ε = ρ d ρ g is equal to 0.01, to order of magnitude, so that the influence of the duston the gas dynamics can often be neglected. On the other hand, the results of simulations indicate that thedust particles may be concentrated in specific regions of the disk (e.g., in spiral arms [33], in the inner part ofthe disk [43], in self-gravitating gaseous clumps [10]), increasing the density ratios in these locations to values ε ≈ To approximate the interphase interactions in the equations of motion in (1)-(2), we linearized the drag termon the relative velocity between the gas and particles, with the stopping time t stop ( ρ g , u − v ) computed usingknown quantities from the previous step and the relative velocity of the gas and particles computed from thesubsequent step. A simple description of this approach has the form: v n +1 − v n τ + ( v n · ∇ ) v n = − ∇ p n ρ n g + g n − ρ n d ρ n g v n +1 − u n +1 t n stop + f n g ρ n g , (15) u n +1 − u n τ + ( u n · ∇ ) u n = g n + v n +1 − u n +1 t n stop + f n d ρ n d . (16)We solved the entire system (1)-(2) using a two-stage scheme based on the operator splitting method withrespect to physical processes (as in our earlier studies [40, 37]). We solved the continuity equation and theequation of motion using the finite-difference method described in detail for a pure gas in [36]. In the first stageof the scheme, the advective terms and transport of mass and momentum are computed: ∂ρ g ∂t + ∇ ( ρ g v ) = 0 , ρ g (cid:20) ∂v∂t + ( v · ∇ ) v (cid:21) = 0 , (17) ∂ρ d ∂t + ∇ ( ρ d u ) = 0 , ρ d (cid:20) ∂u∂t + ( u · ∇ ) u (cid:21) = 0 . (18)We computed the spatial derivatives in each stage using the piecewise-parabolic method of [12]. In this stage,we determined the quantities ρ n +1 / , ρ n +1 / , v n +1 / , u n +1 / .In the next stage, we computed the influence of the forces of pressure, drag, gravitation, and other effectsusing the updated densities and velocities of the gas from the first stage: v n +1 − v n +1 / τ = − ∇ p n +1 / ρ n +1 / + g n +1 / − ρ n +1 / ρ n +1 / v n +1 − u n +1 t n +1 / + f n +1 / ρ n +1 / , (19)10 n +1 − u n +1 / τ = g n +1 / + v n +1 − u n +1 t n +1 / + f n +1 / ρ n +1 / . (20)Thanks to the fact that v n +1 and u n +1 can be expressed explicitly using (15), (16), the computational costsrequired for the semi-implicit approximation of friction are similar to those for an explicit approximation. In Section (3.3.1), we will consider the degree of numerical overdissipation of the scheme (15), (16) during thecomputation of intense interphase interaction using a time step τ , determined by the Courant condition, andnot t stop . In Section 3.3.2 we will compare the accuracy of solutions obtained using the firstorder scheme (15),(16) and a fourth-order scheme to simulate the motion of single dust grain in a gaseous disk whose angularvelocity is slightly sub-Keplerian. In the one-dimensional equations (1) and (2), we set g = 0, S g = S d = 0, f g = f d = 0, f D = ρ d ( v − u ) t stop andwrote the energy equation: ρ g (cid:18) ∂e∂t + v ∂e∂x (cid:19) = − p ∂v∂x , (21)where e is the internal energy of the gas, which is related to the pressure as p = ρ g e ( γ − , (22)where γ is the adiabatic index. We made this system dimensionless using the length l ∗ , mass m ∗ , and time t ∗ units, so that the derivative dimensional quantities have the form ρ ∗ = m ∗ l ∗ , v ∗ = l ∗ t ∗ , p ∗ = m ∗ t ∗ l ∗ , e ∗ = l ∗ t ∗ . Wesolved the system in the dimensionless variables.For the resulting system (1), (2), (21), (22), we specified a no-flow condition at the boundaries of an interval.We specified the velocities to be zero at the initial time, and specified the jumps in the gas pressure and in thedensities of the gas and dust to be (cid:20) ρ g ρ ∗ , pp ∗ , ee ∗ , ρ d ρ ∗ (cid:21) l = [1 , , . , , (cid:20) ρ g ρ ∗ , pp ∗ , ee ∗ , ρ d ρ ∗ (cid:21) r = [0 . , . , , . ,γ l = γ r = 7 / . We specified λ = 5 · − l ∗ ρ ∗ ρ g , ρ s = 2 . ρ ∗ and considered two particle sizes: small grains with a = 5 · − l ∗ and large grains with a = 10 − l ∗ .It was shown in [40] that the semi-implicit, first-order operator splitting scheme (15), (16) preserves theasymptotic solutions in computations of intense interphase interactions with a constant drag coefficient. Thissame conclusion follows from Fig.. 7 for the case when C D varies in space and time and t stop = t stop ( (cid:107) u − v (cid:107) ).Figure 7 presents simulations with the Henderson drag coefficient (10)–(13) for a grain size a = 5 · − l ∗ .We adopted the results for our simulations with detailed spatial resolution with N = 400 cells per segment anda time step for which the Courant parameter is CFL = 0 .
005 as a reference. It is clear that the computationsof the gas and dust velocities and densities with N = 200 , CFL = 0 . ρ d /ρ g = 1, we also carried out simulations for the same problemwith the dust content relative to the gas content enhanced by a factor of 10 and 1000. We confirmed that thelevel of numerical overdissipation of the solution remained insignificant for these cases in computations usingthe asymptotic-preserving scheme.An example of a semi-implicit drag term approximation that does not preserve the asymptotic behavior is v n +1 − v n τ + ( v n · ∇ ) v n = − ∇ p n ρ n g + g n − ρ n d ρ n g v n +1 − u n t n stop + f n g ρ n g , (23)11igure 7: Solution for a Sod shock tube for the scheme (15), (16), which preserves the asymptotic behavior, forKn ≥ a = 5 · − l ∗ , the Epstein regime, Henderson drag coefficient (10)–(13)).Figure 8: Solution for a Sod shock tube for the scheme (23), (24), which does not preserve the asymptoticbehavior, with various time steps for Kn ≥ a = 5 · − l ∗ , Epstein regime, Henderson dragcoefficient (10)–(13)). 12igure 9: Solution for a Sod shock tube for the scheme (15), (16), which preserves the asymptotic behavior,with various steps in space and time, for Kn (cid:28) u n +1 − u n τ + ( u n · ∇ ) u n = g n + v n +1 − u n +1 t n stop + f n d ρ n d . (24)The solution of this problem obtained using the scheme (23), (24) is presented in Fig. 8. This figure showsthat, when CFL is varied from 0.5 to 0.00005, the computed wave-propagation speeds vary significantly, withthe accuracy of the resulting solution being satisfactory only when τ < . t stop .It follows from the upper panel of Fig. 7 that the Knudsen number for the small grains exceeds unity.This means that, according to the regime classification (8), the interaction of the gas and dust occurs in theEpstein free-molecular-flow regime. The computation of t n stop was carried out using the general formula (14)for linear and non-linear drag term. Note that the influence of the numerical resolution on the results is visibleonly in the lower right panel of Fig. 7, which presents the quantity Re, which depends linearly on (cid:107) u − v (cid:107) . Itis known, however (see, e.g., [22]), that the relative velocity of the gas and dust can be ignored in the case t stop max( u, v ) (cid:28) l ∗ , and the solution for a gas–dust medium can be obtained from the solution for a gaseousmedium by replacing c s by c ∗ s = c s (cid:18) ρ d ρ g (cid:19) − / . (25)This means that the solution of the stationary problem does not depend on the drag coefficient and t stop , and isdetermined by the ratio of the dust and gas densities. Therefore, the appreciable deviations of Re for differentspatial resolutions presented in the lower right panel of Fig. 7, do not lead to substantial differences in thequantities depicted in the other panels. This conclusion is supported by the right panels of Fig. 10, which showthe velocities of the gas and dust obtained for small grains with N = 200 and CFL = 0 . a = 10 − l ∗ . It follows from the right panels of Fig. 9 that, for the specifiedgrain size, the gas flows around the particles like a continuous medium (Kn (cid:28)
1) with the Stokes regime,transition regime, and Newton regime being realized. The simulation results for N = 200 , CFL = 0 . The equation of motion of a particle in the field of a massive central body in polar coordinates is presented, forexample, in [26]. If gas drag acts on the particle together with gravitational and the centrifugal force, accordingto (4), these equations acquire the form: drdt = v r ,dv r dt = v ϕ r − GMr − C D ρ g ( u r − v r ) || u − v || aρ s ,dv ϕ dt = − v r v ϕ r − C D ρ g ( u ϕ − v ϕ ) || u − v || aρ s , (26)where r is the orbital radius of the particle, ( v r , v ϕ ) the radial and azimuthal velocities of the particle, ( u r , u ϕ )the velocity of the gas, M the mass of the central body, and G the gravitational constant. We took thedistributions of the density, temperature, and velocity in the axially symmetric gaseous disk to be given by themodel [33] (a detailed description of the model is presented in the Appendix B), and to display the followingdependences on the radius: ρ g ∼ r − . , c s ∼ r − . , u r = 0. Let us supposed that the central body has a massof 0 . M (cid:12) and the disk extends from 1 au to 100 au and has a mass of 0 . M (cid:12) . Let the material density ofsolid bodies be ρ s = 2 . − .At the initial time, we specified the particle to have an orbital radius r = 10 au and a velocity v ϕ =0 . V K ( r ), where V K ( r ) = (cid:114) GMr , v r = − − au · Ω( r = 1 au), where Ω( r ) = V K ( r ) r . We assumed that thesize of particles drifting in the disk grows as a = a (cid:18) rr (cid:19) − , with the particles a distance r from the centralbody having sizes a = 134 .
64 cm (Toy Model 1) or a = 13 .
464 cm (Toy Model 2). We simulated the motionof a particle toward the central body until its orbital radius reached 1 au. We used the semi-implicit first-orderscheme r n +1 − r n τ = v nr ,v n +1 r − v nr τ = ( v nϕ ) r n − GM ( r n ) − C n D ρ n g ( u n +1 r − v n +1 r ) || u n − v n || a n ρ s ,v n +1 ϕ − v nϕ τ = − v nr v nϕ r n − C n D ρ n g ( u n +1 ϕ − v n +1 ϕ ) || u n − v n || a n ρ s , (27)14able 1: Errors in the values of variables in the system (27), found using two different methods: an explicitfourth-order Runge–Kutta method and the first-order semi-implicit method, relative to the reference solutionModel—Drag coefficient Weidenschilling (1977) Henderson (1976)Toy Model 1 ∆ r/r = 1 . × − ∆ r/r = 7 . × − ∆ v r /v r = 1 . × − ∆ v r /v r = 8 . × − ∆ v ϕ /v ϕ = 4 . × − ∆ v ϕ /v ϕ = 3 . × − Toy Model 2 ∆ r/r = 4 . × − ∆ r/r = 1 . × − ∆ v r /v r = 9 . × − ∆ v r /v r = 4 . × − ∆ v ϕ /v ϕ = 1 . × − ∆ v ϕ /v ϕ = 1 . × − with a time step of τ = 0 . − ( r = 1 au). This value of τ means that a particle that moves with Keplerianvelocity in a circular orbit of radius 1 au makes a full orbit in 628 steps. This time step will be obtained from theCourant condition when the dynamics of a gaseous disk is modeled in cylindrical coordinates with a resolutionof 256 cells in azimuth and with Courant number CFL = 0 . C D , Stokes number St = t stop Ω, and radial component of the particle velocity v r , relativeto the Keplerian velocity via radius are presented in Fig. 11. for Toy Models 1 and 2. The red curve showsthe results obtained with the Henderson coefficient and the multicolored curve the results obtained with thestandard astrophysical drag coefficient. The upper panels of Fig. 11 show that, in Toy Model 1, the standardastrophysical drag coefficient undergoes a discontinuity of the first type in the first ten years. The discontinuitypoint coincides with an extremum of the radial velocity. The radial velocity also changes from growing todecreasing in the first ten years of motion of the particle in Toy Model 1 using the continuous Hendersoncoefficient. The discontinuity points in Toy Models 1 and 2, which are extrema or inflection points of C D , arealso extrema or inflection points of v r .On the whole, apart from the discontinuity in C D the simulations for Toy Model 1 with the different dragcoefficients are qualitatively and quantitatively similar: the times over which the particle drifts from 10 au to 1au are of order 3 × yrs, with the Stokes number reaching St = 10 . at a radius of 1 au. We can see that, overthe time of its motion in the disk, the particle interacts with the gas in the Epstein I, transition III, and NewtonIV regimes in (8), bypassing the Stokes II regime. The lower panels of Fig. 11 show that, with decreasing grainsize in Toy Model 2, the drag coefficient does not change at early times [by virtue of (8), it is determined by Main the Epstein regime], however the Stokes number changes by an order of magnitude compared to Toy Model1. As a result, the time for the particle to drift from 10 au to 1 au is shortened, with the particle successivelypassing through regimes I, II, III, and IV during its interaction with the gas. The computational results for ToyModel 2 with the standard astrophysical drag coefficient and with the Henderson coefficient are qualitativelyand quantitatively similar.As in Toy Model 1, the maximum differences in the radial velocity when using the different drag coefficientsarise in the Epstein regime. This differs from the shock tube problem, where the sensitivity of the solution tochanges in the drag coefficient were more clearly expressed for large bodies in the Stokes, transition, and Newtonregimes. The sensitivity of the system (27) to a change in the coefficient of friction in the Epstein regime isdetermined by the fact that the stationary value of the radial velocity of the particle depends on t stop [30] as v r V K = − η St( r ) + St( r ) − , (28)where η is determined from u ϕ = (1 − η ) V and satisfies 0 < η (cid:28)
1. The upper panels of Fig. 5 show that, in theEpstein regime, the Henderson coefficient differs from the standard astrophysical coefficient by approximatelya factor of three, which, by virtue of (28), leads to a substantial difference in the radial velocity of the particle.Two quantities are present in the second equation of (27) whose difference is many orders of magnitude lessthan their magnitudes ( v ϕ /r and GM /r ). Therefore, in computations of the trajectories of particles whosedynamics are determined by the comparable actions of the centrifugal and inward gravitational forces, high-order methods require less computational resources than first-order methods (here, we are referring to methodsfor the solution of an arbitrary system of ordinary differential equations, and not specialized simplex methodsadapted for specific problems). We investigated how the requirements for the numerical method change incomputations of particle trajectories influenced by additional forces.For this, we computed Toy Model 1 and Toy Model 2 using an explicit fourth-order Runge–Kutta methodwith the same step, τ = 0 . − ( r = 1 au). We determined the values of r , v r and v ϕ at the time when15igure 11: Simulations of the trajectory of a growing dust particle in a stationary gaseous disk with mass0.4 M (cid:12) around a central body with mass 0.5 M (cid:12) with various drag coefficients. The red curve shows theHenderson drag coefficient (10)–(13), he multicolored curve the standard astrophysical coefficient (8) with thevarious regimes I–IV indicated. 16igure 12: The left panel shows the relative error in the solution of the system without friction after 40 orbitalperiods [the integrals (29), (30)], obtained using the first- and fourth-order methods in the time step, the centraland right panels the errors in v r in the solution of the system with drag using the first-order semi-implicit methodfor various Stokes numbers. The integration time was two orbital periods for the central panel and 40 orbitalperiods for the right panel.the particle crossed r = 1 au, and further compared these values with those computed using the semi-implicitfirst-order scheme. The results of this comparison are presented in Table 1. In both models and with both dragcoefficients, the accuracy of the computations for the first-order scheme differs from that for the fourth-orderscheme by less than 1%, even though the Stokes number for the large particles exceeds 10 in both models.Among the three variables of the system, the maximum relative errors are displayed by v r . The computations forToy Model 2 using the different methods and drag coefficients differ only insignificantly. However, in Toy Model1, where the standard astrophysical drag coefficient has a discontinuity, the use of the Henderson coefficientmakes it possible to obtain more similar solutions using the first- and fourth-order methods.In order to investigate the differences in the accuracy of the solutions obtained using the first- and fourth-order methods in more detail, we first compared the solutions of the system (27) without drag force. Thissystem of three equations has two first integrals, which express the conservation of angular momentum andenergy: C = rv ϕ , (29) C = v r v ϕ − GMr . (30)At the initial time, we placed the particle at the radius r = 10 au and specified its radial velocity to beclose to zero and its azimuthal velocity to be the Keplerian value. We computed the particle’s trajectory overa time t = 40 × π Ω − ( r ) (40 orbits around the central body) using an explicit Euler method (into which thesemi-implicit, first-order scheme with zero drag is degenerate) and the explicit Runge–Kutta method. At thefinal time, we calculated the values of C and C and compared them to the initial values of these integrals.The dependence of the relative errors in C and C on the time step τ re presented in the left panel of Fig. 12.Both schemes give the stated order of approximation in the solution, with the accuracy of the computationof the angular momentum in the first-order scheme being almost a factor of 100 worse than the accuracy ofthe computation of the energy. Computing C with an accuracy of up to 1% using the first-order schemerequires more than 1000 steps in one orbit of the central body, whereas 10 steps per orbit is sufficient for theRunge–Kutta method.We also studied how many steps per orbit are required by the first-order scheme in order to obtain anaccuracy of 1% when drag force is added to the system (27). We solved the system (27), where the com-ponents of the drag force have the form (cid:18) ( v r − u r )Ω( r )St , ( v ϕ − u ϕ )Ω( r )St (cid:19) with the constant value of St.For St = 10 − , − , − we carried out the integration over a time t = 40 × π Ω − ( r ), and forSt = 10 − , − , , ,
100 over the shorter time t = 2 × π Ω − ( r ). We used the solution of this sameproblem obtained using the fourth-order method with a time resolution of 10 steps per orbit t = 2 π Ω − ( r ),as a reference. The relative errors in v r , computed in the semi-implicit, first-order scheme are presented in17he central and right panels of Fig. 12. For all values of St the semi-implicit scheme demonstrates first-orderin its approximation, but the actual accuracy attained in the solution depends on St. In the computationspresented in the right panel, St (cid:28)
1, and the dynamics of the particle are determined by the gas drag, not bythe balance between centrifugal and gravitational forces. In these computations, the lower the value of St, thehigher the accuracy of the numerical solution. In the computations with St ≥
1, presented in the central panel,the dynamics of the particle are determined by the balance between the centrifugal and gravitational forces,and the actual accuracy in v r does not depend on St ,and is close to the accuracy with which C is calculatedin the computations without friction using the same time steps.On the whole, the data presented in Fig. 12 and Table 1 show the expediency of using the first-order schemefor simulations of the dynamics of particles in a disk at orbital radii greater than 1 au. The sizes of dust grains in circumstellar and protoplanetary disks vary from submicron to several centimeters,and planetesimals have sizes of hundreds of kilometers. Therefore, the exchange of momentum between solidbodies and gas in disks occurs in a variety of regimes, determined by the sizes and velocities of the solid bodies:Epstein, Stokes, and Newton. We have compared different drag coefficients encompassing these regimes. Wehave shown that the standard astrophysical drag coefficient (8) is a discontinuous function for Mach numberscorresponding to the motion of the grains relative to the gas Ma > /
9, but is continuous when Ma < / < / > /
9, it is recommended to use the Henderson dragcoefficient (10)–(13), which is valid and continuous when Ma <
6, Re < × . Each computation of theHenderson coefficient requires of about 100 arithmetic operations.In addition, the need to compute the dynamics of bodies of different sizes in the same way imposes highdemands on numerical methods used for this purpose. We have shown that our semi-implicit, first-orderapproximation scheme, in which the interphase interactions are linearized and the relative velocity is calculatedimplicitly, while the stopping time and other forces, such as the pressure gradient and gravitation are calculatedexplicitly, preserve the asymptotic behavior of the solutions of problems with intense interphase interactions.This means that the scheme is suitable for computing the dynamics of both small bodies that are stronglycoupled to the gas and larger bodies for which the drag force depends non-linearly on the relative velocitybetween the gas and bodies, including the influence of the dispersed phase on the gas dynamics. A APPROACH TO DESCRIBING THE DYNAMICS OF DUSTIN A CIRCUMSTELLAR DISK
The suitability of a hydrodynamical approach to describing a medium is usually determined using two conditions: • The mean free path of the particles making up the medium should be much shorter than the length scaleof the system. We will take the disk scale height H as a length scale of the system.Further, some element of the medium is distinguished whose size is much less than the length scale of the medium,but much larger than the mean free path of the particles. In our case, this is the size of a computational gridcell, δR . • The number of particles in a given volume of the medium δR should be sufficiently large.The first condition is required so that the loss/gain of particles by an element of the medium as a result ofchaotic motions of the particles does not appreciably influence the mean characteristics of the element. Thesecond condition is required for computations using the particle distribution function in the six-dimensionalspace (coordinates and velocities) of the mean characteristics, such as the mean density in a cell, the meanvelocity, and the mean energy.There is no doubt that these conditions are satisfied by the gas in a circumstellar disk. However, we mustinvestigate this question for the dust component. We estimated the mean free path of a solid particle assumingthat the dust subdisk consists of monodisperse, spherical particles. In this case, λ d = 1 √ πa n d = 4 aρ s √ ρ d , (31)18ince, by virtue of their monodisperse nature, n d ( a ) = 3 ρ d πρ s a . (32)We assumed that the volume density of the dust in the equatorial plane of the disk was given by ρ d = Σ d H √ π , (33)where the height of the disk is H = (cid:115) k B Tµ H r GM (34)Let the surface density of the dust and the temperature be power-law functions of the radius, Σ d ( r ) =Σ ( r/ p , T = T ( r/ q . Adopting the characteristic values for a circumstellar disk Σ = 1 g cm − , T =300 K, µ H = 2 . M = 1 M (cid:12) , p = − , q = − .
5, we obtain ρ d = 7 . · − ( r/ − . g cm − . (35)It follows from (31) and (35)that grains smaller than 1 cm in the outer parts of a disk have mean free pathsmuch shorter than δR . It is clear that the dust subdisk can be modeled as a continuous medium in this case. Onthe other hand, these estimates indicate that λ d is close to H in inner regions of the disk. In this case, instead ofthe mean velocity over the volume, we must consider the velocity of the solid phase as a vector plus a dispersion.A possible development of the model could be the use of the Vlasov–Boltzmann equation [35] or moments ofthe Boltzmann equation [44]. The solution of the Vlasov–Boltzmann equation in Lagrange coordinates impliesthe use of discrete particles when modeling the dynamics of the solid phase [5, 50, 48]. Solution for the momentBoltzmann equation can be found in an Euler approach using grid methods [45].The moment Boltzmann equation take into account the possible anisotropic behaviour of the velocity dis-persion of the medium, but computing the components of the velocity dispersion requires the introduction ofadditional equations. When solving the Vlasov–Boltzmann equation using a particle method, the velocity dis-persion varies in a self-consistent way [9], but the difficulty is correctly computing the collision integral. In gasdynamics, this integral is zero for frequent elastic collisions. This may not be the case for the grains.We have concluded that dust in the main part of the disk can be represented as a continuous medium, whilethere may be regimes in the central region of the disk in which solid particles have a velocity dispersion, basedon estimates of the mean free paths of particles, neglecting their interactions with the gas. However, the sameconclusion is supported by estimates of the spatial scale on which exchanges of momentum between the gas anddust occur: λ stop = (cid:107) u − v (cid:107) t stop , (36)where t stop = St / Ω K . According to [8] Eq.(8), the azimuthal velocity of a particle can be estimated as u = v K (cid:18) −
11 + St η (cid:19) . (37)Here, η characterizes the degree of deviation of the velocity of the gas from the Keplerian value, and, for typicaldisk parameters, η ∼ .
01. We then find that λ stop = r St (cid:18) −
11 + St η (cid:19) ∼ r St . (38)It is clear that, when St < .
01, the velocities of grains within a single cell will differ little, while they willdiffer appreciably when St ∼
1. The simulation results [43] indicate that dust that has grown is located in theinner part of the disk.
B MODEL OF A STATIONARY GASEOUS DISK
We calculated the radial distributions of parameters of an axially symmetric, gaseous disk with inner boundary R min and outer boundary R max with mass M disk around a star with mass M using the model of [33], which isbased on the following assumptions: 19igure 13: (Left panel) Radial distributions of the surface density and temperature in a circumstellar disk forToy Models 1 and 2. (Right panel) Radial distribution of η and Toomre parameter for Toy Models 1 and 2. • the surface density of the matter in the disk displays a power-law dependence on the radius with power-lawindex ξ : Σ( r ) = Σ (cid:18) rR (cid:19) − ξ , (39)where Σ R ξ is given by (cid:90) π (cid:90) R max R min Σ( r ) r dr dϕ = M disk , (40) • the density and temperature of the gas in the disk are such that the Toomre parameter varies as Q ( r ) = c s ( r )Ω( r ) πG Σ( r ) = 2 (cid:18) R max r (cid:19) . , (41) • the gaseous disk is in equilibrium, so that u r = 0 , u ϕ ( r ) = (cid:115) GMr + rρ g ( r ) dp ( r ) dr , (42) • the height H and radius r of the disk are related as c s ( r ) V K ( r ) = H ( r ) r , (43) • the gas density at each radius does not vary in the vertical direction, so ρ g ( r ) = Σ( r ) H ( r ) . (44)We assumed ξ = 1 and used the above assumptions to find the distributions c s ( r ), ρ g ( r ), u ϕ ( r ) in the disk,which are required for the solution of (27).It follows from (40) that Σ = M disk (2 − ξ )2 πR ξ ( R − ξ max − R − ξ min ) . (45)It follows from (41) that c s ( r ) = Q ( r ) πG Σ( r )Ω( r ) = πrQ ( r ) G Σ( r ) V K ( r ) , (46)20nd therefore c s ( r ) V K ( r ) = πQ ( r ) G Σ( r ) rV = πQ ( r )Σ( r ) r M . (47)We then obtain from (47) and (43) H ( r ) = πQ ( r )Σ( r ) r M , (48)whence ρ g ( r ) = Σ( r ) H ( r ) = MπQ ( r ) r . (49)Note that, by virtue of (41) ρ g ( r ) = ρ g0 (cid:18) rR (cid:19) − . . (50)We took u ϕ to be related to V K as u ϕ = V + (cid:18) d ln ρ g d ln r (cid:19) c s = V (1 − η ) , (51)where η = c s V (cid:12)(cid:12)(cid:12)(cid:12) d ln ρ g d ln r (cid:12)(cid:12)(cid:12)(cid:12) = 2 . c s V = 2 . π Q ( r )Σ ( r ) r M , (52)so that u ϕ = V K (cid:114) − . π Q ( r )Σ ( r ) r M = V K (cid:115) − . ( r ) ρ ( r ) r (53)The distribution of the surface density and temperature for the disk in Toy Models 1 and 2 is shown inFig. 13. C ACKNOWLEDGMENTS
We thank V.N. Snytnikov and M.I. Osadchii for productive discussions of our results.
D FUNDING
This work was founded by the Russian Science Foundation (grant 17-12-01168).
References [1] V. V. Akimkin, M. S. Kirsanova, Y. N. Pavlyuchenkov, and D. S. Wiebe. Dust dynamics and evolution inexpanding H II regions. I. Radiative drift of neutral and charged grains.
MNRAS , 449:440–450, May 2015.[2] V. V. Akimkin, M. S. Kirsanova, Y. N. Pavlyuchenkov, and D. S. Wiebe. Dust dynamics and evolution inH ii regions - II. Effects of dynamical coupling between dust and gas.
MNRAS , 469:630–638, July 2017.[3] G. Albi, G. Dimarco, and L. Pareschi. Implicit-Explicit multistep methods for hyperbolic systems withmultiscale relaxation. arXiv e-prints , page arXiv:1904.03865, Apr 2019.[4] G.h Bagheri and C. Bonadonna. On the drag of freely falling non-spherical particles.
Powder Technology ,301:526–544, 06 2016.[5] X.-N. Bai and J. M. Stone. Particle-gas Dynamics with Athena: Method and Convergence.
ApJSS ,190:297–310, October 2010.[6] L. Barri`ere-Fouchet, J.-F. Gonzalez, J. R. Murray, R. J. Humble, and S. T. Maddison. Dust distributionin protoplanetary disks. Vertical settling and radial migration.
AAp , 443:185–194, November 2005.217] P. Ben´ıtez-Llambay, L. Krapp, and M. E. Pessah. Asymptotically Stable Numerical Method for MultispeciesMomentum Transfer: Gas and Multifluid Dust Test Suite and Implementation in FARGO3D.
ApJSS ,241:25, April 2019.[8] T. Birnstiel, M. Fang, and A. Johansen. Dust Evolution and the Formation of Planetesimals.
Space ScienceReviews , 205:41–75, December 2016.[9] R. A. Booth and C. J. Clarke. Collision velocity of dust grains in self-gravitating protoplanetary discs.
MNRAS , 458:2676–2693, May 2016.[10] S.-H. Cha and S. Nayakshin. A numerical simulation of a ’Super-Earth’ core delivery from 100 to 8 au.
MNRAS , 415:3319–3334, August 2011.[11] G.-Q. Chen, C. D. Livermore, and T.-P. Liu. Hyperbolic conservation laws with stiff relaxation terms andentropy.
Communications in Pure and Applied Mathematics , 47:787–830, june 1994.[12] P. Colella and P. R. Woodward. The Piecewise Parabolic Method (PPM) for Gas-Dynamical Simulations.
Journal of Computational Physics , 54:174–201, September 1984.[13] N. Cuello, J.-F. Gonzalez, and F. C. Pignatale. Effects of photophoresis on the dust distribution in a 3Dprotoplanetary disc.
MNRAS , 458:2140–2149, May 2016.[14] T. V. Demidova and V. P. Grinin. SPH simulations of structures in protoplanetary disks.
AstronomyLetters , 43:106–119, February 2017.[15] P. S. Epstein. On the Resistance Experienced by Spheres in their Motion through Gases.
Physical Review ,23:710–733, June 1924.[16] T. J. Haworth, J. D. Ilee, D. H. Forgan, S. Facchini, D. J. Price, D. M. Boneberg, R. A. Booth, C. J. Clarke,J.-F. Gonzalez, M. A. Hutchison, I. Kamp, G. Laibe, W. Lyra, F. Meru, S. Mohanty, O. Pani´c, K. Rice,T. Suzuki, R. Teague, C. Walsh, P. Woitke, and Community authors. Grand Challenges in ProtoplanetaryDisc Modelling.
PASA , 33:e053, October 2016.[17] C. B. Henderson. Drag coefficient of spheres in continuum and rarefied flows.
AIAA Journal , 14:707–707,1976.[18] M. Hutchison, D. J. Price, and G. Laibe. MULTIGRAIN: a smoothed particle hydrodynamic algorithmfor multiple small dust grains and gas.
MNRAS , 476:2186–2198, May 2018.[19] S. Ishiki, T. Okamoto, and A. K. Inoue. The effect of radiation pressure on dust distribution inside HIIregions.
MNRAS , 474(2):1935–1943, 2018.[20] S. Jin and C. D. Livermore. Numerical Schemes for Hyperbolic Conservation Laws with Stiff RelaxationTerms.
Journal of Computational Physics , 126:449–467, 1996.[21] A. Johansen and H. Klahr. Dust Diffusion in Protoplanetary Disks by Magnetorotational Turbulence.
ApJ ,634:1353–1371, December 2005.[22] G. Laibe and D. J. Price. Dusty gas with smoothed particle hydrodynamics - I. Algorithm and test suite.
MNRAS , 420:2345–2364, March 2012.[23] G. Laibe and D. J. Price. Dusty gas with smoothed particle hydrodynamics - II. Implicit timestepping andastrophysical drag regimes.
MNRAS , 420:2365–2376, March 2012.[24] G. Laibe and D. J. Price. Dust and gas mixtures with multiple grain species - a one-fluid approach.
MNRAS , 444:1940–1956, October 2014.[25] G. Laibe and D. J. Price. Dusty gas with one fluid.
MNRAS , 440:2136–2146, May 2014.[26] L.D. Landau and E.M. Lifshitz.
Fluid Mechanics. Vol. 6 (2nd ed.).
Butterworth-Heinemann, 1987.[27] P. Lor´en-Aguilar and M. R. Bate. Two-fluid dust and gas mixtures in smoothed particle hydrodynamics:a semi-implicit approach.
MNRAS , 443:927–945, September 2014.2228] P. Lor´en-Aguilar and M. R. Bate. Two-fluid dust and gas mixtures in smoothed particle hydrodynamicsII: an improved semi-implicit approach.
MNRAS , 454:4114–4119, December 2015.[29] J. J. Monaghan and A. Kocharyan. SPH simulation of multi-phase flow.
Computer Physics Communica-tions , 87:225–235, May 1995.[30] Y. Nakagawa, M. Sekiya, and C. Hayashi. Settling and growth of dust particles in a laminar phase of alow-mass solar nebula.
Icarus , 67:375–390, September 1986.[31] R. Pember. Numerical methods for hyperbolic conservation laws with stiff relaxation i. spurious solutions.
SIAM Journal on Applied Mathematics , 53:1293–1330, February 1993.[32] R. F. Probstein and F. Fassio. Dusty Hipersonic Flows.
AIAA Journal , 8:772–779, 1970.[33] W. K. M. Rice, G. Lodato, J. E. Pringle, P. J. Armitage, and I. A. Bonnell. Accelerated planetesimalgrowth in self-gravitating protoplanetary discs.
MNRAS , 355:543–552, December 2004.[34] D.V. Sadin and S.A. Odoev. Comparison of difference scheme with customizable dissipative propertiesand weno scheme in the case of one-dimensional gas and gas-particle dynamics problems.
Scientific andTechnical Journal of Information Technologies, Mechanics and Optics , pages 719–724, 07 2017.[35] V. N. Snytnikov and O. P. Stoyanovskaya. Clump formation due to the gravitational instability of amultiphase medium in a massive protoplanetary disc.
MNRAS , 428:2–12, January 2013.[36] J. M. Stone and M. L. Norman. ZEUS-2D: A radiation magnetohydrodynamics code for astrophysical flowsin two space dimensions. I - The hydrodynamic algorithms and tests.
ApJSS , 80:753–790, June 1992.[37] O. P. Stoyanovskaya, V. V. Akimkin, E. I. Vorobyov, T. A. Glushko, Y. N. Pavlyuchenkov, V. N. Snytnikov,and N. V. Snytnikov. Development and application of fast methods for computing momentum transferbetween gas and dust in supercomputer simulation of planet formation. In
Journal of Physics ConferenceSeries , volume 1103 of
Journal of Physics Conference Series , page 012008, October 2018.[38] O. P. Stoyanovskaya, T. A. Glushko, N. V. Snytnikov, and V. N. Snytnikov. Two-fluid dusty gas in smoothedparticle hydrodynamics: Fast and implicit algorithm for stiff linear drag.
Astronomy and Computing , 25:25–37, October 2018.[39] O. P. Stoyanovskaya, N. V. Snytnikov, and V. N. Snytnikov. Modeling circumstellar disc fragmentation andepisodic protostellar accretion with smoothed particle hydrodynamics in cell.
Astronomy and Computing ,21:1–14, October 2017.[40] O. P. Stoyanovskaya, V. N. Snytnikov, and E. I. Vorobyov. Analysis of methods for computing the trajec-tories of dust particles in a gas-dust circumstellar disk.
Astron.Rep. , 94:1033–1049, December 2017.[41] D V. Sadin. Tvd scheme for stiff problems of wave dynamics of heterogeneous media of nonhyperbolicnonconservative type.
Computational Mathematics and Mathematical Physics , 56:2068–2078, 12 2016.[42] E. I. Vorobyov. Embedded Protostellar Disks Around (Sub-)Solar Protostars. I. Disk Structure and Evo-lution.
ApJ , 723:1294–1307, November 2010.[43] E. I. Vorobyov, V. Akimkin, O. Stoyanovskaya, Y. Pavlyuchenkov, and H. B. Liu. Early evolution of viscousand self-gravitating circumstellar disks with a dust component.
Astronomy and Astrophysics , 614:A98, June2018.[44] E. I. Vorobyov and C. Theis. Boltzmann moment equation approach for the numerical study of anisotropicstellar discs.
MNRAS , 373:197–208, November 2006.[45] E. I. Vorobyov and C. Theis. Shape and orientation of stellar velocity ellipsoids in spiral galaxies.
MNRAS ,383:817–830, January 2008.[46] S. J. Weidenschilling. Aerodynamics of solid bodies in the solar nebula.
MNRAS , 180:57–70, July 1977.[47] F. L. Whipple. On certain aerodynamic processes for asteroids and comets. In A. Elvius, editor,
FromPlasma to Planet , page 211, 1972. 2348] C.-C. Yang and A. Johansen. Integration of Particle-gas Systems with Stiff Mutual Drag Interaction.
ApJS ,224:39, June 2016.[49] Z. Zhu, R. P. Nelson, R. Dong, C. Espaillat, and L. Hartmann. Dust Filtration by Planet-induced GapEdges: Implications for Transitional Disks.
ApJ , 755:6, August 2012.[50] Z. Zhu and J. M. Stone. Dust Trapping by Vortices in Transitional Disks: Evidence for Non-ideal Magne-tohydrodynamic Effects in Protoplanetary Disks.