微织构磁性复合流体抛光的流场仿真建模分析

Flow Field Simulation and Modeling Analysis of Micro-Textured Magnetic Compound Fluid Polishing

摘要: 计算流体力学(Computational Fluid Dynamics)方法因其在流场分析中的广泛应用而被广泛认可。CFD仿真能够提供详细的流场信息，帮助提高材料去除效率和表面质量。因此，本章将利用COMSOL Multiphysics有限元软件对MCF抛光液的运动特性进行研究。材料去除率作为一个关键的性能指标，不仅影响加工效率，还与加工表面质量相关。在材料去除率的理论基础下，可知抛光压力和抛光速度是影响材料去除率的两个因素。传统抛光工件表面相对平整、不存在微结构特征，使得抛光过程中的速度和压力分布较为一致。然而，微织构工件由于其表面存在微结构阵列特征，使得抛光过程中的速度和压力分布呈现出非均匀性，对材料去除率会产生影响。因此，借助有限元仿真对速度和压力分布进行模拟，帮助理解微织构抛光过程中的材料去除行为。通过有限元仿真，分析在不同的抛光间隙和抛光转速条件下，MCF抛光液的流速分布和压力分布以及剪切速率情况。为后续进行MCF微织构抛光试验提供指导与参考依据。

Abstract: The method of Computational Fluid Dynamics (CFD) has garnered widespread recognition due to its extensive application in fluid flow analysis. CFD simulations are capable of providing detailed information about fluid flow fields, thereby facilitating enhanced material removal efficiency and surface quality. Consequently, this chapter employs the COMSOL Multiphysics finite element software to investigate the flow characteristics of MCF (Magnetic Compound Fluid) polishing fluid. Material removal rate, as a pivotal performance indicator, not only influences the efficiency of the machining process but is also intrinsically linked to the surface quality of the workpiece. Based on the theoretical foundation of material removal rate, it is evident that polishing pressure and polishing speed are two factors that affect the material removal rate. Traditional polished workpieces, characterized by relatively flat surfaces devoid of microstructural features, exhibit uniform distributions of speed and pressure during the polishing process. However, micro-textured workpieces, which possess microstructural arrays on their surfaces, display non-uniform distributions of speed and pressure during polishing, thereby affecting the material removal rate. Therefore, finite element simulations are utilized to model the distributions of speed and pressure, aiding in the comprehension of material removal behavior in micro-textured polishing processes. Through finite element simulations, the velocity and pressure distributions of the MCF polishing fluid, as well as the shear rate, are analyzed under different polishing gap and polishing speed conditions. This provides guidance and a reference for subsequent MCF micro-textured polishing experiments.

1. 引言

随着现代精密加工技术的不断发展，对复杂表面和微结构的高精度抛光需求日益增长[1]。磁性复合流体抛光作为一种新兴的柔性加工技术，因其结合了磁流变效应与流体的柔性特性，展现出在微纳结构表面处理中的独特优势[2]。该技术通过磁场调控磁性流体的流动行为，实现对复杂表面的均匀接触和高效材料去除，尤其适用于微织构表面的精密加工。然而，磁性复合流体在抛光过程中的流动特性、微织构表面的相互作用机制仍需深入研究。传统实验方法在揭示复杂流场分布规律和优化工艺参数方面存在局限性，因此，借助数值仿真工具对抛光过程中的流场进行分析成为一种有效的手段。

本文基于COMSOL Multiphysics平台，对微织构表面磁性复合流体抛光过程中的流场分布进行仿真分析[3]。通过建立包含微织构几何特征的三维模型，结合磁流体力学理论，模拟磁性流体在磁场作用下的流动行为，研究微织构几何参数、磁场强度及流体特性对流场分布的影响规律[4]。研究结果不仅为优化磁性复合流体抛光工艺参数提供理论依据，也为微织构表面的高效加工提供了新的技术思路[5]。通过数值仿真与实验验证相结合的方法，本研究旨在为磁性复合流体抛光技术在精密制造领域的应用奠定基础。

2. COMSOL流场仿真分析

2.1. 流场仿真数学模型

本节将采用COMSOL Multiphysics软件对MCF流动流场与磁场的耦合进行建模。MCF是一种新型的功能流体，它既具有流体物质的流动性又具有固体磁性材料的磁性。在无磁场时呈现液体状态，在磁场的作用下发生固态到液态转变，磁性微粒迅速链化，呈现为类固体状态，表现出非牛顿流体的力学特性。MCF抛光研究过程属于不可压缩连续流动，故而需满足流体流动的连续方程式和动量守恒方程式[6]。

流体流动的连续方程

$\frac{\partial \rho }{\partial t}+\frac{\partial \left(\rho u\right)}{\partial x}+\frac{\partial \left(\rho v\right)}{\partial y}+\frac{\partial \left(\rho \omega \right)}{\partial z}=0$ (1)

在上述方程中，ρ代表密度，t代表时间，u代表速度，而矢量u、v、ω分别表示速度矢量u在x、y和z方向上的分量。

动量守恒方程：

$\begin{array}{l}\frac{\partial \left(\rho U\right)}{\partial t}+\left(\frac{\partial }{\partial x}\left(\rho uU\right)+\frac{\partial }{\partial y}\left(\rho vU\right)+\frac{\partial }{\partial z}\left(\rho wU\right)\right)=\left(u\frac{\partial U}{\partial x}+v\frac{\partial U}{\partial y}+w\frac{\partial U}{\partial z}\right)-\frac{\partial p}{\partial x}+S_U\\ \frac{\partial \left(\rho V\right)}{\partial t}+\left(\frac{\partial }{\partial x}\left(\rho uV\right)+\frac{\partial }{\partial y}\left(\rho vV\right)+\frac{\partial }{\partial z}\left(\rho wV\right)\right)=\left(u\frac{\partial V}{\partial x}+v\frac{\partial V}{\partial y}+w\frac{\partial V}{\partial z}\right)-\frac{\partial p}{\partial y}+S_V\\ \frac{\partial \left(\rho W\right)}{\partial t}+\left(\frac{\partial }{\partial x}\left(\rho uW\right)+\frac{\partial }{\partial y}\left(\rho vW\right)+\frac{\partial }{\partial z}\left(\rho wW\right)\right)=\left(u\frac{\partial W}{\partial x}+v\frac{\partial W}{\partial y}+w\frac{\partial W}{\partial z}\right)-\frac{\partial p}{\partial z}+S_W\end{array}$ (2)

在这些方程中：

ρ是流体的密度，U、V、W分别是速度向量$u=\left(U,V,W\right)$ 在x、y、z方向上的分量，p是流体的压力，SU、SV、SW分别是广义源项在x、y、z方向上的分量，包括体积力或其他影响动量的外力，u、v、w是速度向量u的分量。

湍动能标准k方程：

$\frac{\partial u_{i}k}{\partial x_i}-\frac{\partial }{\partial x_i}\left[\left(v+\frac{v_i}{{\sigma }_k}\right)\frac{\partial k}{\partial x_i}\right]=P_r-\epsilon$ (3)

紊动能耗散率ε方程：

$\frac{\partial u_i\epsilon }{\partial x_i}-\frac{\partial }{\partial x_i}\left[\left(v+v_t\right)\frac{\partial \epsilon }{\partial x_i}\right]=\frac{\epsilon }{k}\left(C_{\epsilon 1}{P}_r-C_{\epsilon 2}\epsilon \right)$ (4)

式(4)中，$x_i$ 表示笛卡尔数学坐标系，$u_i$ 表示沿i方向的速度分量，ν表示液体的运动粘度。$P_r$ 表示紊动能生成率，表达式为：

$P_r=v_t\left(\frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i}\right)\frac{\partial u_i}{\partial x_j}$ (5)

$v_t$ 是湍流的涡动粘度(也称为湍流粘性系数)，$u_i$ 和$u_j$ 分别是速度向量u在i和j方向上的分量，$x_i$ 和$x_j$ 是笛卡尔坐标系中的坐标。公式(5)表示了在湍流中，由于速度梯度的存在，流体的粘性导致能量从平均流动转移到湍流涡旋中，从而生成湍动能。公式中的$\left(\partial u_i/\partial x_j+\partial u_j/\partial x_i\right)$ 是速度梯度张量的对称部分，它描述了流体的剪切变形。整个表达式乘以涡动粘度$v_t$ 来计算单位体积内湍动能的生成率。

2.2. COMSOL仿真模型的建立

建立微织构几何参数模型，微织构几何参数如表1所示，微织构几何形状如图1。

$\beta =\frac{{\displaystyle \sum _{j=1}^n{\displaystyle \sum _{k=1}^n{S}_{织构jk}}}}{S_{工件}}=\frac{n^2\cdot \frac{\pi \cdot d^2}{4}}{m^2}=\frac{{\left(\frac{m}{l}\right)}^2\cdot \frac{\pi \cdot d^2}{4}}{m^2}=\frac{\pi }{4}\cdot {\left(\frac{d}{l}\right)}^2$ (6)

β为织构面积率，n为织构数量，m为工件尺寸，j，k表述第j行第k列，l为阵列间距，d为织构直径，面积率仅与直径和阵列间距有关。

本文的微织构工件因其存在周期性，故仿真只分析一个周期单元上的速度场分布和压力分布。如图2所示。

Table 1. Exture parameter table

表1. 织构参数表

织构参数	数值
圆形微织构直径	50 μm
织构深度	50 μm
阵列间距	100 μm
织构面积率	19.635%

Figure 1. Micro-textured array

图1. 微织构阵列

Figure 2. Establishment of texture models

图2. 织构模型的建立

完成模型的绘制后，接下来，我们将从材料库中选择合适的材料属性。而铝合金微织构工件的材料则为铝合金(aluminum)。我们将MCF卧式抛光头的转速设定在300至900转每分钟(rpm)之间。流体的属性被定义为具有10千帕(kPa)的屈服应力和1360千克每立方米(kg/m³)的密度。

并对这些表面的单元进行更细密的划分，以提高有限元仿真的精度，如图3。

Figure 3. Mesh generation

图3. 网格划分

在边界条件方面，我们假设流体是不可压缩的，并在微织构表面上应用无滑移边界条件。卧式抛光头的边界条件被设置为逆时针旋转，其速度与无滑移边界条件相同。此外，我们假设加工对象没有漏磁或渗透等缺陷。在本研究中，我们假设MCF的流动是层流，并选择CFD模块下的“旋转机械，层流”接口来计算流体域和边界条件。

2.3. 微织构抛光速度场分析

为了研究不同工艺参数对MCF流场的影响，需要建立多个模型，将卧式抛光头的转速设置为300 rpm~900 rpm，每间隔为200 rpm，抛光间隙设置为1 mm，的四个模型进行分析。

Figure 4. The influence of polishing speed on velocity field distribution

图4. 不同抛光转速对速度场的影响

从图4中可以看出，当抛光间隙为1 mm时，不同抛光转速下的速度场分布表现出以下特点和变化趋势：抛光转速为300 rpm时：速度场分布较为均匀，流体在微织构区域的流动速度相对较低。速度的最大值出现在微织构的边缘，表明流体在这些区域受到较大的剪切力。流体的速度范围为0到21 × 10⁻³ m/s，标尺显示速度分布较为集中，整体流动较为平稳，抛光转速增加到500 rpm时：速度场分布开始出现局部不均匀性，尤其是在微织构的边缘区域，速度梯度显著增大。速度的最大值仍出现在微织构的边缘，但整体速度范围扩展至0到35 × 10⁻³ m/s，表明流体的流动速度加快，且速度分布的离散性增强。

抛光转速为700 rpm时：速度场分布的不均匀性进一步加剧，流体在微织构区域的流动速度显著增加。速度的最大值仍集中在微织构的边缘，且速度梯度更加明显，表明流体在这些区域的剪切力进一步增大。流体的速度范围仍为0到35 × 10⁻³ m/s，但速度变化更加显著，分布的离散性进一步增强。抛光转速增加到900 rpm时：速度场分布变得非常不均匀，流体在微织构区域的流动速度显著加快，且速度梯度变化较大。速度的最大值仍出现在微织构的边缘，但整体速度范围进一步扩展至0到42 × 10⁻³ m/s，表明流体的流动速度进一步加快。标尺显示速度分布的离散性显著增强，速度变化更加剧烈。如图5所示为抛光转速对速度场大小的影响。

Figure 5. The influence of polishing speed on velocity magnitude

图5. 抛光转速对速度大小的影响

2.4. 微织构抛光压力场分析

从图6中可以看出，当抛光间隙为1 mm时，不同抛光转速下的压力场分布表现出以下特点和变化趋势：

抛光转速为300 rpm时：压力场分布较为均匀，压力值主要集中在微织构区域的边缘。压力范围为1 × 10³ Pa到6 × 10³ Pa，压力分布较为集中，表明流体在低转速下受到的剪切力较小，压力变化较为平稳。抛光转速增加到500 rpm时：压力场分布开始出现局部不均匀性，尤其是在微织构的边缘区域，压力值显著增大。压力范围扩展到1 × 10³ Pa到8 × 10³ Pa，压力分布的离散性增强，表明流体在这些区域受到更大的剪切力，压力梯度开始显现。

抛光转速为700 rpm时：压力场分布的不均匀性进一步加剧，压力值在微织构边缘区域显著集中，且压力范围扩展到1 × 10³ Pa到10 × 10³ Pa。压力分布的离散性进一步增强，表明流体在高转速下受到的剪切力更大，压力梯度更加明显。抛光转速增加到900 rpm时：压力场分布变得非常不均匀，压力值在微织构边缘区域达到最大值(约11 × 10³ Pa)，且压力梯度变化显著。压力范围扩展到1 × 10³ Pa到11 × 10³ Pa，压力分布的离散性显著增强，表明流体在高转速下受到的剪切力进一步增大，压力变化更加剧烈。

Figure 6. Pressure field distribution at different polishing speeds

图6. 不同抛光转速下的压力场分布

Figure 7. The influence of polishing speed on pressure

图7. 抛光转速对压力的影响

如图7为不同转速下的压力场分布折线图。从压力增长趋势来看，随着转速的增大，速度逐渐增大，在转速为900 rpm时，压力达到了最大。

随着抛光转速的增加，压力呈线性增长趋势。压力从低转速100 rpm时的约3 kPa逐渐增加到高转速900 rpm时的约15 kPa。在低转速范围内，压力随转速的增加变化较为显著，曲线斜率较大，表明压力对转速的变化较为敏感。在高转速范围内，压力随转速的增加仍然呈上升趋势，但曲线斜率略有减缓，表明压力的增长速率趋于稳定。

3. 微织构抛光材料去除模型

材料去除数学模型的建立

MCF微织构抛光技术的材料去除原理在本质上和传统抛光轮是一致的，都是机械和化学综合作用的结果，为了研究影响深孔工件表面材料去除率的因素，作出如下假设：

(1) 假设MCF深孔抛光的材料去除主要是机械作用的结果，依靠磨料粒子在接触表面形成一定的压力与相对速度，从而达到对待加工表面的切削作用；

(2) MCF中的磁性颗粒形状为圆形；

(3) 卧式抛光头与待加工工件表面之间的“柔性研磨层”中的链状结构均匀且稳定。

由Preston方程[7]：

$\text{MRR}=\text{KPV}$ (7)

式(7)中，K为Preston系数，由试验测试得出，P为抛光压力，V为抛光速度。磁流变抛光中，抛光压力P是一个较为复杂的参数，抛光过程中主要由流体流动产生的动压力、磁流变液在磁场中产生的磁化压力、磁辅助抛光液的重力造成的压力组成[8]。

$P=P_m+P_d+P_g$ (8)

式(8)中，$P_d$ 为流体动压力，$P_m$ 为磁化压力，$P_g$ 为重力作用下产生的压力。相较于磁化压力$P_m$ ，和流体动压$P_d$ 抛光液重力作用下产生的压力$P_g$ 较小，可以忽略。对于磁辅助抛光液在磁场下产生的磁化压力，其表达式为：

$P_m={\mu }_0{\displaystyle {\int }_0^H{M}_f}\text{d}H$ (9)

式(9)中，${\mu }_0$ 为真空中磁导率，$M_f$ 为磁辅助抛光液的磁化强度，H为外加磁场强度。其中，磁辅助抛光液的磁化强度$M_f$ 为[9]：

$M_f=3\varphi {\mu }_f\frac{{\mu }_p-{\mu }_f}{{\mu }_p+2{\mu }_f}$ (10)

将式(10)代入式(9)中可得到磁化压力，$\varphi$ 为磁辅助抛光液中磁性颗粒的体积比浓度，μ_f为基载液的磁导率，${\mu }_p$ 为磁性颗粒的磁导率，得到$P_m$ ：

$P_m=3\varphi {\mu }_0{\mu }_f\frac{{\mu }_p-{\mu }_f}{{\mu }_p+2{\mu }_f}{\displaystyle {\int }_0^{H}H}\text{d}H$ (11)

最终得到材料去除率MRR为式(12)

$\text{MRR}=K\left(P_d+3\varphi {\mu }_0{\mu }_f\frac{{\mu }_p-{\mu }_f}{{\mu }_p+2{\mu }_f}{\displaystyle {\int }_0^{H}H}\text{d}H\right)V$ (12)

首先计算普雷斯顿系数k，根据试验，由用尺寸为20 × 20 × 3 mm的铝合金工件为实验元件进行磁流变抛光。具体条件为：$h_0=1\text{ mm}$ ，$V=0.314\text{ m/s}$ 、${\eta }_0=14.546\text{ Pa}\cdot \text{s}$ 、${\mu }_0$ 是真空磁导率为$4\pi \times {10}^{-7}\text{H/m}$ ，${\mu }_p=1500{\mu }_0$ ，${\mu }_p$ 是磁性颗粒磁导率，${\mu }_f$ 是基载液的磁导率约等于${\mu }_0$ ，的磁性微粒在磁流变抛光液中所占的体积比是$\varphi =0.33$ ，H为磁场强度，用高斯计测得$H=2.54\times {10}^7\text{A/m}$ ，由公式(11)

$\begin{array}{c}{P}_m=3\varphi {\mu }_0{\mu }_f\frac{{\mu }_p-{\mu }_f}{{\mu }_p+2{\mu }_f}{\displaystyle {\int }_0^{H}H}\text{d}H\\ =3\times 0.33\times {\left(4\pi \times {10}^{-7}\right)}^2\times \frac{1499}{1502}\times \frac{1}{2}\times {\left(2.54\times {10}^7\right)}^2\\ =5.032\times {10}^2\end{array}$ (13)

磁化压力为5.032 × 10² Pa，由于磁化压力仅与磁场强度和抛光液的磁导率有关，且在后续实验中均采用同种抛光液和磁铁，故磁化压力在后续实验中为定值。

流体动压$P_d$ 和速度V均通过仿真得到，带入普雷斯顿方程，通过试验测得铝合金质量前后变化为0.06 mg/min。求解得到k值为4.87 × 10⁻⁴。

即可得到抛光微织构过程中的去除函数[10]，见式(14)

$\text{MRR}=4.87\times {10}^{-4}\left(P_d+5.032\times {10}^2\right)V$ (14)

将抛光转速和仿真做出的流体动压力带入上式(14)得到了微织构抛光的理论去除函数

Figure 8. Material removal rate verification

图8. 材料去除率验证

从图8可以看出，理论材料去除率和实际材料去除率随抛光转速的变化趋势基本一致，表明理论模型能够较好地反映实际抛光过程中的材料去除规律。理论曲线和实际曲线均呈现出随着抛光转速增加，材料去除率逐渐增大的趋势。这表明抛光转速对材料去除率具有显著影响，且理论模型能够有效捕捉这一变化规律。

在低转速范围内，理论曲线与实际曲线的吻合度较高，表明理论模型在低转速下的预测精度较高。此时，材料去除率的变化较为线性，理论值与实际值的偏差较小。在高转速范围内，实际材料去除率略高于理论值，但两者的变化趋势仍保持一致。这可能是由于实际抛光过程中存在一些未被理论模型完全考虑的因素，如流体的湍流效应或微织构的几何复杂性。

4. 结论

本研究通过计算流体力学(CFD)方法，借助COMSOL Multiphysics有限元软件对磁性复合流体(MCF)抛光液在微织构表面抛光过程中的流场分布进行了深入的仿真分析，主要结论如下：

随着抛光转速的增加，速度场和压力场分布逐渐变得不均匀。在低转速下，速度场和压力场分布较为均匀，而在高转速下，速度和压力在微织构边缘区域显著增加，形成明显的梯度变化。

速度场的最大值始终出现在微织构的边缘区域，且随着转速的增加，速度范围从低转速时的0到21 × 10⁻³ m/s扩展到高转速时的0到42 × 10⁻³ m/s。

压力场的最大值也集中在微织构边缘区域，且随着转速的增加，压力范围从低转速时的1 × 10³ Pa到6 × 10³ Pa扩展到高转速时的1 × 10³ Pa到11 × 10³ Pa。

基于Preston方程建立的材料去除模型表明，材料去除率随抛光转速的增加而显著增加，且理论模型与实验结果在低转速范围内吻合度较高。抛光压力和抛光速度是影响材料去除率的关键因素，其中抛光转速对材料去除率的影响尤为显著。

本研究通过数值仿真与实验验证相结合的方法，仿真方法得到的理论材料去除率曲线与实际曲线的吻合度较高，尤其是在低转速区域，验证了仿真方法的可靠性。这一结果表明，通过COMSOL流场仿真和材料去除率建模，可以为后续工艺试验提供可靠的理论指导。揭示了磁性复合流体抛光过程中速度场和压力场的分布规律，以及抛光转速对材料去除率的影响机制，为磁性复合流体抛光技术在精密制造领域的应用奠定了坚实的基础。