SLM技术是利用高能量的激光束照射预先铺覆好的粉末材料，将其直接熔化并固化成形获得具有较好几何精度和力学性能的三维零件。SLM成形过程是：①在基板上预先铺设一定厚度的粉末层；②激光束有选择性地照射粉末层颗粒使其熔化，而后冷却固化；重复步骤①和步骤②进行逐层加工，直至完成整个3D成形件。步骤①为铺粉过程；步骤②为SLM过程。大量研究表明，SLM成形件质量的主要控制因素为[12]：激光工艺参数、粉材配方、粉床质量。其中，粉床质量(通常用堆积密度、均匀性、表面粗糙度等指标定量描述)对粉末层内部热能传递具有显著影响，可以改变成形件的宏观性能。若粉床内部存在空腔或移位缺陷、颗粒呈松散排列状态，则容易形成熔融缺陷(如球化、孔隙或剥蚀)，降低成形件的宏观质量；反之，若粉床内部颗粒排列均匀、紧密、无明显结构缺陷，则更容易获得较高质量的成形件。Ziegelmeier等[13]研究发现，成形件的拉伸强度和断裂韧性随着粉床堆积密度的增加而提高；粉床的表面粗糙度越高，成形件的表面越粗糙。鉴于粉床质量对SLM成形产品几何精度和力学性能的重要影响，需要开展粉床形成机理研究，从细观角度出发，考察铺粉工艺参数、环境因素、粉末颗粒自身物理力学属性等对SLM铺粉过程中颗粒运动行为的作用机制，明确控制粉床质量的主要宏细观参数，在此基础上，提出铺粉设备和工艺参数的优化设计方案。
SLM铺粉过程中，颗粒运动具有随机性和不连续性，导致粉床中颗粒排列具有明显的时变特性，在宏观上表现为粉床具有复杂的流变特性。考虑到大多数连续介质模型及算法(如有限元法)在模拟铺粉过程和粉床细观结构方面存在较大困难，本文提出采用基于离散介质模型的离散单元方法(Discrete Element Method，DEM)。DEM在处理固体颗粒细观运动行为方面具有天然优势，其以单个粒子为研究对象，基于粒子运动控制方程和力-位移准则求解不同时刻粒子的位置、速度、加速度、作用力和力矩等信息，因此可以更加真实地模拟铺粉过程中粉末颗粒间、颗粒与壁面间的细观相互作用。最近，DEM数值仿真方法在SLM铺粉过程中颗粒细观运动行为及粉床流变学特性研究中获得应用。例如，Zohdi[14]采用DEM数值仿真技术模拟了球形颗粒的沉积过程。向召伟等[15]利用DEM算法模拟了单粒径颗粒群、双粒径颗粒群、高斯粒径分布球颗粒系统的随机堆积行为及其铺层和压实过程，并考察了铺粉厚度对粉床质量的影响。Chen等[16]利用三维DEM程序和刚性球颗粒进行铺粉过程数值模拟，并考察范德华力对粉末层状态的影响。Parteli、Poschel[17]和Haeri等[18-19]提出非球形颗粒的建模方法，采用多球近似法建立真实颗粒几何模型，并利用三维DEM模拟了非球形颗粒组成的粉末材料的铺粉过程，定量分析了粉床质量的影响因素，并在此基础上提出了一种铺粉装置的优化设计方案。目前，DEM已经成为铺粉过程细观尺度研究的理想工具。
本文以德国Solution SLM125激光选区熔化成形设备的铺粉装置和模拟月壤为研究对象，建立铺粉仿真三维简化模型，采用DEM进行模拟月壤铺粉过程细观动态行为研究，考察铺粉速度、环境因素等对粉床质量的影响规律，为SBAM-ISRU技术发展提供理论指导和技术支撑。
1 DEM模型及方法 1.1 颗粒几何模型 月壤是在缺乏氧气、液态水、风和生命活动，由大大小小的陨石撞击、宇宙射线和太阳风不断轰击、月表大幅度温差变化导致月岩热胀冷缩破碎等空间风化作用下形成的，而各国研制的模拟月壤(如JSC系列、MLS系列、TYII系列)多以地球火山灰为原料，通过机械破碎、分级处理及矿物成分调配获得，因此，月壤和模拟月壤颗粒形态复杂，从球形到极端棱角形都存在。笔者研究发现，月壤和模拟月壤颗粒形态具有分形特征，在不同尺度下具有一定的统计自相似性[20]，并针对模拟月壤的形状指数和分形特征，提出一个框架性的非球形粒子模型——叠加球模型[21-22]。
根据模拟月壤的SEM照片(见图 1)[19]，采用叠加球模型方法，建立了模拟月壤4种典型颗粒形状的三维几何模型，如图 2所示。在叠加球模型中，每种颗粒形状均由一个基本球单元和若干个附加球单元构成，颗粒内部球单元之间允许重叠且不考虑接触力，颗粒假设为均质三维刚体，在基本球单元质心建立一个局部坐标系，该局部坐标系的3个主轴方向与全局笛卡儿坐标系相同，附加球单元的数量、半径及在局部坐标系的位置则根据真实颗粒形状来确定。
 |
| 图 1 模拟月壤SEM照片[19] Fig. 1 SEM image of lunar regolith simulant[19] |
| 图选项 |
 |
| 图 2 模拟月壤颗粒的简化模型 Fig. 2 Simplified models of lunar regolith simulant particles |
| 图选项 |
1.2 颗粒动力学模型 本文利用DEM分析模拟月壤铺粉中的流变特性时，暂不考虑颗粒二次破碎问题，认为每个颗粒都是均质刚体。由牛顿第二定律，每个颗粒的动力学方程可以表示为[21]
 | (1) |
 | (2) |
 | (3) |
 | (4) |
式中：mp为颗粒质量；xp为颗粒质心坐标；g为重力加速度；Fp, q[c]为颗粒间接触力；Fp, q[a]为颗粒间黏附力；Ip为颗粒转动惯量；ωp为颗粒角速度；Mp为颗粒的合力矩；ρp为颗粒密度；mi[b]、Vi[b]、xi[b]分别为颗粒内部第i个球单元的质量、体积、质心坐标；ΔV为颗粒内部所有球单元重叠部分的体积总和；Nsp为单个颗粒中的球单元总数。
基于刚体动力学理论，颗粒受到的作用力由构成该颗粒的所有球单元所受到的作用力确定。当颗粒-颗粒或颗粒-壁面发生接触时，实际上是颗粒中的某个球单元与相邻颗粒的某个球单元或壁面发生接触。因此，颗粒间相互接触作用可以通过球单元间接触作用模型来表征。目前，常用的颗粒细观接触作用模型有弹簧-阻尼模型、非线性Hertz-Mindlin和Deresiewicz模型等。Renzo和Maio[23]研究证明，当线性弹簧-阻尼接触作用模型参数取值足够准确时，其DEM计算结果与采用非线性Hertz-Mindlin和Deresiewicz接触作用模型获得的DEM结果具有相同的求解精度。考虑到模拟月壤铺粉过程的DEM数值仿真分析将会产生大量颗粒，导致颗粒接触搜索计算和动力学求解所占用的计算资源较大，因此，本文采用一种线性简化的弹簧-阻尼接触作用模型(见图 3)，以提高计算效率。接触作用模型参数将通过三轴压缩试验数据予以校准。
 |
| 图 3 球单元间的接触作用模型 Fig. 3 Contact reaction model between spherical elements |
| 图选项 |
如图 3所示，分属不同颗粒的2个球单元{p, q}发生相互接触作用时，将在其接触位置产生一个相互作用力(F[c])，其可以分解为法向作用力(Fn[c])和切向作用力(Fs[c])。基于线性弹簧-阻尼接触作用模型，球单元间法向作用力由法向弹簧刚度(kn)、法向阻尼系数(cn)及球单元间接触间隙(gs)确定；切向作用力由切向弹簧刚度(ks)、切向阻尼系数(cs)、库伦摩擦系数(μ)确定。相互接触的2个球单元的作用力计算公式如下[20]：
 | (5) |
 | (6) |
 | (7) |
式中：接触间隙gs=(Rp+Rq－||xp－xq||)，R为颗粒半径，x为颗粒质心位置；n为球单元间法向接触单位向量；
在铺粉过程中，模拟月壤颗粒的运动行为不仅与颗粒间接触作用力有关，还受到颗粒间黏附力影响。基于Hamaker理论和Lennard-Jones势，刚性球颗粒间的范德华力可以表示为[24]
 | (8) |
式中：w为黏附功，通常用表面能Δγ表示，Δγ=γp+γq-γpq；z0为两表面间范德华作用平衡间距；z为颗粒间最近点距离；Reff为等效颗粒半径，Reff=(Rp－1+Rq－1)－1。
Fvdw为正值表示斥力，反之为引力。当z=z0时，Fvdw=－2πwReff为最大引力值，被称为拉开力(Pull-Off Force)。另外，对于相同材料颗粒，γp=γq=γ，γpq=0，w=2γ，而γ可以由Hamaker常数计算获得[25]：
 | (9) |
图 4给出了不同粒径月壤颗粒间范德华力随z的变化曲线。其中，z0取值为0.165nm[26]；A取值为4.3×10-20J[27]。如图 4所示，颗粒间范德华力属于近程力，对颗粒间距z非常敏感，随着颗粒间距z的增加，||Fvdw||迅速衰减；同时，||Fvdw||随着颗粒粒径增加呈线性递增。
 |
| 图 4 月壤颗粒间范德华力曲线 Fig. 4 Van der Waals force curves for lunar regolith particles |
| 图选项 |
为了定量表征颗粒间范德华力黏附作用的强弱，本文引入一个无量纲参数——黏结数K，用拉开力(Fpulloff)与重力(Fg)的比值表示[16]，若设2个球颗粒半径相等，则K值可以表示为
 | (10) |
由式(10)，K与颗粒半径Rp的平方呈反比，与颗粒材料黏附功w呈正比，与重力加速度g呈反比。图 5显示了在不同重力场环境下的模拟月壤颗粒的K值随颗粒半径变化曲线。可见，随着颗粒半径减小，颗粒间范德华最大引力与颗粒重力的比值(K值)显著增加，在此情况下，范德华力对颗粒团聚、沉积行为有显著影响，在DEM数值计算中，需要考虑范德华力。此外，由于月面重力加速度仅为地球的1/6，而K值与重力加速度呈反比，因此，月基铺粉时，范德华力影响会增强。
 |
| 图 5 不同重力下月壤颗粒K值曲线 Fig. 5 K curves for lunar regolith particles under lunar and earth gravity |
| 图选项 |
1.3 DEM方法 DEM是一种从离散介质力学理论发展而来的数值方法，其认为介质是由众多离散的细观尺度颗粒(或单元)构成，每个颗粒运动行为服从牛顿第二定律和力-位移准则，介质的宏观变形和演化由颗粒运动和相互位置来描述。DEM采用显式中心差分法求解获得每个单元的运动信息及在运动过程中与外界的相互作用，其是以时间步长为间隔向前推进的计算过程，颗粒的运动定律和力-位移准则交替进行，运动定律计算由不平衡力引起的颗粒运动，力-位移准则计算由颗粒位移引起的接触力。
具体方法如下：利用显式中心差分法求解颗粒动力学方程，计算颗粒在t时刻的平动加速度和转动加速度：
 | (11) |
由式(11)计算出平动加速度和转动加速度，获得(t±Δt/2)时刻的平动速度和转动速度：
 | (12) |
由式(12)计算出平动速度和转动速度，更新颗粒的质心坐标：
 | (13) |
在DEM求解颗粒动力学方程时，Δt取值应确保在一个加载步内，任意一个颗粒的扰动只能传播至最临近的颗粒(不能传播到更远的距离)，颗粒速度和加速度在Δt内保持恒定。另外，显式中心差分法是条件稳定的，Δt应小于一个临界值Δtc，Δtc可由式(14)计算：
 | (14) |
式中：m为颗粒质量；k为弹簧刚度；Δt根据颗粒系统中最小Δtc确定。
2 铺粉数值模拟方法 2.1 铺粉装置建模 本文以德国Solution SLM125激光选区熔化成形设备的铺粉装置为研究对象，根据其工作原理和结构特征，采用刚性壳体单元，建立了包含刮刀组件、送粉盒组件、成形腔组件的三维简化模型，如图 6所示。其中，刮刀组件用于刮平和压实沉积在基板表面的粉末层，为SLM提供高质量的粉床，其由2个圆柱辊和1个圆柱形集粉槽组成；送粉盒组件用于储存粉末颗粒，并在铺粉过程中将粉末颗粒连续不断地输送到基板表面，其由1个送粉盒、输送通道及底部开槽(送粉槽)组成；成形腔组件用于形成粉床和进行SLM成形加工，其由1个基板和成形腔固体壁面组成，通过调节基板高度可以控制单次铺粉的厚度。
 |
| 图 6 铺粉装置全尺度简化模型 Fig. 6 Simplified full-scale model of powder spreading device |
| 图选项 |
SLM125铺粉装置的工作流程为：在单次铺粉开始时刻，成形腔中的基板根据预设的铺粉厚度，沿垂直方向(-Z)下移一定距离；然后，电机驱动刮刀组件和送粉盒组件沿铺粉方向(+X)按照预设的铺粉速度匀速运动，直至刮刀组件完全移出铺粉区域；在铺粉过程中，送粉盒中的粉末颗粒在重力和接触力的作用下，沿着输送通道下落，最后从底部的送粉槽下落到基板表面，随着颗粒连续不断地下落至基板表面，一个沉积粉层逐渐形成；与此同时，安装在送粉槽之后的2个刮刀匀速运动到沉积粉层区域，通过与粉层颗粒的接触作用，对其进行刮平和压实处理，在单次铺粉过程结束后形成一层高质量粉床。SLM125铺粉装置的主要特征尺寸为：基板长度和宽度均为125mm；圆柱刮刀直径(Dblader)为3.5mm；送粉槽长度(Lslot)为3mm、宽度(Wslot)为125mm。
考虑到采用全尺寸模型(见图 6)进行模拟月壤铺粉过程DEM数值模拟时，将产生一个数量庞大的颗粒群系统(颗粒数量达到百万量级)，导致计算周期和计算成本显著增加。为了提高计算效率，节约计算资源，本文提出一种基于周期性边界对称的子模型简化方法。具体思路是：在全局笛卡儿坐标系框架下截取全尺寸模型(见图 6)的1个子计算区域(Ω1)，并在与Y轴正交的Ω1表面设置周期性对称边界条件，在其余4个Ω1表面设置删除边界条件，如图 7所示。在DEM模拟铺粉过程时，如果某个颗粒的质心位置移出“周期性对称边界”，则会立即在另一侧边界自动生成新颗粒；如果某个颗粒的质心位置移出“删除边界”，则会被立刻删除，不再参与后续的DEM计算。此外，本文对送粉盒组件做进一步简化，删除送粉盒和输送通道模型部分，仅保留送粉槽模型。如图 7所示，送粉槽简化为一个底部开槽的矩形容器，在铺粉仿真中，模拟月壤颗粒会自动在送粉槽顶部生成新颗粒流，模拟颗粒从送粉盒经输送通道到达送粉槽的物理现象，颗粒产生通过一个用字子程序实现。通过上述模型简化，可以有效降低铺粉仿真模型中的颗粒数量，显著提高计算效率，同时能准确反映铺粉过程中颗粒的运动行为，使仿真结果具有合理精度。
 |
| 图 7 铺粉仿真DEM模型 Fig. 7 DEM model for powder spreading simulation |
| 图选项 |
2.2 模拟月壤DEM模型参数校准 本文在铺粉过程DEM数值模拟中，假设模拟月壤由4种非球形颗粒(见图 2)构成，并且每种颗粒形状所占比例为25%；颗粒粒径范围取70~120μm，粒度分布与月壤粒径级配曲线相符。由于模拟月壤颗粒的细观模型参数不能直接通过实验确定，因此采用Huang等[21]提出的基于相似理论和三轴压缩实验数值模拟相结合的DEM模型参数校准方法，即建立三轴压缩实验的DEM三维模型(见图 8)，依据三轴压缩实验工况定义载荷和边界条件，通过DEM数值模拟获得模拟月壤颗粒的本构关系曲线并与实验数据对比，再根据基于相似理论推导的Π关系式，对DEM模型参数进行校准，最终获得能够准确表征模拟月壤本构关系的细观参数取值。校准后的模拟月壤DEM模型参数如表 1所示。
 |
| 图 8 三轴压缩实验的DEM模型[21] Fig. 8 DEM model for triaxial compression test[21] |
| 图选项 |
表 1 模拟月壤DEM模型参数 Table 1 DEM model parameters for lunar regolith simulant
| 参数 | 数值 |
| 颗粒密度ρp/(kg·m-3) | 2770 |
| 法向弹簧刚度kn/(N·m-1) | 1×105 |
| 切向弹簧刚度ks/(N·m-1) | 1×105 |
| 法向阻尼系数cn | 0.7 |
| 切向阻尼系数cs | 0.7 |
| 摩擦系数μp | 0.3 |
表选项
2.3 粉床质量评价方法 本文利用堆积密度、平均配位数、表面粗糙度这3个指标，对粉床质量进行定量评价。高质量的粉床应满足：堆积密度和平均配位数较大，表面粗糙度较小。相反，堆积密度和平均配位数计算值越小、表面粗糙度计算值越大，则表示粉床质量越差，这种粉床不利于获得高质量的成形件。
堆积密度是一个宏观参数，定义为固体颗粒体积与总体积之比；颗粒配位数是一个细观参数，定义为单个颗粒与周围颗粒或固壁的接触总数；平均配位数是一个统计学参数，定义为颗粒系统中每个颗粒的配位数之和除以颗粒总数。堆积密度(ρpd)和平均配位数(Nc)的数学表达式为
 | (15) |
 | (16) |
式中：Nci为第i个颗粒的配位数；Np为颗粒总数；Vi为第i个颗粒的体积；V为总体积。
在DEM数值仿真中，单次铺粉形成的粉床，其堆积密度和平均配位数可以通过环刀测试方法计算获得。具体思路是：在粉床区域建立若干个圆柱形环刀取样模型(见图 9)，利用一个用户子程序分别计算每个环刀取样模型内的颗粒总数量、颗粒体积、颗粒接触对总数量、环刀封闭体积，再由式(15)和式(16)求出相应的堆积密度和平均配位数，并将计算结果写入一个列表中；最后，粉床堆积密度和平均配位数定义为所有环刀取样模型计算结果的算数平均值。
 |
| 图 9 粉床上环刀取样模型示意图 Fig. 9 Schematic of cutting-ring sampling model for powder bed |
| 图选项 |
粉床表面粗糙度对成形件质量有重要影响，本文采用轮廓算术平均偏差(Ra)对粉床表面粗糙度进行定量描述。为了获得粉床轮廓上各点高度，此处采用一种射线跟踪技术(Ray-Tracing Technique)[18]，具体方法是：首先，在粉床指定高度生成一组点光源阵列，每个点光源产生一条沿-Z方向的射线；然后，搜寻射线与粉床颗粒表面的相交点，并将该相交点的坐标记录到一个列表中，并标记为一个“轮廓点高度”；最后，利用列表数据和Ra计算公式计算出粉床表面粗糙度。
 | (17) |
式中：Nr为点光源的总数；zi为轮廓点与基准线之间的距离。
3 DEM仿真结果 3.1 地基铺粉过程DEM仿真分析 本文采用三维DEM和SLM铺粉装置简化模型(见图 6)，开展模拟月壤在地球重力环境下的铺粉过程数值模拟研究。为简化计算，铺粉装置(刮刀组件+送粉槽+成形腔)均定义为刚体，采用三维SHELL单元建模；刮刀组件和送粉槽预设恒定的水平速度vx=0.19m/s，即为铺粉速度；成形腔基板预设垂直方向位移δz=-200μm，即为粉床厚度；出粉口与基板垂直距离Δh=3δz；计算区域Ω1宽度W1=1.2mm，长度L1=33.5mm，高度H1= 10.0mm；重力加速度g=9.81m/s2，作用在所有颗粒的质心。
DEM计算结果如下：
接触颗粒组成分子之间的范德华力累积产生颗粒间黏附作用力(颗粒间范德华力)，其对于大颗粒可以忽略，但对小颗粒的接触作用有重要影响。本文通过式(8)将颗粒间范德华力引入颗粒动力学模型中，并考察颗粒间范德华力对地基模拟月壤铺粉形成的粉床质量的影响。如图 10(a)所示，对于地基单次铺粉情况，当忽略颗粒间范德华力作用时，模拟月壤粉材具有较好的流动性，其沉积行为主要由颗粒自身重力控制，颗粒间不发生团聚现象，颗粒-固壁之间物理接触黏结强度为零，非球颗粒间互锁作用、颗粒间摩擦作用时影响颗粒配位数和堆积密度的主要因素。从DEM计算结果可知，在没有范德华力作用时，铺粉形成的模拟月壤粉床内部结构紧密、均匀且连续，没有明显的空穴缺陷，粉床表面粗糙度小。然而，范德华力具有固有性，其在接触颗粒间总是存在。在细颗粒流中，范德华黏附力是影响颗粒间团聚、颗粒-固壁之间形成物理接触的黏附强度的主要因素，对颗粒接触后行为起主导作用。如图 10(b)所示，在范德华力作用下，铺粉形成的模拟月壤粉床内部结构更疏松、局部存在空穴缺陷，造成粉床整体堆积密度、平均配位数降低，表面粗糙度变大。
 |
| 图 10 粉床DEM数值仿真结果 Fig. 10 Powder bed simulated by DEM |
| 图选项 |
采用射线跟踪法，在粉床中心线(图 10黄色虚线)位置上方，沿铺粉方向(X轴)等距生成500个点光源，获得粉床中心线位置的表面形貌，如图 11所示。其中，X轴为经归一化处理的无量纲标量，Y轴为粉床表面上各点与基板之间的相对高度值。如图 11所示，当引入颗粒间范德华力后，粉床表面形貌变得更粗糙，表面粗糙度Ra由10.686μm增加到30.461μm。
 |
| 图 11 粉床表面形貌曲线(Y=0) Fig. 11 Rough-surface profile of powder bed(Y=0) |
| 图选项 |
采用环刀测试方法在粉床上均匀布置9个圆环(见图 9)，通过式(15)和式(16)计算出每个圆环内粉材的堆积密度和平均配位数，结果如图 12和图 13所示。可见，当不考虑颗粒间范德华黏附力时，整个粉床的堆积密度和颗粒配位数的平均值分别为0.588、5.456；当引入颗粒间范德华力时，粉床整体堆积密度和颗粒配位数的平均值均有所降低，分别为0.463、4.993。
 |
| 图 12 粉床堆积密度分布柱状图 Fig. 12 Packing density distribution histogram of powder bed |
| 图选项 |
 |
| 图 13 粉床平均配位数分布柱状图 Fig. 13 Average coordination number distribution histogram of powder bed |
| 图选项 |
从地基铺粉DEM计算结果可知，对于粒径范围在70~120μm的模拟月壤粉末材料，颗粒间范德华力对铺粉形成的粉床质量具有显著影响，其会使粉体颗粒间产生团聚现象，同时增加颗粒与铺粉装置固壁之间形成物理接触的黏附强度，降低粉体材料的流动性，导致所形成的粉床内部结构不均匀、不连续、局部出现空穴缺陷，宏观上表现为粉床的堆积密度小、平均配位数低、表面粗糙度大，粉床质量不佳。因此，对于粒径在70~120μm范围的模拟月壤颗粒系统，颗粒间范德华力在颗粒动力学模型中不能被省略。
3.2 月基铺粉过程DEM仿真分析 针对SBAM-ISRU技术发展需求，在地基铺粉DEM数值仿真基础上，本文进一步开展月基铺粉细观动力学研究。在DEM数值计算中，采用与月基铺粉研究相同的模型参数及铺粉工艺参数，仅调整重力加速度g=1.63m/s2。图 14显示了月基铺粉不同时刻模拟月壤颗粒系统的速度矢量图、范德华力链图、接触力链图。可见，当送粉槽以vx=0.19m/s匀速水平运动时，槽内颗粒在自身重力、颗粒间范德华力、颗粒间接触力共同作用下逐渐沉积到基板表面，颗粒速度主要分成3个区域：在送粉槽内左侧区域，颗粒速度矢量方向近似为水平方向，颗粒速度最大值出现在颗粒-壁面接触界面附近，数值上约等于铺粉速度；在送粉槽内右侧区域，颗粒速度矢量方向近似为垂直方向；在送粉槽出口与基板之间区域，颗粒速度矢量方向与水平存在一个角度。从颗粒速度矢量图可知，处于送粉槽内左侧区域的颗粒不平衡力主要由接触力和范德华黏附力控制，重力影响很小；处于送粉槽内右侧区域的颗粒不平衡力主要由自身重力控制，接触力和范德华黏附力作用较小；处于基板与出粉口之间的颗粒不平衡力由颗粒-颗粒(颗粒-壁面)之间接触力、范德华力、重力共同作用。
 |
| 图 14 月基铺粉DEM仿真结果(g=1.63m/s2) Fig. 14 Lunar-based powder-layering DEM simulation results (g=1.63m/s2) |
| 图选项 |
随着送粉槽水平移动，在基板上逐渐形成具有一定厚度的沉积粉层，但由于颗粒间范德华黏附力及非球颗粒间形成的互锁作用，粉层表面粗糙、颗粒排列不均匀。随后，双圆柱面刮刀逐渐运动到成形腔区域并与沉积粉层表面颗粒发生接触作用，通过刮刀圆柱面对沉积粉层的剪切和压实作用，提高粉层的堆积密度和平均配位数，同时降低粉层的表面粗糙度。
通过环刀测试方法计算获得粉床9个取样点堆积密度和颗粒配位数平均值分别为0.395、3.27；利用射线跟踪技术及式(17)计算出粉床表面粗糙度Ra=48.783μm。模拟月壤地基与月基铺粉DEM计算结果对比显示，月面低重力环境不利于形成高质量粉床，这是因为在月球环境下颗粒自重仅为地球的1/6，颗粒间接触力和范德华力对颗粒运动行为影响增强，范德华力使颗粒团聚同时使颗粒-固壁间形成物理黏结，接触力可能使颗粒产生与重力方向相反的不平衡力，导致铺粉过程中发生扬尘现象。综上所述，月基铺粉形成的粉床具有更低的堆积密度和颗粒配位数，表面也更粗糙，总体质量较差，这将对后续的SLM成形产生不利影响。
3.3 月基铺粉工艺优化方法 从SBAM-ISRU技术应用角度出发，本文提出月基铺粉工艺优化方法研究。传统铺粉工艺研究发现，降低铺粉速度将有利于提高粉床质量，因此，本文分别选取vx=0.05m/s、0.10m/s、0.19m/s开展月基铺粉DEM数值仿真研究，考察铺粉速度对月基铺粉粉床状态的影响规律，并保持其他参数不变。图 15和图 16显示了沉积粉层堆积密度、平均配位数与铺粉速度的关系。可见，铺粉速度越高，沉积粉层的堆积密度和平均配位数越小，说明粉层内部结构疏松、不均匀、不连续。基于环刀测试方法计算出上述3种铺粉速度时粉床整体堆积密度平均值分别为0.473、0.394、0.339；颗粒配位数平均值分别为4.87、3.28、2.495；基于射线跟踪法计算出表面粗糙度Ra分别为30.182、47.463、71.936。图 17和图 18显示了双圆柱面刮刀匀速通过成形腔并完成对沉积粉层的压实和剪切后所形成粉床的堆积密度和平均配位数。同样地，计算出3种铺粉速度时粉床堆积密度平均值分别为0.384、0.390、0.395；配位数平均值分别为4.427、3.945、3.27；表面粗糙度Ra分别为25.464、39.226、48.783。计算结果显示，在其他条件相同情况下，降低铺粉速度可以改善沉积粉层和粉床质量，因此月基铺粉建议选择较低水平速度。
 |
| 图 15 沉积粉层堆积密度分布柱状图 Fig. 15 Packing density distribution histogram of deposited powder layer |
| 图选项 |
 |
| 图 16 沉积粉层平均配位数分布柱状图 Fig. 16 Average coordination number distribution histogram of deposited powder layer |
| 图选项 |
 |
| 图 17 粉床堆积密度分布柱状图 Fig. 17 Packing density distribution histogram of powder bed |
| 图选项 |
 |
| 图 18 粉床平均配位数分布柱状图 Fig. 18 Average coordination number distribution histogram of powder bed |
| 图选项 |
如上所述，降低铺粉速度可以有效减小粉床表面粗糙度、增加颗粒配位数，从而改善月基铺粉的粉床质量，但伴随的不利后果是生产效率降低。从兼顾生产效率和铺粉粉床质量角度出发，本文提出刮刀型面优化设计研究。借鉴Haeri[19]研究成果，本文设计了一种新刮刀结构，刮刀底部型面
设计为超椭圆轮廓，数学表达式为
 | (18) |
式中：超椭圆方程参数取值分别为：a=6Dblader，b=Dblader，n=5.0。与双圆柱刮刀型面相比，超椭圆刮刀型面有效增加了刮刀与粉层颗粒接触面积，减小了铺粉过程中刮刀对粉层颗粒的扰动，避免或降低剪胀作用出现。2种刮刀的月基铺粉DEM仿真结果显示，与双圆柱型面刮刀相比，超椭圆型面刮刀可以有效提高月基铺粉所形成的粉床质量，在相同铺粉工艺条件下，粉床整体堆积密度由0.395增加到0.413，颗粒平均配位数由3.27提高到3.81，表面粗糙度没有明显变化。
4 结论 本文针对基于SBAM-ISRU的月面原位成形技术实验成本过高、地面环境模拟技术难度大等问题，提出一种基于叠加球模型和三维DEM仿真技术的模拟月壤铺粉过程细观数值模拟方法，通过构建合理的数理模型和引入空间环境因素，模拟不同工作条件下粉床形成过程和颗粒细观动态行为，为该技术的理论研究和工程技术发展提供必要的研究工具。结论如下：
1) 采用叠加球模型和三维DEM数值方法可以模拟不同环境和工况条件下的模拟月壤铺粉过程，动态显示铺粉过程中每个颗粒的运动和受力，该仿真方法在处理颗粒运动的随机性、不连续性方面具有优势，可作为实验研究的必要补充。
2) 月面低重力环境对模拟月壤铺粉质量有重要影响，月基铺粉形成的粉床质量相对较差，并且在铺粉过程中会发生扬尘现象。在月面低重力环境下，颗粒间范德华黏附力与非球形颗粒间互锁作用是导致模拟月壤粉材流动性变差的主要原因，范德华力主要引起颗粒间团聚和颗粒与固壁之间产生物理黏结，颗粒间互锁作用相当于增加了颗粒间摩擦力。
3) 铺粉速度是影响月基铺粉沉积粉层状态的重要工艺参数，随着铺粉速度的增加，沉积粉层的堆积密度、颗粒配位数平均值减小，表面粗糙度Ra增加；双圆柱型面刮刀通过压实和剪切沉积粉层可以有效提高堆积密度、颗粒配位数，减小表面粗糙度，铺粉速度越小，粉床的颗粒配位数平均值越大，表面粗糙度越小。因此，在月基铺粉操作时，建议尽量选择较低的铺粉速度，但同时需要兼顾生产效率。
4) 从兼顾生产效率和粉床质量角度出发，本文提出一种超椭圆型面刮刀新结构，通过增加刮刀与粉层颗粒接触面积减小对粉层颗粒的扰动，避免或降低剪胀作用。DEM结果显示，在相同月基铺粉条件下，采用超椭圆型面刮刀获得的粉床状态优于双圆柱型面形成的粉床状态，粉床整体堆积密度由0.395增加到0.413，颗粒平均配位数由3.27提高到3.81，表面粗糙度没有明显变化。建议开展刮刀型面优化设计研究，进一步提高粉床质量，为月基SLM应用技术发展奠定基础。
参考文献
| [1] | SANDERS G B, LARSON W E, PICARD M.Development and demonstration of sustainable surface infrastructure for Moon/Mars exploration[C]//62nd International Astronautical Congress, 2011: JSC-CN-24659. |
| [2] | FATERI M, GEBHARDT A.Experimental investigation of selective laser melting of lunar regolith for in-situ applications[C]//International Mechanical Engineering Congress and Exposition.New York: ASME, 2013: IMECE2013-64334. |
| [3] | HAFLEY R A.Electron beam free-form fabrication in the space environment[C]//Aerospace Sciences Meeting and Exhibit.Reston: AIAA, 2007: 8-11. |
| [4] | PRATER T, WERKHEISER N, LEDBETTER F. 3D printing in zero G technology demonstration mission:Complete experimental results and summary of related material modeling efforts[J]. International Journal of Advanced Manufacturing Technology, 2019, 101(1-4): 391-417. DOI:10.1007/s00170-018-2827-7 |
| [5] | WERKHEISER N.In-space manufacturing: 3D printing in space technology demonstration[C]//National Space & Missile Materials Symposium, 2015: M15-4685. |
| [6] | CECCANTI F, DINI E, KESTELIER X D, et al.3D printing technology for a Moon outpost exploiting lunar soil[C]//61st International Astronautical Congress, 2010: 8812-8820. |
| [7] | CESARETTI G, DINI E, DE KESTELIER X, et al. Building components for an outpost on the lunar soil by means of a novel 3D printing technology[J]. Acta Astronautica, 2014, 93: 430-450. DOI:10.1016/j.actaastro.2013.07.034 |
| [8] | BALLA V K, ROBERSON L B, OCONNOR G W, et al. First demonstration on direct laser fabrication of lunar regolith parts[J]. Rapid Prototyping Journal, 2012, 18(6): 451-457. DOI:10.1108/13552541211271992 |
| [9] | FATERI M, GEBHARDT A. Process parameters development of selective laser melting of lunar regolith for on-site manufacturing applications[J]. International Journal of Applied Ceramic Technology, 2015, 12(1): 46-52. DOI:10.1111/ijac.12326 |
| [10] | GOULAS A, BINNER J G P, HARRIS R A, et al. Assessing extraterrestrial regolith material simulants for in-situ resource utilization based 3D printing[J]. Applied Materials Today, 2017, 6: 54-61. DOI:10.1016/j.apmt.2016.11.004 |
| [11] | GERDESL N, FOKKEN G, LINKE S, et al. Selective laser melting for processing of regolith in support of a lunar base[J]. Journal of Laser Applications, 2018, 30(3): 032018. DOI:10.2351/1.5018576 |
| [12] | GU D, MEINERS W, WISSENBACH K, et al. Laser additive manufacturing of metallic components:Materials, processes and mechanisms[J]. International Materials Reviews, 2012, 57(3): 133-164. DOI:10.1179/1743280411Y.0000000014 |
| [13] | ZIEGELMEIER S, CHRISTOU P, WOLLECKE F, et al. An experimental study into the effects of bulk and flow behaviour of laser sintering polymer powders on resulting part properties[J]. Journal of Materials Processing Technology, 2015, 215(1): 239-250. |
| [14] | ZOHDI T I. Rapid simulation of laser processing of discrete particulate materials[J]. Archives of Computational Methods in Engineering, 2013, 20(4): 309-325. DOI:10.1007/s11831-013-9092-6 |
| [15] | XIANG Z W, YIN M, DENG Z, et al. Simulation of forming process of powder bed for additive manufacturing[J]. Journal of Manufacturing Science and Engineering, 2016, 138(8): 081002-1-081002-9. |
| [16] | CHEN H, WEI Q, WEN S, et al. Flow behavior of powder particles in layering process of selective laser melting:Numerical modeling and experimental verification based on discrete element method[J]. International Journal of Machine Tools and Manufacture, 2017, 123: 146-159. DOI:10.1016/j.ijmachtools.2017.08.004 |
| [17] | PARTELI E J, POSCHEL T. Particle-based simulation of powder application in additive manufacturing[J]. Powder Technology, 2016, 288: 96-102. DOI:10.1016/j.powtec.2015.10.035 |
| [18] | HAERI S, WANG Y, GHITA O, et al. Discrete element simulation and experimental study of powder spreading process in additive manufacturing[J]. Powder Technology, 2017, 306: 45-54. DOI:10.1016/j.powtec.2016.11.002 |
| [19] | HAERI S. Optimisation of blade type spreaders for powder bed preparation in additive manufacturing using DEM simulations[J]. Powder Technology, 2017, 321: 94-104. DOI:10.1016/j.powtec.2017.08.011 |
| [20] | LI W, HUANG Y, CUI Y, et al. Trafficability analysis of lunar mare terrain by means of the discrete element method for wheeled rover locomotion[J]. Journal of Terramechanics, 2010, 47(3): 161-172. DOI:10.1016/j.jterra.2009.09.002 |
| [21] | HUANG Y, LU X, ZHAO R, et al. Three dimensional simulation of lunar dust levitation under the effect of simulated sphere body[J]. Journal of Terramechanics, 2011, 48(4): 297-306. DOI:10.1016/j.jterra.2011.06.005 |
| [22] | HUANG Y, ZHAO R, LI W, et al. Radiative characteristics of nonspherical particles based on a particle superposition model[J]. Journal of Geophysical Research, 2013, 118(20): 762-769. |
| [23] | RENZO A D, MAIO F P. Comparison of contact-force models for the simulation of collisions in DEM-based granular flow codes[J]. Chemical Engineering Science, 2004, 59(3): 525-541. DOI:10.1016/j.ces.2003.09.037 |
| [24] | JOHNSON K L, GREENWOOD J A. An adhesion map for the contact of elastic spheres[J]. Journal of Colloid and Interface Science, 1997, 192(2): 326-333. DOI:10.1006/jcis.1997.4984 |
| [25] | GOTZINGER M, PEUKERT W. Dispersive forces of particle-surface interactions:Direct AFM measurements and modelling[J]. Powder Technology, 2003, 130(1-3): 102-109. DOI:10.1016/S0032-5910(02)00234-6 |
| [26] | ISRAELACHVILI J N. Intermolecular and surface forces[M]. 2nd ed.Salt Lake City: Academic Press, 1992. |
| [27] | POPEL S I, ZELENYI L M, GOLUB A P, et al. Lunar dust and dusty plasmas:Recent developments, advances, and unsolved problems[J]. Planetary and Space Science, 2018, 156: 71-84. DOI:10.1016/j.pss.2018.02.010 |
