用于运动锥形阀芯小间隙二维流场计算的动网格更新方法与流程

本发明属于计算流体动力学技术领域，尤其是涉及一种用于锥形阀芯减压阀流场计算的动网格更新方法。

背景技术：

减压阀利用流体通过狭窄截面时产生的节流效应进行减压，通过阀体内阀芯的开度来调节流体的流量和压力，在能源、化工、机械等领域有广泛的使用。尤其是随着氢能等清洁能源的大规模推广，为了提升单位体积下的储能密度，气体的储存压力也越来越高，对减压阀的性能也提出了更为严格的要求。

锥形阀芯是常用的减压阀阀芯形式，阀芯固定在弹簧上，当高压气体加载到阀芯上时，阀芯在气体压力的推动下离开阀体壁面，并形成狭窄的间隙，高压气体通过间隙时产生节流效应，导致压力降低，温度增加。由于阀芯同时受到气体压力和弹簧力的作用，根据牛顿力学定律，阀芯会发生振动，进而导致气体流量和压力的波动，这在实际的服役环境中会影响能源的利用效率，甚至导致泄露、爆炸等更为严重的后果。对减压阀内部流场的研究有助于探究减压阀的减压规律，控制阀芯的振动和流场状态，因此具有重要意义。

计算流体力学方法是研究减压阀内部的流场的重要手段。一般采用有限体积法将流体域划分成很多的离散单元，称之为网格。网格可分为结构化网格和非结构化网格，采用结构化网格可以节省计算资源，提高计算效率和精度。由于阀芯的运动是影响流场状态的关键因素，因此在计算时必须要考虑在内，而阀芯运动时意味着流场的边界发生了改变，原有的结构化网格就必须进行更新。

本方法尤其适用于高压差小流量的减压阀，这类减压阀在稳态下，锥形阀芯与阀体之间的间隙很小，约为其1％，在进行瞬态计算时，fluent现有的动网格更新方法(即光顺法、动态层法和网格重构法)无法获得较好的网格质量，甚至产生负体积网格，导致计算终止。为了能够借助fluent软件进行此类减压阀内部流场计算，急需开发一种新的适用于锥形阀芯减压阀内部流场计算的动网格更新方法。

技术实现要素：

本发明提供了一种用于减压阀内部流场计算的动网格更新方法，在阀芯产生较大位移时仍能保持良好的网格质量，使得流场计算能够顺利进行。

本发明的技术方案如下：

一种用于运动锥形阀芯小间隙二维流场计算的动网格更新方法，所述锥形阀芯位于减压阀中，该方法包括：

s1：将待模拟的三维锥形阀芯前后流道简化为轴对称的二维模型，并对二维模型进行结构化网格划分；

s2：用fluent软件建立锥形阀芯前后流道的流场瞬态计算模型，阀芯在流体力、摩擦力、弹簧弹力和重力的作用下会沿着轴线方向作直线运动；

s3：根据动力学理论，建立阀芯的运动方程；

s4：在fluent软件中进行二维瞬态流场计算，计算当前时间步中作用在阀芯两端的流体力；

s5：在流场下一时间步的瞬态流场计算开始前，在udf中利用计算得到的作用在阀芯两端的流体力，再根据阀芯的运动方程，计算阀芯在当前时间步中发生的振动位移；

s6：根据s5中的振动位移计算结果，更新阀芯整个流道的流体域中相关网格节点的坐标值，并以更新后的流体域网格为基础进行下一时间步的瞬态流场计算。

在上述技术方案基础上，本发明的各步骤可以采用以下优选方式实现。

优选的，步骤s1中，二维模型的整个流道流体域的网格划分均采用四边形结构化网格，划分步骤为：

s11：在二维模型中，将整个流道的流体域划分为入口直管道区、出口直管道区、间隙区、入口过渡区和出口过渡区5个子流体域；其中入口直管道区、出口直管道区、间隙区的子流体域均为矩形，而入口过渡区和出口过渡区的子流体域均为截角三角形；

s12：以间隙区子流体域作为基础块，然后将该基础块中与阀芯边界平行的一组对边分别向两个过渡区进行延伸，直至达到两个过渡区的子流体域边界，将每个过渡区分别划分为两个块；

s13：然后再以所述对边与每个过渡区流体域边界的两个交点为起点，分别沿垂直于该起点所在的直管道区子流体域边界的方向继续延伸，直至达到直管道区另一侧的子流体域边界，将入口直管道区划分为两个块，而出口直管道区划分为三个块；

s14：再将每个过渡区中的三角形块进行y型2d块切分，然后以切分后新生成的且不位于整个流道的流体域边界上的每个顶点为起点，分别沿垂直于该起点所在的块边界的方向继续延伸，直至达到整个流道的流体域边界，所述块边界不包含y型2d块切分过程中新产生的三条分割线，完成整个流道的流体域中块的划分；

s15：根据整个流道的流体域中划分的所有块，对每个区域进行网格划分，其中间隙区的网格节点进行网格加密处理，形成加密网格，同时两个过渡区和两个直管道区中也对应形成加密网格；间隙区的网格节点在垂直边界线和沿边界线两个方向上均保持等距离分布。

优选的，步骤s1中用icem软件对二维模型进行结构化网格划分。

优选的，步骤s3中，阀芯为刚体，且只有沿着轴线方向的一个运动自由度，阀芯的运动方程包括流体力、摩擦力、重力及弹簧弹力，其中流体力由fluent软件计算得到的压力数据进行积分求得。

优选的，步骤s4中，利用fluent软件计算作用在阀芯两端的流体力，计算结果保存在fluent软件的udf接口中。

优选的，步骤s5中，在fluent软件的udf接口中利用s4中计算得到的作用在阀芯两端的流体力，根据s3中建立的阀芯运动方程计算当前时间步阀芯的位移，同时将当前时间步结束时的状态向量保存在另一个数据文件中，作为下次计算的初始状态向量。

优选的，其特征在于，步骤s6中，整个流道的流体域中相关网格节点的坐标值更新方法具体为：

s61：针对间隙区进行网格移动，网格移动方法为：将区域中距离阀芯最远的网格线上的节点保持不动，并以这条网格线为参考线，以其余网格线上节点到参考线的距离从近到远，依次将各条网格线上的节点记为第1、2、3、…、n层，同一层中的各节点到参考线的距离相等；根据s5中计算得到的阀芯振动位移，对于间隙区中的任一节点，该节点需要移动的位移计算公式为：

按照每个节点需要移动的位移，更新各节点的坐标，实现流体域的网格移动；

s62：针对两个过渡区贴近锥形阀芯的加密网格，也按照s61所述的网格移动方法计算节点位移，两个过渡区中加密网格之外的其余网格节点不移动；

s63：针对出口直管道区中的加密网格，也按照s61所述的网格移动方法计算节点位移；出口直管道区中除加密网格以外的其余网格节点，进行适应性调整。

进一步的，所述的s63中，适应性调整方法为：

以出口直管道区内具有最大移动距离的一层节点为第n层节点，距离其最近的一层节点为第n-1层节点，依次类推，取第0层为参考线，该参考线上的节点保持不动，第1到n-1层的节点位移根据s61中的位移计算公式计算每个节点需要移动的位移，并更新节点的坐标。

优选的，所述的锥形阀芯位于高压差小流量的减压阀中。

进一步的，所述减压阀中的压差为70mpa，阀内流量为1g/s，阀芯处的通道间隙为0.01mm。

本发明的另一目的在于提供一种减压阀内部流场模拟方法，其在模拟过程中可利用前述任一方案所述的动网格更新方法更新模型网格。

与现有技术相比，本发明具有以下有益效果：

1、本发明采用完全的结构化网格划分方法，能有效的节省计算资源，提高计算的效率和精度。

2、本发明提出的动网格更新方法，能保证生成网格的质量，在边界移动过程中，新生成的网格保持着较高的网格质量，不会产生网格畸变，由此保证计算的顺利执行。

附图说明

图1为本发明实施例中的锥形阀芯减压阀模型。

图2为计算步骤流程图。

图3为本发明实施例中抽象出的二维轴对称模型，中心线为对称轴。

图4为本发明实施例中二维轴对称模型的结构化网格划分图。

图5为本发明实施例中网格划分时的等价拓扑模型。

图6为本发明实施例中初始的结构化网格划分及局部网格放大图。

图7为本发明实施例中边界分别向左和向右移动后，图6中区域①的局部网格放大图。

图8为本发明实施例中边界分别向左和向右移动后，图6中区域②的局部网格放大图。

图9为本发明实施例中边界分别向左和向右移动后，图6中区域③的局部网格放大图。

图10为本发明实施例中边界分别向左和向右移动后，图6中区域④的局部网格放大图。

如图11为计算过程流程图。

具体实施方式

下面结合附图和实施例对本发明做进一步详细描述，需要指出的是，以下所述实施例旨在便于对本发明的理解，而对其不起任何限定作用。如图1所示是本实施例中所采用的减压阀模型，包括阀体1、阀芯2和弹簧3。

高压气体从阀体1的入口进入，锥形阀芯2在气体压力的推动下离开阀体壁面，并形成狭窄的间隙，高压气体通过间隙时产生节流效应，导致压力降低，温度增加。

本发明针对有锥形阀芯，狭窄间隙的高压差小流量减压阀流场计算，提供了一种流场边界及动网格的更新方法，计算步骤如图2所示。该用于运动锥形阀芯小间隙二维流场计算的动网格更新方法，具体步骤如下：

s1：将待模拟的三维锥形阀芯前后流道简化为轴对称的二维模型，阀芯处的通道间隙尺寸相对于入口直径而言极其狭小，约为入口直径的1％。图3为流场的二维轴对称模型，加粗的部分代表阀芯边界线，计算时，这部分会发生水平方向的移动，图中阀芯与阀体的间隙仅为0.01mm，而入口直径为1mm，是间隙尺寸的一百倍。

二维模型构建完毕后，需要对二维模型进行结构化网格划分。在本实施例中，二维模型的整个流道流体域的网格划分均采用四边形结构化网格，划分步骤为：

s11：在二维模型中，根据几何特征的不同，可以将整个流道的流体域划分为入口直管道区、出口直管道区、间隙区、入口过渡区和出口过渡区5个子流体域。其中入口直管道区、出口直管道区、间隙区的子流体域均为矩形，而入口过渡区和出口过渡区的子流体域均为截角三角形(三角形靠近间隙区的一角被截取，实际是一个四边形)。流体按如下顺序流经各个区域的顺序为：入口直管道区、入口过渡区、间隙区、出口过渡区、出口直管道区。

s12：以间隙区子流体域作为基础块，由于间隙区子流体域是一个矩形，因此这个块的边界也是矩形，具有两组对边，一组对边是与阀芯边界平行的，而另一组边界是与阀芯边界垂直的。然后，将该基础块中与阀芯边界平行的一组对边分别向两个过渡区(即入口过渡区、出口过渡区)进行延伸，直至达到两个过渡区的子流体域边界。由此，将入口过渡区和出口过渡区中，每个过渡区实际上被划分为两个块，一个块呈三角形，另一个块呈梯形，梯形块靠近阀芯，而三角形块远离阀芯。总体而言，梯形块的形式与间隙区子流体域的基础块是基本一致的，因此能够保证后续划分网格时两者的网格点能够出现较为均匀的分层，便于后续的网格点位置更新。

s13：在前述的每个过渡区中，延伸的对边与该过渡区流体域边界会出现两个交点，因此需要以两个交点为起点，分别沿垂直于该起点所在的直管道区子流体域边界的方向继续延伸，直至达到直管道区另一侧的子流体域边界。由于从两个交点分别延伸出一条分割线，因此将入口直管道区划分成了两个块(其中一条分割线与入口直管道区流体域边界重合)，而出口直管道区划分成了三个块。

s14：再针对入口过渡区和出口过渡区，将每个过渡区中的三角形块进行y型2d块切分，也就是将二维的三角形块通过y型切分分成三个四边形的块。每个三角形块经过切分后，在其三条边上会分别出现1个新生成的顶点(也就是分割线与边界的交点)，其中一个顶点位于整个流道的流体域边界(也就是模型的边界)上，另外两个顶点没有位于整个流道的流体域边界上，也就是处于整个流道内部。然后，以切分后新生成的且不位于整个流道的流体域边界上的每个顶点为起点，分别沿垂直于该起点所在的块边界的方向继续延伸，其中此处所说的块边界不包含y型2d块切分过程中新产生的三条分割线。具体来说，此处，每个过渡区中，有两个顶点需要向相应的直管道区延伸，其中一个顶点位于过渡区和直管道区的交界线上，该顶点沿着垂直该过渡区和直管道区交界线的方向朝直管道区延伸，直至达到整个流道的流体域边界；另一个顶点位于过渡区的三角形块和梯形块交界线上，该顶点沿着垂直该三角形块和梯形块交界线的方向梯形块延伸，直至达到整个流道的流体域边界。当完成上述划分工作后，即完成了整个流道的流体域中块的划分。

当然，上述的划分过程是以构建的二维模型为基础说明的，具体的划分方法可以根据实际调整，以达到该划分结果为准。在本实施例中，网格划分通过icem软件进行。图4为模型的结构化网格划分图。而图5为网格划分时所采用的等价拓扑模型。为了获得更好的网格质量，在图4的三角形过渡区域采用y型2d块切分方式，然后在icem软件中，将图5中的l型2d块切分方式通过点和线的关联操作，即可转化为图4的y型2d块切分。

s15：当整个流道的流体域中各块划分完毕后，根据整个流道的流体域中划分的所有块，即可通过软件设置相应的网格点密度，自动对每个区域进行网格划分。其中需要注意的是，由于在间隙区本身空间较小，因此为了保证模拟的准确性，其空间网格应当相对于其他的位置进行加密。

但是，由于间隙区的块与过渡区中的两个梯形块、两个直管道区中与梯形块邻边相接的矩形块之间，具有延伸性，因此当间隙区的网格节点进行网格加密处理，形成加密网格后，两个过渡区和两个直管道区中也对应形成加密网格。在划分过程中，需要保证间隙区的网格节点在垂直阀芯边界线和沿阀芯边界线两个方向上均保持等距离分布。最终划分出的网格如图6所示，其中间隙区网格宽度约为直管道区网格宽度的十分之一，并且这些加密网格沿着锥形阀芯表面延伸，直到网格边界。但锥形阀芯表面处的加密网格已经不是严格的矩形，而是梯形，梯形网格的高度与间隙区矩形网格宽度保持一致。

s2：用fluent软件建立锥形阀芯前后流道的流场瞬态计算模型，阀芯在流体力、摩擦力、弹簧弹力和重力的作用下会沿着轴线方向作直线运动。

s3：根据动力学理论，建立阀芯的运动方程。

阀芯为刚体，且只有沿着轴线方向的一个运动自由度。阀芯同时受到流体压力、摩擦力、弹簧弹力和重力的作用，其合力由式(1)计算得到：

f＝fp+f+kx+mg(i)

其中fp为流体作用在阀芯上的压力，f为摩擦力，根据不同的阀门结构和密封类型确定，k为弹簧的刚度系数，m为阀芯质量，g为重力加速度。

其中fp由式(2)计算得到。

其中p为流体的压力，从数据文件中读入，平行于阀芯的运动方向，此处选取x轴正方向，为阀芯表面的微分，其方向为垂直于该处的阀芯表面，指向流体域外。

根据牛顿第二定律可以由阀芯的受力和质量计算出位移的二阶导数，见式(3)，进而可以求得阀芯速度和下一时间步的移动距离，见式(4)和式(5)。

其中流体压力、弹簧弹力均与阀芯位移有关，尤其是流体压力受到多种因素的影响，无法简单地表示为关于x的函数，故无法直接求得x随时间变化的函数，需要在fluent中进行离散化的数值计算。由fluent软件计算得到的压力数据进行积分求得流体力。

s4：在fluent软件中进行二维瞬态流场计算，计算当前时间步中作用在阀芯两端的流体力，并将其保存在fluent软件的udf接口的自编程序中。

s5：在流场下一时间步的瞬态流场计算开始前，先在fluent软件udf接口的自编程序中，利用s4中计算得到的作用在阀芯两端的流体力，再根据阀芯的运动方程，计算阀芯在当前时间步中发生的阀芯振动位移，同时将当前时间步结束时的状态向量保存在另一个数据文件中，作为下次计算的初始状态向量。

s6：根据s5中的振动位移计算结果，通过udf接口的自编程序更新阀芯整个流道的流体域中相关网格节点的坐标值，并以更新后的流体域网格为基础进行下一时间步的瞬态流场计算。

整个流道的流体域中相关网格节点的坐标值更新方法具体为：

s61：针对间隙区进行网格移动，网格移动方法为：将间隙区中距离阀芯最远的网格线上的节点保持不动，并以这条网格线为参考线。步骤s1中已说明，该间隙区内网格在垂直边界线方向上等距离分布，因此可以根据节点到参考线的距离将节点分为不同的层，第1层节点距离参考线最近(但不为0)，之后为第2层、第3层以此类推，第n层节点，也就是阀芯边界线处的节点，距离参考线最远，n即为节点的层数，每一层上节点到参考线的距离相等。间隙区中除参考线外的其余网格线，以其余网格线上节点到参考线的距离从近到远，依次将各条网格线上的节点记为第1、2、3、…、n层，同一层中的各节点到参考线的距离相等。

根据s5中计算得到的阀芯振动位移，对于间隙区中的任一节点，首先计算出节点到参考线的距离，然后确定该节点在第几层，则该节点需要移动的位移计算公式为：

按照以上公式，第n层的节点，也就是阀芯界面处节点的位移和计算得到的阀芯位移相同，其它节点根据其距离阀芯的位置逐渐减少。按照每个节点需要移动的位移，更新各节点的坐标，实现流体域的网格移动；

s62：针对两个过渡区贴近锥形阀芯的梯形块的加密网格，也按照s61中的网格移动方法计算节点位移，这里的参考线与(6-1)中的参考线在一条直线上，将计算得到的节点位移加到该节点对应的坐标上，即可实现过渡区加密网格区域的网格移动。两个过渡区中除加密网格之外的其余网格节点不移动。

s63：结构化网格的划分方式，决定了加密网格区会一直延伸到边界，所以在两个直管道区也存在加密网格。此处的节点位移按照如下原则计算：如果参考线延伸方向与阀芯移动方向平行，节点无位移；如果参考线延伸方向与阀芯移动方向垂直，仍按照s61中的网格移动方法计算节点位移，此时的参考线为s61中的参考线在此处延伸后的直线，第n层节点仍然具有最大位移，该层节点所在网格线即为阀芯边界线在此处延伸后的直线。具体而言，在本实施例中，入口直管道区中的网格无位移，而出口直管道区中的网格有位移，因此需要对出口直管道区中的网格进行移动距离计算。

因此，针对出口直管道区中的加密网格，也按照s61中的网格移动方法计算节点位移。二出口直管道区中除加密网格以外的其余网格节点，可进行适应性调整，以免出现网格变形，调整原则是：越靠近出口直管道区内具有最大移动距离的一层节点的网格节点，移动距离越大但不超过该最大移动距离。在本实施例中，以出口直管道区内具有最大移动距离的一层节点为第n层节点，距离其最近的一层节点为第n-1层节点(即图6中第n层节点右侧的相邻层节点)，依次类推直到第0层，取第0层为参考线，该参考线上的节点保持不动，第1到n-1层的节点位移根据s61中的位移计算公式计算每个节点需要移动的位移，并更新节点的坐标。

由此，每当完成网格节点的更新后，即可根据s4步骤进行下一时间步的二维瞬态流场计算得到流体力，然后再根据s5步骤计算得到的阀芯振动位移，再按照s6继续更新网格节点，不断循环直至完成所有时间步的模拟。

按照上述的网格更新流程，其结果如下：图6为初始的结构化网格划分及局部网格放大图，而图7到图10分别为本发明实施例中边界分别向左和向右移动后，图6中区域①②③④的局部网格放大图。可以看出，在边界移动过程中，新生成的网格保持着较高的网格质量。网格更新过程如图11所示。

以上所述的实施例对本发明的技术方案和有益效果进行了详细说明，应理解的是以上所述仅为本发明的具体实施例，并不用于限制本发明，凡在本发明的原则范围内所做的任何修改、补充和等同替换，均应包含在本发明的保护范围之内。