一种两相流泵的自适应网格动态加权误差评估方法

摘要

本发明提供了一种两相流泵的自适应网格动态加权误差评估方法，包括如下步骤：(1)生成模型，初始网格划分并进行CFD数值计算，(2)计算速度梯度值和压力梯度值，(3)计算各单元的梯度值和梯度差值，并归一化处理，(4)误差评估；本发明所述的误差评估方法，将单元的速度梯度和压力梯度结合在一起对单元误差进行评判，从而解决了现有评判方法无法同时考虑速度和压力梯度的不足；在给定误差阈值本发明并不是给定一个固定值，而是在归一化处理后，每次递减误差阈值，从而加快了自适应网格的加密过程；根据本发明方案，可以实现对各种类型网格的进行误差评估，且评估效率高。

权利要求

1.一种两相流泵的自适应网格动态加权误差评估方法，其特征在于，包括如下步骤：(1)通过对两相流泵进行外特性试验，获得两相流泵的扬程和效率数据，生成模型泵三维模型并进行初始网格划分；

(2)进行CFD数值计算；

(3)根据步骤(2)所述的CFD数值计算结果获得的单元i的速度梯度值 和压力梯度值 所述单元i的梯度值为 其中λ1为速度梯度的权重系数，λ2为压力梯度的权重系数，λ1+λ2＝1；

(4)所述单元i的梯度差值为 将单元i的梯度差值进行归一化处理，得归一化梯度差值 其中 为所有单元梯度差值的最大值；

(5)对所述 大于给定的误差阈值的单元进行标定并加密，所述给定的误差阈值在0～1之间；

(6)重复步骤(2)～(5)，且步骤(5)每次给定的误差阈值逐渐减小，直至得到的扬程和效率与步骤(1)所述的扬程和效率相比满足精度要求为止。

2.如权利要求1所述的两相流泵的自适应网格动态加权误差评估方法，其特征在于，步骤(3)所述λ1>λ2。

3.如权利要求1或2所述的两相流泵的自适应网格动态加权误差评估方法，其特征在于，所述单元i的压力梯度值为：其中 分别为所述单元i在X、Y、Z方向上的压力梯度。

4.如权利要求1或2所述的两相流泵的自适应网格动态加权误差评估方法，其特征在于，所述单元i的速度梯度值为：其中 分别为所述单元i在X、Y、Z方向上的速度梯度。

说明书

一种两相流泵的自适应网格动态加权误差评估方法

技术领域

[0001] 本发明涉及两相流泵自适应网格生成领域，特别是涉及一种自适应误差评估方法。

背景技术

[0002] 自适应计算技术是具有某种最优性质的反馈方法，其基本思想在于：以常规数值计算为基础，自适应网格生成为核心，依据计算结果获取误差分布，自动加密误差较大区域，进而优化网格直至误差满足要求，其循环流程如图1所示。自适应网格生成包括误差评估、自适应网格加密和网格优化。误差评估是对网格误差做出有效的估计，并标记出误差较大区域；自适应网格加密则根据标记的信息进行加密；自适应加密后难以保证网格单元几何特性，易出现引起计算发散的畸变单元(体积为零或为负单元)，必须通过网格优化来改善网格单元的几何形状。

[0003] 误差评估是自适应网格生成的核心和关键，直接影响着自适应计算的效率和可靠性。现有方法大都采用相邻点之间伴随变量梯度的变化率来估计误差，认为梯度变化大的区域对流场计算精度的影响也大。研究表明对于压力梯度变化剧烈的地方进行相应的自适应加密能够较好地捕捉尾流涡旋。但对于两相流泵这类几何形状复杂的旋转机械而言，局部细化内部流动不稳定区域(如失速、流动分离等)，压力梯度并不能很好地衡量误差，即压力梯度变化较大并不代表流场计算误差也大，这样有可能造成在某些不需要加密的区域出现过加密的情况。因此，合理地选择流场伴随变量可有效避免因错误判断导致的网格过加密。对于某些流动问题，伴随变量的选择很明确，例如，在冲击流区域压力梯度是一个较优的标准，而对于剪切流问题，速度梯度的评估效果要好于压力梯度。基于不能表征多流动现象的单变量评估，因而易造成评估失效。

[0004] 此外，自适应网格生成需要通过比较误差估计值和阈值来确定需要细化的区域(称为误差单元标记)。固定阈值标记方式是目前较为常用的一种误差单元标记策略，该策略在整个自适应计算过程中保持阈值固定不变，直到所有单元误差都小于该值时停止计算。然而采用固定阈值标记方式进行自适应计算后，仍有少量误差大于给定阈值的单元，为满足误差要求需多次迭代计算，故算法效率较低。

发明内容

[0005] 针对现有技术中存在不足，本发明提供了一种两相流泵的自适应网格动态加权误差评估方法，采用逐渐减小阈值的标记方式，即每次加密，计算单元的数量随着阈值的减小而减少，这样有效避免了只针对某些区域的反复加密计算，因而减少了迭代的次数，能够正确和完整的评估出误差较大单元，从而能够指导加密工作的顺利进行。

[0006] 本发明是通过以下技术手段实现上述技术目的的。

[0007] 一种两相流泵的自适应网格动态加权误差评估方法，包括如下步骤：

[0008] (1)通过对两相流泵进行外特性试验，获得两相流泵的扬程和效率数据，生成模型泵三维模型并进行初始网格划分；

[0009] (2)进行CFD数值计算；

[0010] (3)根据步骤(2)所述的CFD数值计算结果获得的单元i的速度梯度值 和压力梯度值 所述单元i的梯度值为 其中λ1为速度梯度的权
重系数，λ2为压力梯度的权重系数，λ1+λ2＝1；

[0011] (4)所述单元i的梯度差值为 将单元i的梯度差值进行归一化处理，得归一化梯度差值 其中 为所有单元梯度差值的最大值；

[0012] (5)对所述 大于给定的误差阈值的单元进行标定并加密，所述给定的误差阈值在0～1之间；

[0013] (6)重复步骤(2)～(5)，且步骤(5)每次给定的误差阈值逐渐减小，直至得到的扬程和效率与步骤(1)所述的扬程和效率相比满足精度要求为止。

[0014] 进一步，步骤(3)所述λ1>λ2。

[0015] 对于两相流泵而言，相比压力梯度速度梯度表征流场的误差更为有效(尤其在不稳定流动区域)，因此本发明在建立速度和压力梯度误差评估函数时，速度梯度的权重系数大于压力梯度。

[0016] 在上述方案中，所述单元i的压力梯度为：

[0017]

[0018] 其中 分别为所述单元i在X、Y、Z方向上的压力梯度。

[0019] 在上述方案中，所述单元i的速度梯度为：

[0020] 其中 分别为所述单元i在X、Y、Z方向上的速度梯度。

[0021] 单元i的速度梯度与压力梯度计算方法除本发明中提供的方法以外，在现有技术中还有其他计算方法，在此不一一列出。

[0022] 本发明的有益效果：

[0023] (1)针对现有评估方法都采用单一变量评估策略的弊端，本发明采用速度和压力梯度相结合的方式进行评估，故能够更好地反应对计算精度影响较大的单元位置，从而解决了现有评判方法无法同时考虑速度和压力梯度的不足。

[0024] (2)为了解决不同算例情况无法准确给定误差阈值的情况，本发明采用在计算各个单元的梯度差值时采用归一化处理方式，使得所有单元梯度差值都在0～1之间，从而加快了自适应网格的加密过程。

[0025] (3)此外通过动态调整误差阈值的方式对梯度差值大于给定值得单元进行标定并加密，即在第一次给一个较大的误差阈值(接近最大值1)，在随后的迭代过程中的误差阈值逐渐减小直到两相流泵的自适应计算的扬程和效率结果与试验相比满足给定精度为止，进一步提高自适应加密的过程的效率和有效性，有效地避免了不同算例无法精确给定误差阈值的问题。

[0026] (4)根据本发明方案，可以实现对各种类型网格的进行误差评估，且评估效率高。

附图说明

[0027] 图1为本发明所述自适应网格生成过程流程图。

具体实施方式

[0028] 下面结合附图以及具体实施例对本发明作进一步的说明，但本发明的保护范围并不限于此。

[0029] 本发明所述的速度梯度和压力梯度均为单元重心上的(单元可以是四面体也可以是六面体)。

[0030] 如图1为本发明所述自适应网格生成过程流程图。

[0031] 步骤一，搭建两相流泵外特性试验台，泵的扬程有进出口压力表测量得到；采用电测法测量泵的功率，外特性结果得到的最优工况点流量为34.48m3/h扬程为4.76，效率为57.3％，根据两相流泵的水力模型采用Pro/E进行三维建模，并采用ICEM生成初始网格。

[0032] 步骤二，计算采用CFX，标准k-ε湍流模型、SIMPLE算法，非耦合隐式方案进行求解。

[0033] 进口边界条件：采用速度进口，由质量守恒定律和无旋假设确定进口轴向速度。

[0034] 出口边界条件：自由出流，假定出口边界处流动已充分发展，出口区域距离回流区较远。

[0035] 壁面条件：固体壁面采用边壁无滑移条件。

[0036] 扬程：

[0037] 效率：

[0038] 式中：Pout表示泵出口的压力，Pin表示泵进口压力，ρ表示泵内液体的密度，Q表示泵的流量，H表示泵的扬程，g为重力加速度；M′为叶片工作面、背面和前后盖板内、外表面的力矩之和；ω表示泵的角速度；η'为包含容积损失、圆盘损失后泵全流场计算域预测效率值；轴承和密封的损失取3％；因此预测泵的效率η＝η′×(1-3％)。

[0039] 步骤三，根据计算结果获得各个单元的速度梯度值 和压力梯度值并按照公式 计算各个单元的值，其中λ1为速度梯度的权重系
数，λ2为压力梯度的权重系数，λ1+λ2＝1，i为单元的标号，λ1推荐值为0.6，λ2推荐值为0.4；

[0040] 1、单元i的压力梯度值 计算方法如下：

[0041] (1)在X方向上单元i的压力梯度：

[0042] Pi+1为标号i+1单元的压力，Pi-1为标号i-1单元的压力，Δxi+1，i为单元i+1与单元i在X方向的距离(单元重心之间X方向上的距离)，Δxi，i-1为单元i与单元i-1在X方向的距离(单元重心之间X方向上的距离)；

[0043] 所述单元i在Y、Z方向上的压力梯度与所述单元i在X方向上的压力梯度计算方法一致；

[0044] (2)单元i的压力梯度值为：

[0045] 其中 分别为所述单元i在X、Y、Z方向上的压力梯度。

[0046] 2、单元i的速度梯度值 计算方法如下：

[0047] (1)本发明定义X，Y，Z方向的速度分别为u，v，w，则X方向上单元i的速度梯度计算方法如下：

[0048]

[0049] 其中ui+1为标号i+1单元X方向的速度，ui-1为标号i-1单元X方向的速度，Δxi+1,i为单元i+1与单元i在X方向的距离(单元重心之间X方向上的距离)，Δxi,i-1为单元i与单元i-1在X方向的距离(单元重心之间X方向上的距离)。

[0050] (2)X方向的速度梯度大小为：

[0051]

[0052] 对于其他方向(如Y，Z)的计算方法和X方向计算方法一致；

[0053] (3)在单元i的速度梯度大小为：

[0054] 其中 分别为所述单元i在X、Y、Z方向上的速度梯度。

[0055] 3、边界上单元的压力梯度计算方法：

[0056] 左边界上的单元(例如i＝0)，在X方向上单元i的压力梯度计算方法如下：

[0057] 对于其他方向的计算方法和X方向计算方法一致，因此，在单元0的压力梯度大小为：

[0058]

[0059] 右边界上的单元(例如i＝imax，imax表示最大单元数)，在X方向上单元imax的压力梯度计算方法如下， 对于其他方向的计算方法和X方向计算方法一致，因此，在单元imax点的压力梯度大小为：

[0060]

[0061] 4、边界上单元的速度梯度计算方法：

[0062] 左边界上的单元(例如i＝0)，在X方向上单元0的速度梯度计算方法如下：

[0063]

[0064] X方向的速度梯度大小为：

[0065]

[0066] 对于其他方向(如Y，Z)的计算方法和X方向计算方法一致，因此，在单元0的速度梯度大小为：

[0067]

[0068] 右边界上的单元(例如i＝imax，imax表示最大单元数)，在X方向上单元imax的速度梯度计算方法如下：

[0069]

[0070] X方向的速度梯度大小为：

[0071]

[0072] 对于其他方向(如Y，Z)的计算方法和X方向计算方法一致，因此，在单元imax的速度梯度大小为：

[0073]

[0074] 步骤四：根据公式 计算各个单元的梯度差值，并对各个单元的梯度差值进行归一化处理，即 其中 为所有单元梯度差值的最大值；

[0075] 步骤五：根据给定的误差阈值对大于该值的单元进行标定，如果每次都给定一个固定值会增加整个自适应过程的时间，故在第一次自适应过程中给一个较大的误差阈值推荐采用0.95，即标定所有单元的 的值大于0.95，而在随后的计算迭代过程中不断减小该值；

[0076] 步骤六：对所有标定的单元进行加密，即在标定单元各边的中点处添加一个加密节点，在此基础上接着继续进行CFD计算；

[0077] 步骤七：重复步骤二到步骤六，且步骤五中给定的误差阈值逐渐减小(如第二次给定0.9，第三次给定0.88，以此类推)，直到CFD计算的得到的扬程和效率与试验误差小于2％为止(推荐计算误差为2％)，可以结束自适应迭代过程。

[0078] 本实施例在初始网格为86万，自适应加密后得到的网格数为183万后的自适应计算结果与试验的外特性结果对比，扬程误差为1.2％，效率误差为1.6％，计算时间为6小时26分钟；而如果初始网格为185万，计算的结果为，扬程误差为4.2％，效率误差为5.1％，4小时14分钟；而网格数达到398万，才能得到与自适应网格相同的计算精度，且计算时间为16小时43分钟，可见本专利提出的方法能够较好的提高计算效率和计算精度。

[0079] 所述实施例为本发明的优选的实施方式，但本发明并不限于上述实施方式，在不背离本发明的实质内容的情况下，本领域技术人员能够做出的任何显而易见的改进、替换或变型均属于本发明的保护范围。