CN116070550A

CN116070550A - 一种改进的基于时间解析piv的重构流场压力场方法

Info

Publication number: CN116070550A
Application number: CN202310209958.0A
Authority: CN
Inventors: 陈朗生; 叶青青; 陈奕宏; 邵雪明
Original assignee: Zhejiang University ZJU
Current assignee: Zhejiang University ZJU
Priority date: 2023-03-07
Filing date: 2023-03-07
Publication date: 2023-05-05
Anticipated expiration: 2043-03-07
Also published as: CN116070550B

Abstract

本发明公开了一种改进的基于时间解析PIV的重构流场压力场方法，属于图像识别和处理领域。方法包括：由时间解析PIV获取的速度场信息出发，通过本征正交分解重构去除测量噪声，完成速度场无散修正；通过不可压流体NS方程获取全场压力梯度，进而构造泊松方程；根据流场湍流度设置边界条件，给定时间连续迭代初值，使用SOR算法按多路径循环迭代求解泊松方程求得流场空间压力；最后向固体表面积分，最终获得流、固全场压力信息。本发明通过引入降阶重构、改进泊松方程求解、向固体表面积分，改善了传统压力场重构受PIV速度测量噪声影响降低压力求解准确度、求解泊松方程耗时、缺乏固体表面压力信息的缺点。

Description

一种改进的基于时间解析PIV的重构流场压力场方法

技术领域

本发明涉及粒子图像识别和处理领域，具体涉及一种改进的基于时间解析PIV(粒子图像测速法，Particle Image Velocimetry)的重构流场压力场方法。

背景技术

压力分布作为流场的重要动力学参数，能反映流场的受力情况，压力的非定常脉动还是造成空化、气动噪声等现象的重要因素。因此实现流场中全局压力信息测量具有十分重要意义。

但是传统的压力测量方法如接触式传感器只能以接触式方式测量局部有限点并且测压点的布置会对流场形成干扰作用。压敏漆 (Pressure Sensitive Paint，PSP)技术却只能用于空气中的固体壁面压力测量，这限制了流场全场压力分布信息的获取。PIV测量技术是在流动显示基础上基于图像处理发展起来的非侵入、全局、瞬时流场速度测量技术，被广泛运用于科学研究与工业测量中。有学者提出将PSP与PIV技术结合同步采集流场速度与压力信息。如CN107655517A公开了一种基于压敏粒子光强测量的空间流体速度压力同步测量系统，可以同步获得流场中流体的压力与速度信息，其原理是将PSP系统的压敏漆吸附在PIV系统的示踪粒子上，制作出压敏粒子。通过在流场中撒布压敏粒子，并使用高速CCD相机以双曝光模式采集压敏粒子图像，后通过对比两次曝光中压敏粒子的位置变化得到流场速度信息，通过光强比-标准压力曲线将压敏粒子的光强信息转化为压力信息从而同步得到流场的压力与速度信息。但是该方法的压敏粒子制作过程繁琐、成本高昂，此外压敏材料依赖氧气含量变化的特性使其在水中无法使用，并且压敏漆的吸附会减弱示踪粒子的光散射特性而降低粒子图像质量从而使流场速度计算误差增大。

不同于将PSP与PIV结合以同步测量流场速度与压力信息，基于时间解析PIV的压力场重构直接基于时间解析PIV采集的高频速度信息，通过Navier-Stokes（N-S）方程推导出流场的压力信息。无需对PIV测试硬件改进，通过后续压力场重构算法，便可在额外获取流场压力信息。近年来，基于时间解析PIV的压力场重构技术的提出，使得能通过实验方式同步得到全流场的速度与压力信息，为后续的流场动力学分析提供了数据基础。

目前通过PIV测量重构流场压力场的方法主要有直接积分法和求解泊松方程法，前者在积分过程容易出现累计误差，仅适用于小雷诺数、压力梯度平缓的流动。如CN201710220822.4公开了一种基于粒子图像测速(PIV)技术的可压缩流体压力场的计算方法，但该文献重点关注可压缩流场(尤其是含有激波的超声速复杂流场)压力场的计算，克服了现有的可压缩流体压力场计算的局限性，但未关注流场中固体壁面的压力信息的获取；且在方程求解步骤中，仅着重于方程的离散建立与边界条件设置，并未考虑如何基于大量PIV数据快速求解多个时刻的压力方程。另外CN201710829729.3公开了一种基于粒子图像测速(PIV)技术的超声速翼型的测力方法，其根据PIV测得的超声速流场的速度场，直接使用MacCormack方法计算对应的密度场和压力场，并未对PIV获取的速度场进行去噪处理，PIV的测量噪声计算过程中会传递到压力场中，从而降低所预测压力场的准确度。可见，现有基于PIV的流场压力场重构技术存在受PIV测量噪声影响，测量速度场不满足无散条件而导致流场压力场求解准确性降低；基于大量PIV数据求解多时刻压力场时间成本高；同时也不能获得流场中固体壁面的压力信息。

因此，目前基于PIV重构压力场技术在解决PIV测量误差影响、快速求解大量时刻流场压力和流场中固壁面压力信息获取方面都有待优化。

发明内容

基于背景中所阐述的问题，本发明的目的在于提供一种改进的基于时间解析PIV的重构流场压力场方法。本发明能解决现有技术高压力梯度下受PIV测量噪声影响较大，求解泊松方程耗时长，无法同步获得固壁面和空间压力的问题。本发明为基于二维速度场实现高精度时间解析压力场重构提取提供了一种优化解决方案。

本发明公开了一种改进的基于时间解析PIV的重构流场压力场方法，具体包括：

1）基于高速CMOS相机成像原理，由时间解析PIV测得流场中的高频速度场信息；

2) 经快照本征正交分解方法(Proper Orthogonal Decomposition，POD)重构速度场，得到去除测量噪声的速度场，完成速度场的无散修正；

3) 基于POD无散修正速度场，通过欧拉法求得流场压力梯度并对压力梯度表达式取散度，并构建基于速度空间导数的泊松方程；

4）设置边界条件，给定时间连续迭代初值，使用逐次超松弛(Successive OverRelaxation，SOR)迭代算法，进行多路径循环迭代求解，求得流场速度网格节点上的压力信息；

5) 基于求得的流场速度网格上的压力与压力梯度，向固壁面积分，求得固壁面压力。

作为本发明的优选方案，步骤2）中的快照POD重构速度场，包括如下步骤：综合流场x、y方向速度分量

、

进行POD分解，按照模态阶数从低到高顺序，在含能高的模态中，确定结构明显，有物理意义的模态用于重构速度场，进而去除速度测量噪声，完成速度场无散修正。

优选的，本发明的步骤3)为基于POD无散修正的速度场，通过欧拉法求得流场压力梯度，构建并求解泊松方程完成流场内部的压力场重构。其中，通过欧拉法求得流场压力梯度为：将流场物质导数分为当地加速度项、对流加速度项与粘性项，分别将POD修正速度场在时间与空间上离散差分，并带入欧拉法的压力梯度表达式求得流场压力梯度

。

步骤3）中的构建基于速度空间导数的泊松方程为：对流场压力梯度

表达式取散度,导出仅含速度空间导数的泊松方程表达式；通过对POD无散修正速度场进行空间差分运算，求得压力拉普拉斯算子

，构建基于速度空间导数的泊松方程。

作为本发明的优选方案，步骤4）所述多路径循环迭代具体为：迭代求解泊松方程的过程中以多次迭代为一个循环周期，每个周期内按不同路径遍历流场中的速度网格节点来更新流场网格节点的压力值；其中每次迭代的终止节点均为下一次迭代的起始节点，并且该周期的最后一次迭代的终止节点需回到该周期的首次迭代的起始节点；以此去除求解泊松方程时因迭代方向单一性引起的累积误差。

本发明与现有技术对比所具有的主要有益效果包括：

1)本发明采用了POD重构技术，因此克服了目前基于重构PIV流场压力场受PIV测量噪声影响无散条件不满足、边界条件稳定性问题，从而重构出更准确的压力场。

2)本发明采用了时间连续初值设置以及SOR迭代格式，因此克服了目前基于PIV重构流场压力场在求解泊松方程中耗时长的困难，能快速重构出多时刻流场压力场。

3）本发明采用了多路径循环迭代方向，因此降低了目前基于PIV重构流场压力场在求解泊松方程中因迭代方向单一而产生的误差，使重构的压力场更准确。

4)本发明通过向固壁面积分获取了固壁面边界的压力，因此提供了目前基于PIV重构流场压力场技术所缺乏的固壁面信息，有利于后续对流场中固体进行振动、声学分析。

附图说明

图1为本发明实施例提供的改进的基于时间解析PIV的重构流场压力场方法的流程图。

图2为本发明实施例提供的时间解析PIV实验布置图。

图3为本发明实施例提供的POD重构前后流场压力梯度对比图。

图4为本发明实施例提供的时间连续迭代与普通迭代压力收敛趋势对比图。

图5为本发明实施例提供的瞬时压力云图。

图6为本发明实施例提供的固体表面压系数力分布图。

图1中：2-1、流场测量视窗，2-2、PIV光学镜组，2-3、高频双腔激光器，2-4、同步器，2-5、数据采集分析系统，2-6、CMOS高速相机。

具体实施方式

这里结合附图内容，以实例来说明本发明方法的具体实施方式。这里需要说明的是，以下仅为本发明方法的一个应用实例，该方法并不限于处理本实例中涉及中的圆柱扰流的PIV压力场重构。

下面通过具体实施例对本发明的方法技术做进一步说明，需要指出的是，以下所述实施例旨在便于对本发明的理解，而对其不起任何限制作用：

如图1所示，本实施例提供了一种改进的基于时间解析PIV的重构流场压力场方法，本发明所述重构流场压力场包括获取流场速度网格节点上的压力信息和固璧表面的压力信息；本发明方法主要包括以下步骤：

步骤1)：基于高速CMOS相机成像原理，由时间解析PIV测得流场中的高频速度信息；

步骤2)：经快照POD重构速度场，去除测量噪声，完成速度场无散修正；

步骤3)：基于POD无散修正速度场，通过欧拉法求得流场压力梯度并构建基于速度空间导数的泊松方程；

步骤4)：设置边界条件，给定时间连续迭代初值，使用SOR迭代算法按多路径循环迭代求解泊松方程，得到流场速度网格节点上的压力信息；

步骤5)：基于求得的流场速度网格上的压力与压力梯度，向固壁面积分，求得固壁面压力。

在本发明的一个具体实施例中，所述的步骤1)采用高速CMOS相机对粒子图像进行高频采集，测试区域由高速激光（Nd：YLF）照射，时间分辨率最高为4kHz。通过物理标定、互相关计算获得高时空分辨率速度场信息。

具体的，根据实验待测流场工况搭建相应的实验平台，其总体布置如图2所示，实验平台包括流场测量视窗2-1、PIV光学镜组 2-2、高频双腔激光器2-3、同步器2-4、数据采集分析系统2-5、CMOS高速相机2-6。

本实施例的流场为水中的圆柱绕流，PIV测试区域如图3所示。

高频双腔激光器2-3发射出的点光源由PIV光学镜组2-2调制成片光源，照亮流场测量视窗2-1，通过同步器2-4将高频双腔激光器2-3和CMOS高速相机2-6同步，数据采集分析系统2-5向同步器 2-4 发出指令开始采集粒子图像，采集到的粒子图像经数据采集分析系统2-5分析处理得到流场的速度信息。

本实施例中，该步骤获得的流场速度信息的形式为速度场矩阵。

本发明的步骤2)的目的是去除PIV测量噪声，得到修正（去除测量噪声）的速度场。在本发明的一个具体实施例中，按如下方法进行实施：

由时间解析PIV获取速度场矩阵后求得流向速度脉动矩阵

与法向速度脉动矩阵

。

将各个时刻的速度脉动矩阵重组为列向量，并将流向方向速度脉动向量

放在上方，法向方向速度脉动向量

放在下方串联成一个综合的速度脉动列向量

，按照时间顺序从左到右将综合的速度脉动列向量组装成总速度脉动矩阵：

；

其中，上标t1~tN分别指的是1~N时刻，下标1~M指的是测量域中的1~M个空间节点。

包含了流场的非定常湍流结构信息，对其进行快照POD分解可以获取流场各个模态

，而各个模态也可以反映出相应能量对应的流动结构；保证重构总能量大于85%，按模态阶数从低到高的顺序，依次选取结构明显与流动结构对应的模态，进行POD重构滤去PIV测量噪声，完成速度场无散修正。

本实施例中选取前20阶模态按下式进行速度场重构。

其中

表示重构的第k时刻速度场，

表示PIV所测得的平均速度场，n表示参与速度场重构模态的阶数，

表示时间系数。

速度场作为PIV压力场重构的重要输入，经无散修正后带入压力梯度计算以及泊松方程构建，使得整个压力求解过程建立在更符合二维不可压流动的物理状态上，提高了压力求解的精准度。

具体的，在本发明的一个具体实施例中，按如下方法实施步骤3）：

欧拉法计算得到流场压力梯度

的分量形式为：

；

其中

为流场压力，

为流场流体密度，

为流向速度，

为法向速度，

为时间，

为流体动力粘度。

注意到流场压力梯度分量中包含速度时间导数，速度空间导数及速度空间二阶导数。而流场压力梯度作为后续求解泊松方程的第二类边界条件，其求解精度影响压力场求解的准确性，因此其中包含的时间导数及空间导数均按照五点中心差分格式求解，以提高精度。在求流场t时刻x方向速度分量的时间导数

时（

求解一致），先将速度分量

离散到t-2、t-1、t、t+1、t+2共5个时刻，继而基于中心差分格式求得时间导数：

此外，针对测量的起始与终止时段，无法将流场速度向前或向后将速度离散到t-2、t-1时刻，则向后、向前差分构造五点差分格式。

进一步地，速度空间导数

、

、

与

参考时间导数的求解，先在流场空间相应维度（x、y方向）上离散，并按五点中心差分格式进行求解（对于PIV测量边界以及圆柱壁面边界等无法向前或向后离散的空间节点，参考时间导数的起始与终止时段，按向后、向前差分格式求解）。最后，速度空间二阶导数

、

、

与

可以看作空间导数的导数，对已经求得的速度空间导数再在流场空间上进行一次离散，按照相同的五点差分格式计算即可求得。

本实施例为二维不可压流动，经POD无散修正后，速度场满足无散条件，泊松方程可以舍去散度相关项，使得压力求解得到简化。速度时间导数与空间二阶导数消失，仅保留速度空间导数。速度空间导数在前序步骤欧拉法求解流场压力梯度时已经按照五点差分格式求得，带入便求得流场压力的拉普拉斯算子

。图3为基于POD重构前后速度场求得的

对比图，可以看出通过POD无散修正后流场的特征结构更突出，由于测量误差引起的破碎的细小结构被剔除，以此减小了PIV测量噪声对后续压力场重构造成的影响。

完成步骤3)后，开始进行步骤4)中的求解泊松方程，首先需确定方程的边界条件。通过湍流度计算，远离固体壁面且湍流度低的PIV视窗边界为第一类边界，第一类边界节点不参与泊松方程迭代，通过拓展伯努利方程直接给出该节点压力值；湍流度大的PIV视窗边界与固体边界面为第二类边界，给定边界条件为该边界的流场压力梯度值，即由欧拉法计算得到的流场压力梯度值；在泊松方程的每一次求解迭代过程中第二类边界节点的压力由附近节点积分求得。

本实施例在平行于来流的上下测量区域边界处湍流度仅有1%。因此本实施例的边界条件设置为：平行于来流的上下测量区域边界为第一类边界，按照伯努利方程给定该边界的压力值；在垂直于来流的左右测量区域边界以及临近圆柱的节点为第二类边界，给定由步骤3)中由欧拉法求得的流场压力梯度

。

给定边界条件后按照下式所示的SOR迭代算法，计算速度网格的非边界节点的压力值；

；

其中，

为迭代次数，

为松弛因子，

为节点步长，

、

为流场速度网格节点编号。

对于一般的泊松方程求解，

的值在1.2-1.6范围内都有较好的加速迭代收敛的效果；优选地，本实施例

设置为1.3。

考虑到时间解析PIV的高频采样率，两个相邻时刻时间间隔仅为1/4000s，其对应流场流动结构、压力分布相近。通过将上一时刻流场压力的求解结果设置为下一时刻泊松方程求解初值，使求解初值更加靠近方程的解，可以达到加快方程求解收敛速度、稳定其收敛过程的效果。

具体地，泊松方程求解的初值

按时间连续迭代初值设置：

（a）在初始时刻

，求解初值设置为来流压力：

；

（b）在非初始时刻

，求解初值设置为前一时刻

泊松方程的解：

。

按时间连续迭代初值设置与按普通来流压力迭代初值设置，非初始时刻流场中某点压力求解值随迭代次数变化如图4所示：通过图4对比可以看出，本发明使用的时间连续迭代初值设置，在经过5000次迭代后便完成了收敛；而按普通来流压力迭代初值设置需20000次才能达到相同收敛水平，并且迭代过程中压力值的振荡程度大。由此可以看出使用本发明提出的时间连续迭代初值设置使泊松方程求解收敛更快，更稳定，这极大地降低了求解时间。

同时为了避免求解泊松方程时因为迭代方向单一性导致累积误差，提高压力求解精度，本发明通过多路径循环迭代求解泊松方程。本实施例多路径循环迭代的具体实现方式为：在迭代求解泊松方程过程中以四次迭代为一个循环周期；选取流场速度网格的四个边界点为迭代路径起点，每个周期内的每次迭代求解按不同路径从起点开始遍历流场中的速度网格节点来更新流场节点的压力值，其中每次迭代的终止节点回到起点之外的另一个角点，在第四次迭代完成时，其终止节点回到该周期首次迭代的起始节点。

经过迭代，使泊松方程的解收敛，获得流场的压力信息。

本实施例重构的瞬时压力云图如图5所示，从图5可以看出，圆柱绕流的卡门涡街中交替出现的涡结构与圆柱后方两侧间隙出现的的负压区一一对应，这符合卡门涡街的压力分布特性，此外圆柱头部的压力系数也近似为1，也满足圆柱绕流的驻点压力特性。

由于本实施例的流场为水中的圆柱绕流，因此，本实施例的步骤5）为向圆柱壁面积分，获取圆柱壁面压力信息。

首先通过原始粒子图像，确定圆柱圆心在图像中对应的像素坐标，进而转化为流场空间中的几何坐标。

本实施例以5°为步长，按圆柱半径绕圆心一周，构造了72个圆柱壁面边界节点。基于圆心几何坐标计算出圆柱壁面边界节点对应的几何坐标，以此将PIV测量的速度网格与圆柱几何边界耦合。

以临近圆柱几何边界的速度网格节点为积分起点,向圆柱几何边界节点积分求取圆柱壁面压力，积分表达式为：

；

其中，

为圆柱几何边界点的压力值，

、

分别为临近圆柱几何边界的第

个速度网格节点上的压力值和压力梯度值，

为速度网格节点与圆柱几何边界的位移差值，

为临近圆柱几何边界的速度网格节点个数。

优选地，本实施例

，即选择离圆柱边界节点最近的7个点向其积分；

对所有圆柱边界节点都完成积分后便获得圆柱壁面压力信息。

对圆柱壁面压力按照下式无量纲化得到压力系数

；

，

圆柱表面压力系数

分布如图6所示，与传统用压力传感器测得的压力系数曲线的走势吻合良好，60°/270°左右都存在一个极大值，这也对应了圆柱绕流的分离区，0°对应的驻点也满足

。

至此流场内部与固壁面边界的压力全部求出。由实施例可见，本发明通过引入POD重构、改进泊松方程求解、向固体表面积分，改善了传统压力场重构受PIV速度测量噪声影响大、求解泊松方程耗时、缺乏固体表面压力信息的缺点。本发明采用了时间连续初值设置以及SOR迭代格式，克服了目前基于PIV重构流场压力场在求解泊松方程中耗时长的困难，能快速重构出多时刻流场压力场。本发明采用了多路径循环迭代方向，降低了目前基于PIV重构流场压力场在求解泊松方程中因迭代方向单一而产生的误差，使重构的压力场更准确。本发明通过向固壁面积分获取了固壁面边界的压力，提供了目前基于PIV重构流场压力场技术所缺乏的固壁面信息，有利于后续对流场中固体进行振动、声学分析。

以上所述，仅为本发明较佳实施例而已，并非用于限定本发明的保护范围，发明并不限于上述实施方式中的具体细节，在本发明的技术构思范围内，可以对本发明的技术方案进行多种等同变换，这些等同变换均属于本发明的保护范围。

Claims

1.一种改进的基于时间解析PIV重构流场压力场的方法，其特征在于，所述重构流场压力场包括获取流场速度网格节点上的压力信息和固璧表面的压力信息；所述方法包括如下步骤：

1) 基于高速CMOS相机成像原理，由时间解析PIV测得流场中的高频速度信息；

2) 经快照POD重构速度场，去除速度场测量噪声，完成速度场无散修正；

3) 基于POD无散修正速度场，通过欧拉法求得流场压力梯度并构建基于速度空间导数的泊松方程；

4) 设置边界条件，给定时间连续迭代初值，使用SOR迭代算法，通过多路径循环迭代求解泊松方程，得到流场速度网格节点上的压力信息；

2.根据权利要求1所述的改进的基于时间解析PIV重构流场压力场的方法，其特征在于，所述步骤2）中的快照POD重构速度场，包括如下步骤：

综合流场x、y方向速度分量

、

3.根据权利要求1所述的改进的基于时间解析PIV重构流场压力场的方法，其特征在于，所述步骤3）中的通过欧拉法求得流场压力梯度为：将流场物质导数分为当地加速度项、对流加速度项与粘性项，分别将POD修正速度场在时间与空间上离散差分，并带入欧拉法的压力梯度表达式求得流场压力梯度

。

4.根据权利要求1所述的改进的基于时间解析PIV重构流场压力场的方法，其特征在于，所述步骤3）中的构建基于速度空间导数的泊松方程为：对流场压力梯度

,构建基于速度空间导数的泊松方程。

5.根据权利要求1所述的改进的基于时间解析PIV重构流场压力场的方法，其特征在于，所述步骤4）的边界条件，具体为：通过湍流度计算，远离固体壁面且湍流度低的PIV视窗边界为第一类边界，通过拓展伯努利方程直接给出该节点压力值；

湍流度大的PIV视窗边界与固体边界面为第二类边界，给定边界条件为该边界的压力梯度值，即由欧拉法计算得到的流场压力梯度

；在泊松方程的每一次求解迭代过程中第二类边界节点的压力由附近节点积分求得。

6.根据权利要求1所述的改进的基于时间解析PIV重构流场压力场的方法，其特征在于，所述步骤4）的时间连续迭代初值，具体为：按照时间连续给定方程求解初值

，

a) 在初始时刻

，泊松方程求解的求解初值

设置为来流压力

；

b) 在非初始时刻

，泊松方程求解的求解初值

，都设置为前一时刻

泊松方程的压力求解结果

；

其中，

为N时刻泊松方程求解过程中，迭代次数

的压力场矩阵，

为N-1时刻泊松方程求解完成后得到的解，即N-1时刻的流场压力场矩阵。

7.根据权利要求1所述的改进的基于时间解析PIV重构流场压力场的方法，其特征在于，所述步骤4）的SOR迭代算法，具体为：使用下式所示的逐次超松弛SOR迭代格式求解泊松方程；

；

；

其中，

是泊松方程求解过程中P在流场内点

上的第

次迭代的值，

为点

的压力拉普拉斯算子，i、j分别为流场空间y、x方向节点计数指标，

为松弛因子，

为节点步长，m、n分别为流场速度网格的行数与列数。

8.根据权利要求1所述的改进的基于时间解析PIV重构流场压力场的方法，其特征在于，所述步骤4）的通过多路径循环迭代求解泊松方程，具体为：迭代求解泊松方程的过程中以多次迭代为一个循环周期，每个周期内按不同路径遍历流场中的速度网格节点来更新流场网格节点的压力值；其中每次迭代的终止节点均为下一次迭代的起始节点，并且该周期的最后一次迭代的终止节点需回到该周期的首次迭代的起始节点；以此去除求解泊松方程时因迭代方向单一性引起的累积误差。

9.根据权利要求1所述的改进的基于时间解析PIV重构流场压力场的方法，其特征在于，所述步骤4）中的向固壁面积分，求得固壁面压力，具体为：

将固壁面几何边界节点与流场速度网格节点耦合，以与固壁面边界临近的流场速度网格点为积分起点，向固壁面几何边界节点积分求得压力，如下式所示：

；

其中，

为固壁面边界节点上的压力值，

、

分别为临近固体表面边界的第

个速度网格节点上的压力值和压力梯度值，

为速度网格节点与物理几何边界的位移差值，

为临近固壁面边界的速度网格节点个数。