摘要：为改善液压滑阀开启过程中阀芯的工作性能，基于流体- 固体耦合理论和动网格技术，建立滑阀开启过程三维流体- 固体耦合数值分析模型，对阀芯在驱动力、弹簧力及液动力作用下的动作过程进行了模拟分析。针对阀芯受力分析结果，改进滑阀内部流道结构；以改进后阀芯槽口和凸台的4个结构参数为优化对象，利用最小二乘拟合和反向传播神经网络构建最大冯米塞斯应力与液动力峰值的目标函数，借助遗传算法确定了槽口和凸台的4个结构参数最优值。研究结果表明：优化后阀芯所受液动力和最大冯米塞斯应力的峰值分别减小了16.3%和22.0%；优化设计阀芯的结构参数可明显提高滑阀开启性能。
关键词：流体力学； 液压滑阀； 流体- 固体耦合； 动网格； 遗传算法
0 引言
滑阀结构被液压系统的流量、方向和压力控制阀广泛采用，其开启过程中流场不稳定，会导致阀芯控制精度不高。阀芯开启过程的动力学特征直接影响液压滑阀的控制性能，对于液压系统有着较大的影响。考虑流体- 固体耦合的研究可追溯到19世纪初的机翼气动弹性问题[1]。目前，流体- 固体耦合的研究主要集中于水轮机、冷却器盖[2-4]等，针对液压阀的研究并不多。而对于液压阀阀芯的动态仿真大多是直接给定阀芯的开口度或阀芯的运动速度，并未考虑阀芯开启过程的瞬态受力情况和流体- 固体耦合的影响[5-6]。因此，有必要对液压滑阀的阀芯及其开启过程进行重点研究，分析其瞬态受力情况并加以优化。针对上述问题，借助流体- 固体耦合理论和动网格技术，以流体力学软件CFD-ACE+为平台，对液压滑阀开启过程进行流体- 固体耦合特性分析；针对分析结果，改进液压滑阀的内部流道结构并进行优化；建立优化后的模型，对优化结果进行仿真验证，并比较优化前后液压滑阀的开启性能。
1 几何模型和网格划分
1.1 几何模型
液压滑阀通常由多个阀腔组成，基于相似性原理，用其中的一个阀腔就能说明滑阀内部流体的流动状态。借鉴文献[7-9]的简化方法，归纳出具有共性的液压滑阀作为研究对象，其结构如图1所示。
图1 液压滑阀结构
Fig.1 Structure diagram of hydraulic spool valve
在该滑阀阀芯的开启过程中，液压油的流动具有对称性，为提高计算效率，取流动区域和阀芯的一半作为研究对象，建立用于流体- 固体耦合分析的三维几何模型，如图2所示。
图2 液压滑阀几何模型
Fig.2 Geometric model of hydraulic spool valve
1.2 网格划分
流体- 固体耦合会导致网格变形较大，出现三维瞬态流场解析发散的问题，单一的网格变形方法很难同时具有较大的网格变形能力和得到较高的变形质量。为此，本文将弹簧法和超限插值法相结合，解决结构网格变形问题。基于网格分块的思想，将计算区域分成16个网格块；采用弹簧连接不同网格块的角点，通过弹簧驱动各角点的变形；对于网格块内部节点，采用超限插值法将角点位移通过插值函数分布到网格块内部节点[10]；网格划分采用前处理软件GEOM实现；为更好地捕捉阀口处的流场变化，对阀口区域网格进行细化。在大量试算的基础上，最终确定网格总数为95 252个。其中，流域网格51 029个，阀芯网格44 223个，如图3所示。
图3 网格划分
Fig.3 Grid division
2 控制方程
2.1 动网格条件下的流场控制方程
在网格运动过程中，流场内控制体将发生改变，一般采用积分形式来表示诸守恒方程：
(1)
式中：ρ为液体密度；uf为流体速度；V为控制体体积；p为控制体压力；φ为通用变量，代表uf、V、p等求解变量；∂V为控制体体积V的边界；A为控制体截面积向量；ug为动网格变形速度；Γ为广义扩散系数，代表μ+μt、μ+μt/σk、μ+μt/σε，μ为动力黏度，μt为湍动黏度，k为湍动能，ε为湍流耗散率，σk是与k对应的Prandtl数，σε是与ε对应的Prandtl数；
φ为通用变量梯度；Sφ为φ的源项。
2.2 流体- 固体耦合控制方程
在流体- 固体耦合交界面处应力：
(2)
式中：ΔFi为在i方向上的力(N)；ΔAi为在i方向的受力面积(m2)。
力平衡方程
(3)
式中：σxx、σyy、σzz为应力分量；τxy、τxz、τyx、τyz、τzx、τzy为应变分量。
几何方程
(4)
式中：ξxx、ξyy、ξzz为应变分量；u、v、w为x、y、z方向的速度分量；γxy、γyx、γzy、γyz、γxz、γzx为切应变分量。
本构方程
(5)
式中：λ为泊松比；E为杨氏模量(GPa)。
3 边界条件与参数选取
3.1 流域参数
入口压力取6.00 MPa，出口压力取0.15 MPa；液压油型号VG46，密度ρ=885 kg/m3，运动黏度ν=46 mm2/s，由于液压油压缩性很小，视为不可压缩流体；弹簧刚度80 N/mm；滑阀阀芯两端驱动压差1 MPa；采用k-ε湍流模型。
3.2 阀芯参数
阀芯材料为45号钢。密度ρs=7 850 kg/m3；杨氏模量Es=210 GPa；泊松比λs=0.3.
4 仿真结果分析与优化
4.1 阀芯运动情况分析
阀芯的运动状态由驱动力、弹簧力及液动力共同决定，仿真得到的阀芯运动情况如图4所示。图中点1为阀芯开启初期(t=0.42 ms)，点2为阀芯位移最大时(t=1.80 ms)，点3为阀芯位移达到稳态时(t=47.00 ms)。
图4 阀芯开启过程位移变化
Fig.4 Displacement change of spool during opening
从图4可知，液压滑阀开启过程中，阀芯位移存在一定的波动，最大位移达到1.49 mm，稳态位移为0.78 mm. 阀芯位移的波动会引起流量的波动，从而影响执行元件的运动状态；位移的波动也会产生振动和噪声，影响流场及阀芯瞬态受力情况。
液压滑阀阀芯稳态位移和稳态液动力的理论计算公式[11]为
(6)
式中：D为流道直径，D=16 mm；xv为阀芯稳态位移(mm)；W为阀口的面积梯度(mm2/mm)；kt为弹簧刚度，kt=80 N/mm；F为滑阀阀芯所受的驱动力(N)；Fs为阀芯所受的稳态液动力(N)；Δp为滑阀进出口压差，Δp=5.85 MPa；Δpq为驱动压差，Δpq=1.00 MPa.
联立(6)式，可以解得阀芯稳态位移xv=0.71 mm，稳态液动力Fs=-43.73 N.
仿真得到的阀芯稳态位移xvf=0.78 mm，稳态液动力Fsf=-47.96 N. 比较理论计算与仿真计算结果，阀芯稳态位移和稳态液动力的误差分别为9.9%和9.7%，可以认为仿真计算结果与理论计算结果基本吻合，证明了仿真的可靠性。误差的产生主要是因为在理论计算中稳态液动力采用经验公式求解，没有考虑内部结构不同及流场微观流动带来的影响。
4.2 流体- 固体耦合分析
图5 不同时刻流体- 固体耦合分析
Fig.5 Fluid-structure interaction analysis at different moments
针对阀芯开启初期(t=0.42 ms，图4中点1)，最大位移(t=1.8 ms，图4中点2)，稳态位移(t=47 ms，图4中点3)，三维流体- 固体耦合仿真结果如图5所示，图中给出了液压滑阀流域的压力场、压力等值线分布以及液压滑阀阀芯所受的冯米塞斯应力分布。其中，冯米塞斯应力是根据第四强度理论得到的一种当量应力，考虑了第1、第2和第3主应力。冯米塞斯应力值越大的地方材料变形就越大，往往也是越容易受到破坏的地方，因而可以用来对疲劳、破坏等进行评价[12]。
从图5可以看出：
1)在流场的A区域和B区域会产生负压，负压区域面积随着阀芯位移的增大而变大。流场负压区是引起气蚀的主要原因，气蚀会导致阀芯产生气蚀破坏；
2)在液压滑阀开启初期，阀芯应力集中于阀芯的1、2、3和4区域且四周对称分布。当阀口位移较大后，应力集中于3和4区域。
阀芯在开启过程中所受最大冯米塞斯应力的大小和位置随时间变化情况如图6所示。由图6可看出，最大冯米塞斯应力主要集中在1、2、3、4区域，在阀芯开启瞬间达到峰值。
图6 不同时刻最大冯米塞斯应力分析
Fig.6 Analysis of maximum von Mises stresses at different moments
4.3 液压滑阀结构优化分析
为改善滑阀开启过程中阀芯受力情况，需减小阀芯所受的液动力(包括稳态液动力与瞬态液动力)和因流体- 固体耦合作用所产生的最大冯米塞斯应力。现改进滑阀内部流道结构，在阀芯上增设凸台和槽口，新结构如图7所示。
图7 新结构滑阀
Fig.7 New structure of spool valve
图8分别给出了原结构、仅有凸台、仅有槽口和新结构阀芯在开启过程中最大冯米塞斯应力和液动力的变化情况。其中，液动力是通过求解阀芯各壁面轴向力的合力得到[13]，为瞬态液动力和稳态液动力之和，本文将其作为其中一个优化指标。从图8中可看出，阀芯凸台结构能有效地减小液动力，槽口则对最大冯米塞斯应力减小较为明显，而同时增设槽口和凸台后，最大冯米塞斯应力和液动力均减小。
图8 结构变化对阀芯受力影响
Fig.8 Effect of structure change on force of spool
经大量试算，选取如图7所示的L1、L2、L3、L4的尺寸为优化对象，以滑阀开启过程中不同瞬时阀芯所受最大冯米塞斯应力和液动力的峰值为目标函数，建立33个滑阀试验模型，进行三维流体- 固体耦合仿真分析，结果如表1所示。为了便于表达，实验因子L1、L2、L3、L4分别用x1、x2、x3、x4表示，评价指标σmax、Fsmax分别用y1、y2表示。其中，1.0 mm≤x1≤5.0 mm，0.1 mm≤x2≤2.0 mm，0.1 mm≤x3≤2.0 mm，0.1 mm≤x4≤2.0 mm.
表1 液压滑阀开启过程实验设计表
Tab.1 Test design of opening process of spool valve
优化阀芯结构要求最大冯米塞斯应力和液动力均最小，取(7)式作为优化目标函数：
y=wvfv+wtft，
(7)
式中：fv、ft分别为最大冯米塞斯应力和液动力峰值的隶属度(无因次量)；wv、wt为对应权重，均取0.5. 隶属度函数[14]按(8)式计算：
(8)
式中：si、so和sw分别为表1中某一样本的指标值、该指标的最优值和最差值，i=1，2，…，33.
为取得每个指标的最优值和最差值，采用2阶函数分别对最大冯米塞斯应力和液动力峰值的仿真结果进行最小二乘拟合，拟合值和期望值如图9所示。拟合方程分别为
图9 最大冯米塞斯应力和液动力峰值最小二乘拟合情况
Fig.9 Least square fitting of peak value of maximum von Mises stress and transient flow force
拟合方程y1和y2评价标准为误差决定系数[15]，检验误差决定系数为
(9)
式中：N为检验点的数量；y*(i)为拟合值；y(i)为仿真值；
为仿真的平均值。经计算，拟合R2分别为0.948和0.964，这表明响应模型与仿真模型基本一致，可以作为最大冯米塞斯应力和液动力峰值关于各设计因子的数学模型。
利用遗传算法[16]对拟合出的两个方程求解最值。在处理约束时，采用罚函数的方法，将其转化为无约束优化问题。罚函数具体操作是：在约束条件内，适应度值高；约束条件外，适应度值低，容易被淘汰。得出最大冯米塞斯应力峰值的最优值和最差值分别为68 MPa、130 MPa，液动力峰值的最优值和最差值分别为-68 N、-101 N. 根据(7)式计算得到目标函数(见表1)。
由于最小二乘拟合各设计因子与目标函数之间的关系时效果不太理想，采用反向传播(BP)神经网络构建其响应面模型。选取表1中27组数据作为训练数据，其余6组与预测结果进行对比分析误差，拟合情况如图10所示。
图10 目标函数BP神经网络模型预测结果
Fig.10 Prediction of BP neural network for objective function
从图10可看出，预测误差绝对值在8%以内，因此构建的BP神经网络预报精度完全可以接受。再次借助遗传算法寻求目标函数的最优值，寻优适应度收敛过程曲线如图11所示。最终得到优化结果如表2所示。
图11 遗传算法适应度收敛曲线
Fig.11 Fitness convergence curves of GA
表2 阀芯优化结果
Tab.2 Optimimized results of spool
根据寻优结果，选取L1、L2、L3、L4的尺寸分别为3.1 mm、0.9 mm、1.6 mm、0.9 mm，建立优化后的三维流体- 固体耦合模型，得出相同设定条件下阀芯所受的最大冯米塞斯应力和液动力变化情况，与优化前的结果对比，如图12所示。原结构、优化设计及仿真结果对比如表3所示。
图12 滑阀开启过程优化前后对比情况
Fig.12 Comparison of results during opening before and after optimization
表3 阀芯优化设计验证对比
Tab.3 Verification and comparison of optimization design for spool
由图12和表3可知，优化后的阀芯在开启过程中不同瞬时受到最大冯米塞斯应力的峰值为94.7 MPa，减小了22.0%，液动力峰值为83.2 N，减小了16.3%，与优化设计结果较为接近，但是均比优化前的值要小，且在整个开启过程中最大冯米塞斯应力变化较为平稳，明显改善了滑阀开启性能。
5 结论
1)液压滑阀开启过程中流域会产生负压区，负压区主要集中在阀口处靠近阀芯和阀座的两个区域，且随阀芯位移的增加而扩散；
2)液压滑阀开启过程中不同位置阀芯所受应力是不断变化的，且在阀芯开启初期最大冯米塞斯应力达到峰值；
3)液压滑阀开启过程中，阀芯所受最大冯米塞斯应力值主要集中于4个截面突变区域。改进阀芯结构，可减小这些区域的最大冯米塞斯应力，同时也能减小液动力；
4)内部流道结构的优化设计值与仿真结果基本吻合，证明仿真分析方法和优化方法具有可行性。
