论文部分内容阅读
【摘 要】在探月工程中，月球探测器的软着陆过程中必须考虑到缓冲发动机羽流激起月壤尘埃的问题，因为扬起的月壤很可能对探测器设备产生多种不良的影响。本文以缓冲发动机羽流与月壤所形成的气固两相流场为主要研究对象，介绍了如何应用FLUENT-EDEM两相耦合模块对其进行模拟仿真。以一个简化模型作为实例，阐述了该耦合方法相较于以往的两相流仿真所具有的优势。
【关键词】两相流;月壤;耦合;EDEM
引言
实现月球探测器的软着陆是我国探月工程第二阶段的关键技术，而火箭发动机反向喷气是一种常见的着陆缓冲技术。但在探测器软着陆过程中缓冲发动机羽流激起的月壤可能会衰减导航敏感器件的探测信号，甚至导致仪器无法正常工作，激起的月尘还会造成的视觉模糊，严重影响探测器的安全及准确着落，因此探月过程中必须考虑羽流中的月壤对探测器所能造成的影响。对气固两相流的研究于探月工程具有不可忽视的积极意义，本文将要讨论CFD-EDEM耦合方法在气固两相流研究中的应用。
1 概述
气固两相流的理论分析相较于单相流的分析困难得多，而描述两相流的通用微分方程组至今尚未建立。大量理论工作主要采用以下两类简化模型[1]：
1）均相模型。将两相介质看成是一种混合得非常均匀的混合物，假定处理单相流动的概念和方法仍然适用于两相流，但需对它的物理性质及传递性质作合理的假定;
2）分相模型。认为单相流的概念和方法可分别用于两相系统的各个相，同时考虑两相之间的相互作用。
目前，对如何处理气固两相流在数值计算上主要采用四种方法：一是将某相看成连续的，根据连续理论导出欧拉基本方程，称为欧拉法;二是将某相视为不连续的离散相，对每个质点进行拉格朗日追踪，称为拉格朗日法。将上述的这两种方法综合起来，气固两相流的研究方法有欧拉-欧拉、欧拉-拉格朗日、拉格朗日-拉格朗日三种方法。第四种方法就是基于分子运动理论的直接模擬蒙特卡洛法与计算流体力学方法相结合对气固两相流数值计算的方法。
而本文所要探讨的CFD-EDEM耦合方法相较于上述四种方法，则具有更加独到的优势，我们将在下文作展开说明。
2 CFD-EDEM耦合方法
2.1 EDEM软件简介
仿真软件EDEM的主要功能是仿真、分析和观察颗粒流的运动规律，它是全球第一个被专门设计为使用离散元模型技术作为基础的通用CAE软件。
EDEM的Creater（前处理模块）可以快速、简便的为散货颗粒进行参数化建模。
Simulator（求解器）基于先进的离散元方法，结合了经典的碰撞模型，在拉格朗日框架下，离散求解系统的运动学方程、动力学方程和本构方程。从而获得每个颗粒及壁面单元体的速度、受力等参数。
Analyst是EDEM的后处理模块，含有很多工具对计算结果进行加工处理，方便地获得颗粒的速度、能量分布，设备的受力情况等，并可利用统计方法对各物理量间的相关性进行分析。
2.2 Fluent-EDEM耦合方案
EDEM不仅可以单独进行颗粒力学模拟，也可以和通用CFD软件Fluent耦合进行流体-颗粒系统的模拟。其主要特点如下：
1）EDEM-Fluent耦合计算时，EDEM的耦合模块内置于Fluent软件中，两者无缝连接;2）EDEM能够自动读取Fluent中的网格;3）EDEM可以使用Fluent中离散相模型和欧拉多项流模型进行多相流耦合;4）采用EDEM直接对颗粒运动进行模拟;5）完全双向动力耦合。
Fluent-EDEM可以模拟绝大多数的流体-颗粒系统，特别是在以下方面，更是具有不可代替的优势：
?颗粒尺度分布范围广;非球形颗粒;流体中充满颗粒，颗粒的体积分数大;颗粒凝聚和分离;机械设备上颗粒沉积和黏结;可以直观化地观察计算过程中颗粒的运动情况。
Fluent-EDEM耦合模拟流程如下：先由Fluent将某一时间点的流场计算收敛，将流场信息通过曳力模型转化为EDEM中作用在颗粒上的流体曳力，EDEM计算每个颗粒所受的外力，并由此更新颗粒的位置、速度等信息，最后这些颗粒属性以动量汇的形式加到CFD计算中，从而影响流场。
2.3 CFD-EDEM耦合理论
2.3.1 CFD-EDEM升力模型
2.3.1.1 Saffman升力模型
在边界层，剪切层内，或者是流体发生剪切时，会出现较高的速度梯度。当颗粒处于高速流体当中，较高的速度梯度造成了颗粒表面的压力差，从而产生了升力。Saffman升力可用下式表达：
2.3.1.2 Magnus升力模型
Magnus升力是由于颗粒自旋而产生的升力，球体的阻力系数与雷诺数较高的情况下的Magnus升力的关联可以表示为下式：
这个模型即使对于雷诺数高达2000时对于颗粒仍然适用 [2]。
2.3.1.3 Fluid-Induced Torque（流激诱导转矩）
在以上耦合的模型中，没有考虑由颗粒自旋产生的流体剪切力的影响。显然，若在高速流体中这个力不能忽略。根据Dennis等人的研究[3]，有以下公式：
通过这样的处理，这一扩展模型可用于颗粒雷诺数在32到1000的情况。
2.3.2 CFD-EDEM阻力模型
Fluent-EDEM耦合采用一种改进的自由流阻力来计算作用在颗粒球形上的力，对于颗粒表面为非球面的情况，先用对颗粒的球型包围盒进行阻力的计算。这种处理方法仅适合处理两种境况下的颗粒：（1）那些大小相同或者是体积小于流体网格单元的颗粒;（2）周围的流体参数（速度、密度和粘度等）的变化基本为常量的颗粒。阻系数CD，取决于雷诺数Re： 式中：r 是流体密度，h 为流体的粘度，L 是颗粒球的直径，v 是颗粒与流体间的相对速度，a 是CFD网格单元的自由体积。除了阻力，颗粒的固有浮力也必须考虑。浮力的计算公式是：
一般对于Lagrangian耦合来说，这种改进的阻力法则可以用来考虑颗粒载荷的影响;而在Eulerian耦合中，颗粒载荷需通过结合体积分数的方法来考虑这一影响。
3 FLUENT-EDEM耦合实例
3.1 计算模型建立
气固两相流数值计算过程中，颗粒在气流中场中的运动是从静止开始计算的。为了简化计算，同时使流场设置不至于过大，分析探测器发动机距月表较近时（发动机喷管出口垂直月表向下）羽流与月面月尘（月壤）的相互作用。
据此，将模型的参数设置为：基本流场为以5mm为半径的底面、15mm为高的圆柱，圆柱的上底面作为喷流的入口，下底面作为月面。月壤颗粒密度设为3100
，泊松比为0.42，弹性模量为10MPa，环境压强为0.01Pa;喷管的出口（亦即流场的入口）速度为1000m/s（因为喷管出口压力在0.02MPa左右，而环境压强仅为0.01Pa，喷流在喷管内处为极度膨胀不足的状态，所以达到这样大的速度也是合理的）。重力加速度为1.62
。
边界条件设置如下：流场入口为压力入口边界条件;出口则为压力出口边界条件;在固体壁面上，气相流场采用无滑移固壁边界条件，颗粒与壁面之间采用Hertz Mindlin（No Slip）碰撞模型。
3.2 数值计算方法
采用FLUENT与EDEM軟件耦合求解气固两相流场，气流流场采用FLUENT软件，选用不可压缩流体、
紊流模型计算;在EDEM软件中采用Lagrangian坐标方法分析颗粒在弯管中运动和分布规律。首先由FLUENT得到原始气流场结果，并将流场参数传递给EDEM;由EDEM中DEM模块计算在该流场下颗粒受力的大小，然后通过内部迭代计算颗粒下一时间步长运动参数，同时更新颗粒运动轨迹;最后将EDEM计算结果返回FLUENT进行下一时间步长迭代。
采用Lagrangian方法计算颗粒运动轨迹，由Rayleigh数控制迭代步长。Rayleigh数由下面公式给出：
式中，R是颗粒半径，ρ是颗粒密度，G是剪切模量，υ是泊松比。则颗粒运动的迭代时间步长为：
FLUENT时间步长应为EDEM软件时间步长的1/10～1/100，由此来确定流场的网格尺度，故FLUENT和EDEM网格不同、时间尺度不同[4]。
3.3计算结果及分析
经过耦合计算之后，流场显示收敛，保存计算结果。通过EDEM自带的Analyst模块，我们可以导出进行加工处理后的数据，再对其进行相关的分析讨论。
计算流场的速度云图如下所示：
月壤颗粒的各项特性随仿真时间的推进的变化图如下：
为了加快仿真时间，在EDEM中的颗粒工厂只加入了1000个月壤颗粒进行计算，由上图可知，仿真只进行到了0.3s左右，颗粒已经几乎全部被吹出了仿真区域，而且被吹起的月壤的平均速度大多能够达到二三十米每秒，说明激起的扬尘足以对探测器上所携带的仪器产生干扰，会对探测器的安全及准确降落造成影响，与Apollo15的录像测量数据相吻合。
计算完成后，我们可以创建图表来记录指定的时间步长或者一定时间段内的仿真情况。以下便是不同时间点下一定速度范围内各月壤颗粒数目统计的直方图，从图中我们可以非常直观地得到仿真过程中某一时间点下的不同速度的颗粒分布状况，从而对我们的后续研究及观察提供有利的参考依据。
以下为各个时间点下流场中不同位置的颗粒的数目分布图，据此我们可以得出流场内颗粒的分布密度，从而判定是否对探测器上的设备造成影响。
4 结论
相较于以往的两相流计算仿真，FLUENT-EDEM耦合方法不需要像使用UDF添加颗粒相需要进行编程那么麻烦，而且EDEM由于是采用表面网格来描述边界表面，使其可以与FLUENT的流体网格的边界表面元素的点对点耦合;加上耦合方案为双向动力耦合，令仿真效果更加真实可信。由于可以直接在软件中手动设置颗粒的相关参数，可以很方便地进行各项属性的更改，而且EDEM也提供了非球体颗粒的设置，可以让仿真更加贴近实际情况。此外，EDEM中还加入了先进的静电模型来仿真带电颗粒，这与月壤因太阳风的吹扫而带电的情况相吻合。
参考文献：
[1]倪晋仁，李振山.风沙两相流理论及其应用[M].北京：科学出版社，2006.
[2]Oesterlé，B.and Bui Dinh，T.，Experiments on the lift of a spinning sphere in a range of intermediate Reynolds numbers，Experiments in Fluids，25，pp.16-22，1998.
[3]Zhou，H.，Flamant，G.，Gauthier，D.，DEM-LES simulation of coal combustion in a bubbling fluidized bed.Part II：coal combustion at the particle level，Chemical Engineering Science，59 pp.4205-4215，2004.
[4]许峰，王彤，傅耀.弯管内小尺度颗粒运动的数值模拟与分析.四川兵工学报.2010.02
【关键词】两相流;月壤;耦合;EDEM
引言
实现月球探测器的软着陆是我国探月工程第二阶段的关键技术，而火箭发动机反向喷气是一种常见的着陆缓冲技术。但在探测器软着陆过程中缓冲发动机羽流激起的月壤可能会衰减导航敏感器件的探测信号，甚至导致仪器无法正常工作，激起的月尘还会造成的视觉模糊，严重影响探测器的安全及准确着落，因此探月过程中必须考虑羽流中的月壤对探测器所能造成的影响。对气固两相流的研究于探月工程具有不可忽视的积极意义，本文将要讨论CFD-EDEM耦合方法在气固两相流研究中的应用。
1 概述
气固两相流的理论分析相较于单相流的分析困难得多，而描述两相流的通用微分方程组至今尚未建立。大量理论工作主要采用以下两类简化模型[1]：
1）均相模型。将两相介质看成是一种混合得非常均匀的混合物，假定处理单相流动的概念和方法仍然适用于两相流，但需对它的物理性质及传递性质作合理的假定;
2）分相模型。认为单相流的概念和方法可分别用于两相系统的各个相，同时考虑两相之间的相互作用。
目前，对如何处理气固两相流在数值计算上主要采用四种方法：一是将某相看成连续的，根据连续理论导出欧拉基本方程，称为欧拉法;二是将某相视为不连续的离散相，对每个质点进行拉格朗日追踪，称为拉格朗日法。将上述的这两种方法综合起来，气固两相流的研究方法有欧拉-欧拉、欧拉-拉格朗日、拉格朗日-拉格朗日三种方法。第四种方法就是基于分子运动理论的直接模擬蒙特卡洛法与计算流体力学方法相结合对气固两相流数值计算的方法。
而本文所要探讨的CFD-EDEM耦合方法相较于上述四种方法，则具有更加独到的优势，我们将在下文作展开说明。
2 CFD-EDEM耦合方法
2.1 EDEM软件简介
仿真软件EDEM的主要功能是仿真、分析和观察颗粒流的运动规律，它是全球第一个被专门设计为使用离散元模型技术作为基础的通用CAE软件。
EDEM的Creater（前处理模块）可以快速、简便的为散货颗粒进行参数化建模。
Simulator（求解器）基于先进的离散元方法，结合了经典的碰撞模型，在拉格朗日框架下，离散求解系统的运动学方程、动力学方程和本构方程。从而获得每个颗粒及壁面单元体的速度、受力等参数。
Analyst是EDEM的后处理模块，含有很多工具对计算结果进行加工处理，方便地获得颗粒的速度、能量分布，设备的受力情况等，并可利用统计方法对各物理量间的相关性进行分析。
2.2 Fluent-EDEM耦合方案
EDEM不仅可以单独进行颗粒力学模拟，也可以和通用CFD软件Fluent耦合进行流体-颗粒系统的模拟。其主要特点如下：
1）EDEM-Fluent耦合计算时，EDEM的耦合模块内置于Fluent软件中，两者无缝连接;2）EDEM能够自动读取Fluent中的网格;3）EDEM可以使用Fluent中离散相模型和欧拉多项流模型进行多相流耦合;4）采用EDEM直接对颗粒运动进行模拟;5）完全双向动力耦合。
Fluent-EDEM可以模拟绝大多数的流体-颗粒系统，特别是在以下方面，更是具有不可代替的优势：
?颗粒尺度分布范围广;非球形颗粒;流体中充满颗粒，颗粒的体积分数大;颗粒凝聚和分离;机械设备上颗粒沉积和黏结;可以直观化地观察计算过程中颗粒的运动情况。
Fluent-EDEM耦合模拟流程如下：先由Fluent将某一时间点的流场计算收敛，将流场信息通过曳力模型转化为EDEM中作用在颗粒上的流体曳力，EDEM计算每个颗粒所受的外力，并由此更新颗粒的位置、速度等信息，最后这些颗粒属性以动量汇的形式加到CFD计算中，从而影响流场。
2.3 CFD-EDEM耦合理论
2.3.1 CFD-EDEM升力模型
2.3.1.1 Saffman升力模型
在边界层，剪切层内，或者是流体发生剪切时，会出现较高的速度梯度。当颗粒处于高速流体当中，较高的速度梯度造成了颗粒表面的压力差，从而产生了升力。Saffman升力可用下式表达：
2.3.1.2 Magnus升力模型
Magnus升力是由于颗粒自旋而产生的升力，球体的阻力系数与雷诺数较高的情况下的Magnus升力的关联可以表示为下式：
这个模型即使对于雷诺数高达2000时对于颗粒仍然适用 [2]。
2.3.1.3 Fluid-Induced Torque（流激诱导转矩）
在以上耦合的模型中，没有考虑由颗粒自旋产生的流体剪切力的影响。显然，若在高速流体中这个力不能忽略。根据Dennis等人的研究[3]，有以下公式：
通过这样的处理，这一扩展模型可用于颗粒雷诺数在32到1000的情况。
2.3.2 CFD-EDEM阻力模型
Fluent-EDEM耦合采用一种改进的自由流阻力来计算作用在颗粒球形上的力，对于颗粒表面为非球面的情况，先用对颗粒的球型包围盒进行阻力的计算。这种处理方法仅适合处理两种境况下的颗粒：（1）那些大小相同或者是体积小于流体网格单元的颗粒;（2）周围的流体参数（速度、密度和粘度等）的变化基本为常量的颗粒。阻系数CD，取决于雷诺数Re： 式中：r 是流体密度，h 为流体的粘度，L 是颗粒球的直径，v 是颗粒与流体间的相对速度，a 是CFD网格单元的自由体积。除了阻力，颗粒的固有浮力也必须考虑。浮力的计算公式是：
一般对于Lagrangian耦合来说，这种改进的阻力法则可以用来考虑颗粒载荷的影响;而在Eulerian耦合中，颗粒载荷需通过结合体积分数的方法来考虑这一影响。
3 FLUENT-EDEM耦合实例
3.1 计算模型建立
气固两相流数值计算过程中，颗粒在气流中场中的运动是从静止开始计算的。为了简化计算，同时使流场设置不至于过大，分析探测器发动机距月表较近时（发动机喷管出口垂直月表向下）羽流与月面月尘（月壤）的相互作用。
据此，将模型的参数设置为：基本流场为以5mm为半径的底面、15mm为高的圆柱，圆柱的上底面作为喷流的入口，下底面作为月面。月壤颗粒密度设为3100
，泊松比为0.42，弹性模量为10MPa，环境压强为0.01Pa;喷管的出口（亦即流场的入口）速度为1000m/s（因为喷管出口压力在0.02MPa左右，而环境压强仅为0.01Pa，喷流在喷管内处为极度膨胀不足的状态，所以达到这样大的速度也是合理的）。重力加速度为1.62。
边界条件设置如下：流场入口为压力入口边界条件;出口则为压力出口边界条件;在固体壁面上，气相流场采用无滑移固壁边界条件，颗粒与壁面之间采用Hertz Mindlin（No Slip）碰撞模型。
3.2 数值计算方法
采用FLUENT与EDEM軟件耦合求解气固两相流场，气流流场采用FLUENT软件，选用不可压缩流体、
紊流模型计算;在EDEM软件中采用Lagrangian坐标方法分析颗粒在弯管中运动和分布规律。首先由FLUENT得到原始气流场结果，并将流场参数传递给EDEM;由EDEM中DEM模块计算在该流场下颗粒受力的大小，然后通过内部迭代计算颗粒下一时间步长运动参数，同时更新颗粒运动轨迹;最后将EDEM计算结果返回FLUENT进行下一时间步长迭代。
采用Lagrangian方法计算颗粒运动轨迹，由Rayleigh数控制迭代步长。Rayleigh数由下面公式给出：
式中，R是颗粒半径，ρ是颗粒密度，G是剪切模量，υ是泊松比。则颗粒运动的迭代时间步长为：
FLUENT时间步长应为EDEM软件时间步长的1/10～1/100，由此来确定流场的网格尺度，故FLUENT和EDEM网格不同、时间尺度不同[4]。
3.3计算结果及分析
经过耦合计算之后，流场显示收敛，保存计算结果。通过EDEM自带的Analyst模块，我们可以导出进行加工处理后的数据，再对其进行相关的分析讨论。
计算流场的速度云图如下所示：
月壤颗粒的各项特性随仿真时间的推进的变化图如下：
为了加快仿真时间，在EDEM中的颗粒工厂只加入了1000个月壤颗粒进行计算，由上图可知，仿真只进行到了0.3s左右，颗粒已经几乎全部被吹出了仿真区域，而且被吹起的月壤的平均速度大多能够达到二三十米每秒，说明激起的扬尘足以对探测器上所携带的仪器产生干扰，会对探测器的安全及准确降落造成影响，与Apollo15的录像测量数据相吻合。
计算完成后，我们可以创建图表来记录指定的时间步长或者一定时间段内的仿真情况。以下便是不同时间点下一定速度范围内各月壤颗粒数目统计的直方图，从图中我们可以非常直观地得到仿真过程中某一时间点下的不同速度的颗粒分布状况，从而对我们的后续研究及观察提供有利的参考依据。
以下为各个时间点下流场中不同位置的颗粒的数目分布图，据此我们可以得出流场内颗粒的分布密度，从而判定是否对探测器上的设备造成影响。
4 结论
相较于以往的两相流计算仿真，FLUENT-EDEM耦合方法不需要像使用UDF添加颗粒相需要进行编程那么麻烦，而且EDEM由于是采用表面网格来描述边界表面，使其可以与FLUENT的流体网格的边界表面元素的点对点耦合;加上耦合方案为双向动力耦合，令仿真效果更加真实可信。由于可以直接在软件中手动设置颗粒的相关参数，可以很方便地进行各项属性的更改，而且EDEM也提供了非球体颗粒的设置，可以让仿真更加贴近实际情况。此外，EDEM中还加入了先进的静电模型来仿真带电颗粒，这与月壤因太阳风的吹扫而带电的情况相吻合。
参考文献：
[1]倪晋仁，李振山.风沙两相流理论及其应用[M].北京：科学出版社，2006.
[2]Oesterlé，B.and Bui Dinh，T.，Experiments on the lift of a spinning sphere in a range of intermediate Reynolds numbers，Experiments in Fluids，25，pp.16-22，1998.
[3]Zhou，H.，Flamant，G.，Gauthier，D.，DEM-LES simulation of coal combustion in a bubbling fluidized bed.Part II：coal combustion at the particle level，Chemical Engineering Science，59 pp.4205-4215，2004.
[4]许峰，王彤，傅耀.弯管内小尺度颗粒运动的数值模拟与分析.四川兵工学报.2010.02