基于流体体积模型模具微细流道磨料流精密加工控制方法

摘要

一种基于流体体积模型模具微细流道磨料流精密加工控制方法，包括以下步骤：(1)基于流体体积模型建立磨料流运动的数学模型；(2)磨料流的求解过程，具体包括：2.1)始条件和边界条件；2.2)网格的划分；2.3)建立离散方程；2.4)离散初始条件和边界条件；2.5)给定求解控制参数；2.6)求解离散方程；2.7)断解的收敛：对于稳态问题的解，或是瞬态问题在某个特定时间步上的解，通过多次迭代得到；对于瞬态问题，若采用显式格式进行时间域上的积分得到；在求解值达到指定精度后，结束迭代过程。本发明能够实现磨料流形的观察、成本低、速度快。

说明书

技术领域

本发明涉及模具精密光整加工领域，更具体的说，是涉及塑料模具凹模和凸模中细小结构精密加工的分析方法。

背景技术

随着科技的进步和产品的精密化，模具应用领域对注塑件表面对应的模具面以及模具流道的粗糙度要求越来越高。为了消除模具表面所残留的机械加工痕迹，光整加工技术成为必要的工艺环节，占整个模具制造时间的50％以上，现有方法一股需借助工具接触或靠近待加工表面进行加工。但面向模具制造中所涉及细小结构的精密光整加工技术的研究，对于尺寸小或几何形态特殊的表面难以使用工具进行接触式光整加工，这一问题目前尚无有效方法解决。

磨料流加(abrasive flow machining，AFM)是一种较新的表面加工方法，可以提高表面质量，去毛刺，抛光、倒圆并可去除由电火花加工或激光加工后的再铸层。但是，现有的磨粒流加工方法一股不能直接应用于模具结构化表面的光整加工或尚不能获得满意的表面粗糙度。其中磁流变光整、磁射流光整和电流变液光整等方法仍须借助工具。磨粒水射流光整加工由于射流的直线特性使得它难以均匀地冲刷到结构化表面的各个部位。挤压珩磨使用强黏性磨粒流强力挤压进行加工会破坏结构化表面，另外其磨粒流的运动形式使表面纹理呈规律性，限制了其用于镜面级表面粗糙度的超精密加工。

发明内容

为了克服已有微细流道磨料流精密加工控制方法的难以进行磨料流形的观察、成本高、速度快的不足，本发明提供一种能够实现磨料流形的观察、成本低、速度快的基于流体体积模型模具微细流道磨料流精密加工控制方法。

本发明解决其技术问题所采用的技术方案是：

一种基于流体体积模型模具微细流道磨料流精密加工控制方法，所述模具微细流道磨料流精密加工控制方法包括以下步骤：

(1)、建立磨料流运动的数学模型，控制方程组如下所示：

1.1)体积分数方程：跟踪相与相之间的界面是通过求解一相或多相的体积分数的连续方程来完成的，对第q相方程如下：

$\frac{{\partial \alpha }_q}{\partial t}+v_q▿{\alpha }_q=\frac{S_{{\alpha }_q}}{{\rho }_q}+\frac{1}{{\rho }_q}\underset{p=1}{\overset{2}{\Sigma }}({\stackrel{·}{m}}_{\mathrm{pq}}-{\stackrel{·}{m}}_{\mathrm{qp}})---\left(1\right)$

式中：α_q是第q相的体积分数；v_q是q相的速度；ρ_q是第q相的物理密度；是α_q的源项；是p相到q相得质量输送，是q相到p相得质量输送；主相体积分数的计算基于如下的约束：$\underset{q=1}{\overset{2}{\Sigma }}{\alpha }_q=1;$

1.2)属性方程：如果第二相的体积分数被跟踪，那么每一单元中的密度由式：

ρ＝α₂ρ₂+α₁ρ₁                                    (2)

1.3)动量方程：

$\frac{\partial }{\partial t}\left(\rho \overrightarrow{v}\right)+▿·\left(\rho \overrightarrow{v}\overrightarrow{v}\right)=-▿p+▿·\left[\mu (▿\overrightarrow{v}+▿{\overrightarrow{v}}^T)\right]+\rho \overrightarrow{g}+\overrightarrow{F}---\left(3\right)$

式中：ρ为体积分数平均密度；p是流体微元体积上的压力(静压)；μ为物质的粘性；为重力加速度；为体积力；

1.4)能量方程：

能量方程是各相共享的，表示如下：

$\frac{\partial }{\partial t}\left(\mathrm{\rho E}\right)+▿\left[\overrightarrow{v}(\mathrm{\rho E}+p)\right]=▿·(k_{\mathrm{eff}}▿T)+S_h---\left(4\right)$

VOF模型处理能量E和温度T，作为质量平均变量：

$E=\frac{\underset{q=1}{\overset{n}{\Sigma }}{\alpha }_q{\rho }_q{E}_q}{\underset{q=1}{\overset{n}{\Sigma }}{\alpha }_q{\rho }_q}---\left(5\right)$

每一相的E_q是基于该相的比热和共享温度；

1.5)附加的标量方程：

采用湍流流动加工，采用了SST湍流模型，湍动能方程k为：

$\frac{\partial \mathrm{\rho k}}{\partial t}+\mathrm{div}[\mathrm{\rho vk}-](\mu +\frac{{\mu }_t}{{\rho }_{k1}}▿k)=G-{\beta }^{\prime }\mathrm{\rho k\omega }---\left(6\right)$

湍流频率ω方程为：

$\frac{\partial \mathrm{\rho \omega }}{\partial t}+\mathrm{div}[\mathrm{\rho v\omega }-(\mu +\frac{{\mu }_t}{{\sigma }_{\omega 1}})▿\omega ]={\alpha }_1\omega \frac{G}{k}-{\beta }_1{\mathrm{\rho \omega }}^2---\left(7\right)$

使用选择器

$\bar{G}=\min (G,c_{\lim }ϵ)c_{\mathrm{lin}}=10---\left(8\right)$

将k-ε湍流模型进行变形：

$\frac{\partial \mathrm{\rho k}}{\partial t}+\mathrm{div}[\mathrm{\rho vk}-(\mu +\frac{{\mu }_t}{{\sigma }_{k2}})▿k]=G-{\beta }^{\prime }\mathrm{\rho k\omega }---\left(9\right)$

$\frac{\partial \mathrm{\rho \omega }}{\partial t}+\mathrm{div}[\mathrm{\rho v\omega }-(\mu +\frac{{\mu }_t}{{\sigma }_{\omega 2}})▿\omega ]={\alpha }_2\omega \frac{G}{k}-{\beta }_2{\mathrm{\rho \omega }}^2+\frac{2\rho {\sigma }_{\omega 2}}{\omega }▿k▿\omega ---\left(10\right)$

将基本的k-ω模型、变形后的k-ε湍流模型分别乘以函数F₁和1-F₁，得到SST湍流模型：

Φ₃＝F₁Φ₁+(1-F₁)Φ₂                                        (11)

式中Φ₁为基本的k-ω模型；

Φ₂为变形的k-ε湍流模型；

$F_1=\tanh \left({\arg }_1^4\right)$

${\arg }_1=\min [\max \left(\frac{\sqrt{k}}{{\beta }^{\prime }\omega y^{,}}\frac{500\upsilon }{y^2\omega }\right),\frac{4\rho {\sigma }_{\omega 2}k}{{\mathrm{CD}}_{\mathrm{k\omega }}{y}^2}]$

${\mu }_t=\frac{\rho a_{1}k}{\max (a_1\omega ,{\mathrm{SF}}_2)}$

$F_2=\tanh \left({\arg }_2^2\right)$

${\arg }_2=\max \left(\frac{2\sqrt{k}}{{\beta }^{\prime }\omega y^{,}}\frac{500\upsilon }{y^2\omega }\right)$

y为近壁节点最近距离；

υ为运动粘度

常数及经验系数如下

β′＝0.09    α₁＝5/9     σ_k1＝2

β₁＝0.075    α₂＝0.44    σ_ω1＝2

β₂＝0.0828   σ_k2＝1      σ_ω2＝0.856；

(2)、磨料流的求解过程，具体包括：

2.1)始条件和边界条件；

2.2)、网格的划分：不同的计算区域采用结构化网格和非结构化网格相结合的方法来划分网格；

2.3)、建立离散方程，通过数值的方法把计算区域内有限数量位置上的因变量当作未知量来处理，从而建立一组关于未知量的代数方程，然后通过求解代数方程组来得到节点值，而计算区域内的其他位置上的值则根据节点位置上的值来确定；

2.4)、离散初始条件和边界条件：如在静止壁面上速度为0，针对所生成的网格，将连续型的初始条件和边界条件转化为特定节点上的值；

2.5)、给定求解控制参数：在离散空间上建立离散化的代数方程组，并施加离散化的初始条件和边界条件，给定流体的物理参数和湍流模型的经验系数、迭代计算的控制精度、瞬态问题的时间步长和输出频率；

2.6)、求解离散方程：对生成的具有定解条件的代数方程组进行求解；

2.7)、断解的收敛：对于稳态问题的解，或是瞬态问题在某个特定时间步上的解，通过多次迭代得到；对于瞬态问题，若采用显式格式进行时间域上的积分得到；在求解值达到指定精度后，结束迭代过程。

本发明的技术构思为：将计算流体力学中基于流体体积模型(Volume Of Fluid，VOF)和SST湍流模型相结合的的方法引入到模具细小结构精密加工中，能够实现流道内磨料流流形的观察。该方法不仅具有成本低、速度快等较显著的优点，甚至还能解决由于实验技术所限难以进行测量的问题。

采用流体体积模型和SST湍流模型相结合的的方法引入到模具细小结构精密加工领域，主要由模型的建立过程、求解过程和后期分析过程。

所述的模型的建立过程是基于磨粒流的二相流基本理论和连续介质理论，因没有热交换发生，则直接将连续方程与动量方程作为控制方程使用。由于磨料流的流动大多是处于湍流范围，一股情况下，需要增加湍流方程。因此建立起了基于流体体积模型和SST湍流模型相结合的的控制方程。建立控制方程，是求解任何问题前都必须首先进行的。

所述的求解过程，作为一个数值求解方案；该过程主要包括初始条件和边界条件的确定、网格的划分、建立离散方程、离散初始条件和边界条件、给定求解控制参数、求解离散方程、判断解的收敛。

所述的后期分析过程，是显示和输出计算结果。通过上述求解过程得出了各计算节点上的解后，需要通过适当的手段将整个计算域上的结果表示出来。这时，我们可采用线值图、矢量图、等值线图、流线图、云图等方式对计算结果进行表示。

与现有技术相比，本发明的有益效果是：1)该方法具有成本低、速度快等较显著的优点，还能解决由于实验技术所限难以进行测量的问题。2)该方法能实现单向流动循环加工，提高磨料流的加工效率。3)该方法具有重复性和经济性，只需更改参数就可以完成所要求的目的。

附图说明

图1是数学模型框图；

图2是求解过程示意图；

图3是具体实例模型示意图；

图4是网格划分示意图；

图5是具体实例出口径向速度示意图；

具体实施方式

下面结合附图对本发明作进一步描述。

参照图1～图5，一种基于流体体积模型模具微细流道磨料流精密加工控制方法，包括以下步骤：

(1)、建立磨料流运动的数学模型

磨料流运动没有正确完善的数学模型，仿真模拟就缺少了理论依据，基于磨粒流的二相流基本理论和连续介质理论，因没有热交换发生，则直接将连续方程与动量方程作为控制方程使用。由于磨料流的流动大多是处于湍流范围，一股情况下，需要增加湍流方程。因此建立起了基于流体体积模型和SST湍流模型相结合的的控制方程(见附图1)，VOF模型是用体积率函数表示流体自由面的位置和流体所占的体积，其方法占内存小，是一种简单而有效的方法，VOF所涉及的相关流体控制方程有：体积分数方程；属性方程；动量方程；能量方程；附加的标量方程。SST模型虽然与标准模型、RNG模型一样都采用了各向同性湍流的假定，但它通过适当修正湍流模型可用于计算近壁区(壁函数)的粘性内层，因而能适合各种压力梯度变化的流场模拟；综合了k-ε和k-ω湍流模型在边界层内外计算的优点，能够准确及时预测分离的特性。具体的说就是要建立反映问题各个变量之间的关系方程以及相应的初始条件，这个就是数值模拟的出发点。具体控制方程组如下所示：

1.1)体积分数方程

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

$\frac{{\partial \alpha }_q}{\partial t}+v_q▿{\alpha }_q=\frac{S_{{\alpha }_q}}{{\rho }_q}+\frac{1}{{\rho }_q}\underset{p=1}{\overset{2}{\Sigma }}({\stackrel{·}{m}}_{\mathrm{pq}}-{\stackrel{·}{m}}_{\mathrm{qp}})---\left(1\right)$

式中：α_q是第q相的体积分数；v_q是q相的速度；ρ_q是第q相的物理密度；是α_q的源项；是p相到q相得质量输送，是q相到p相得质量输送。体积分数方程(the volume fraction equation)不是为主相求解，主相体积分数的计算基于如下的约束：$\underset{q=1}{\overset{2}{\Sigma }}{\alpha }_q=1.$

1.2)属性方程

出现在输运方程中的属性(Properties)是由于存在于每一控制体积中的分相决定的，在两相流系统中，如果用下标1和2表示，如果第二相的体积分数被跟踪，那么每一单元中的密度由式：

ρ＝α₂ρ₂+α₁ρ₁                                        (2)

1.3)动量方程

$\frac{\partial }{\partial t}\left(\rho \overrightarrow{v}\right)+▿·\left(\rho \overrightarrow{v}\overrightarrow{v}\right)=-▿p+▿·\left[\mu (▿\overrightarrow{v}+▿{\overrightarrow{v}}^T)\right]+\rho \overrightarrow{g}+\overrightarrow{F}---\left(3\right)$

式中：ρ为体积分数平均密度；p是流体微元体积上的压力(静压)；μ为物质的粘性；为重力加速度；为体积力。

1.4)能量方程

能量方程(the energy equation)也是各相共享的，表示如下：

$\frac{\partial }{\partial t}\left(\mathrm{\rho E}\right)+▿·\left[\overrightarrow{v}(\mathrm{\rho E}+p)\right]=▿·(k_{\mathrm{eff}}▿T)+S_h---\left(4\right)$

VOF模型处理能量E和温度T，作为质量平均变量：

$E=\frac{\underset{q=1}{\overset{n}{\Sigma }}{\alpha }_q{\rho }_q{E}_q}{\underset{q=1}{\overset{n}{\Sigma }}{\alpha }_q{\rho }_q}---\left(5\right)$

这里对每一相的E_q是基于该相的比热和共享温度

1.5)附加的标量方程

依赖于需要求解的问题，求解时需要涉及附加的标量方程(addition scalarequation)。由于采用湍流流动加工，采用了SST湍流模型，

湍动能方程k为

$\frac{\partial \mathrm{\rho k}}{\partial t}+\mathrm{div}[\mathrm{\rho vk}-](\mu +\frac{{\mu }_t}{{\rho }_{k1}}▿k)=G-{\beta }^{\prime }\mathrm{\rho k\omega }---\left(6\right)$

湍流频率ω方程为

$\frac{\partial \mathrm{\rho \omega }}{\partial t}+\mathrm{div}[\mathrm{\rho v\omega }-(\mu +\frac{{\mu }_t}{{\sigma }_{\omega 1}})▿\omega ]={\alpha }_1\omega \frac{G}{k}-{\beta }_1{\mathrm{\rho \omega }}^2---\left(7\right)$

为了避免在流域滞止区湍动能的错误计算，使用了选择器

$\bar{G}=\min (G,c_{\lim }ϵ)c_{\mathrm{lin}}=10---\left(8\right)$

为综合k-ω模型在近壁区模拟和k-ε湍流模型在外部区域计算的优点，将k-ε湍流模型进行变形

$\frac{\partial \mathrm{\rho k}}{\partial t}+\mathrm{div}[\mathrm{\rho vk}-(\mu +\frac{{\mu }_t}{{\sigma }_{k2}})▿k]=G-{\beta }^{\prime }\mathrm{\rho k\omega }---\left(9\right)$

$\frac{\partial \mathrm{\rho \omega }}{\partial t}+\mathrm{div}[\mathrm{\rho v\omega }-(\mu +\frac{{\mu }_t}{{\sigma }_{\omega 2}})▿\omega ]={\alpha }_2\omega \frac{G}{k}-{\beta }_2{\mathrm{\rho \omega }}^2+\frac{2\rho {\sigma }_{\omega 2}}{\omega }▿k▿\omega$

(10)

将基本的k-ω模型、变形后的k-ε湍流模型分别乘以函数F₁和1-F₁，得到SST湍流模型

Φ₃＝F₁Φ₁+(1-F₁)Φ₂

(11)

式中Φ₁-基本的k-ω模型

Φ₂-变形的k-ε湍流模型

$F_1=\tanh \left({\arg }_1^4\right)$

${\arg }_1=\min [\max \left(\frac{\sqrt{k}}{{\beta }^{\prime }\omega y^{,}}\frac{500\upsilon }{y^2\omega }\right),\frac{4\rho {\sigma }_{\omega 2}k}{{\mathrm{CD}}_{\mathrm{k\omega }}{y}^2}]$

${\mu }_t=\frac{\rho a_{1}k}{\max (a_1\omega ,{\mathrm{SF}}_2)}$

$F_2=\tanh \left({\arg }_2^2\right)$

${\arg }_2=\max \left(\frac{2\sqrt{k}}{{\beta }^{\prime }\omega y^{,}}\frac{500\upsilon }{y^2\omega }\right)$

y-近壁节点最近距离

υ-运动粘度

常数及经验系数如下

β′＝0.09    α₁＝5/9     σ_k1＝2

β₁＝0.075    α₂＝0.44    σ_ω1＝2

β₂＝0.0828   σ_k2＝1      σ_ω2＝0.856

(2)、磨料流的求解过程

作为一个数值求解方案；该过程主要包括初始条件和边界条件的确定、网格的划分、建立离散方程、离散初始条件和边界条件、给定求解控制参数、求解离散方程、判断解的收敛。

初始条件和边界条件，是控制方程有确定解的前提，控制方程与相应的初始条件、边界条件的组合构成对一个物理过程完整的数学描述。初始条件是所研究对象在过程开始时刻各个求解变量的空间分布情况。对于瞬态问题，必须给定初始条件。对于稳态问题，不需要初始条件。由于本发明是利用了湍流的壁面效应，是一个瞬态的问题，因此需要给定初始条件。边界条件是在求解区域的边界上所求解的变量或其导数随地点和时间的变化规律。对于任何问题，都需要给定边界条件。

网格的划分，采用数值方法求解控制方程时，都是想办法将控制方程在空间区域上进行离散，然后求解得到的离散方程组。要想在空间域上离散控制方程，必须使用网格。不同的问题采用不同数值解法时，所需要的网格形式是有一定区别的，但生成网格的方法基本是一致的。目前，网格分结构网格和非结构网格两大类。简单地讲，结构网格在空间上比较规范，如对一个四边形区域，网格往往是成行成列分布的，行线和列线比较明显。而对非结构网格在空间分布上没有明显的行线和列线。对于二维问题，常用的网格单元有三角形和四边形等形式；对于三维问题，常用的网格单元有四面体、六面体、三棱体等形式。在整个计算域上，网格通过节点联系在一起。针对本发明，在不同的计算区域采用了结构化网格和非结构化网格相结合的方法来划分网格，使网格的划分更加贴近工程实际。

建立离散方程，针对求解域内所建立的偏微分方程的复杂性、困难性，因此，就需要通过数值的方法把计算区域内有限数量位置上的因变量当作未知量来处理，从而建立一组关于这些未知量的代数方程，然后通过求解代数方程组来得到这些节点值，而计算区域内的其他位置上的值则根据节点位置上的值来确定。这样，用变量的离散分布近似解来代替定解问题精确的连续数据，可以预料，当网格点很密时，离散方程的解将趋近于相应微分方程的精确解。根据离散原理的不同，离散的方法有有限差分法、有限元法、有限体积法。有限差分法是将求解域划分为差分网格，用有限个网格节点代替连续的求解域，只考虑网格点上的数值而不考虑其他网格点之间如何变化，有限元法是基于变分原理和加权余量法，把计算区域划分为有限个互不重叠的单元，必须假定值符合网格点之间的变化规律，并将其作为近似解。有限体积法是基于因变量在有限大学的控制体积中的守恒原理，要求因变量的积分守恒对任意一组控制体积都得到满足，在整个计算区域内自然也就得到满足，但有限体积法在寻求控制体积的积分时，必须假定值在网格点之间的分布。因此针对不同的求解区域，采对微分方程中不同的项采取不同的插值函数。

离散初始条件和边界条件，前面所给定的初始条件和边界条件是连续性的，如在静止壁面上速度为0，现在需要针对所生成的网格，将连续型的初始条件和边界条件转化为特定节点上的值，如静止壁面上共有90个节点，则这些节点上的速度值应均设为0。

给定求解控制参数，在离散空间上建立了离散化的代数方程组，并施加离散化的初始条件和边界条件后，还需要给定流体的物理参数和湍流模型的经验系数等。此外，还要给定迭代计算的控制精度、瞬态问题的时间步长和输出频率等。

求解离散方程，在进行了上述设置后，生成了具有定解条件的代数方程组。对于这些方程组，数学上已有相应的解法，如线性方程组可采用Gauss消去法或Gauss-Seidel迭代法求解，而对非线性方程组，可采用Newton-Raphson方法。

判断解的收敛，对于稳态问题的解，或是瞬态问题在某个特定时间步上的解，往往要通过多次迭代才能得到。有时，因网格形式或网格大小、对流项的离散插值格式等原因，可能导致解的发散。对于瞬态问题，若采用显式格式进行时间域上的积分，当时间步长过大时，也可能造成解的振荡或发散。因此，在迭代过程中，要对解的收敛性随时进行监视，并在系统达到指定精度后，结束迭代过程(见附图2)。

应用基于流体体积模型模拟流道中的磨料流流动情况，以“之”字形微型管道为研究对象，对其内部的软性磨粒流流场进行数值模拟。具体操作过程如下(见附图3)：网格划分采用了结构化和非结构化相结合的方法(见附图4)，上端为磨料流的入口处，下端为磨料流的出口处，通过特定的夹具和接口构成封闭流道，使磨料流顺利流过模具细小表面(见附图5)显示了出口处径向速度示意图，从该图中可以清晰的看出磨料流在径向方向上速度值的分布。

最后，还需要注意的是，以上列举的仅是本发明的一个具体实施例子。显然，本发明不限于以上实施例，还可以有许多变形。本领域的普通技术人员能从本发明公开的内容直接导出或联想到的所有变形，均应认为是本发明的保护范围。