0 引 言
随着近年来对太阳系外缘天体研究的不断深入，冰质小天体逐渐成为天文学和行星科学领域研究的热点。冰质小天体泛指太阳系中具有显著挥发分冰含量，且其活动主要由冰升华驱动的固态天体。这类天体通常具有低密度、高孔隙度的冰-尘埃混合结构，分布在太阳系较冷的区域，如小行星带外围、“柯伊伯”带（Kuiper Belt）、“奥尔特”云（Oort Cloud）等。主要包括短周期彗星（如67P/Churyumov-Gerasimenko，以下简称“67P”）、主带彗星（如彗星133P（133P/Elst-Pizarro）），以及外太阳系的小型冰质天体（如“柯伊伯”带天体中的较小成员）。其中，彗星因其活跃的挥发分气体与尘埃的喷发活动（形成可见的彗发与彗尾），成为研究太阳系原始物质演化的关键对象。
因此，本文选用67P作为研究对象，以探究冰质小天体近核区域的尘埃喷发动力学机制。该彗星轨道周期为6.45年，自转周期12.4 h，并于2021年11月2日抵达近日点。作为“罗塞塔”任务的主要探测目标，67P积累了包括高分辨率形貌、气体成分及尘埃分布在内的多维度观测数据，为近核区域尘气耦合研究提供了不可替代的实测依据。
当彗星运行至近日点附近时，太阳辐射加热彗核表面，导致水冰等挥发分物质升华。升华气体向太空高速膨胀，速度可达数百米每秒[1]。彗核表面的尘埃粒子在这些高速挥发分气体的拖曳作用下脱离表面，向外喷发，最终形成尘埃喷流。
针对此过程中尘气耦合的喷发动力学机制，现有数值模型主要分为两类：一类基于粒子动力学方法，另一类则基于流体力学方法。使用粒子动力学方法的研究人员通常基于直接模拟蒙特卡洛（Direct Simulation Monte Carlo，DSMC）方法对描述近核流场的玻尔兹曼方程进行求解[2]。与流体力学方法相比，粒子动力学方法的主要优势在于可适用于任何程度的非平衡流场与稀薄气体流场[2]。彗发动力学的DSMC数值模拟始于1988年Combi等[3]的奠基性工作。2008年，Tenishev等[4]通过耦合彗核边界层热物理过程与彗发的光化学过程，首次建立了涵盖彗核表面至其外$1 \times {10^6}{\text{ km}}$ 区域的全域模型，为后续研究提供范式基础。2016年，Lai等[5]针对67P构建了气体-尘埃耦合流场，定量揭示了近日点附近尘埃自南半球向北半球的跨区域输运机制。第二类模型通常与使用DSMC获得的计算结果进行比对，以验证其适用性，进而论证在一定情况下，流体力学方法可视为计算效率更高的简化算法。2017年，Shou等[6]开发了基于BATS-R-US程序的多流体尘埃模型，能够模拟瞬态现象，与DSMC方法具有很好的一致性，并对67P的尘埃观测结果进行了部分解释。
本文结构如下：第1节建立三维尘埃粒子动力学模型，针对全球尘气喷发（气体均匀包裹彗核）与局域喷发（限定小面积喷发源）两种模式，分别构建了包含自转效应的控制方程；第2节基于67P的物性参数，开展全球与局域喷发的数值仿真，系统分析了尘埃的驱离阈值、运动轨迹及未逃逸粒子落点的经纬度分布规律；第3节归纳全文主要结论，阐明喷发尺度与自转效应对近核区尘埃动力学的调控机制。
1 三维条件下的尘埃动力学建模
首先，本文假设67P为一个半径${r_0} = {\text{2 000 m}}$ 的均匀球体，并考虑彗星自转对尘埃粒子运动产生的影响。值得注意的是，67P实际为不规则双叶形结构，此简化假设遵循2021年Fink等[7]模型。为了简化计算，本文设立一固连坐标系，该坐标系随着彗核以相同角速度自转。坐标系原点位于彗核中心，令z轴与自转轴重合，指向北极；x轴位于彗星的赤道平面，并指向本初子午线的方向；y轴位于彗星的赤道平面，并指向90°东经方向。
显然，此固连坐标系是一个非惯性系。本文假设尘埃粒子在这个非惯性系下仅受到四个力的作用，分别为气体拖曳力$ {{\boldsymbol{F}}_{{\text{gas}}}} $ 、彗核引力$ {{\boldsymbol{F}}_{{\text{gravity}}}} $ ，以及由于彗核自转所产生的离心力$ {{\boldsymbol{F}}_{{\text{centrifugal}}}} $ 与科氏力$ {{\boldsymbol{F}}_{{\text{coriolis}}}} $ 。其运动方程由牛顿第二定律给出：
式中，$ {m_{\mathrm{d}}} $ 为尘埃粒子的质量，$ {{{\boldsymbol{a}}}_{\mathrm{d}}} $ 为尘埃粒子的加速度。下面将分别推导全球喷发与局域喷发场景中尘埃粒子所受各力的具体表达式。
1.1 三维全球喷发尘埃动力学建模
在三维全球尘气喷发模型中，假设升华产生的气体可以均匀地包裹整个彗核，并向外进行各向同性膨胀。图1 中未画出离心力与科氏力。其中，蓝色圆代表彗核，红色点代表尘埃粒子，覆盖在彗核外边的灰色圆环表示彗核尘气喷发时的气体环境。因此，尘埃粒子所受到的气体拖曳力的方向如下：将尘埃粒子与彗核中心进行连线，尘埃粒子所受到的气体拖曳力的方向沿着这条直线并指向背离彗核的方向，如图1 中$ {{\boldsymbol{F}}_{{\text{gas}}}} $ 所示，大小为：
式中，$ {\sigma _{\mathrm{d}}} $ 为尘埃粒子的横截面积，$ {\rho _{\mathrm{g}}} $ 为气体密度，$ {C_{\mathrm{d}}} $ 为气体阻力系数。本文取$ {C_{\mathrm{d}}} = 2 $ ，与Fink等[7]模型保持一致。$ {{\boldsymbol{v}}_{\mathrm{g}}} $ 为气体速度。$ {{\boldsymbol{v}}_{{\text{dg}}}} $ 比较特殊，为尘埃粒子的速度$ {{\boldsymbol{v}}_{\mathrm{d}}} $ 与粒子此时所处位置气体速度$ {{\boldsymbol{v}}_{\mathrm{g}}} $ 的点积，意为尘埃粒子在此地气体速度方向上的速度分量大小。
因此，气体拖拽力${{\boldsymbol{F}}_{{\text{gas}}}}$ 可表示为：
式中，$ x、y、z $ 为尘埃粒子的位置坐标。
然后，计算引力${{\boldsymbol{F}}_{{\text{gravity}}}}$ ，引力的大小由万有引力公式可以得出：
式中，${m_{\mathrm{d}}}$ 为尘埃粒子的质量，$G$ 为万有引力常数，${M_{\mathrm{n}}}$ 为彗核质量，$ r $ 为尘埃粒子到彗核中心的距离。
引力的方向指向彗核中心，因此引力${{\boldsymbol{F}}_{{\text{gravity}}}}$ 可表示为：
在此非惯性系下，粒子除了受到上述提到的气体拖曳力${{\boldsymbol{F}}_{{\text{gas}}}}$ 和彗核引力${{\boldsymbol{F}}_{{\text{gravity}}}}$ 的作用，还会受到离心力${{\boldsymbol{F}}_{{\text{centrifugal}}}}$ 与科氏力${{\boldsymbol{F}}_{{\text{coriolis}}}}$ 的作用，离心力${{\boldsymbol{F}}_{{\text{centrifugal}}}}$ 的计算表达式为：
式中，$ \omega $ 为67P自转的角速度。$ {{\boldsymbol{R}}} $ 表示由彗核自转轴指向尘埃粒子的距离向量，将其分解为分量形式，可得：
科氏力${{\boldsymbol{F}}_{{\text{coriolis}}}}$ 的计算表达式如下：
式中，$ {m_{\mathrm{d}}} $ 为尘埃粒子的质量，$ {{\boldsymbol{\omega}} } $ 为67P自转的角速度矢量，$ {{{\boldsymbol{v}}}^{'}} $ 为粒子相对于非惯性系的速度矢量，将向量$ {{\boldsymbol{\omega}} } $ 与$ {{{\boldsymbol{v}}}^{'}} $ 展开为分量形式，计算可得：
将式(3)、式(5)、式(7)、式(9)代入式(1)，可得：
那么，尘埃粒子的加速度可表示为：
另外，有：
式中，$ {z_{{\text{gn}}}} $ 为表面气体产率，意为在单位时间单位面积下，由于气体喷发所产生的气体分子数；$ {m_{{\mathrm{gas}}}} $ 为气体分子的分子量；$ {m_{\mathrm{u}}} $ 为原子质量，${m_{\mathrm{u}}} = 1.66 \times {10^{ - 27}}{\text{ kg}}$ 。
将式(12)和式(13)代入式(11)，那么其中等号右边第一项可化简为[7]：
那么，式(11)可化简为：
至此，完成了三维全球尘气喷发下的尘埃动力学模型的构建。
1.2 三维局域喷发尘埃动力学建模
在三维局域喷发模型中，尘埃粒子所受到的彗核引力${{\boldsymbol{F}}_{{\text{gravity}}}}$ 、离心力${{\boldsymbol{F}}_{{\text{centrifugal}}}}$ 和科氏力${{\boldsymbol{F}}_{{\text{coriolis}}}}$ 的大小与方向的计算均与1.1节相同。然而，由于喷流空间受限，气体拖曳力的计算公式需进行修正。
在此情况下，本文假设喷发区域为圆形（直径为${d_0}$ ），且喷出的气体不再可以均匀地包裹整个彗核，而是形成一个气体锥（如图2 所示）。这意味着仅彗核中的局部区域产生喷发，喷发面积只占彗核表面积的一小部分。当尘埃粒子位于锥面内时，粒子便会受到气体拖曳力作用；而当粒子不断运动从而位于锥面外时，则不再受到气体拖曳力影响。图2 为二维示意图，展示了这一过程。需注意图2 中未画出离心力与科氏力。
从图2 中可以看出，喷出气体形成的锥面延伸至彗核内部，汇聚于一点a，形成一个完整的圆锥。该圆锥的轴线与直线aO（即彗核中心O与气体锥顶点a的连线）重合。在此模型中，尘埃粒子所受到的气体力不再与尘埃粒子与彗核中心的连线Op重合，而是与粒子位置p与圆锥顶点a的连线ap重合，方向指向背离彗核的方向。
为便于计算，假设此圆形喷发区域中心点（即圆心）的纬度为M°，经度为N°，圆锥的锥角为53°（此角度仅为方便计算，实际中可假设为任意可能的度数）。据此，可以计算出圆锥顶点a距离彗核表面的距离:
从而可以得到点a在此非惯性坐标系下的坐标：
式中，$ {r_0} $ 为彗核半径。
通过上述描述，可以求得气体拖曳力${{\boldsymbol{F}}_{{\text{gas}}}}$ 的方向应为尘埃粒子的坐标减去点a坐标这一向量的单位向量，即：
气体力大小的计算表达式与式(2)相同：
假设尘埃粒子半径为s，密度为$ {\rho _{\mathrm{g}}} $ 的刚性球体，气体阻力系数$ {C_{\mathrm{d}}} = 2 $ ，那么化简可得：
在此模型下，表面气体产率$ {z}_{\text{gn}}\left(r\right) $ 的计算表达式为：
式中，$ {z_{{\text{gn}}}}\left( {{r_0}} \right) $ 为位于彗核表面位置的表面气体产率，$ d_0^{} $ 为喷发源的直径，D为尘埃粒子此时所处圆锥底面的直径，具体计算表达式为：
将式(1)、式(5)、式(7)、式(9)、式(17)、式(18)、式(20)~(22)组合起来，便可以得到三维条件下的局域尘气喷发尘埃加速度模型。
2 三维喷发下的尘埃动力学行为分析
2.1 三维全球尘气喷发下的尘埃动力学行为分析
2.1.1 可驱离尘埃粒子临界半径的计算
在三维全球喷发模型中，首先计算可驱离尘埃粒子的临界半径，即可克服引力脱离彗核表面的最大粒子尺寸。由式(15)可知，当尘埃粒子位于彗核表面，还未开始运动（即初速度为0）时，科氏力为0。因此，本文在计算可驱离尘埃粒子临界半径时，仅考虑气体拖曳力${{\boldsymbol{F}}_{{\text{gas}}}}$ 、彗核引力${{\boldsymbol{F}}_{{\text{gravity}}}}$ 和离心力${{\boldsymbol{F}}_{{\text{centrifugal}}}}$ 的平衡关系。
在运动初始阶段，气体拖曳力${{\boldsymbol{F}}_{{\text{gas}}}}$ 与离心力${{\boldsymbol{F}}_{{\text{centrifugal}}}}$ 均指向背离彗核的方向，而彗核引力${{\boldsymbol{F}}_{{\text{gravity}}}}$ 的方向朝向彗核中心。因此，在计算可驱离尘埃粒子临界半径时，气体拖曳力${{\boldsymbol{F}}_{{\text{gas}}}}$ 和离心力${{\boldsymbol{F}}_{{\text{centrifugal}}}}$ 被视为使得尘埃粒子运动的力，而彗核引力${{\boldsymbol{F}}_{{\text{gravity}}}}$ 作为阻碍尘埃粒子运动的力。
由式(6)可知，当尘埃粒子位于彗核赤道时，尘埃粒子所受到的离心力${{\boldsymbol{F}}_{{\text{centrifugal}}}}$ 最大，其大小为：
因此，位于赤道上的尘气喷发可以抬起全彗核最大尺寸的尘埃粒子，其粒子半径的计算表达式为：
三维全球尘气喷发模型的计算参数见表1 。通过计算，可解得${s_{\max }} = 12.29{\text{ mm}}$ 。但在赤道以外的地区，随着纬度升高，离心力随之衰减，其提供的抬升辅助减弱，因此可驱离尘埃粒子的临界半径逐渐减小。在两极附近，离心力趋近于零，此时临界半径达到全局最小值。经过计算，发现在南、北两极，可驱离尘埃粒子的临界半径${s_{\max }} = 9.36{\text{ mm}}$ 。
表 1 三维全球尘气喷发模型的计算参数Table 1 Parameters for the 3D global dust-gas eruption model |
| 计算参数 | 数值 |
| 气体分子量 | |
| 喷出气体速度 | |
| 全球气体产率 | |
| 彗核半径 | |
| 尘埃粒子密度 | |
| 万有引力常数 | |
| 彗核质量 | |
| 彗核自转角速度 |
2.1.2 尘埃尺寸对运动轨迹的影响
由此，可根据观测到的67P的全球气体产率$Q(t) = 5 \times {10^{27}}{\text{ molecules/s}}$ [7]，计算其表面气体产率，将上述参数代入式(24)，使用4/5阶自适应步长龙格-库塔（Runge-Kutta）法（Matlab ode45求解器），进行迭代求解。初始条件设为粒子静止于赤道（0°N，0°E），以简化分析。该模型同样适用于彗核表面的任意初始位置。
由第2.1.1节可知，在0°N，可驱离尘埃粒子的最大半径不能超过$ 12.29\;\mathrm{m}\mathrm{m} $ 。因此，本文取尘埃粒子半径分别为$ 4.5\text{ mm}、\text{ 8}\text{.0 mm}、\text{ }11.\text{5 mm}、\text{ 15}\text{.0 mm} $ ，在这四种情况下，探究尘埃尺寸对运动轨迹的影响。
由图3 可以看出，在全球尘埃喷发模型下，粒子尺寸对其动力学行为起着重要作用：对于半径为4.5 mm的小尺寸粒子，因其质量小，气体拖曳力起主导作用，使其不断加速，最终获得足够速度逃逸；对于半径为8.0 mm的中等尺寸粒子，其运动处于临界状态，气体拖曳力、彗核引力和离心力接近平衡，导致粒子动能与势能反复转换，从而形成了复杂的振荡轨迹，但最终因动能不足，无法完全克服引力势垒而回落；对于半径为11.5 mm的大尺寸粒子，因其质量大，气体拖曳力难以有效加速，其运动很快由引力主导，因迹近似于简单抛物线，迅速回落至表面；对于半径为15.0 mm的粒子，其尺寸已远超临界驱离半径，气体拖曳力甚至无法克服彗核引力将粒子抬离表面，因此未产生有效位移（图3 中未显示其轨迹）。
2.1.3 尘埃尺寸对落点的影响
由图3 可以看出，即使最终回落的粒子，其弹道轨迹也存在显著差异。对于中等尺寸的粒子，其运动轨迹较为复杂；而大尺寸粒子在引力主导下，轨迹相对简单。鉴于上述差异，为有效归纳落点分布规律，本文选取了运动轨迹均由引力主导的9.5 mm、10.0 mm和11.0 mm粒子进行落点分析，以聚焦于尺寸效应本身，避免复杂轨迹引入的干扰。
从图4 可以看出，尽管不同尺寸尘埃粒子的落点曲线存在一定差异，但整体趋势一致：起始位置位于北半球的尘埃粒子，其落点纬度均低于起始纬度，即落点偏南；反之，起始位置位于南半球的尘埃粒子，其落点纬度均高于起始纬度，即落点偏北。该现象主要是由于尘埃粒子受到离心力的作用，在本文坐标系中，离心力的方向始终垂直于彗核的自转轴。对于从北半球向外运动的粒子，离心力可分解为一个背离粒子与彗核连线方向的分力和垂直于此方向且指向偏南的分力，该分力使得粒子轨迹向南偏转（相对于速度方向），导致其落点偏南；反之，从南半球出发的粒子，离心力的对应分力使其向北偏转，导致其落点偏北。
由图4 还可以看出：在高纬度地区，不同尺寸尘埃粒子的落点曲线均呈现斜率为1的线性关系（即落点纬度等于起始纬度）。这表明，在这些高纬度区域，尘埃粒子所受的气体拖曳力与离心力不足以克服重力将其抬升，导致粒子未被有效抬升而直接“落”在起始位置附近。图5 进一步佐证了这一结论：在各条曲线对应的高纬度区间（即落点纬度等于起始纬度的区间），粒子的落点经度始终保持为0°，与其初始经度一致。此外，根据第2.1.1节可驱离尘埃粒子临界半径的计算（南、北两极处最大可被驱离半径约为9.36 mm），可以推断在南、北两极附近，图5 中所示尺寸（9.5 mm、10.0 mm、11.0 mm）的尘埃粒子半径均明显大于当地的最大可被驱离半径，因此无法被有效抬升，这与图4 中观测到的现象完全吻合。
本文推测，这一异常行为是小尺寸尘埃粒子所受多种力非线性耦合的体现。对于尺寸较大的粒子（如11.0 mm），其运动主要受彗核引力影响，气体拖曳力与惯性力的影响相对较弱，因此其轨迹和落点表现出更简单、一致的规律性，即经度偏移随纬度升高而单调减小。而对于小尺寸尘埃粒子在低纬度区域，随着纬度升高，其运动轨迹的南、北向偏移量更大（见图4 ），运动的时间更长，科氏力的累积效应也更显著，从而暂时克服了彗核引力的主导作用，表现为经度偏移量的短暂增加。
2.2 三维局域尘气喷发下的尘埃粒子动力学行为分析
2.2.1 可驱离尘埃粒子临界半径的计算
与第2.1节相同，为了探究三维局域尘气喷发下的尘埃粒子动力学行为，本文首先对可驱离尘埃粒子的临界半径进行了计算。仍以赤道位置的尘气喷发为例，该喷发可以抬起全彗核最大尺寸的尘埃粒子。其粒子半径的计算表达式为：
基于2008年Tenishev等[4]的数据，在计算$ {z}_{\mathrm{g}\mathrm{n}}\left({r}_{0}\right) $ 过程中，选取了日心距1.29 AU、太阳天顶角0°条件下67P的水汽喷发通量。根据其数据图估算，此时水通量约为$ 3.5\times {10}^{20}\text{molecules}/({\text{m}}^{2}\cdot \text{s}) $ ，其余参数与第2.1.1节保持一致（见表1 ）。由此解得${s_{\max }} = 43.426{\text{ mm}}$ 。但在赤道以外的地区，尤其是靠近两极的地区，由于离心力的减小，可驱离尘埃粒子临界半径会小于43.236 mm。同样，在南、北两极，经计算得可驱离尘埃粒子临界半径为${s_{\max }} = 32.936{\text{ mm}}$ 。由此可见，在本文的参数设置下，局域喷发的抬升能力显著强于全球喷发。
2.2.2 尘埃尺寸对运动轨迹的影响
同样地，基于2015年Gulkis等[8]对67P在3.5 AU处的观测数据（气体喷发速度$ {\boldsymbol{v}}_{{\mathrm{g}}}=680\;\mathrm{m}/\mathrm{s} $ ），本文仍使用$ {\boldsymbol{v}}_{{\mathrm{g}}}=700\;\mathrm{m}/\mathrm{s} $ 作为仿真基准值，并假设67P喷发的气体分子均为水分子，其物性参数与2021年Fink等[7]保持一致，具体可见表1 。同时，本文依然将初始条件设为：粒子静止于赤道（0°N，0°E）（实际上可将尘埃粒子置于任意位置），仅在计算关键参数$ {z}_{\mathrm{g}\mathrm{n}}\left({r}_{0}\right) $ 时，采用了2008年Tenishev等[4]提供的日心距1.29 AU、太阳天顶角0°条件下67P的水汽喷发通量数据，根据其数据图估算，此时水通量${z_{{\text{gn}}}}({r_0})$ 约为$ 3.5\times {10}^{20}\;\text{molecules}/({\text{m}}^{2}\cdot \text{s}) $ 。为了对比，本文仍取尘埃粒子半径为4.5 mm、8.0 mm、11.5 mm、15.0 mm，探究在这四种情况下，尘埃尺寸对运动轨迹的影响。
2.2.3 尘埃尺寸对落点的影响
3 结 语
本文将2021年Fink等[7]提出的简化模型拓展至三维情形，并同步整合了彗核自转效应（离心力与科氏力）及喷发尺度（全局性与区域性）的影响，建立了三维全球尘气喷发与三维局域尘气喷发（限定喷发源面积）的尘埃粒子动力学模型。基于67P的物理参数开展数值模拟，系统探究了尘埃粒子的驱离阈值、运动轨迹及落点分布规律。主要结论如下：
（1）在考虑彗核自转效应后，尘埃粒子临界驱离半径随纬度升高而减小，并计算了相应模型下临界半径的最大值与最小值。
（2）尘埃粒子的运动轨迹受其尺寸和初始位置的影响显著。小尺寸粒子易获得逃逸速度；中等尺寸粒子轨迹复杂多变，最终回落彗核；大尺寸粒子主要受引力主导，轨迹相对简单。
（3）对于尘埃粒子的落点分布，研究发现，全球喷发的落点偏移规律显著：起始于北半球的尘埃粒子落点普遍偏南（纬度降低），起始于南半球的尘埃粒子落点普遍偏北（纬度升高）。而局域喷发与喷发源面积显著相关，若喷发源面积很小，其纬度偏移量很小。落点经度偏移量与粒子尺寸、初始纬度及喷发模式相关，通常情况下，粒子尺寸越小、初始纬度越低，经度偏移越大。值得注意的是，在全球喷发模型下的低纬度区域，本文观测到了可能源于自转与气体拖曳非线性耦合的经度偏移异常现象，其内在物理机制仍有待后续深入研究。
综上所述，本文构建的简化三维动力学模型有效揭示了喷发尺度与彗核自转效应对近核区尘埃动力学的调控机制，为深入理解冰质小天体近核区尘埃环境的时空演化特征提供了理论工具，并为评估深空探测器在飞越或伴飞此类天体时遭遇尘埃撞击的风险奠定了基础。