推导尘埃漂移与气体的响应

选一个转动环片与停止时间

局部坐标 $x$ 指向径向外侧， $y$ 沿绕行方向， $z$ 指向垂直方向。环片以角频率 $\Omega$ 转动，速度以局部 Kepler 剪切为参考，因此参考半径上自由绕行的粒子，其相对速度为零。气体速度记为 $\mathbf u$ ，尘埃速度记为 $\mathbf v$ 。

固定线性阻力让尘埃加速度为 $-(\mathbf v-\mathbf u)/t_{\rm stop}$ 。定义 $\mathrm{St}=\Omega t_{\rm stop}$ 与局部体密度比 $\epsilon=\rho_{\rm d}/\rho_{\rm g}$ 。相同 Stokes 数在不同气体条件下，可以对应不同粒径。气体承受相反的单位体积力，除以自己的密度后，阻力加速度多出 $\epsilon$ 因子。

对应阅读：Youdin & Johansen 2007, Linear Evolution and Numerical Methods, v1 · Benítez-Llambay, Krapp & Pessah 2019, FARGO3D multifluid tests, v3

只给气体加压强作用

在均匀径向压强梯度下，把向外加速度写成 $2\eta v_{\rm K}\Omega$ 。正的 $\eta$ 表示压强向外降低，气体转得比 Kepler 速度慢。忽略平衡速度的径向导数、垂直引力与粘性。

径向、方位向各写一条稳态受力式，共四条。系数二与二分之一，来自科里奥利作用和 Kepler 背景剪切的共同影响。只有气体方程出现压强项。

0=2v_y-\frac{v_x-u_x}{\mathrm{St}},\qquad0=-\frac{v_x}{2}-\frac{v_y-u_y}{\mathrm{St}}

尘埃平衡转动与阻力；时间已用角频率倒数归一。

0=2u_y+2\eta v_{\rm K}+\frac{\epsilon(v_x-u_x)}{\mathrm{St}},\qquad0=-\frac{u_x}{2}+\frac{\epsilon(v_y-u_y)}{\mathrm{St}}

气体另有径向压强支撑。将气体方程加上密度比乘尘埃方程，阻力项相抵。

对应阅读：Youdin & Goodman 2005, Streaming Instabilities, arXiv v1 · Youdin & Johansen 2007, Linear Evolution and Numerical Methods, v1

把耦合平衡解出来

将气体方位式与 $\epsilon$ 倍的尘埃方位式相加，得到 $u_x=-\epsilon v_x$ ；径向加权和则给出 $u_y+\epsilon v_y=-\eta v_{\rm K}$ 。

尘埃径向式给出 $v_y=(1+\epsilon)v_x/(2\mathrm{St})$ ，方位式给出 $u_y=v_y+\mathrm{St}v_x/2$ 。一起代回加权径向关系， $v_x$ 前面的系数变成 $[(1+\epsilon)^2+\mathrm{St}^2]/(2\mathrm{St})$ ，于是得到下面的速度。

Q=(1+\epsilon)^2+\mathrm{St}^2,\quad v_x=-\frac{2\mathrm{St}}Q\eta v_{\rm K},\quad u_x=\frac{2\epsilon\mathrm{St}}Q\eta v_{\rm K}

正压强支撑参数下，尘埃向内，气体向外；均匀平衡中，质量加权径向流量相抵。

v_y=-\frac{1+\epsilon}{Q}\eta v_{\rm K},\qquad u_y=-\frac{1+\epsilon+\mathrm{St}^2}{Q}\eta v_{\rm K}

这些小偏差加上 Kepler 绕行，才得到完整轨道速度。

对应阅读：Youdin & Johansen 2007, Linear Evolution and Numerical Methods, v1

完整算例：漂移有多快

取示例轨道半径 $1.50\times10^{11}\,\mathrm m$ 、 $v_{\rm K}=3.00\times10^4\,\mathrm{m\,s^{-1}}$ 与 $\eta=10^{-3}$ 。反馈可忽略、 $\mathrm{St}=1$ 时，径向速度为 $-30.0\,\mathrm{m\,s^{-1}}$ 。半径除以速率，局部漂移时标约为158年。这是在该半径的尺度估算，没有积分穿过整个变化中的盘。

保持停止时间，改取 $\epsilon=1$ ，此时 $Q=5$ ，尘埃以 $12.0\,\mathrm{m\,s^{-1}}$ 向内，气体以相同速率向外。两者密度相同，径向净质量流量为零。气体响应了阻力，尘埃的漂移也慢下来。

让初始气压峰自行演化

有限陷阱有空间梯度，仅用局部稳态速度不能让它演化。两种物质分别满足连续性 $\partial_t\rho+\partial_x(\rho v_x)=0$ ，同时输运动量，气体采用 $P=c_s^2\rho_{\rm g}$ ，再交换阻力。在径向实验中，长度以 $H=c_s/\Omega$ 、速度以 $c_s$ 为单位；后面的线性波则采用 $\eta r$ 与 $\eta v_{\rm K}$ ，两套尺度要分开。

气体初态在周期域 $L=8H$ 上正比于 $\exp[A\cos(2\pi(x-x_{\rm peak})/L)]$ ，再归一到平均密度一。尘埃最初均匀，速度从局部耦合漂移平衡给出。空间梯度随后会让两种物质继续演化，气压峰也可变形。

下面的实际固定结果采用128格，预热八个角频率倒数。场景也读取同一份气尘场；尘埃峰来自计算的聚集。

八个角频率倒数后的实际合成有限体积结果：128周期格，气体声速为一，域宽八个气体标高，初始对数密度振幅0.6，尘埃初态均匀。freeze-radial.ts 保存两种物质的密度、动量场，以及分辨率和时间步长比较。
情形	Stokes 数	初始平均密度比	尘埃峰均比	气体最大径向速率 · 声速单位
默认陷阱	$0.300$	$0.300$	$1.74$	$0.00508$
紧耦合	$0.0300$	$0.300$	$1.07$	$0.000986$
快速漂移	$1.00$	$0.300$	$2.97$	$0.00467$
强反馈	$0.300$	$3.00$	$1.15$	$0.0121$

对应阅读：Youdin & Johansen 2007, Linear Evolution and Numerical Methods, v1 · Benítez-Llambay, Krapp & Pessah 2019, FARGO3D multifluid tests, v3

把物理反馈与数值平滑分开

有限体积法把每个面的通量从一格减去、向相邻格加回，包括周期边界，因此分别守恒两种物质的质量。局部阻力步保持质心速度，令相对速度按 $\exp[-(1+\epsilon)\Delta t/t_{\rm stop}]$ 衰减。只考虑阻力时，相对动能会下降。

\frac{dE_{\rm rel}}{dt}=-\frac{\rho_{\rm d}}{t_{\rm stop}}|\mathbf v-\mathbf u|^2,\quad E_{\rm rel}=\frac{\rho_{\rm g}\rho_{\rm d}}{2(\rho_{\rm g}+\rho_{\rm d})}|\mathbf v-\mathbf u|^2

阻力耗散相对运动能量，守恒混合物动量。演变的剪切与气压可提供运动，等温模型则带走热量。

对应阅读：Benítez-Llambay, Krapp & Pessah 2019, FARGO3D multifluid tests, v3 · Youdin & Johansen 2007, Linear Evolution and Numerical Methods, v1

检查已经分辨到的尺度

径向求解器采用一阶空间 Rusanov 通量、二阶 SSP Runge–Kutta 输运，以及对称的精确阻力/周转源项步。空间数值扩散会平滑陡结构。默认时刻，128格与1024格参考在整个周期域上的相对 L1 尘埃密度差异为百分之2.92，峰高低约百分之6.1。这些是实测的分辨率差异。

无压流体还要求每个位置只有一个尘埃速度，不能表示粒子多流交叉或稠密焦散。有限的早期实验适合讨论捕获与气体反馈，后续非线性集中需要粒子法或更完整的动力学处理。径向方程也缺少线性流动不稳定所需的垂直模式；下一章把这个方向补上。