基于VOF的模具结构化表面软性磨粒流数值模拟

计时鸣 唐 波 谭大鹏

浙江工业大学特种装备制造与先进加工技术教育部／浙江省重点实验室，杭州，310014

0 引言

为了消除模具表面所残留的机械加工痕迹，光整加工技术成为必要的工艺环节，光整加工时间占整个模具制造时间的50%以上，现有方法一般需借助工具接触或靠近待加工表面进行加工[1－2]。在光学元件加工领域，将凹槽、棱柱、棱镜阵列等立体结构的表面称为结构化表面[3－4]，借鉴此称呼，将模具中沟、槽、孔、棱柱、棱锥、窄缝等复杂异型面统一称为结构化表面。涉及这些结构化表面精密光整加工技术的研究却比较薄弱。

磨粒流加工是一种较新的表面加工方法，该加工方法可通过去毛刺、抛光、倒圆角并去除由电火花加工或激光加工后的再铸层来提高表面质量[5－7]。但是，现有的磨粒流加工方法一般不能直接应用于模具结构化表面的光整加工。针对模具制造中结构化表面光整加工的技术难题，笔者提出了一种基于软性磨粒流（即弱黏性或无黏性液固两相磨粒流）的模具结构化表面无工具精密光整加工新方法，即利用软性磨粒流在约束流道内的湍流壁面效应实现对封闭流道内模具结构化表面的微力微量切削，完成对结构化表面无工具化精密光整加工。

计算流体力学（computational fluid dynamics，CFD）方法克服了传统理论的方法与实验方法的弱点。本文基于流体体积（volume of fluid，VOF）模型和标准k－ε模型相结合的方法，以之字形微型流道为例，揭示其内部的流动规律，为应用软性磨粒流实现模具结构化表面精密加工提供理论依据。数值模拟方法不仅具有成本低、速度快等较显著的优点，甚至还能解决由于实验技术所限难以进行测量的问题[8]。

1 软性磨粒流的数学模型

VOF模型是基于界面模拟的现代运动界面追踪方法之一，模型适用于两种或多种互不穿透流体间界面的跟踪计算。模型对每一相引入体积分数变量，通过求解每一控制单元内体积分数值确定相间界面。

1.1 控制方程

设多相流在计算网格内的第q相流体的体积分数为αq。在每个控制体积内，所有项的体积分数之和为1。当αq＝1，则说明该第q相流体单元全部为指定相流体所占据；若αq＝0，则说明该单元全部为另一相流体所占据，相对于另一相流体则称为空单元；当αq＜1时，则说明该单元包含了第q相流体和另一相或者其他多项流体的界面。基于αq的局部值、适当的属性和变量在一定范围内分配给每一控制体积。流体函数的输运方程和流体运动控制方程一起构成VOF模型的控制方程[9]。

1.1.1 体积分数方程

在VOF模型中，跟踪相与另一相之间的界面是通过求解一相或多相的体积分数的连续方程来完成的，对第q相（可取值为1，2），方程如下：

式中，αq为第q相的体积分数；νq为q相的速度；ρq为第q相的物理密度；Sαq为αq的源项；为p相到q相的质量输送；为q相到p相的质量输送。

1.1.2 动量方程

动量方程是任何流动系统必须满足的基本定律，流体体积模型的动量方程取决于通过属性ρ和μ的所有相的体积分数。动量方程如下：

式中，ρ为软性磨粒流的密度，ρ＝α2ρ2＋α1ρ1；α1为液体相的体积分数；α2为固体颗粒相的体积分数；ρ1为液体相的密度；ρ2为固体颗粒相的密度；p为流体微元体积上的压力（静压）；μ为软性磨粒流的黏度；g为重力加速度；F为体积力。

1.1.3 附加的标量方程

为了能消除表面纹理的方向性，须使软性磨粒流的运动处于湍流状态，湍流状态下的软性磨粒流运动方向具有随机性，故采用了标准k－ε的湍流模型[10-11]。

湍动能k输运方程如下：

式中，t为时间；xi、xj分别为x、y、z三个方向的坐标；ui为x、y、z三个方向的速度分量；μt为湍流黏度系数，μt＝ρ为由于平均速度梯度引起的湍动能k的产生项；σk为湍动能k 的湍 流普朗 特数，σk＝1.0[12-13]。

由于软性磨粒流的不可压缩，且不考虑软性磨粒流的浮力、自定义的源项，则由浮力引起的湍动能k的产生项Gb＝0，可压湍流中脉动扩张的贡献项YM＝0，用户定义的源项Sk＝0。

耗散率ε输运方程如下：

式中，uj为x、y、z三个方向的速度分量；σε为耗散率ε的湍流普朗特数，σε＝1.3[12-13]。

1.2 边界条件

为了揭示之字形微型流道内软性磨粒流流动规律，以FLUENT 6.3软件为计算平台，采用VOF和标准k－ε模型相结合的方法对其进行了数值模拟。模拟中采用的边界条件是速度入口边界条件和压力出口边界条件。

模具结构化表面区域在求解过程中考虑了边壁粗糙影响的壁面函数[14]，方程为

式中，uP为近壁面网格点P的速度；u＊为摩阻流速；yP为网格点P到模具结构化表面的距离；κ为卡门常数，κ＝0.42；E为表征粗糙率的参数，对光滑壁面取E为9.8；ΔB为模具结构化表面的粗糙修正函数，其为模具结构表面的粗糙常数，取值在0到1之为模具结构化表面平均粗糙高度。

计算中取粗糙常数Cks为0.5，模具结构化表面平均粗糙高度为0.003。

2 数值模拟及其结果分析

2.1 计算模型及初始条件

为了能使熔融的塑料顺利光滑地进入模具内部，采用了圆形的微型流道，流道的结构示意图如图1所示，图1左上方为软性磨粒流的入口处，右下方为软性磨粒流的出口处。结合模具微型流道的实际情况确定了之字形微型流道结构尺寸如表1所示。

图1 之字形微型流道的结构示意图

表1 之字形微型流道结构尺寸mm

为了研究之字形微型流道内部流场的流动特性，对不同直径的微型流道和同尺寸不同速度的微型流道其内的软性磨粒流的流动形态进行分析，由于流道尺度特征和结构简单，单精度求解器就能很好地满足计算精度要求，采用了2D单精度非稳态和3D单精度非稳态隐式分离求解器进行求解。模拟中使用的入口速度参数和与此对应的雷诺数如表2所示，流体产生湍流的基本条件是雷诺数达到一个临界值（一般雷诺数Re＜2000为层流状态，Re＞4000为湍流状态，Re＝2000～4000为过渡状态），因此表2所示的软性磨粒流处于湍流状态。

表2 5组不同的速度以及对应的雷诺数Re

模拟中使用的其他仿真工艺和物性参数如下：模具材料为模具钢；磨粒为碳化硅；磨粒与流体（46号机械润滑油）体积比为1∶9；磨粒流的动力黏度为2.064Pa·s。由于磨粒是粉末状且其颗粒足够细，因此假设磨粒为连续的流体相。将46号机械润滑油考虑为第一相，磨粒碳化硅为第二相，求解方程组时，磨粒体积分数在10%以下，可以忽略颗粒与颗粒之间的相互作用、颗粒体积分数对连续相的影响。

2.2 网格划分

网格生成技术是CFD的重要内容，网格质量对CFD计算精度和计算效率有重要影响，因此，有必要对网格生成方式给予足够的关注。软性磨粒流为不可压缩的流体且流道是圆形，软性磨粒流在沿直径截面的运动规律分布是一致，采用了二维轴对称几何结构模拟，计算区域和网格划分由前处理器软件GAMBIT来完成。之字形微型流道壁面处采取边界层网格，使网格的划分更加合理。具体网格划分如图2所示，网格数量为2620个。

为了直观地展现软性磨粒流在流道内部的流动轨迹，对不同结构尺寸的流道采用三维数值模拟。由于之字形微型流道几何结构分布不均匀，因此分段划分流道的网格同时采用边界层网格技术，使网格更贴近实际。直径为4mm微型流道的网格划分示意图如图3所示，其网格数量为22 026个。

图2 之字形微型流道平面网格划分示意图

图3 之字形微型流道三维网格划分示意图

2.3 仿真结果与讨论

根据材料去除的类型，一般分为脆性去除和塑性去除两大类[15－16]。脆性去除时，材料的去除是通过裂纹的形成及扩展引起材料脆性破裂或材料压碎的方式来实现的，塑性去除一般是有剪切作用造成的。由于软性磨粒流抛光属于低速加工[17]，因此磨粒流流过模具结构化表面时只发生塑性去除。

软性磨粒流的流型受两相流的组分、边界条件、流速等因素影响，因此将这些参数作为调控因子，实现对软性磨粒流的控制。改变流速和直径大小来模拟微型流道内软性磨粒流的流动情况，通过比较内部流场的湍流强度、湍流耗散率以及磨粒流的速度分布等影响湍流流动的参数，为合理选取模具结构化表面软性磨粒流精密加工的参数提供依据。

不同直径流道内软性磨粒流的轴向速度轨迹如图4、图5所示，之字形微型流道内部磨粒流各处的速度分布不均匀，中心处的速度比近壁面处的速度高，越靠近壁面磨粒流所受的阻碍就越大，同时各磨粒流流层存在摩擦切应力，速度较快的磨粒流通过切应力拖动相临的低速磨粒流向前运动，而速度低的磨粒流会阻止高速磨粒流向前运动；随着流道尺寸的增大，使得单位时间内磨粒流与壁面相互作用的有效磨粒数减小，边界层磨粒与壁面的剪切能力下降，加工效率将会降低。

图4 直径2mm时内部轴向速度分布图

在之字形微型流道内部软性磨粒流的加工过程中，径向速度具有重要的作用。不仅是因为径向速度促使磨粒靠近壁面，而且径向速度对表面的粗糙度有很重要的影响。图6、图7给出了不同直径出口处径向速度的分布图，从图6、图7可以看出，径向速度在不同入口速度时分布是不一致的，同一截面上向两边的减小程度也不同，当完全靠近壁面时，径向速度变为零。由于模具结构化表面的限制和流体的附着力的作用，近壁区其流速很小，速度梯度大，脉动消失，在距离模具结构化表面稍远处壁面对磨粒流质点的影响减弱。因此，为了保证模具结构化表面加工的均匀性，根据加工情况调整软性磨粒流的出入口方向，磨粒流从入口加工一段时间后，换成从出口方向作为入口进行加工，这样有利于模具结构化表面的表面形貌取得一致。

图5 直径4mm时内部轴向速度分布图

图6 直径2mm时出口速度径向分布图

图7 直径4mm时出口速度径向分布图

图8～图10分别给出了之字形微型流道入口直径为2mm时不同入口速度其内部的轴向速度、湍流耗散率和湍动能比较图。

图8所示的之字形微型流道不同入口速度的流道内部轴向速度比较表明：随着入流速度的增加，流道内软性磨粒流的速度也随之增加，磨粒流与壁面碰撞的频率也随之增大；当速度大于10m／s时，速度在两个拐角处的波动幅度较大，是由于在拐角处磨粒流与壁面以及磨粒流与磨粒流之间质量、动量和能量发生了剧烈的变化和交换。

图8 流道内部流场速度比较图

图9 流道内部流场湍流动能比较图

图10 流道内部流场湍流耗散率比较图

图9所示的之字形微型流道湍流动能比较表明：随着入口速度的增大，流道内部软性磨粒流的湍动能变化幅度也随之增大，湍动能的波动比速度波动幅度大；速度梯度变化大的地方湍动能变化也大；当入口湍动能为0.4m2／s2时，磨粒流的全程运动过程中，湍动能的波峰值与波谷值之间的差值变化较小，能提供较平稳的加工条件，有利于模具结构化表面取得一致的形状精度。

图10所示的之字形微型流道湍流耗散率比较表明：随着入流速度的增加，湍流耗散率也随之增加，湍流耗散率的波动比速度波动幅度大；由于流道中存在两个拐角，从图10可以看出第一个峰值出现位置到第二个峰值出现位置之间湍流耗散率变化较大，入口处到第一个峰值位置以及第二个峰值位置到出口处湍流变化较平稳；由于湍流耗散率与涡旋尺度存在一一对应关系，耗散率越小，涡旋尺度越大，当入口湍流耗散率为19.8 m2／s3时，湍流耗散率较小，变化范围小，涡流尺度就大，从入口处到出现第一个峰值的位置，从第一个峰值位置变化到第二个峰值的位置以及第二个峰值到出口变化的较为平稳，有利于工件取得较一致的表面形貌。

3 加工试验

搭建了软性磨粒流加工试验平台，如图11所示。平台主要由电气控制柜、溶液箱、隔膜泵、试验工作台四大部分及一整套管路阀门及附属部件组成。磨粒流存储在溶液箱中，通过隔膜泵驱动对待加工试样进行加工，然后通过管道返回溶液箱，这样既保证了磨粒流加工的连续性，也减小了磨粒流用量，具有较好的经济性。

图11 “软性”磨粒流精密光整加工试验平台

利用上面提到软性磨粒流加工实验平台上采用模具钢试件，试验过程中由搅拌棒全程搅拌，保证接近于悬浮液，磨料和液体混合良好。根据以上的仿真结果，取流道直径为2mm、入口速度为10m／s、湍动能为0.4m2／s2和湍流耗散率为19.8m2／s3。

用美国维易科精密仪器有限公司的Veeco NT9000光学轮廓仪检测工件的粗糙度，加工前后工件微观粗糙度对比情况如图12所示，加工前微观粗糙度为93.28nm；加工后微观粗糙度为42.24nm。

4 结论

（1）采用VOF模型和标准k－ε模型相结合的计算方法，应用Fluent软件对之字形微型流道的软性磨粒流流场进行了数值模拟，数值模拟结果为磨粒流加工过程中的工艺参数选择提供了参考依据。

图12 加工前后工件微观粗糙度对比

（2）采用VOF多相流模型进行的软性磨软性磨粒流数值模拟很好地揭示了之字形微型流道内的软性磨软性磨粒流的流动情况。数值模拟可以形象地再现流体的流动情景，甚至还能解决由于实验技术所限难以进行测量的问题。

（3）模拟结果表明，随着入口流速的提高，管道内的磨粒流的平均速度也随之增大，有利于近壁区磨粒流与流道壁面的相互作用，从而提高磨粒流的加工效率。同时模具流道的形状和结构对磨粒流加工也有重要影响。当微型流道直径为2mm、入口速度为10m／s、湍动能为0.4m2／s2、湍流耗散率为19.8m2／s3时，与其他几组参数工艺参数相比：微型流道内轴向、切向速度分布及湍动能分布等更加均匀，能够得到较高的加工效率和表面纹理形态。

（4）通过试验验证了数值模拟的结果，数值模拟的结果能指导软性磨粒流加工参数的合理选取，数值模拟手段为深入研究软性磨粒流的基本规律提供了一种理论工具。

[1]Shiou F J，Ciou H S.Ultra－precision Surface Finish of the Hardened Stainless Mold Steel Using Vibration－assisted Ball Polishing Process[J].International Journal of Machine Tools and Manufacture，2008，48（7／8）：721-732.

[2]Gianpaolo S，Roberto M，Gianmaria C.A Surface Rough－ness Predictive Model in DeterministicPolishing of Ground Glass Moulds[J].International Journal of Machine Tools and Manufacture，2009，49（1）：1－7.

[3]Brinksmeier E.Polishing of Structured Molds[J].CIRP Annals－Manufacturing Technology，2004，53（1）：247－250.

[4]Brinksmeier E.Finishing of Structured Surfaces by Abrasive Polishing[J].Precision Engineering，2006，30（3）：325－336.

[5]Rhoads L J.Abrasive Flow Machining：a Case Study[J].Journal of Materials Processing Technology，1991，28（2）：107－116.

[6]Rhoads L J.Abrasive Flow Machining[J].Manufacturing Engineering，1988，101：75－78.

[7]Jain R K，Jain V K，Dixit P M.Modling of Material Removal and Surface Roughness in Abrasive Flow Machining Process[J].International Journal of Machine Tools ＆ Manufacture，1999，39（12）：1903－1923.

[8]Patankar S V.Numerical Heat Transfer and Fluid Flow[M].New York：Megraw－Hil，1980.

[9]江帆，黄鹏.Fluent高级应用与实例分析[M].北京：清华大学出版社，2008.

[10]Verseeg H K，Malalasekera W.An Introduction to Computationl Fluid Dynamics：the Finite Volume Method[M].New York：Wiley，1995.

[11]陶文铨.数值传热学[M].2版.西安：西安交通大学出版社，2001.

[12]Annaland M V S，Deen N G，Kuipers J A M.Multi－level Computational Fluid Dynamics Models for the Description of Particle Mixing and Granulation in Fluidized Beds[J].Handbook of Powder Technology，2007，11：1071－1107.

[13]Du Wei，Bao Xiaojun，Xu Jian，et al.Computational Fluid Dynamics（CFD）Modeling of Spouted Bed：Influence of Frictional Stress，Maximum Packing Limit and Coefficient of Restitution of Particles[J].Chemical Engineering Science，2006，61（14）：4558－4570.

[14]李志勤，李洪，李嘉，等.溢流丁坝附近自由水面的实验研究与数值模拟[J].水利学报，2003，34（8）：53－57.

[15]Booij S M，Van B H，Braat J J M.Nanometer Deep Shaping with Fluid Jet Polishing[J].Optical Engineering，2002，41（8）：1926－1931.

[16]方惠，郭培基，余景池.液体喷射抛光材料去除机理的研究[J].光学技术，2004，30（2）：248－250.

[17]李兆泽，李圣怡，戴一帆，等.磨料射流抛光中各工艺参数对材料去除率及抛光区形貌的影响[J].中国机械工程，2008，19（21）：2532－2535.