CN106598912A - 一种基于cfd‑dem耦合模型的磨粒流场分析方法 - Google Patents

Abstract

本发明公开了一种基于耦合模型的磨粒流场分析方法，包括以下步骤：(1)建立并求解面向磨粒流加工的流体相控制方程组；(2)根据步骤(1)的求解结果，通过DEM理论求解磨粒运动方程；(3)根据磨粒运动方程，更新磨粒在流道中的位置，并求解新的体积分数以及磨粒受到的作用力；(4)结合新的流体体积分数以及磨粒受到的作用力，将其代入到步骤(1)中开始新一轮的计算；依照上述计算过程，不断循环更新每个磨粒的速度和位置，进而演化整个磨粒流加工系统，得到磨粒流运动规律。本发明将计算流体力学模型与离散元理论模型相结合的方法引入到磨粒流加工领域，能够实现对磨粒运动的精确分析，进而实现对磨粒流加工工艺参数的准确调控。

Description

一种基于CFD-DEM耦合模型的磨粒流场分析方法

技术领域

本发明涉及精密光整加工技术领域，更具体的说，尤其涉及一种基于CFD-DEM耦合模型的磨粒流场分析方法。

背景技术

精密光整加工技术是现代工业发展的重要支撑技术。其中，基于流体的磨粒流加工技术由于使用了柔性性良好的流体，使其在某些加工场合中具有无可比拟的优势。传统的磨粒流加工总体上可以分为两类，一类是使用强粘度流体进行挤压加工的传统挤压磨粒流加工技术，如挤压珩磨；另一类则是使用弱粘度或无粘度流体的磨粒流加工技术，如磨粒射流加工、弹性发射抛光、浮法抛光、化学机械抛光以及磁流变抛光等。磨粒流加工技术从本质来上均是通过流体驱动磨粒颗粒，使得磨粒颗粒不断冲击作用于工件的待加工表面，以此实现工件的光整加工。虽然磨粒冲击碰撞工件是一个基本的加工过程，但却是一个高度非线性行为。不同的流体物理特性(如流体相粘度、密度)、不同的颗粒尺寸大小以及不同的流道结构，均会影响磨粒运动。如何精确分析磨粒运动，从而准确得出磨粒流场运动规律，以实现能够对磨粒流加工工艺参数的精密调控，是目前磨粒流加工领域中亟待解决的技术难题。

目前有诸多科技人员研发了多种基于计算流体力学(Computational FluidDynamics,CFD)模型分析的调控方法。总体分为两类：离散相(Discrete Phase Model,DPM)模型以及欧拉-欧拉双流体模型。但基于该两类模型的调控方法均存在一定缺陷，其中采用DPM模型的调控方法将磨粒视为无体积质点，无法处理磨粒与工件碰撞运动，只能得出在流体驱动下的磨粒运动过程，精度受到限制；而采用欧拉-欧拉双流体模型的调控方法将磨粒相视为连续流体相，故借助该方法能得出工件表面大致的磨粒分布状态，从而估计工件表面加工强弱区域，但是无法得出磨粒运动过程。因此，研发一种能真正意义上考虑磨粒的物理存在，同时能考虑到磨粒-工件碰撞运动的分析方法，是一件非常有实际意义的工作。因此本发明通过将离散元模型(Discrete Element Method,DEM)引入到CFD模型中，通过数值求解技术，建立一种基于CFD-DEM耦合模型的磨粒流加工控制方法。

发明内容

本发明的目的在于克服现有磨粒流加工中难以对磨粒流运动进行精确分析，进而难以精确调控加工工艺参数的不足，提供了一种将离散元模型引入到CFD模型中，通过数值求解技术建立的能够有效分析整个流场内磨粒流运动，进而实现对磨粒流加工精确控制的基于耦合模型的磨粒流场分析方法。

本发明通过以下技术方案来实现上述目的：一种基于CFD-DEM耦合模型的磨粒流场分析方法，主要包括如下步骤：

(1)建立并求解面向磨粒流加工的流体相控制方程组：在CFD-DEM耦合策略分析中，磨粒采用DEM理论建模并视为真正的离散磨粒存在，流体相和磨粒之间不存在质量传递；假设由磨粒占据后，流体网格单元中流体相剩余的体积分数为α_l，磨粒受到的作用力为磨粒受到的平均作用力为S，由此得到CFD-DEM耦合框架下的流体相控制方程组为：

其中，ρ_l与分别为流体相的密度及速度；

为液相应力-应变张量，表达式如下：

其中，μ_l以及λ_l为流体相剪切粘度和体积粘度；为流体相湍动能强度；磨粒受到的平均作用力S表达式如下：

其中，V_l为流体单元网格体积；

(2)根据步骤(1)的求解结果，通过DEM理论求解磨粒运动方程：磨粒运动求解采用DEM理论求解，在DEM理论求解框架下，每一个磨粒运动在Lagrange坐标系下通过牛顿第二运动定律求解，磨粒运动求解方程为：

其中，与分别代表碰撞接触力以及重力，为流体对磨粒的作用力，为磨粒的重力，表达式为：

其中m_p为磨粒质量；

磨粒运动方程中磨粒各项受力求解如下：

2.1)流体对磨粒的作用力求解：

当磨粒处于高速流体当中，较高的速度梯度会造成磨粒表面的压力差，因此，必须考虑Saffman升力对磨粒的作用。流体对磨粒的驱动作用可用下式表示：

式中，为流体对磨粒的曳力，为流体对磨粒的压力梯度力，为流体对磨粒的浮力，为流体对磨粒的Saffman升力；

在磨粒流加工中，存在高浓度磨粒区域，因此必须考虑空隙率对曳力的影响，曳力由下式计算：

其中，χ＝3.7-0.65exp[-(1.5-logRe_p)²/2]，

式中，ρ_p为磨粒密度，为磨粒速度以及d_p为磨粒直径，C_d为曳力系数，Re_p为磨粒雷诺数；

曳力系数C_d通过下式计算：

磨粒雷诺数Re_p通过下式计算：

其中μ_l为流体相动力粘性系数；

流体对磨粒的压力梯度力由下式计算：

流体对磨粒的浮力由下式计算：

流体对磨粒的Saffman升力由下式计算：

V_p为磨粒的体积；

2.2)磨粒的碰撞接触力求解：

由于在流体驱动下，磨粒不断冲击工件，依靠磨粒与工件的碰撞接触力实现工件的光整加工，因此采用软球模型精确计算碰撞接触力，并采用弹性系数及阻尼系数来量化弹簧、阻尼器、滑动器的作用；接触力可以分解为法向接触力和切向接触力，假设弹性系数及阻尼系数分别为k及η，则法向接触力可由下式计算：

式中，α为磨粒接触时的法向重叠量；为磨粒i相对于磨粒j的速度， 为磨粒i球心位置到磨粒j球心位置的单位矢量；和分别为磨粒i和磨粒j的速度；

k_n是磨粒的法向弹性系数，其计算公式为：

η_n是磨粒的法向阻尼系数，其计算公式为：

式中，E和υ是磨粒材料的弹性模量和泊松比，a是磨粒半径；

切向接触力可由下式计算：

式中，k_t和η_t是切向弹性系数和切向阻尼系数；δ是接触点的切向位移；

是接触点的滑移速度，可由下式计算：

式中，a_i和a_j分别是磨粒i和磨粒j的半径；和分别是磨粒i和磨粒j的角速度；

k_t是磨粒的切向弹性系数，由下式计算：

η_t是磨粒的切向阻尼系数，由下式计算：

综合考虑法向力和切向力的作用，磨粒i所受到的合力为：

磨粒i所受到的力矩为：

磨粒流加工过程中，多个磨粒可能同时发生接触，此时，磨粒i所受到的合力为：

磨粒i所受到的合力矩为：

(3)根据步骤(2)中求得的磨粒运动方程，更新磨粒在流道中的位置，并求解新的体积分数α_p以及磨粒受到的作用力

α_p计算过程如下：首先在每一个磨粒周围采用边界盒将其包围起来；之后在每一个边界盒内规则的抽取样本点，如果样本点位于磨粒表面内，则将其保存；最后，对每一个样本点进行检验，确认样本点是否位于网格单元内；因此，某个网格单元内的固体体积分数就是网格单元内的样本点数量在所有样本点数量中所占的比例，由下式计算：

其中，V_particle为磨粒体积，n_c为网格单元中样本点数量，N为边界盒内样本点总数量；

(4)结合步骤(3)中新的流体体积分数α_l以及磨粒受到的作用力将流体体积分数α_l以及磨粒受到的作用力代入到步骤(1)中开始新一轮的计算；依照步骤(1)、步骤(2)和步骤(3)的计算过程，不断的循环更新每个磨粒的速度和位置，进而确定性的演化整个磨粒流加工系统，得到磨粒流运动规律。

本发明的技术构思为：将计算流体力学模型与离散元理论模型相结合的方法引入到磨粒流加工领域，能够实现对磨粒运动的精确分析，进而实现对磨粒流加工工艺参数的准确调控。

本发明的有益效果在于：

1)离散磨粒相通过DEM理论建模，视为真正的固体颗粒存在，可以直观的获得磨粒运动轨迹。

2)采用软球模型求解磨粒与工件的碰撞作用过程，可以较为精确的获得磨粒与工件碰撞作用力以及碰撞接触分布等影响加工效率与均匀性的关键参数。

3)该方法具有成本低、精度高等显著优点，并可以较为容易获得某些通过实验方法不容易获得的磨粒运动规律。

附图说明

图1是本发明CFD-DEM耦合计算过程图。

图2是本发明实施例中三入口面约束磨粒流加工装置主视图。

图3是本发明实施例中三入口面约束磨粒流加工装置俯视图。

图4是本发明流道的网格划分及边界条件示意图。

图5是本发明稳定运动状态下磨粒速度分布图。

图6是本发明工件表面磨粒碰撞分布图。

图中，1-磨粒流入口、2-非约束空间、3-磨粒流、4-加工工具、5-工件、6-约束空间、7-夹具；8-速度入口边界、9-自由出流边界、10-壁面边界。

具体实施方式

下面结合某一具体实例对本发明进行详细说明。

如图1-6所示，具体实施例采用一种三入口面约束磨粒流加工装置对工件2进行加工，圆形加工工具沿圆周方向开有三个磨粒流入口1，磨粒流3在水泵的驱动作用下从磨粒流入口1进入到圆形加工工具中部的非约束空间2；在充分混合后，磨粒流3以更高的流速进入到约束空间6内，进而对安装在圆形加工工具底部的工件2进行光整加工。采用本发明分析该加工方法中工件表面处的磨粒流场，具体过程如下：

(1)建立流道有限元模型并设定边界条件

根据三入口面约束磨粒流加工装置的结构，建立相应流道模型，并采用非结构化网格划分，计算区域共有71857个四面体网格单元，如图4所示；流道入口采用速度入口边界条件，出口采用自由出流条件，其余均为壁面条件；

(2)设定基本求解参数

流体采用10#航空液压油，密度900kg/m3，动力粘度0.008kg/m·s；磨粒选用碳化硅(SiC)，密度3170kg/m3，直径0.05mm，磨粒体积分数为0.1；压力速度耦合处理使用压力关联方程的半隐式算法，保证收敛率；压力离散插值方式采用标准方式(Standard)；采用基于一阶迎风格式的有限体积法进行动量、湍动能以及湍流耗散率的离散。流体相求解采用非稳态的方式，时间步设为0.0000121s。

(3)建立模型并开始求解

3.1)建立并求解流体相控制方程组：在CFD-DEM耦合策略分析中，磨粒采用DEM理论建模并视为真正的离散磨粒存在，流体相和磨粒之间不存在质量传递；假设由磨粒占据后，流体网格单元中流体相剩余的体积分数为α_l，磨粒受到的作用力为磨粒受到的平均作用力为S，由此得到CFD-DEM耦合框架下的流体相控制方程组为：

其中，ρ_l与分别为流体相的密度及速度；

为液相应力-应变张量，表达式如下：

其中，μ_l以及λ_l为流体相剪切粘度和体积粘度；为流体相湍动能强度；磨粒受到的平均作用力S表达式如下：

其中，V_l为流体单元网格体积；

3.2)根据步骤(1)的求解结果，通过DEM理论求解磨粒运动方程：磨粒运动求解采用DEM理论求解，在DEM理论求解框架下，每一个磨粒运动在Lagrange坐标系下通过牛顿第二运动定律求解，磨粒运动求解方程为：

其中，与分别代表碰撞接触力以及重力，为流体对磨粒的作用力，为磨粒的重力，表达式为：

其中m_p为磨粒质量；

磨粒运动方程中磨粒各项受力求解如下：

a)流体对磨粒的作用力求解：

当磨粒处于高速流体当中，较高的速度梯度会造成磨粒表面的压力差，因此，必须考虑Saffman升力对磨粒的作用。流体对磨粒的驱动作用可用下式表示：

式中，为流体对磨粒的曳力，为流体对磨粒的压力梯度力，为流体对磨粒的浮力，为流体对磨粒的Saffman升力；

在磨粒流加工中，存在高浓度磨粒区域，因此必须考虑空隙率对曳力的影响，曳力由下式计算：

其中，χ＝3.7-0.65exp[-(1.5-logRe_p)²/2]，

式中，ρ_p为磨粒密度，为磨粒速度以及d_p为磨粒直径，C_d为曳力系数，Re_p为磨粒雷诺数；

曳力系数C_d通过下式计算：

磨粒雷诺数Re_p通过下式计算：

其中μ_l为流体相动力粘性系数；

流体对磨粒的压力梯度力由下式计算：

流体对磨粒的浮力由下式计算：

流体对磨粒的Saffman升力由下式计算：

V_p为磨粒的体积；

b)磨粒的碰撞接触力求解：

由于在流体驱动下，磨粒不断冲击工件，依靠磨粒与工件的碰撞接触力实现工件的光整加工，因此采用软球模型精确计算碰撞接触力，并采用弹性系数及阻尼系数来量化弹簧、阻尼器、滑动器的作用；接触力可以分解为法向接触力和切向接触力，假设弹性系数及阻尼系数分别为k及η，则法向接触力可由下式计算：

式中，α为磨粒接触时的法向重叠量；为磨粒i相对于磨粒j的速度， 为磨粒i球心位置到磨粒j球心位置的单位矢量；和分别为磨粒i和磨粒j的速度；

k_n是磨粒的法向弹性系数，其计算公式为：

η_n是磨粒的法向阻尼系数，其计算公式为：

式中，E和υ是磨粒材料的弹性模量和泊松比，a是磨粒半径；

切向接触力可由下式计算：

式中，k_t和η_t是切向弹性系数和切向阻尼系数；δ是接触点的切向位移；

是接触点的滑移速度，可由下式计算：

式中，a_i和a_j分别是磨粒i和磨粒j的半径；和分别是磨粒i和磨粒j的角速度；

k_t是磨粒的切向弹性系数，由下式计算：

η_t是磨粒的切向阻尼系数，由下式计算：

综合考虑法向力和切向力的作用，磨粒i所受到的合力为：

磨粒i所受到的力矩为：

磨粒流加工过程中，多个磨粒可能同时发生接触，此时，磨粒i所受到的合力为：

磨粒i所受到的合力矩为：

3.3)根据步骤(2)中求得的磨粒运动方程，更新磨粒在流道中的位置，并求解新的体积分数α_p以及磨粒受到的作用力

α_p计算过程如下：首先在每一个磨粒周围采用边界盒将其包围起来；之后在每一个边界盒内规则的抽取样本点，如果样本点位于磨粒表面内，则将其保存；最后，对每一个样本点进行检验，确认样本点是否位于网格单元内；因此，某个网格单元内的固体体积分数就是网格单元内的样本点数量在所有样本点数量中所占的比例，由下式计算：

其中，V_particle为磨粒体积，n_c为网格单元中样本点数量，N为边界盒内样本点总数量；

3.4)结合步骤(3)中新的流体体积分数α_l以及磨粒受到的作用力将流体体积分数α_l以及磨粒受到的作用力代入到步骤3.1)中开始新一轮的计算；依照步骤3.1)、步骤3.2)和步骤3.3)的计算过程，不断的循环更新每个磨粒的速度和位置，进而确定性的演化整个磨粒流加工系统，得到磨粒流运动规律。

(4)结果分析

计算收敛后，即可对磨粒流运动规律进行分析。附图5所示为达到稳定运动状态后磨粒速度分布图。从该图中可以看出，磨粒运动经历了一个先加速，后碰撞减速，再加速的过程，最高速度可以达到90m/s以上。附图6所示为使用该分析方法得出的磨粒对工件表面的碰撞落点分布图。从图中可以得出，该加工方法不能保证加工的均匀性，需要在加工中周期性旋转工件以实现整体工件表面的加工均匀。

上述实施例只是本发明的较佳实施例，并不是对本发明技术方案的限制，只要是不经过创造性劳动即可在上述实施例的基础上实现的技术方案，均应视为落入本发明专利的权利保护范围内。

Claims (1)

1.一种基于CFD-DEM耦合模型的磨粒流场分析方法，其特征在于：主要包括如下步骤：

(1)建立并求解面向磨粒流加工的流体相控制方程组：在CFD-DEM耦合策略分析中，磨粒采用DEM理论建模并视为真正的离散磨粒存在，流体相和磨粒之间不存在质量传递；假设由磨粒占据后，流体网格单元中流体相剩余的体积分数为α_l，磨粒受到的作用力为磨粒受到的平均作用力为S，由此得到CFD-DEM耦合框架下的流体相控制方程组为：

其中，ρ_l与分别为流体相的密度及速度；

为液相应力-应变张量，表达式如下：

其中，μ_l以及λ_l为流体相剪切粘度和体积粘度；为流体相湍动能强度；磨粒受到的平均作用力S表达式如下：

其中，V_l为流体单元网格体积；

(2)根据步骤(1)的求解结果，通过DEM理论求解磨粒运动方程：磨粒运动求解采用DEM理论求解，在DEM理论求解框架下，每一个磨粒运动在Lagrange坐标系下通过牛顿第二运动定律求解，磨粒运动求解方程为：

其中，与分别代表碰撞接触力以及重力，为流体对磨粒的作用力，为磨粒的重力，表达式为：

其中m_p为磨粒质量；

磨粒运动方程中磨粒各项受力求解如下：

2.1)流体对磨粒的作用力求解：

当磨粒处于高速流体当中，较高的速度梯度会造成磨粒表面的压力差，因此，必须考虑Saffman升力对磨粒的作用；流体对磨粒的驱动作用可用下式表示：

式中，为流体对磨粒的曳力，为流体对磨粒的压力梯度力，为流体对磨粒的浮力，为流体对磨粒的Saffman升力；

在磨粒流加工中，存在高浓度磨粒区域，因此必须考虑空隙率对曳力的影响，曳力由下式计算：

其中，χ＝3.7-0.65exp[-(1.5-logRe_p)²/2]，

式中，ρ_p为磨粒密度，为磨粒速度以及d_p为磨粒直径，C_d为曳力系数，Re_p为磨粒雷诺数；

曳力系数C_d通过下式计算：

$C_d={(0.63+\frac{4.8}{{\mathrm{Re}}_p^{0.5}})}^2,$

磨粒雷诺数Re_p通过下式计算：

其中μ_l为流体相动力粘性系数；

流体对磨粒的压力梯度力由下式计算：

流体对磨粒的浮力由下式计算：

流体对磨粒的Saffman升力由下式计算：

V_p为磨粒的体积；

2.2)磨粒的碰撞接触力求解：

由于在流体驱动下，磨粒不断冲击工件，依靠磨粒与工件的碰撞接触力实现工件的光整加工，因此采用软球模型精确计算碰撞接触力，并采用弹性系数及阻尼系数来量化弹簧、阻尼器、滑动器的作用；接触力可以分解为法向接触力和切向接触力，假设弹性系数及阻尼系数分别为k及η，则法向接触力可由下式计算：

式中，α为磨粒接触时的法向重叠量；为磨粒i相对于磨粒j的速度， 为磨粒i球心位置到磨粒j球心位置的单位矢量；和分别为磨粒i和磨粒j的速度；

k_n是磨粒的法向弹性系数，其计算公式为：

$k_n=\frac{4}{3}{(\frac{1-{\mathrm{\υ}}_i^2}{E_i}+\frac{1-{\mathrm{\υ}}_j^2}{E_j})}^{-1}{\left(\frac{a_i+a_j}{a_i{a}_j}\right)}^{-1/2},$

η_n是磨粒的法向阻尼系数，其计算公式为：

${\mathrm{\η}}_n=2\sqrt{{\mathrm{mk}}_n},$

式中，E和υ是磨粒材料的弹性模量和泊松比，a是磨粒半径；

切向接触力可由下式计算：

式中，k_t和η_t是切向弹性系数和切向阻尼系数；δ是接触点的切向位移；

是接触点的滑移速度，可由下式计算：

式中，a_i和a_j分别是磨粒i和磨粒j的半径；和分别是磨粒i和磨粒j的角速度；

k_t是磨粒的切向弹性系数，由下式计算：

$k_t=8{\mathrm{\α}}^{1/2}{(\frac{1-{\mathrm{\υ}}_i^2}{G_i}+\frac{1-{\mathrm{\υ}}_j^2}{G_j})}^{-1}{\left(\frac{a_i+a_j}{a_i{a}_j}\right)}^{-1/2}$

η_t是磨粒的切向阻尼系数，由下式计算：

${\mathrm{\η}}_t=2\sqrt{{\mathrm{mk}}_t},$

综合考虑法向力和切向力的作用，磨粒i所受到的合力为：

磨粒i所受到的力矩为：

磨粒流加工过程中，多个磨粒可能同时发生接触，此时，磨粒i所受到的合力为：

磨粒i所受到的合力矩为：

(3)根据步骤(2)中求得的磨粒运动方程，更新磨粒在流道中的位置，并求解新的体积分数α_p以及磨粒受到的作用力

α_p计算过程如下：首先在每一个磨粒周围采用边界盒将其包围起来；之后在每一个边界盒内规则的抽取样本点，如果样本点位于磨粒表面内，则将其保存；最后，对每一个样本点进行检验，确认样本点是否位于网格单元内；因此，某个网格单元内的固体体积分数就是网格单元内的样本点数量在所有样本点数量中所占的比例，由下式计算：

${\mathrm{\α}}_p=1-{\mathrm{\α}}_l=\underset{particles}{\Σ}\frac{n_c}{N}{V}_{particle},$

其中，V_particle为磨粒体积，n_c为网格单元中样本点数量，N为边界盒内样本点总数量；

(4)结合步骤(3)中新的流体体积分数α_l以及磨粒受到的作用力将流体体积分数α_l以及磨粒受到的作用力代入到步骤(1)中开始新一轮的计算；依照步骤(1)、步骤(2)和步骤(3)的计算过程，不断的循环更新每个磨粒的速度和位置，进而确定性的演化整个磨粒流加工系统，得到磨粒流运动规律。