CN103729564B - 一种基于粒子图像测速技术的压力场计算方法和装置 - Google Patents

Abstract

本发明公开了一种基于粒子图像测速（PIV）技术的压力场计算方法，该方法包括：计算得到初始压力梯度矢量场后，计算所述初始压力梯度矢量场的修正矢量场；依据所述初始压力梯度矢量场和修正矢量场计算得到修正后的压力梯度矢量场；对所述修正后的压力梯度矢量场进行积分，得到压力场。本发明还同时公开了一种实现所述方法的装置，从而提高压力场的计算精度，且可简化压力场积分算法的复杂度，提高压力积分效率。

Description

一种基于粒子图像测速技术的压力场计算方法和装置

技术领域

本发明涉及激光测速技术领域中的粒子图像测速（Particle ImageVelocimetry，简称PIV）技术，尤其涉及一种基于PIV技术的压力场计算方法和装置。

背景技术

PIV技术是一种现代激光测速技术，主要运用于流场速度测量，通过追踪示踪粒子在流场中的运动来得到流场速度场。PIV技术通过片光源或者体光源能实现二维或者三维速度场测量，通过高速相机成像还能实现时间解析的PIV实验测量，并利用流体力学控制方程，即N-S方程重构得到和速度场耦合的流场压力场。在压力场重构过程中主要涉及两个技术环节，即：流场压力梯度矢量场的计算和对压力梯度的积分。

目前，所述流场压力梯度矢量场的计算和压力积分方法通常采用多积分路径方法。在其实现过程中，需要首先对速度场进行测量，但是实验测量得到的速度场往往存在误差，而且该误差会在由时间解析的速度场求解压力梯度矢量场的过程中被引入，并进一步影响通过对压力梯度矢量场进行积分得到的压力场。由于压力积分方法比较复杂，很难分析出速度场中的误差如何传播并影响到最后的压力场。可见，现有流场压力场的计算精度并不高。此外，所述多积分路径的压力积分方法需要对流场某点的压力采用大量不同路径的积分来消除随机误差，因此压力积分的工作量大、效率低。

发明内容

为解决现有技术存在的问题，本发明实施例提供一种基于PIV技术的压力场计算方法和装置。

本发明实施例提供了一种基于PIV技术的压力场计算方法，该方法包括：

计算得到初始压力梯度矢量场后，计算所述初始压力梯度矢量场的修正矢量场；依据所述初始压力梯度矢量场和修正矢量场计算得到修正后的压力梯度矢量场；对所述修正后的压力梯度矢量场进行积分，得到压力场。

其中，计算所述初始压力梯度矢量场的修正矢量场的方法，包括：

设置与所述初始压力梯度矢量场对应的修正矢量场，并形成由初始压力梯度矢量场和对应的修正矢量场构成的中心差分格式的无旋方程；

通过拉格朗日乘数法求解所述修正矢量场的最小二范数，依据所述最小二范数和所述中心差分格式的无旋方程形成拉格朗日函数，通过计算所述拉格朗日函数极值问题得到初始压力梯度矢量场的修正矢量场。

优选的，如果流场为三维结构，设流场在（x，y，z）三个坐标方向上的网格节点数分别为N_x,N_y,N_z，则所述初始压力梯度矢量场和对应的修正矢量场构成的中心差分格式的无旋方程为：

$\left\{\begin{array}{cc}\frac{G_{i,j+1,k}^{z^{*}}-G_{i,j-1,k}^{z^{*}}}{2\mathrm{\Δy}}-\frac{G_{i,j,k+1}^{y^{*}}-G_{i,j,k-1}^{y^{*}}}{2\mathrm{\Δz}}-R_{\mathrm{ijk}}^{x0}=0& 1\left(a\right)\\ \frac{G_{i,j,k+1}^{x^{*}}-G_{i,j,k-1}^{x^{*}}}{2\mathrm{\Δz}}-\frac{G_{i+1,j,k}^{z^{*}}-G_{i-1,j,k}^{z^{*}}}{2\mathrm{\Δx}}-R_{\mathrm{ijk}}^{y0}=0& 1\left(b\right)\\ \frac{G_{i+1,j,k}^{y^{*}}-G_{i-1,j,k}^{y^{*}}}{2\mathrm{\Δx}}-\frac{G_{i,j+1,k}^{x^{*}}-G_{i,j-1,k}^{x^{*}}}{2\mathrm{\Δy}}-R_{\mathrm{ijk}}^{z0}=0& 1\left(c\right)\end{array}\right.$

其中， $R_{\mathrm{ijk}}^{x0}=\frac{G_{i,j+1,k}^{z0}-G_{i,j-1,k}^{z0}}{2\mathrm{\Δy}}-\frac{G_{i,j,k+1}^{y0}-G_{i,j,k-1}^{y0}}{2\mathrm{\Δz}}$ 、 $R_{\mathrm{ijk}}^{y0}=\frac{G_{i,j,k+1}^{x0}-G_{i,j,k-1}^{x0}}{2\mathrm{\Δz}}-\frac{G_{i+1,j,k}^{z0}-G_{i-1,j,k}^{z0}}{2\mathrm{\Δx}}$ 、 $R_{\mathrm{ijk}}^{z0}=\frac{G_{i+1,j,k}^{y0}-G_{i-1,j,k}^{y0}}{2\mathrm{\Δx}}-\frac{G_{i,j+1,k}^{x0}-G_{i,j-1,k}^{x0}}{2\mathrm{\Δy}},$ 所述和分别是在直角坐标系（x，y，z）三个坐标轴方向上的分量，所述i、j、k为三个坐标方向上网格节点的索引，其取值范围分别为（1，N_x），（1，N_y）和（1，N_z），所述表示初始压力梯度矢量场，在三个坐标轴方向上的分量用上标分别表示为 所述表示修正矢量场，在三个坐标轴方向上的分量用上标分别表示为 ${\stackrel{\‾}{G}}^{*}=(G^{x^{*}},G^{y^{*}},G^{z^{*}}).$

其中，所述修正矢量场的最小二范数为：

$F={\left|\right|{\stackrel{\‾}{G}}_{\mathrm{ijk}}^{*}\left|\right|}_2=\sqrt{\underset{i,j,k}{\mathrm{\Σ}}[{\left(G_{\mathrm{ijk}}^{x^{*}}\right)}^2+{\left(G_{\mathrm{ijk}}^{y^{*}}\right)}^2+{\left(G_{\mathrm{ijk}}^{z^{*}}\right)}^2]},$ 所述F是的最小二范数。

其中，所述拉格朗日函数为：

$L=F+\underset{i,j,k}{\mathrm{\Σ}}{\mathrm{\λ}}_{\mathrm{ijk}}^x{M}_{\mathrm{ijk}}^x+\underset{i,j,k}{\mathrm{\Σ}}{\mathrm{\λ}}_{\mathrm{ijk}}^y{M}_{\mathrm{ijk}}^y+\underset{i,j,k}{\mathrm{\Σ}}{\mathrm{\λ}}_{\mathrm{ijk}}^z{M}_{\mathrm{ijk}}^z,$ 其中，所述分别为所述无旋方程1(a)、1(b)、1(c)等号的左边项，分别为与之对应的拉格朗日乘子。

其中，所述修正后的压力梯度矢量场的计算方法为：其中，所述为修正后的压力梯度矢量场，所述为初始压力梯度矢量场，所述为修正矢量场。

其中，所述对修正后的压力梯度矢量场进行积分的方法包括：相互正交的积分路径法。

其中，所述对修正后的压力梯度矢量场进行积分的方法包括：单积分路径法或多积分路径法。

本发明实施例还提供了一种基于粒PIV技术的压力场计算装置，该装置包括：矢量场计算模块、修正模块和积分模块；其中，

所述矢量场计算模块，用于计算得到初始压力梯度矢量场后，计算所述初始压力梯度矢量场的修正矢量场；

所述修正模块，用于依据所述初始压力梯度矢量场和修正矢量场计算得到修正后的压力梯度矢量场；

所述积分模块，用于对所述修正后的压力梯度矢量场进行积分，得到压力场。

其中，所述矢量场计算模块计算初始压力梯度矢量场的修正矢量场，为：

设置与所述初始压力梯度矢量场对应的修正矢量场，并形成由初始压力梯度矢量场和对应的修正矢量场构成的中心差分格式的无旋方程；

通过拉格朗日乘数法求解所述修正矢量场的最小二范数，依据所述最小二范数和所述中心差分格式的无旋方程形成拉格朗日函数，通过计算所述拉格朗日函数极值问题得到初始压力梯度矢量场的修正矢量场。

本发明提供的基于PIV技术的压力场计算方法和装置，得到初始压力梯度矢量场后，计算所述初始压力梯度矢量场的修正矢量场；依据所述初始压力梯度矢量场和所述初始压力梯度矢量场的修正矢量场计算得到修正后的压力梯度矢量场；对所述修正后的压力梯度矢量场进行积分，得到压力场。可见，本发明实施例使压力梯度矢量场在离散形式下满足旋度为零，一定程度上减小了速度场中的误差对压力场重构的影响，提高了压力场重构的精度。而且，对修正后的压力梯度矢量场进行积分，积分结果和积分路径无关，只需要非常少的积分路线（理论上只要一条积分路线就能精确还原压力）来消除离散格式的数值误差，减少了积分工作量，大大提高了积分效率。

附图说明

在附图（其不一定是按比例绘制的）中，相似的附图标记可在不同的视图中描述相似的部件。具有不同字母后缀的相似附图标记可表示相似部件的不同示例。附图以示例而非限制的方式大体示出了本文中所讨论的各个实施例。所有带x，y坐标的附图，其物理量均为无量纲量，其中x，y坐标分别代表流向和展向，分别用测量区域长度和宽度无量纲化，速度用自由来流速度无量纲化，矢量场用矢量模的最大值无量纲化，压力场用远前方来流的动压无量纲化。

图1为本发明实施例所述基于PIV技术的压力场计算方法实现流程示意图；

图2为本发明实施例所述修正矢量场的计算方法流程示意图；

图3-a为本发明实施例所述应用场景中t₁时刻实测流场速度矢量图；

图3-b为本发明实施例所述应用场景中t₂时刻实测流场速度矢量图；

图4-a为本发明实施例所述应用场景中实测流场初始矢量场的矢量图；

图4-b为本发明实施例所述应用场景中修正后的矢量场的矢量图；

图5-a为本发明实施例所述应用场景中由初始矢量场得到的压力等值线图；

图5-b为本发明实施例所述应用场景中由修正后的矢量场得到的压力等值线图；

图6为本发明实施例所述修正矢量场的计算装置的结构示意图。

具体实施方式

根据三维不可压缩流体的N-S方程（如下）可知，压力是势函数：

$\▿P=-\mathrm{\ρ}\frac{D\stackrel{\‾}{U}}{\mathrm{Dt}}+\mathrm{\μ}{\▿}^2\stackrel{\‾}{U}+\stackrel{\‾}{f}$

其中，是哈密顿算子，P是压力，ρ是流体的密度，是速度矢量，μ是动力粘性系数，是体积力，是压力梯度，是拉普拉斯算子，是速度的调和函数。

记N-S方程等号右边三项的和为矢量场即则有矢量场满足如下无旋条件：

$\▿\×\▿P=\▿\×\stackrel{\‾}{G}=\stackrel{\‾}{0}$

其中，表示取旋度，即：对矢量场取旋度，代表零矢量。

因为实验测量存在误差，由时间解析的速度场计算得到的矢量场不能严格满足无旋，即因此，本发明实施例提出了一种对PIV速度场计算得到的压力梯度矢量场进行修正，使所述压力梯度矢量场在数值离散形式下严格满足无旋条件，并对其进行正交积分路径法积分得到压力场的方法。

本发明的实施例中：计算得到初始压力梯度矢量场后，计算所述初始压力梯度矢量场的修正矢量场；依据所述初始压力梯度矢量场和所述初始压力梯度矢量场的修正矢量场计算得到修正后的压力梯度矢量场；对所述修正后的压力梯度矢量场进行积分，得到压力场。从而提高压力场的计算精度，且可简化压力场积分算法的复杂度，提高压力积分效率。

下面结合附图及具体实施例对本发明作进一步详细说明。

图1为本发明实施例所述基于PIV技术的压力场计算方法实现流程示意图，如图1所示，包括：

步骤101：计算初始压力梯度矢量场；

这里，可利用PIV技术获得实验流场的时间解析的速度场其时间分辨率为△t，记两个连续的时刻t₁和t₂的速度场分别为和对于任意一个实测流场，在基于PIV的压力场计算过程中，上述速度场和和所对应的时刻t₁和t₂都为已知量。根据上述已知的时间解析的速度场在拉格朗日坐标系下可求得其中间时刻的实质导数，通过空间差分方程可求出中间时刻对应的粘性项，最终由实质导数项、粘性项和体积力项可获得初始压力梯度矢量场，记为此步骤的方法可通过现有技术实现，此处不再详述。

步骤102：计算所述初始压力梯度矢量场的修正矢量场；

如图2所示，该步骤包括：

步骤1021：设置与所述初始压力梯度矢量场对应的修正矢量场，并形成由初始压力梯度矢量场和对应的修正矢量场构成的中心差分格式的无旋方程；

这里，因为初始压力梯度矢量场通常情况下不满足无旋条件，因此需对其进行修正。初始矢量场沿直角坐标系（x，y，z）三个坐标轴方向上的分量用上标分别表示为设修正矢量场为其沿直角坐标系（x，y，z）三个坐标轴方向上的分量用上标分别表示为 初始矢量场和修正矢量场满足因此，流场初始矢量场和修正矢量场的中心差分格式（在流场边界上需要采用单侧差分格式）的无旋方程为：

$\left\{\begin{array}{cc}\frac{G_{i,j+1,k}^{z^{*}}-G_{i,j-1,k}^{z^{*}}}{2\mathrm{\Δy}}-\frac{G_{i,j,k+1}^{y^{*}}-G_{i,j,k-1}^{y^{*}}}{2\mathrm{\Δz}}-R_{\mathrm{ijk}}^{x0}=0& 1\left(a\right)\\ \frac{G_{i,j,k+1}^{x^{*}}-G_{i,j,k-1}^{x^{*}}}{2\mathrm{\Δz}}-\frac{G_{i+1,j,k}^{z^{*}}-G_{i-1,j,k}^{z^{*}}}{2\mathrm{\Δx}}-R_{\mathrm{ijk}}^{y0}=0& 1\left(b\right)\\ \frac{G_{i+1,j,k}^{y^{*}}-G_{i-1,j,k}^{y^{*}}}{2\mathrm{\Δx}}-\frac{G_{i,j+1,k}^{x^{*}}-G_{i,j-1,k}^{x^{*}}}{2\mathrm{\Δy}}-R_{\mathrm{ijk}}^{z0}=0& 1\left(c\right)\end{array}\right.$

其中， $R_{\mathrm{ijk}}^{x0}=\frac{G_{i,j+1,k}^{z0}-G_{i,j-1,k}^{z0}}{2\mathrm{\Δy}}-\frac{G_{i,j,k+1}^{y0}-G_{i,j,k-1}^{y0}}{2\mathrm{\Δz}}$ 、 $R_{\mathrm{ijk}}^{y0}=\frac{G_{i,j,k+1}^{x0}-G_{i,j,k-1}^{x0}}{2\mathrm{\Δz}}-\frac{G_{i+1,j,k}^{z0}-G_{i-1,j,k}^{z0}}{2\mathrm{\Δx}}$ 、 $R_{\mathrm{ijk}}^{z0}=\frac{G_{i+1,j,k}^{y0}-G_{i-1,j,k}^{y0}}{2\mathrm{\Δx}}-\frac{G_{i,j+1,k}^{x0}-G_{i,j-1,k}^{x0}}{2\mathrm{\Δy}},$ 所述和分别是在直角坐标系（x，y，z）三个坐标轴方向上的分量，代表了初始场相对于满足无旋条件的理论场的偏差。假设流场在（x，y，z）三个坐标方向上的网格节点数分别为N_x,N_y,N_z，则所述i、j、k为三个坐标方向上网格节点的索引，其取值范围分别为（1，N_x），（1，N_y）和（1，N_z）。

步骤1022：通过拉格朗日乘数法求解所述修正矢量场的最小二范数，依据所述最小二范数和所述中心差分格式的无旋方程形成拉格朗日函数，通过计算所述拉格朗日函数极值问题得到初始压力梯度矢量场的修正矢量场；

具体的，通过拉格朗日乘数法求解所述修正矢量场在全流场的最小二范数，即对所述全流场进行最小程度的修正，其二范数设为：

$F={\left|\right|{\stackrel{\‾}{G}}_{\mathrm{ijk}}^{*}\left|\right|}_2=\sqrt{\underset{i,j,k}{\mathrm{\Σ}}[{\left(G_{\mathrm{ijk}}^{x^{*}}\right)}^2+{\left(G_{\mathrm{ijk}}^{y^{*}}\right)}^2+{\left(G_{\mathrm{ijk}}^{z^{*}}\right)}^2]},$

形成拉格朗日函数：

$L=F+\underset{i,j,k}{\mathrm{\Σ}}{\mathrm{\λ}}_{\mathrm{ijk}}^x{M}_{\mathrm{ijk}}^x+\underset{i,j,k}{\mathrm{\Σ}}{\mathrm{\λ}}_{\mathrm{ijk}}^y{M}_{\mathrm{ijk}}^y+\underset{i,j,k}{\mathrm{\Σ}}{\mathrm{\λ}}_{\mathrm{ijk}}^z{M}_{\mathrm{ijk}}^z$

其中，F是的最小二范数，分别为上述差分格式的无旋方程1(a)、1(b)、1(c)等号的左边项，分别为与之对应的拉格朗日乘子。

求解拉格朗日函数极值问题，即求解方程组：

$\left\{\begin{array}{cc}\frac{\∂L}{\∂G_{\mathrm{ijk}}^{x^{*}}}=0,& \frac{\∂L}{\∂{\mathrm{\λ}}_{\mathrm{ijk}}^x}=0\\ \frac{\∂L}{\∂G_{\mathrm{ijk}}^{y^{*}}}=0,& \frac{\∂L}{\∂{\mathrm{\λ}}_{\mathrm{ijk}}^y}=0\\ \frac{\∂L}{\∂G_{\mathrm{ijk}}^{z^{*}}}=0,& \frac{\∂L}{\∂{\mathrm{\λ}}_{\mathrm{ijk}}^z}=0\end{array}\right.$

通过上述方程组求解获得全场的修正矢量场

步骤103：依据所述初始压力梯度矢量场和所述初始压力梯度矢量场的修正矢量场计算得到修正后的压力梯度矢量场；

这里，利用步骤102中求解得到的修正矢量场对初始矢量场进行修正，即：求所述与的差值，得到修正后的矢量场

${\stackrel{\‾}{G}}^n={\stackrel{\‾}{G}}^0-{\stackrel{\‾}{G}}^{*}$

即为修正后的差分离散格式的严格满足无旋条件的矢量场，且是满足无旋条件的最近似于初始矢量场的矢量场。

步骤104：对所述修正后的压力梯度矢量场进行积分，得到压力场；

具体的，将步骤103中得到的修正后的矢量场通过方程：进行压力场P的积分。

这里，由于本发明实施例修正后的压力梯度场严格满足无旋条件，因此可以采用一种相互正交的积分路径法，对三维问题用三条两两相互正交的积分路径，而不需采用现有技术中提到的多积分路径方法来进行积分。

以三维问题为例，若流场压力的节点位置用（i，j，k）表示，则三条积分路径分别为从（1，j，k）沿i增大的方向积到（i，j，k）、从（i，1，k）沿j增大的方向积到（i，j，k）和从（i，j，1）沿k增大的方向积到（i，j，k），将从流场边界出发的通过三条相互正交的积分路径积分到（i，j，k）节点的压力平均值作为该节点处的压力值。整个流场的参考压力点可以设置在流场的任何位置，通常情况在节点（1，1，1）处给定一个参考压力并作为积分过程的起始点。对于二维问题与三维问题类似，只需要在二维流场平面内通过“┘”形的正交双积分路径进行压力场积分。这种相互正交的积分路径法较现有方法大大简化了压力场重构的积分过程，减少了压力积分的工作量。

此外，除了正交正积分路径法，修正后的压力梯度场也可以仅采用单积分路径法来进行压力重构，或者多积分路径法，即类似于背景技术提到的多积分路径方法。

可见，本发明实施例使压力梯度矢量场在离散形式下满足旋度为零，一定程度上减小了速度场中的误差对压力场重构的影响，提高了压力场重构的精度。而且，对修正后的压力梯度矢量场进行积分，积分结果和积分路径无关，只需要非常少的积分路线（理论上只要一条积分路线就能精确还原压力）来消除离散格式的数值误差，减少了积分工作量，大大提高了积分效率。

下面结合具体应用场景对本发明方法进行描述。

这里以二维圆柱绕流流场为例进行说明。所述圆柱直径为20mm，二维自由流场的流速度是70mm/s，流场速度场时间分辨率△t为0.01s，密度为水的密度1g/cm³。

对流动进行简化，忽略了粘性力和体积力，不可压缩流体的N-S方程如下：

$\▿P=-\mathrm{\ρ}\frac{D\stackrel{\‾}{U}}{\mathrm{Dt}}$

其中，是哈密顿算子，P是压力，ρ是流体的密度，是速度矢量，为压力梯度。

记N-S方程等号右边一项为矢量场即则有矢量场满足如下无旋条件：

$\▿\×\stackrel{\‾}{G}=\stackrel{\‾}{0}$

其中，表示取旋度，即：对矢量场取旋度。

对于二维问题，所有矢量和微分算子都只有x、y两个方向上的分量。具体二维圆柱绕流流场压力重构步骤如下：

步骤一：利用PIV技术获得实验流场的时间解析的速度场已知时间分辨率为△t=0.01s，记两个连续的时刻t₁和t₂的已知速度场分别为和如图3所示。在拉格朗日坐标系下求得其中间时刻的实质导数，从而得出初始矢量场如图4-a所示。

步骤二：记初始矢量场沿直角坐标系（x，y）两个坐标轴方向上的分量用上标分别表示为记修正矢量场为其沿直角坐标系（x，y）两个坐标轴方向上的分量用上标分别表示为将初始矢量场和修正矢量场代入中心差分格式（在边界上用单侧差分格式进行离散）的二维无旋方程：

$\frac{G_{i+1,j}^{y^{*}}-G_{i-1,j}^{y^{*}}}{2\mathrm{\Δx}}-\frac{G_{i,j+1}^{x^{*}}-G_{i,j-1}^{x^{*}}}{2\mathrm{\Δy}}-R_{\mathrm{ij}}^0=0,$

其中，所述 $R_{\mathrm{ij}}^0=\frac{G_{i+1,j}^{y0}-G_{i-1,j}^{y0}}{2\mathrm{\Δx}}-\frac{G_{i,j+1}^{x0}-G_{i,j-1}^{x0}}{2\mathrm{\Δy}}$ 代表残差。

步骤三：通过拉格朗日乘数法求解修正矢量场的最小二范数，其二范数

设为： $F={\left|\right|{\stackrel{\‾}{G}}_{\mathrm{ij}}^{*}\left|\right|}_2=\sqrt{\underset{i,j}{\mathrm{\Σ}}[{\left(G_{\mathrm{ij}}^{x^{*}}\right)}^2+{\left(G_{\mathrm{ij}}^{y^{*}}\right)}^2]}$

形成拉格朗日函数：

$L=F+\underset{i,j}{\mathrm{\Σ}}{\mathrm{\λ}}_{\mathrm{ij}}{M}_{\mathrm{ij}}$

其中，F是的最小二范数，M_ij为上述二维无旋方程中等号的左边项，λ_ij为与之对应的拉格朗日乘子。

求解拉格朗日函数极值问题得到修正矢量场即求解方程组：

$\frac{\∂L}{\∂G_{\mathrm{ij}}^{*}}=0,\frac{\∂L}{\∂{\mathrm{\λ}}_{\mathrm{ij}}}=0.$

步骤四：用求解得到的对初始矢量场进行修正，得到修正后的矢量场

${\stackrel{\‾}{G}}^n={\stackrel{\‾}{G}}^0-{\stackrel{\‾}{G}}^{*}$

即为最终要得到的满足无旋条件的矢量场，如图4-b所示。比较图4-a和图4-b可以看出，有涡旋位置的矢量场经过无旋修正由之前杂乱的状态，得到了很大改善，变成接近于中心对称的指向涡核中心的矢量场。

步骤五：将步骤四中得到的矢量场根据公式进行积分。通过相互正交的“┘”形的双积分路径求出压力场P。

为便于比较，此处还对未进行无旋修正的初始场进行了压力积分。图5-a为矢量场压力积分等值线图，图5-b为矢量场压力积分等值线图。比较两图可以看出，经过无旋修正，低压中心更加突出，低压区域的压力分布更接近圆柱绕流尾迹流场旋涡的流动特性。

本发明实施例还提供了一种基于PIV技术的压力场计算装置，如图6所示，该装置可位于计算机内部，包括：矢量场计算模块601、修正模块602和积分模块603；其中，

所述矢量场计算模块601，用于计算得到初始压力梯度矢量场，并计算所述初始压力梯度矢量场的修正矢量场；

所述修正模块602，用于依据所述初始压力梯度矢量场和修正矢量场计算得到修正后的压力梯度矢量场；

所述积分模块603，用于对所述修正后的压力梯度矢量场进行积分，得到压力场。

其中，所述矢量场计算模块601计算所述初始压力梯度矢量场，为：

通过PIV技术获得的实验流场时间解析的速度场，在拉格朗日坐标系下求得初始压力梯度矢量场；

其中，所述矢量场计算模块601计算初始压力梯度矢量场的修正矢量场，为：

设置与所述初始压力梯度矢量场对应的修正矢量场，并形成由初始压力梯度矢量场和对应的修正矢量场构成的中心差分格式的无旋方程；

通过拉格朗日乘数法求解所述修正矢量场的最小二范数，并依据所述最小二范数和所述中心差分格式的无旋方程形成拉格朗日函数，通过计算所述拉格朗日函数极值问题得到初始压力梯度矢量场的修正矢量场。

本发明实施例中，所述矢量场计算模块601和修正模块602可通过计算机中的中央处理器（Central Processing Unit，CPU）、数字信号处理器（Digital Signal Processor，DSP）或可编程逻辑阵列（Field－Programmable Gate Array，FPGA）实现，所述积分模块603可通过积分器实现。

本领域内的技术人员应明白，本发明的实施例可提供为方法、装置、或计算机程序产品。因此，本发明可采用硬件实施例、软件实施例、或结合软件和硬件方面的实施例的形式。而且，本发明可采用在一个或多个其中包含有计算机可用程序代码的计算机可用存储介质（包括但不限于磁盘存储器和光学存储器等）上实施的计算机程序产品的形式。

本发明是参照根据本发明实施例的方法、设备（装置）、和计算机程序产品的流程图和/或方框图来描述的。应理解可由计算机程序指令实现流程图和/或方框图中的每一流程和/或方框、以及流程图和/或方框图中的流程和/或方框的结合。可提供这些计算机程序指令到通用计算机、专用计算机、嵌入式处理机或其他可编程数据处理设备的处理器以产生一个机器，使得通过计算机或其他可编程数据处理设备的处理器执行的指令产生用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的装置。

这些计算机程序指令也可存储在能引导计算机或其他可编程数据处理设备以特定方式工作的计算机可读存储器中，使得存储在该计算机可读存储器中的指令产生包括指令装置的制造品，该指令装置实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能。

这些计算机程序指令也可装载到计算机或其他可编程数据处理设备上，使得在计算机或其他可编程设备上执行一系列操作步骤以产生计算机实现的处理，从而在计算机或其他可编程设备上执行的指令提供用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的步骤。

以上所述，仅为本发明的较佳实施例而已，并非用于限定本发明的保护范围。

Claims (9)

1.一种基于粒子图像测速PIV技术的压力场计算方法，其特征在于，该方法包括：

计算得到初始压力梯度矢量场后，计算所述初始压力梯度矢量场的修正矢量场；依据所述初始压力梯度矢量场和修正矢量场计算得到修正后的压力梯度矢量场；对所述修正后的压力梯度矢量场进行积分，得到压力场；

其中，所述修正后的压力梯度矢量场在离散形式下满足无旋条件，采用单积分路径法或相互正交的积分路径法对修正后的压力梯度矢量场进行积分。

2.根据权利要求1所述的方法，其特征在于，计算所述初始压力梯度矢量场的修正矢量场的方法，包括：

设置与所述初始压力梯度矢量场对应的修正矢量场，并形成由初始压力梯度矢量场和对应的修正矢量场构成的中心差分格式的无旋方程；

通过拉格朗日乘数法求解所述修正矢量场的最小二范数，依据所述最小二范数和所述中心差分格式的无旋方程形成拉格朗日函数，通过计算所述拉格朗日函数极值问题得到初始压力梯度矢量场的修正矢量场。

3.根据权利要求2所述的方法，其特征在于，如果流场为三维结构，设流场在(x，y，z)三个坐标方向上的网格节点数分别为N_x,N_y,N_z，则所述初始压力梯度矢量场和对应的修正矢量场构成的中心差分格式的无旋方程为：

$\left\{\begin{array}{cc}\frac{G_{i,j+1,k}^{z*}-G_{i,j-1,k}^{z^{*}}}{2\mathrm{\Δ}y}-\frac{G_{i,j,k+1}^{y*}-G_{i,j,k-1}^{y^{*}}}{2\mathrm{\Δ}z}-R_{ijk}^{x0}=0& 1\left(a\right)\\ \frac{G_{i,j,k+1}^{x*}-G_{i,j,k-1}^{x^{*}}}{2\mathrm{\Δ}z}-\frac{G_{i+1,j,k}^{z*}-G_{i-1,j,k}^{z^{*}}}{2\mathrm{\Δ}x}-R_{ijk}^{y0}=0& 1\left(b\right)\\ \frac{G_{i+1,j,k}^{y*}-G_{i-1,j,k}^{y*}}{2\mathrm{\Δ}x}-\frac{G_{i,j+1,k}^{x*}-G_{i,j-1,k}^{x^{*}}}{2\mathrm{\Δ}y}-R_{ijk}^{z0}=0& 1\left(c\right)\end{array}\right.$

其中， 所述和分别是在直角坐标系(x，y，z)三个坐标轴方向上的分量，所述i、j、k为三个坐标方向上网格节点的索引，其取值范围分别为(1，N_x)，(1，N_y)和(1，N_z)，所述表示初始压力梯度矢量场，在三个坐标轴方向上的分量用上标分别表示为(G^x0，G^y0，G^z0)，所述表示修正矢量场，在三个坐标轴方向上的分量用上标分别表示为

4.根据权利要求3所述的方法，其特征在于，所述修正矢量场的最小二范数为：所述F是的最小二范数。

5.根据权利要求4所述的方法，其特征在于，所述拉格朗日函数为：

其中，所述分别为所述无旋方程1(a)、1(b)、1(c)等号的左边项，分别为与之对应的拉格朗日乘子。

6.根据权利要求1-5中任一项所述的方法，其特征在于，所述修正后的压力梯度矢量场的计算方法为：其中，所述为修正后的压力梯度矢量场，所述为初始压力梯度矢量场，所述为修正矢量场。

7.根据权利要求1-5中任一项所述的方法，其特征在于，所述对修正后的压力梯度矢量场进行积分的方法包括：多积分路径法。

8.一种基于粒子图像测速PIV技术的压力场计算装置，其特征在于，该装置包括：矢量场计算模块、修正模块和积分模块；其中，

所述矢量场计算模块，用于计算得到初始压力梯度矢量场后，计算所述初始压力梯度矢量场的修正矢量场；

所述修正模块，用于依据所述初始压力梯度矢量场和修正矢量场计算得到修正后的压力梯度矢量场；

所述积分模块，用于对所述修正后的压力梯度矢量场进行积分，得到压力场；

其中，所述修正后的压力梯度矢量场在离散形式下满足无旋条件，所述积分模块采用单积分路径法或相互正交的积分路径法对修正后的压力梯度矢量场进行积分。

9.根据权利要求8所述的装置，其特征在于，所述矢量场计算模块计算初始压力梯度矢量场的修正矢量场，为：

设置与所述初始压力梯度矢量场对应的修正矢量场，并形成由初始压力梯度矢量场和对应的修正矢量场构成的中心差分格式的无旋方程；

通过拉格朗日乘数法求解所述修正矢量场的最小二范数，依据所述最小二范数和所述中心差分格式的无旋方程形成拉格朗日函数，通过计算所述拉格朗日函数极值问题得到初始压力梯度矢量场的修正矢量场。