基于湍流弥散模型的批量薄壁孔件抛光数值模拟与分析①

单晓杭 汪佳成 张 利 叶必卿

(浙江工业大学机械工程学院 杭州 310032）

0 引言

医学诊疗技术进步的关键推动因素是医疗器械的发展,我国医疗器械行业已经进入快速发展阶段。近年来,医疗器械市场规模增长率一直保持在20%左右,远超全球市场增长率,并且早在2014 年我国已经成为医疗器械的全球第二大市场[1]。但是我国医疗器械行业的发展目前主要依靠低端产品,高精零部件大多还是要依赖进口。

临床医疗器械和植入类医疗器械通常具有结构紧凑、加工精度要求高、材料要求苛刻等不利于加工制造的特点。在医疗器械中微小孔是最常见的结构,不仅用于连接,甚至有的器械本身就是细长的小孔,如内窥镜导光口、导尿管等。国外目前通过电加工技术来制造导尿管、内视镜测量件以及缝合器械等临床医疗器械[2],但是电解加工存在诸多缺点,如加工精度与加工稳定性较差,电极工具的设计比较复杂,电解加工需要的设备繁多,加工成本高,环境污染严重等。我国通过引入国外高精密数控机床来加工微创外科手术器械小孔[3],但是通过机加工来制造微创器械仍无法脱离手工制造,工艺链中仍需通过磨抛和钳工来提高精度,导致加工效率低下,并且精度误差大。对于体外医疗器械,如助听器、假肢等医疗模型和医疗器械,3D 打印技术如今应用非常广泛[4]。3D 打印技术具有节省材料、可以实现复杂曲面的制造等优点,但是3D 成型加工周期长,无法规模化生产,并且由于打印材料受到限制导致成型后的器材小孔无法进行后续的加工来保证精度要求。

金属医疗器械的表面质量决定了治疗效果的好坏,尤其是植入人体的金属医疗器械,其表面质量直接关系到其疗效和并发症的发生[5]。钛合金材料具有强度高、抗腐蚀性强、生物相容性优越等特点,被视为可植入体的首选材料[6],但是钛合金工艺性能比较差、表面硬度高的特性为加工带来了困难。目前对于医疗器械中的小孔加工主要分为传统加工方法和非传统加工方法,传统加工方法主要是钻孔、冲孔、磨孔等典型机加工方案,由于医疗器械通常采用钛合金作为原材料,钛合金的热传导性比较差,在切削过程中切削热引起的月牙洼易导致加工误差[7]。非传统加工主要包括超声波、离子束、激光等方法,但是这些方式均存在设备费用与技术水平要求高、难以大规模推广[8]等问题,因此在控制成本的前提下实现微小孔批量化精密加工是当前医疗器械行业的燃眉之急。

针对复杂曲面和小孔工件加工的研究,文献[9-12]提出了一种软性磨粒流的加工方式,由于磨粒流自身的流动性,被加工的零件几乎没有几何尺寸与零件类型限制。磨粒在流场的作用下可以对工件表面进行微量切削以实现光整加工。文献[13-16]证实了磨粒流抛光小孔可以有效降低工件的表面粗糙度,但是目前关于这方面的研究仅仅停留在如何抛光单个工件的层面上,从实现单个工件的抛光到批量抛光,难点在于需要专门的设备以及专门的夹具来保障所有工件抛光均匀性的同时还要考虑到加工效率。对此本文借助计算流体力学方法、材料冲蚀模型以及正交实验思想模拟夹具几何尺寸的变化对每一个工件的抛光效果的影响,选取对于单个工件抛光效果最好、对于整体抛光均匀性最好的最优方案,研究结果可为今后薄壁类零件的批量化抛光提供一定的方案指导。

1 材料去除量分析

软性磨粒流是一种低黏度固液两相流,具有良好的流动性。磨粒在连续相的作用下经过孔壁时做随机运动,磨粒颗粒可视为微型刀具,颗粒会对所加工零件表面的凸起部分产生碰撞和剪切的作用,可以实现孔壁表面纹理的无序化光整加工。

Preston 方程[17]被广泛应用于磨削加工。在Preston 方程中,磨粒对于工件的材料去除量可以用磨料的速度和磨料对壁面的压力以及比例常数Kp来表示。

式中,Δz表示磨削去除量;v表示磨粒在壁面处的相对速度,P表示磨粒在壁面处的压力;Kp为Preston 常数,该常数与磨粒和工件本身的物理性质以及磨粒与壁面作用时的因素(如磨粒与壁面的碰撞角度)有关。

根据式(1)可以得出,零件的材料去除量与磨料在壁面处的相对速度、压力大小、加工时间、切削系数成正比。决定切削系数的因素主要有磨粒的直径、磨粒与流体的混合比、磨粒的弹性、工件的硬度和表面粗糙度,在同一个模型中,同样的加工条件下决定抛光效果是否均匀的主要影响因素是流体磨料对工件的壁面压力与壁滑速度[18]。因此要使一批工件的抛光效果均匀,只需要保证近壁面处磨粒的相对速度与压力乘积在各个工件上保持均匀即可。本文据此理论来分析一批工件的磨粒流抛光均匀性问题,观测指标是所有工件的速度与压力乘积的均方差和平均值的大小以及均方差W1和平均值W2的比值α。

式中,W1与W2分别表示所有工件内孔的磨粒在近壁面处相对速度与压强乘积值,α用来表示一批工件材料去除量的离散程度。

2 磨粒流抛光数学模型

2.1 液相湍流模型

液相湍流方程采用雷诺时均方程,雷诺时均法对N-S 方程进行平均化处理,抹平所有尺度脉动,计算出流动的平均量。大量的工程实践应用表明,该模型可以较好地预测管流、通道流动、喷管内的流动等,适用于大雷诺数的场景,与本文研究的问题相符合。雷诺时均连续方程为

动量方程为

式中,ui表示略去平均符号的雷诺平均速度分量,ρ为密度,p为压强,为脉动速度,σij为应力张量分量。

k -ε模型假设湍流粘性和湍动能以及耗散率有关,标准k -ε方程形式如式(5)所示。

式中k和ε分别为湍动能和湍流耗散率,Pt为湍动能生成项,ut为湍流粘性系数,模型常数通常为Cζ1=1.44,Cζ2=1.92,σζ=1.3,σk=1.0,Cu=0.09。

2.2 磨粒受力分析

磨粒在两相流受力主要分为体积力和表面力,在本文中只考虑重力和曳力。根据动力学方程推导出单颗磨料力学方程:

(1)曳力

在连续相-离散相的两相流流动过程中,当磨粒与流体介质速度存在差异时就会产生相互作用力。

曳力的表达式为

2.3 磨料冲蚀模型

冲蚀是指工件受到流动粒子冲击时表面出现破坏的一类磨损现象。冲蚀现象是由冲击速度、磨粒直径与硬度、冲击角等多种因素决定,目前业界对此没有统一的说法。文献[19]做了大量的冲蚀实验,最终提出了可以用于不同材质的材料冲蚀模型,即Oka 模型。Oka 模型的冲蚀率表达式为

式中,E90为冲击角为直角时的参考冲蚀率,f(γ)为冲击角函数,γ为磨粒冲击角度,V为磨粒冲击速度,Vref为磨粒与壁面间的相对速度,k2为速度指数,d为磨粒直径,dref为用户指定的参考磨粒直径,k3为直径系数。

式中,n1与n2为冲击角函数常数,Hv为工件材料维氏硬度。

根据Oka 模型,冲蚀模型中n1、n2、k2和k3分别可表示为

Oka 模型用于估算固体颗粒碰撞而引起的侵蚀破坏,所以不用考虑工件的几何形状,只需确定工件与磨粒颗粒的材料参数即可确定模型参数,因此本文颗粒冲蚀模型选用Oka 模型来模拟粒子对壁面的冲蚀效果。本文选定工件材料为Ti-6Al-4V 钛合金、磨粒材质为碳化硅,因此可确定Oka 冲蚀模型的参数如表1 所示。

表1 冲蚀模型

2.4 湍流弥散模型

通过上文可以得出在流场中磨粒的受力主要来自曳力,由式(7)可得,曳力的主要影响因素为流体的速度,而在雷诺时均法中抹平了湍流场中的脉动,只体现了平均值,并且雷诺方程是各向同性模型,即湍流量在给定位置的方向都相同,这一模型简化了求解,但是在一定程度上对于磨粒在湍流场中的运动状态难以还原,同时Oka 模型的冲蚀率预测的准确性关键在于磨料的速度矢量与真实值是否接近,因此本文提出一种基于线性同余法和Tanaka 碰撞模型[20]及Nanbu 算法的湍流弥散模型。因为脉动值的变化本身是随时间以一种无规则的方式在变化,虽然目前通过计算机软件还无法生成真正意义上的随机数,但是具有类似于随机数的统计特征,并且线性同余法是生成随机数的理想模型之一。线性同余法的递归表达式如下:

式中,A为乘子,B为增量,m为模,当n为1 时Xn-1为初始值,初始值可以任意指定一个奇数,根据递推关系求得数值后再取Yn=Xn/m便可以得到均匀随机数序列。为了使生成的随机序列周期比较大,那么则需要保证m的值尽可能大,同时为了提高随机性,取m为系统函数的当前毫秒值,取A为湍流的时均速度,取增量B为数值计算过程的时间t。

因为湍流的脉动量是在时均量的上下波动,因此还需要用生成的随机值Yn与时均值做取余运算,同时为了实现各向异性的效果还需要加上随机向量:

式中,U为脉动速度矢量,u为时均速度,I为随机向量。因此曳力的计算公式就用此湍流弥散模型中重新构建的速度U来代替,即:

由于在雷诺时均场中速度是时均量,因此粒子的自然碰撞行为也被忽略了,本文借助Tanaka 粒子碰撞模型以及Nanbu 算法来模拟粒子间碰撞。Tanaka 假设在某体积为V的网格内任意两颗粒的碰撞概率为

式中,Dp为颗粒直径,为两颗粒的相对速度,nσ为网格的颗粒浓度。

由式(14)可知,追踪颗粒i与颗粒1,2,3,…,j,…,N碰撞概率的总和为

在一个时间步长Δt内上述概率满足:

如图1 所示,Nanbu 算法将单位长度划分为N+1 个区间,每一个区间的长度为粒子i与网格内另一个粒子碰撞的概率,然后取在0～1 上均匀分布的随机数列中的一个随机值R,如果R∈[Pi,1] 则发生碰撞,否则不碰撞。

图1 Nanbu 算法示意图

为了直观体现湍流弥散方法对粒子在湍流场中无序运动的还原,本文设置了如下情形:在直径为10 mm 的圆管中设置入口流速为10 m/s,在入口处每隔固定时间释放2000 个粒子,颗粒在曳力的作用下运动,分别模拟了未应用湍流弥散方法和使用湍流弥散方法时颗粒的运动过程,如图2 所示。

图2 不同模型的粒子轨迹

图2 是在同一时刻下粒子所处的位置,可以从图2(a)中看到,在曳力的作用下所有的粒子速度相差无几,运动轨迹也相同,基本上是在同一时间进入管道,从同一时间离开管道,并且忽略了粒子碰撞,很多粒子几乎重叠在了一起,这显然不符合自然现象;但图1(b)中相同时间相同位置放入的粒子在湍流弥散方法的作用下运动轨迹发生了较大的改变,变得杂乱、无序,且粒子与粒子之间会发生碰撞。

为了验证湍流弥散模型对材料冲蚀的有效性,本文设置了Vieira 实验[21]工况下的仿真实验并与其实测的材料去除量进行比对,Vieira 实验参数如表2 所示。

表2 Vieira 实验工况

图3 和图4 是在第1 组实验工况下的两种模型的冲蚀率图,对比图3 和图4 可以明显看出,湍流弥散模型由于通过随机扰动模拟了粒子的脉动速度以及考虑了粒子的碰撞效果,冲蚀率要高于无湍流弥散模型,V 形冲蚀痕迹也得到了显著增强且更接近Vieira 实测值。并且通过表3 以及图5 对比11 组工况下的仿真结果,可以明显看出使用湍流弥散模型由于考虑了粒子速度的脉动量以及粒子相互碰撞的效应,对于预测材料去除量的误差有了明显改善,平均误差下降了7.77%。

图3 无湍流弥散模型冲蚀率图

图4 湍流弥散模型冲蚀率图

表3 仿真模型误差对比

图5 两种模型误差量

3 研究思路

本文中的抛光对象为孔径为3 mm、壁厚为1 mm、长度为15 mm 的薄壁细孔,该孔径与孔长度尺寸是医疗器械中典型的结构,如心脏起搏器中的连接器、内窥镜组件等,其孔径表面粗糙度一般需要达到Ra0.8 的精度要求。文献[22]已经可以实现单个工件磨粒流抛光后将表面粗糙度降低至0.5 μm,因此本文只要实现批量抛光单个工件之间材料去除量均匀即可。本文设计了如图6 和图7 所示的结构作为夹具。其中区域Ⅰ代表入口区域,区域Ⅱ代表缓冲区域,区域Ⅳ为工件实际加工区域。夹具的入口区域Ⅰ作为磨粒流的入口直接与泵联通,区域Ⅱ的下方是放置装夹工件的夹具。由于批量化抛光的需求,因此为了能够一次性加工比较多的工件,区域Ⅱ的直径会相对比较大,而区域Ⅰ的直径比区域Ⅱ小得多,从区域Ⅰ到区域Ⅱ因为流道直径的突变就会造成流场各处速度与压强分布不均匀,从而导致抛光区域Ⅳ中的工件抛光可能存在不均匀的现象,如图8 所示。若区域Ⅰ与区域Ⅱ的直径一致,则不存在沿着圆周分布每个工件抛光不均匀的问题。但是区域Ⅰ与区域Ⅱ的直径如果一致将会带来泵的选型困难和耗能过大的问题,因为磨粒流抛光要求流场中的速度、压力不能太小。区域Ⅲ属于过渡区域,因为从区域Ⅱ到区域Ⅲ存在缩口结构也会导致流场在区域Ⅲ中分布不均匀,所以这一区域不作为实际加工区域。

图6 夹具三维结构示意图

图7 夹具结构示意图

图8 湍流动能分布图

磨粒流从区域Ⅱ流入到工件内孔是一个从大口径流道突变到小口径流道的过程,因而就出现孔壁面处流场特征量分布紊乱,不利于抛光均匀性,所以需要区域Ⅲ来引流,消除这一段不均匀的流道对抛光均匀性的影响。这一段区域的流道直径与工件的孔径大小相同但并不作为实际加工区域,磨粒流流过这一段区域以后就会达到速度与压强分布相对均匀的状态。在此模型中影响抛光效果均匀性的因素为入口直径、入口流速、入口区域长度、缓冲区域直径和缓冲区域长度。

从理论上来说,本文提出的研究思路可以适用于任意数量的批量工件抛光问题,对于数值模拟没有任何影响。但是从实际的加工效果考虑,当工件的数量增加时,夹具的几何尺寸必然增大,这就造成泵的选型困难问题,因此为了保证后续实际加工的抛光效果,根据已有的实验设备和磨粒流抛光的研究经验综合确定夹具内流道直径最大处为120 mm以下,相应的工件装夹数量在70 个以下。

基于以上分析可以得知,流场中速度与压强分布的不均匀性将会导致工件受到的切削力在流场中随着空间位置的改变而变化,最终导致材料去除量存在显著的差异。因此为了实现均匀抛光的效果,需要研究夹具的几何尺寸的变化对工件的抛光均匀性产生怎样的影响。本文采用正交实验法设计了16 组对比实验,通过上文中提及的湍流弥散模型以及雷诺时均方程来计算湍流场中粒子速度并采用Oka 冲蚀模型来计算磨粒的冲蚀效果。结合式(2)计算每一次实验的离散度、每个工件速度与压强乘积的平均值与均方差值,然后结合统计学原理综合平衡法进行分析,首先分别对每个指标进行单因素分析,然后把每个指标的分析结果进行综合平衡得出结论。

4 数值模拟

4.1 计算对象与网格划分

使用三维建模软件建立流道的实体模型,然后通过Workbench-Meshing 划分流体网格。所有流道网格种类均为四面体网格,四面体网格比六面体网格具有更为均匀的疏密程度,因此能够更好地模拟零件表面的形态。入口区域与缓冲区域的网格比加工区域粗糙一些,网格单元大小为0.7 mm,而加工区域因为要计算每个工件壁面的压强与速度的乘积值,网格质量要求比入口区域与缓冲区域要高,为0.28 mm。网格顶点单元数为156,边单元数为3565,边界单元数为81 328,总单元数为1 756 012。然后通过COMSOL 软件进行数值模拟,确定不同参数对抛光均匀性的影响,流道结构模型如图9 所示,工件分布如图10 所示。

图9 流道模型图

图10 工件分布图

本文为每个因素设置4 个水平,表4 为正交实验因素水平表。设计L16(45)的正交表如表5所示。

表4 因素水平表

表5 正交实验表

4.2 数值计算初始条件设置

运用欧拉-拉格朗日方法对加工流场与磨粒运动进行数值计算,采用湍流弥散模型来修正磨粒在流体中的运动状态。本文模型的液相为水,磨粒颗粒为碳化硅,本模型具有如下特点:(1)连续相为恒温不可压缩流体,离散相的物理特性均为常数,不受流场的影响而变化;(2)所有颗粒的几何形状都按直径大小相同的球形来处理;(3)磨粒在流场中只考虑重力与曳力的影响,不考虑浮力、附加质量力、Maguns 力、Saffaman 力、压力梯度力和Basset 力等[23]。

湍流模型选择RANS 标准k -ε模型计算连续相流场运动,使用Oka 冲蚀模型模拟磨粒对工件表面的冲蚀效果。设置入口速度作为边界条件,出口条件为抑制回流,磨粒体积分数为10%。根据以上仿真参数以及模型结构得到的仿真结果如表6 所示。

4.3 仿真结果分析与讨论

从仿真结果中可以看出,当夹具的几何尺寸、流场的初始条件发生改变时,抛光的均匀性变化特别显著,可以明显看出当因素水平发生变化时会显著影响工件的抛光效率和抛光均匀性。

分析本次正交实验的结果,最终分别得到缓冲区域直径、缓冲区域长度、入口区域长度、入口流速、入口直径对速度与压力乘积的平均值和均方差影响的主次因素,从有限的组合中找到最好的实验参数水平组合。将表6 中的16 组数据导入到Minitab 软件,通过数理统计中极差分析法和综合平衡法分析实验的因素水平对抛光均匀性的影响,首先通过单指标分析每一个因子对评价指标影响的主次。将每一个因子中同一个水平所对应的实验响应综合记为Ki(i=1,2,3),极差Rj=[max(Ki) -min(Ki)](i=1,2,3)越大说明当前因子对观测指标的影响程度就越大,最后通过交互作用分析出最优的组合从而得出结论。

表6 仿真结果

4.3.1 参数对抛光后所有工件的速度与压力乘积的均方差影响

从表7 极差影响可以得出参数对批量抛光工件速度与压力乘积的均方差影响的主次顺序为:入口直径＞缓冲区域直径＞入口流速＞缓冲区域长度＞入口区域长度。通过均方差响应图(图11)可知:(1)缓冲区域直径的变化会对均方差值的大小产生显著改变,从直径70 mm 增大到80 mm 均方差变大,增大到80 mm至90 mm之间均方差变小,并且影响更显著,往后再持续增大缓冲区域直径均方差的变化量变小;(2)随着缓冲区域长度的增大,总体上来看,均方差的值呈变小的趋势,只有在缓冲区域长度为50 mm 至60 mm 之间时缓冲区域长度增加,均方差的值局部变大;(3)入口区域的长度变大,均方差的值会变小;(4)随着入口流速的增加,均方差的值会变大,但是这并不意味着离散程度就增加,因为入口流速增加意味着流场中的特征量也发生了变化,所以离散程度不能仅凭均方差的改变来说明;(5)入口直径变大,均方差随之也变大。但是由于速度不变的情况下改变入口直径相当于增加流量,导致湍流情况产生变化,所以这不能说明抛光均匀性变差。

图11 均方差主效应图

表7 均方差均值响应值(望小)

4.3.2 参数对抛光后所有工件的速度与压力乘积的平均值影响

从表8 可以得出,参数对批量抛光工件效率影响的主次顺序为:入口直径＞入口流速＞缓冲区域直径＞入口区域长度＞缓冲区域长度。通过分析均值响应表8 以及图12 可以得知:(1)缓冲区域直径在70 mm 到80 mm 之间时平均值变大,持续增大到80 mm以上平均值开始变小,缓冲区域直径持续增大会导致抛光效率下降,因为流场强度不变,但是流场空间变大导致流场各处的湍动能下降;(2)缓冲区域长度以60 mm 为拐点,左侧呈负相关,右侧呈正相关;(3)入口区域长度对平均值的影响呈三段折线关系,但是影响并不显著;(4)随着入口流速的变大均值变大,但是在入口流速值为6 m/s 到7 m/s之间对平均值的影响比较显著,大于7 m/s 之后影响就很小,说明持续增大入口流速并不能一直提高抛光效率;(5)入口直径对平均值的影响也是呈正相关的趋势,并且影响较为显著。

表8 速度与压强乘积平均值均值响应值(望大)

图12 速度与压强乘积均值响应图

4.3.3 参数对抛光后所有工件的速度与压力乘积的均匀程度的影响

从表9 可以得出参数对批量抛光工件离散程度影响的主次顺序为:入口直径＞缓冲区域直径＞入口流速＞缓冲区域长度＞入口区域长度,并且入口直径是显著影响离散度的因素。

表9 抛光离散度均值响应表(望小)

以抛光后所有工件的近壁处的磨粒相对速度与压强的乘积值的均方差与平均值作为衡量抛光效果的两个指标,均方差与抛光效果呈负相关,均值与抛光效果呈正相关,定义其比值α为离散度。根据统计学标准化效应t统计量来检验效应为0 的原假设,绘制出参考线来指示哪些效应在统计意义上显著,从而绘制出Pareto 图,如图13 所示。可以分析出影响显著的交互作用,综合图14 的单因素分析以及交互作用图15,选用以下参数为最优组合:缓冲区域直径80 mm,缓冲区域长度40 mm,入口区域长度40 mm,入口流速7.5 m/s,入口直径40 mm。

图13 离散度的Pareto 图

图14 离散度响应图

图15 离散度交互作用图

5 结论验证与分析

根据上文分析出的参数最优组合,对该组参数进行仿真分析,结果表明所有工件的PV 值乘积的均方差和均值分别为7.08 MPa·m/s 和561.5 MPa·m/s,离散度为1.26%,整体而言所有工件的抛光均匀性较好,同时每个工件的PV 平均值比较大,抛光效率良好。各个工件之间的材料去除量相差不大,实现了抛光均匀的要求。图16 是整个夹具中PV 值的分布图,可以看到在经过缓冲区的缓冲后,流场的PV 值已经达到了比较均匀的程度,只有在圆心处局部高于其他位置,并且相差并不大。而且经过引流区后,在实际的加工区域PV 值基本一致,说明了在当前夹具的几何尺寸参数下抛光均匀效果有明显改善。从图17 冲蚀速率图可以看出,整体冲蚀速率良好,且切削量较大处都位于工件的孔壁表面,达到了良好的抛光效果。

图16 流场PV 值分布图

图17 冲蚀速率

以下从单个工件抛光的整体均匀性来加以验证,图18 是分布在夹具中沿着径向的一批工件孔壁PV 值分布图,图19 是所有工件沿着任意母线的PV值线图。从分布图18 中可以看出,在引流区域与缓冲区域交界处由于从大口径流道突然缩口到小口径流道导致引流区域的流场不均匀,而在加工区域PV值则大致是均匀的。从图19 可以得出,由于流道缩口以及沿程压力损失,在工件入口处压力较大,随着流体的流动压力逐渐变小,压力的分布趋势与速度相反,但是在经过引流区的长度以后,PV 值沿着轴向达到了比较均匀的程度,实现了均匀抛光批量工件的目标。

图18 工件轴向PV 值分布图

图19 工件PV 值线图

6 结论

本文提出一种较为有效的还原粒子在湍流场中的运动以及计算颗粒冲蚀效果的湍流弥散模型,仿真结果说明了湍流弥散模型对于还原粒子在流场中的无序运动以及预测材料冲蚀量的有效性。通过仿真实验设计了一种可以实现给定工件几何尺寸批量抛光的夹具,仿真结果得出了夹具几何参数对于抛光工件均匀性影响的主次,同时得到了当前工件批量抛光时夹具的最佳几何尺寸以及入口速度。