基于SPH的管道内液压油可视化仿真效果分析

　　0 引言
　　 液压油普遍应用于生产生活的各种机械中，但是由于油体被管壁遮挡，不能实时看到管内的液压油流动状况。然而，液压教学、液压安全检测、液压系统生产运行以及相关液压系统培训等都需要液压油管内流动的视觉再现，因此，对管内液压油可视化模拟具有现实意义。
　　 液压油是多种添加剂与矿物油混合的一种液体，它粘性很强、密度大、高可压缩性，除此之外，液压油的粘性还随着温度的变化而变化，通常称为液压油的粘温性。外界温度越高，液压油的体积就会越大，其在管道内所受的压力也就越大。液压油是一种流体，所以流体的模拟方法也可应用在液压油的模拟中。目前，液压油的模拟主要有基于参数建模的方法和基于物理的方法。参数建模方法一般运用波动理论函数来模拟小范围内流体波动的效果，不能真实地描述现实世界中液压油在管道内流动的情况，而且控制比较困难。时培成^[1]等应用数字建模仿真模拟推导出的液压油密度及油液中空气含量随液压油液压力与温度变化的函数关系。文献[2]通过数学建模，对不同流速比和相同流速比但雷诺数不同的汇流管道进行了模拟。基于物理的模拟方法是基于立体运动的物理学规律对油体运动形态进行模拟，从而使液压油流动形态的模拟更具真实性;基于物理的方法有欧拉法和拉格朗日法两种。欧拉法计算量较小、方法简单，但是因为它以流体运动的整个区域作为观察对象，所以必须对整个空间进行处理，控制起来比较困难。相对而言，拉格朗日法实施比较容易，表达较清晰，不需要对流体的整个区域进行处理，在保证物理质量守恒的条件下，容易控制。虽然拉格朗日法重建平滑运动界面比较难，但其计算量也会随着系统粒子数量的增加而快速增加，随着计算机硬件的发展，这一缺点的影响也在下降，人们仍然普遍采用拉格朗日法模拟流体。    SPH（Smoothed Particle Hydrodynamics）是目前应用最多的拉格朗日方法之一。Losasso^[3]等结合SPH和Particle level set 的方法模拟了流体。随着硬件的发展，采用GPU来加速模拟的整个计算成为又一有效手段。文献[4]基于GPU采用SPH方法较成功地模拟了管道内液体流动和液体出管道口后的漩涡、溅射等现象，证明了基于GPU的计算速度的高效性。但是，它只对管道内粒子赋予基本属性，构建了简单的模型，且不能提供真实的视觉效果。
　　 在液压系统中，经常遇到多个支流混合为一个管道的现象，管道交叉处漩涡现象是支流混合的普遍现象，使用以上方法并不能达到粘温特性的视觉效果、较低计算量以及较容易控制的目的。本文提出的适用于液压油模拟的Stoke非定常平行剪切流模拟方法可以达到容易控制的目的。结合物理模拟方法的优点，引入人工粘性，设置粘温系数，对其参数进行定量设置，降低了计算量，提高了粘温特性的视觉效果。
　　1 Stoke非定常平行剪切流数学原理
　　 在拉格朗日粒子法中，SPH模型是对Navier-Stokes方程进行求解的^[5]。运用SPH模型求解，其核心原理可以理解为在空间中离散位置的流体粒子所携带的状态量，通过对最近相邻粒子的状态量值进行累加求和得到，这些粒子的状态量包括粒子本身的质量、粒子在管道内所受的压强力、粒子相互间的粘性力、外部作用力、位置、颜色、透明度等，其粒子运动方程为：
　　dνdt=-1ρP+μ2vρ+Fρ（1）
　　 其中，-1ρP为流体指点所受的压应力，μ2vρ为流体指点所受的切应力，μ为流体粘性系数， Fρ表示流体指点所受到的外场力质量。
　　1.1 管道液压油粒子模型建立
　　管道内液压油可以看作由很多不同形状、大小、位置、颜色、透明度、运动速度、运动方向和生命周期等属性的粒子组成。这些粒子是随着时间动态变化的，在这些变化过程中，有些粒子会消亡，有些粒子属性会改变。液压油粒子状态量值的积分形式表示为：F（x）=∫Ωf（x′）W（x-x′，h）dx′ （2）
　　 式（2）以离散求和的形式表示为：
　　Fi（x）=∑j≠iFjmjρjW（xi-xj，h）（3）
　　 式中：j：遍历所有流体粒子，mj：粒子j的质量，ρj：粒子j的密度，Fi为位于xi的场量，Ω为包含x的积分体积，h为光滑半径，W（x-x′，h）为具有一定作用半径的光滑核函数。针对不同的受力状态量，所选择的光滑函数是不同的。 密度核函数为：
　　W（r，h）poly=kpoly（h2-r2）3，0≤r≤h
　　0，otherwise
　　kpoly为固定系数，在三维模拟中 kpoly=31564πh9。
　　粒子压强力的计算通常取计算粒子压强的平均值，单个粒子压力可用理想气体状态方程计算，其光滑核函数为
　　W（r，h）spiky=kspiky（h-r）3，0≤r≤h
　　0，otherwise
　　其中kspiky=15πh6。
　　考虑到液压油的强粘性，本文借鉴XSPH^[6]方法，在求解通用SPH模型时引入了人工粘性，降低了方程中使用的动力学系数，有助于在较大步长下保证流体方程的稳定性和收敛性 [7]。为了凸显液压油的粘温特性，在粒子受力中增加了粘温影响因子μ，它是受温度影响形成的粘温指数。μ=μ₀·TT₀n
　　其中，T₀，μ₀分别为参考温度和相应的粘度，n为随温度不同而设定的不同指数。一般而言，
　　n=12，T>3000
　　1，300 　　8ρ，90≤T≤300
　　本文结合XSPH得到液压油粒子粘性力方程：
　　Fxsphi=∑j≠i2mjvj-viρi+ρjWviscosity
　　其中
　　W（r，h）viscosity=kviscosity（-r32h3+r2h2+h2r-1）
　　，0≤r≤h
　　0，otherwise ，kviscosity=152πh3。
　　而粘性力与粘温指数的关系为Fxsph=μuy，根据温度的设定即可求出粘温指数进而求出粘性力。这里将粘性力与粘温指数的求解与液压油的温度关联起来，可以有效地表现出液压油的粘温特性。
　　1.2 管道粒子边界处理
　　运用虚粒子来处理粒子边界条件和解决粒子缺陷问题是有效方案之一。在Monaghan^[8]的整个流体粒子模拟研究中，在边界粒子i影响领域h内，所有相邻粒子被分为两大类：实粒子和虚粒子。实粒子是所有与粒子i相邻的内部粒子;虚粒子是所有与粒子i相邻的边界 粒子。该方法的主要思想是：在固定边界上分布一组虚粒子，用于对邻近边界的粒子强排斥，从而阻止边界粒子的非物理性穿透。
　　本文引入Monaghan虚粒子思想，将管道壁看成由一层固定的虚粒子组成的粒子壁，内部实粒子一旦进入虚粒子的作用域就执行碰撞检测，根据实粒子与虚粒子之间的作用距离赋值作用力大小，距离越近作用力越大。实粒子与虚粒子碰撞后的反弹速度大小及方向与赋予的虚粒子刚性作用力和虚粒子与实粒子的距离有关。虚粒子和实粒子在管道内流动情况如图1所示，空圆点代表虚粒子，实圆点代表实粒子。
　　设定管道边界液压油粒子虚粒子与实粒子的作用力主要是碰撞力，此处的粘性力忽略不计。
　　1.3 管道汇流漩涡模拟数学原理
　　管道内液压油汇流是单相流体在管道内的混合流动，即液压油经由管道射入充满同种油液的流动液体中，这种汇流不可避免地会产生漩涡现象，即湍流现象。湍流产生的一般规律如下：在粘性底层中，在平面上有顺流向的高速区和低速区相间的带状流动结构，低速带向下游流动中其头部会缓缓上抬，与壁面形成横向漩涡，而横向漩涡在流场作用下会上升，且在向下游的流动中会变形成为马蹄形涡，马蹄涡头部的上举最终形成底部低速流体向上层高速流体区域喷射，然后出现外层高速流体向下游流体俯冲清扫。流向速度分布曲线上出现拐点，增加了流动的不稳定性，促成层流向湍流的转变。漩涡的仿真模拟是基于粒子进行的，本文在使用SPH方法的同时在流体粒子的作用力中引入漩涡力，如图2所示。漩涡力使流体粒子在两条支路交汇处做圆环运动，形成漩涡视觉效果。
　　图1 管道内粒子分布
　　图2 漩涡力
　　本文主要借用κ-ε模式来模拟双管道交叉处产生的漩涡现象。κ-ε模式^[8]结合涡粘模式即确定合理湍流脉动的特征长度l来求解。湍流脉动长度尺寸可以由湍动能κ和湍动能耗散率ε来估计。模型方程为：κκt+κxj=2vT2+xj[（v+vTσk）κxj]-ε （4）
　　 ε模型方程为：
　　εt+εxj=Cε1εκ（2vT2）+
　　xj[（v+vTσε）εxj]-Cε2ε2κ（5）
　　 νT涡粘系数方程为
　　νT=Cμκ2ε（6）
　　 涡粘假设的平均运动方程为：
　　<μi>t+<μj><μi>xj
　　=-1ρ<Π>xi+xj[（v+vT）<μi>xj]+（7）
　　 其中，
　　<Π>=<ρ>+2k3（8）
　　 平均运动的连续方程为：
　　<μi>xi=0（9）
　　 平均运动的应变率张量为：
　　=12（<μi>xj+<μj>xi）（10）
　　 uj为湍流流动的平均速度，Cε1、Cε2、δε为经验常数。由式（4）～（10）可以构成湍流封闭模式方程组来求解。
　　Stoke非定常平行剪切流模拟方法的数学模型是：运用粒子的动态属性，设置粒子的初始值，并根据数学模型规定粒子的运动轨迹，使粒子沿轨迹进行运动来模拟液压油的流动特性。在粒子的生命周期中一般有以下几个步骤：①新粒子产生，并被赋予属性;②根据粒子的动态属性对粒子进行移动和变换并赋予粒子新属性;③判断粒子生命值;④若粒子生命值已为0，则删除该死亡粒子;⑤若粒子生命值不为0，则绘制该粒子并在屏幕上显示出来。
　2 实验结果与分析
　　在本实验中运用流体模拟软件建立Stoke非定常剪切流液压油模拟数学模型，进行不同角度双管交叉处漩涡的模拟，具体实验步骤如下：
　　（1）建立三维管道物理模型。首先运用Gamebit软件进行三维物理模型的创建，在创建时划分模型混合型网格，设定模型进出口位置，针对交叉管道模型设定中间和左侧的端口为液压油入口，右侧为液压油出口。
　　（2）将模型导入仿真软件中，进行收敛性检测。以90度交叉管道为例，将Gambit软件导出的文件导入fluent中，在计算出网格数之后进行收敛性检测。图3所示为90度交叉管道三维模型的收敛性检测图，从图中可以看出该模型在迭代260次时就能很好地收敛。这样的收敛性可以使仿真结果更符合实际，也更加准确。
　　（3）初始化数学模型参数。将进行收敛性检测的模型初始化为数学参数，并将这些参数属性赋给粒子，这里设定：n=1，μ₀=1.031，T₀=40℃。在湍流封闭模式方程组中设定Cu=0.09，Cε1=1.4，Cε2=1.90，δε=1.30。
　　（4）设定仿真条件。根据液压油的一般工作温度，设定一个液压油在交叉管道内冷却的特定工作环境。实验时，假设两入口端的温度不同，左侧管道入口液压油温度为800℃，中间管道入口液压油温度为400℃，左侧液压油流速为20m/s，中间入口处液压油流速也是20m/s。
　　 （5）运行仿真模型，显示交叉管道内液压油的流动形态。图4所示为双管90度交叉管道内液压油模拟图。不同颜色区域表示液压油在管道内粘性作用力产生的压强不同。整个管道内的压强是以颜色的过渡表示的，不同的颜色表示不同的压强。红色表示压强最高，蓝色表示压强最低（本文只能黑白表示）。我们可以根据颜色的划分轻易地分辨出液压油的粘温层次。
　　 根据以上实验步骤，建立60度交叉管道和30度交叉管道三维模型，进行60度双管交叉管道内液压油粘温性仿真模拟实验，得到如图5和图6所示的液压油流动形态模拟结果。
　　图3 双管90度交叉模型收敛检测
　　图4 双管90度交叉管道内液压油模拟
　　图5 双管60度交叉管道内液压油模拟
　　图6 双管30度交叉管道内液压油模拟
　　 从上述实验可以清楚地看到该方法较真实地反映出液压油体在管道内流动的形态，液压油温粘性的层次能够清楚地表现出来。
　　3 结语
　　 本文以改进的SPH方法为基础，结合湍流标准模型，提出了一种基于SPH的管道内液压油流动性模拟方法。该方法是对现有方法的扩展，能够较逼真地展示出管道内液压油粘温的层次视觉效果以及双管道不同角度交叉处的漩涡形态。在此基础上可以进一步研究改善液压油的粘性计算模型，使其粘温特性得到更充分的展示。
　　参考文献：
　　＼[1＼] 时培成，王幼民，王立涛.液压油液数字建模与仿真[J].农业机械学报，2007，38（12）：148-176.
　　.沈阳：东北大学，2007.
　　.IEEE Transactions on Visualization and Computer Graphics，2008，14（4）：794-804.
　　[4] 陈英建. SPH实现管道流体实时可视化及在GPU上的优化[D].秦皇岛：燕山大学，2012
　　[5] 徐士彪，张晓鹏. 交互式水滴效果模拟[J].计算机辅助设计与图形学报，2013，25（8）：1159-1168 .
　　.Journal of Computational  Physics，1994，110（2）：399-406.
　　.In：Proceedings of Pacific Graphics 2004，Seoul，2004：247-256.
　　[8] G R LIU ， M B LIU . 光滑离子流体动力学——一种无网格粒子法[M].韩旭，杨刚，强洪夫，译.长沙：湖南大学出版社，2005.