一种基于流体力学连续方程的速度场快速修正方法及装置

摘要

本发明公开了一种基于流体力学连续方程的速度场快速修正方法及装置，速度场快速修正装置，包括连续方程离散形式模块、速度场散度计算模块、求导算子分解模块、速度场重构模块，速度场快速修正方法通过对一维求导算子的矩阵的特征值分解，实现了对修正速度场的直接计算，得到的速度场在离散形式下处处满足连续方程，且是对原始速度场的最保守的修正，本发明能大大提高计算效率，降低了内存的消耗，使得修正方法具有更大的适用性。

说明书

技术领域

本发明涉及一种基于流体力学连续方程的速度场快速修正方法及装置，属于流体力学速度测量技术领域。

背景技术

图像粒子测速技术(Particle Image Velocimetry，PIV)是实验流体力学领域最重要的定量测量手段之一。该技术通过分析示踪粒子在相机上的成像，获得测量区域内的速度场。随着层析粒子图像测速(Tomographic Particle Image Velocimetry，TPIV)等测量技术的发展，三维体空间内全部三个速度分量(Three-dimensional Three-component,3D3C)的获取越来越方便。3D3C速度场包含所有的速度梯度分量，允许使用流体力学的基本方程，对速度场测量精度进行检验。对于不可压缩流体，流动满足连续方程(设u,v,w是三个速度分量；表示速度场的散度)，因此速度场的散度应该处处为零。然而实验测量的结果由于存在误差，所以不能严格满足全场散度为零，有不容忽视的散度残差存在。这些残差不仅会影响流场显示的结果，造成统计量的偏差，还会传播到由速度场重构的压力场中。基于连续方程速度场修正的任务是消除原始速度场散度的残差，或者说从原始测量速度场中重构出一个新的满足连续方程的速度场。前人的研究结果表明，通过对速度场的修正可以减小测量误差10％-20％。

目前，对于速度场散度残差的修正目前已有诸多方法，包括求解泊松方程、无散基函数重构和数学优化方法。一种基于非线性优化问题的最小二乘法修正方法是最简便的方法，并被数值与实验测试证明能够有效地改善湍流统计量。这种最小二乘法修正方法以连续方程为约束条件，目标函数为修正量的二范数的优化修正方法，旨在寻求满足连续方程的最小修正。基于优化方法的修正方案，可以保证对原始速度场进行最小的修正，从而使修正速度场忠于原始实验数据。然而，目前求解这种优化问题的方法既耗内存，效率又低，缺乏实用性，对计算机硬件要求较高，很能难用于实际的速度场修正。其他的非优化的修正方法不能保证是对原始速度场的最小修正，因此可能造成对速度场的过度修正，从而引入额外的误差。

发明内容

本发明的目的是为了解决上述优化问题，提出一种快速、精确的速度场修正方法及装置，降低了对计算机硬件的要求，使速度场修正更容易实现。

本发明的一种基于流体力学连续方程的速度场快速修正方法，通过对一维求导算子的矩阵的特征值分解，实现了对修正速度场的直接计算，得到的速度场在离散形式下处处满足连续方程，且是对原始速度场的最保守的修正，本发明能大大提高计算效率，降低了内存的消耗，使得修正方法具有更大的适用性，具体包括以下几个步骤：

步骤101：确定连续方程的离散形式；

步骤102：计算原始速度场的散度；

步骤103：对矩阵D_nxD_nx^T,D_nyD_nx^T,D_nzD_nx^T进行特征值分解；

步骤104：完成对速度场的修正。

本发明的一种基于流体力学连续方程的速度场快速修正装置，包括连续方程离散形式模块、速度场散度计算模块、求导算子分解模块、速度场重构模块；

连续方程离散形式模块通过对速度场网格坐标进行分析，得到连续方程的离散形式，将散度算子通过一维求导算子表示，同时得到一维求导算子矩阵的形式。

速度场散度计算模块通过连续方程的离散形式，计算原始速度场散度的残差。

求导算子分解模块通过对一维求导算子进行特征值分解，得到该算子的特征值与特征向量。

速度场重构模块根据原始速度场散度的残差，利用算子的特征值与特征向量重构出满足连续方程的速度场，即修正速度场。该模块可以通过计算机的并行计算实现，大大提高了计算效率。

本发明的优点在于：

(1)本发明方法能够有效的消除测量误差造成的速度场散度的残差，修正后的速度场处处满足离散形式的连续性方程；

(2)本发明方法是对原始速度场的最小修正，修正速度场高度忠于原始实验数据；

(3)本发明方法通过对散度算子的降维分解，实现对速度场的快速修正；方法简单、直接，容易程序实现；

(4)本发明方法对内存占用较低，允许对三维大空间内的速度场快速修正；

(5)本发明方法可以通过矩阵运算实现，便于更快速的并行计算；

(6)本发明提供的速度场修正方法及装置，能够实现对实验测量的三维速度场的快速修正，修正速度场忠实于原始实验数据，又满足流体力学的控制方程之一：连续方程，从而能更好地反应真实流场的特征。

附图说明

图1速度场修正流程图；

图2速度场修正效果图。

图3本发明实施例装置的结构示意图。

具体实施方式

下面将结合附图和实施例对本发明作进一步的详细说明。

实施例一、

本发明的一种基于流体力学连续方程的速度场快速修正方法，流程如图1所示，具体实施步骤如下：

步骤101：确定连续方程的离散形式；

不可压缩流动的TPIV测量结果是分布在nx×ny×nz规则结构网格结点上的3D3C瞬时速度场，Δx,Δy,Δz是x,y,z三个方向的网格间距。设u,v,w是三个速度分量，流体力学连续方程的形式为该方程涉及速度场的沿着三个坐标方向的偏导数。为了方便后面的运算，首先给出连续方程的离散形式。假设u在网格结点(i,j,k)上的取值为u^ijk，相应的取值为在网格内部结点采用中心差分格式，在边界上采用前向或者后向差分格式，相应的偏导数计算公式得到

$\left(\begin{array}{c}\frac{\partial u}{\partial x}{|}^{1,jk}=\frac{1}{\Delta x}(u^{2,jk}-u^{1,jk})\\ \frac{\partial u}{\partial x}{|}^{ijk}=\frac{1}{2\Delta x}(u^{i+1,jk}-u^{i-1,jk}),(1<i<nx)\\ \frac{\partial u}{\partial x}{|}^{nx,jk}=\frac{1}{\Delta x}(u^{nx,jk}-u^{nx-1,jk})\end{array}\right)---\left(2\right)$

可以简写为

$\frac{\partial u}{\partial x}{|}^{ijk}=\underset{i^{\prime }=1}{\overset{nx}{\Sigma }}{D_{nx}}^{{\mathrm{ii}}^{\prime }}{u}^{i^{\prime }jk}---\left(3\right)$

其中，∑是求和符号，i'是遍历的指标索引，D_nx^ii'是下列矩阵D_nx的i行i'列元素

求导算子矩阵D_nx与速度场的偏导数计算有关，另外两项偏导数也可以类似地写出：

$\frac{\partial v}{\partial y}{|}^{ijk}=\underset{j^{\prime }=1}{\overset{ny}{\Sigma }}{D_{ny}}^{{\mathrm{jj}}^{\prime }}{v}^{{\mathrm{ij}}^{\prime }k}---\left(5\right)$

$\frac{\partial w}{\partial z}{|}^{ijk}=\underset{k^{\prime }=1}{\overset{nz}{\Sigma }}{D_{nz}}^{{\mathrm{kk}}^{\prime }}{w}^{{\mathrm{ijk}}^{\prime }}---\left(6\right)$

其中，j',k'是遍历的指标索引，分别为在网格结点(i,j,k)上的取值，D_ny^jj'D_nx^kk'分别是下面矩阵的j(k)行j'(k')列元素

于是连续方程的离散形式为

$\underset{i^{\prime }=1}{\overset{nx}{\Sigma }}{D_{nx}}^{{\mathrm{ii}}^{\prime }}{u}^{i^{\prime }jk}+\underset{j^{\prime }=1}{\overset{ny}{\Sigma }}{D_{ny}}^{{\mathrm{jj}}^{\prime }}{v}^{{\mathrm{ij}}^{\prime }k}+\underset{k^{\prime }=1}{\overset{nz}{\Sigma }}{D_{nz}}^{{\mathrm{kk}}^{\prime }}{w}^{{\mathrm{ijk}}^{\prime }}=0---\left(9\right)$

最后指出，D_nx,D_ny,D_nz与速度场的偏导数计算有关，这里称为求导算子矩阵。采用不同的离散格式，求导算子矩阵的具体形式也各异，例如采用更高阶精度的中心格式或者迎风格式等。

步骤102：计算原始速度场的散度。

速度场散度通过公式计算得到，S_o^ijk是原始速度场散度在结点(i,j,k)处的数值，利用步骤101给出的D_nx^ii'D_ny^jj'D_nx^kk'可以得到速度场散度；

${S_o}^{ijk}=\underset{i^{\prime }=1}{\overset{nx}{\Sigma }}{D_{nx}}^{{\mathrm{ii}}^{\prime }}{u_o}^{i^{\prime }jk}+\underset{j^{\prime }=1}{\overset{ny}{\Sigma }}{D_{ny}}^{{\mathrm{jj}}^{\prime }}{v_o}^{{\mathrm{ij}}^{\prime }k}+\underset{k^{\prime }=1}{\overset{nz}{\Sigma }}{D_{nz}}^{{\mathrm{kk}}^{\prime }}{w_o}^{{\mathrm{ijk}}^{\prime }}---\left(10\right)$

其中：u_o,v_o,w_o分别是实验测量的原始速度场u_o的三个速度分量，下标o表示原始实验测量的结果。

速度场散度的计算方式同样可以采用不同的离散形式。速度场散度的大小反应实验测量误差的大小。本发明的任务是消除散度的残差，因此首先需要对原始速度场散度的残差进行计算。

步骤103：对矩阵D_nxD_nx^T,D_nyD_nx^T,D_nzD_nx^T进行特征值分解。

$\left(\begin{array}{c}{D}_{nx}{D_{nx}}^T={\Phi }_{nx}{\Lambda }_{nx}{{\Phi }_{nx}}^T\\ D_{ny}{D_{ny}}^T={\Phi }_{ny}{\Lambda }_{ny}{{\Phi }_{ny}}^T\\ D_{nz}{D_{nz}}^T={\Phi }_{nz}{\Lambda }_{nz}{{\Phi }_{nz}}^T\end{array}\right)---\left(11\right)$

其中Φ_nx,Φ_ny,Φ_nz分别是D_nxD_nx^T,D_nyD_nx^T,D_nzD_nx^T的特征向量矩阵，

设相应的Φ_nx^mn,Φ_ny^mn,Φ_nz^mn分别表示Φ_nx,Φ_ny,Φ_nz的m行n列元素。Λ_nx,Λ_ny,Λ_nz是特征值矩阵(对角阵)。Λ_nx^l,Λ_ny^l,Λ_nz^l分别表示Λ_nx,Λ_ny,Λ_nz的第l个对角元。

通过对D_nxD_nx^T,D_nyD_nx^T,D_nzD_nx^T三个矩阵进行分解，得到了正交的特征向量矩阵。利用特征向量矩阵可以方便快速地对速度场进行修正。这三个矩阵维度分别为nx×nx，ny×ny，nz×nz，远远小于速度场的维度nx×ny×nz，因此可以很快地计算特征值与特征向量。而且，由于三个矩阵均是实对称阵，它们的所有特征值均为实数，所有特征向量正交，这为后面的计算提供了方便。

步骤104：依次按照如下公式完成对速度场的修正：

Γ^ijk＝Λ_nxⁱ+Λ_ny^j+Λ_nz^k>

${\mu }^{ijk}=\underset{l,m,n=1}{\overset{nx,ny,nz}{\Sigma }}\left[{{\Phi }_{nx}}^{il}{{\Phi }_{ny}}^{jm}{{\Phi }_{nz}}^{kn}{\left({\Gamma }^{lmm}\right)}^{-1}\underset{i^{\prime },j^{\prime },k^{\prime }=1}{\overset{nx,ny,nz}{\Sigma }}\left({{\Phi }_{nx}}^{i^{\prime }l}{{\Phi }_{ny}}^{j^{\prime }m}{{\Phi }_{nz}}^{k^{\prime }n}{S_o}^{i^{\prime }{j}^{\prime }{k}^{\prime }}\right)\right]---\left(13\right)$

$\left(\begin{array}{c}{u_c}^{ijk}={u_o}^{ijk}-\underset{i^{\prime }=1}{\overset{nx}{\Sigma }}{D_{nx}}^{i^{\prime }i}{\mu }^{i^{\prime }jk}\\ {v_c}^{ijk}={v_o}^{ijk}-\underset{j^{\prime }=1}{\overset{ny}{\Sigma }}{D_{ny}}^{j^{\prime }j}{\mu }^{{\mathrm{ij}}^{\prime }k}\\ {w_c}^{ijk}={w_o}^{ijk}-\underset{k^{\prime }=1}{\overset{nz}{\Sigma }}{D_{nz}}^{k^{\prime }k}{\mu }^{{\mathrm{ijk}}^{\prime }}\end{array}\right)---\left(14\right)$

其中，Φ_nx^mn,Φ_ny^mn,Φ_nz^mn与Λ_nx^l,Λ_ny^l,Λ_nz^l为求导算子的特征值与特征向量，由步骤103得到。S_o^i'j'k'是原始实验速度场散度，由步骤102得到。(i,j,k)(i',j',k')(l,m,n)均为三个方向的坐标索引，取值范围均为(1～nx,1～ny,1～nz)。Γ^ijk与μ^ijk是计算的中间变量，与速度场散度S_o^ijk类似。另外，Γ^ijk一般会出现一个零元素。为了方便后面的公式计算，令Γ^ijk所有零元素均改为1。实际上，Γ^ijk零元素改成任意非零数字均不影响修正结果。可以证明，对Γ^ijk零元素修改后，按照上述公式计算的速度场是满足连续方程的最小修正的结果。u_c^ijk,v_c^ijk,w_c^ijk是修正速度场，即本发明的最终目标。

该步骤运用了步骤103矩阵分解的结果，将对速度场的修正过程通过明确的公式计算来实现。这是本发明得以快速实现的重要原因。另一方面，上述计算公式均是对原始速度场的线性运算，可以通过矩阵乘法运算来实现，这非常利于并行运算。

本发明得到的修正速度场u_c^ijk,v_c^ijk,w_c^ijk，消除了原始速度场散度的残差，对测量误差具有抑制作用，提供了一个更能反映真实流场结构的速度场，为后面对流动机理的研究提供了更好的素材。

实施例二、

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

下面速度场修正的实施是针对如图2a所示的施加高斯噪声的瞬时速度场进行的。也可以采用不同人工噪音源进行测试。

在直接数值模拟(Direct Numerical Simulation，DNS)的一个100×100×16速度场上施加高斯随机噪声。Δx＝Δy＝Δz＝1。高斯随机噪声的均值μ为0，方差σ是速度方差的0.1倍。高斯随机噪声的概率密度函数f(x)为

$f\left(x\right)=\frac{1}{\sqrt{2\pi }\sigma }\exp (-\frac{{(x-\mu )}^2}{2{\sigma }^2})---\left(15\right)$

施加噪声的速度场用以模拟TPIV得到的原始速度场U_o,V_o,W_o。图2a展示了施加噪声后速度场U_o分量云图。下面对该速度场进行修正。

步骤一：

根据连续方程的离散形式，生成三个求导算子矩阵D_nx,D_ny,D_nz。采用中心差分格式，边界上采用前向或后向差分格式。D_nx,D_ny,D_nz的形式如下

步骤二：

根据步骤102的公式(10)计算速度场的散度S_o^ijk。

${S_o}^{ijk}=\underset{i^{\prime }=1}{\overset{nx}{\Sigma }}{D_{nx}}^{{\mathrm{ii}}^{\prime }}{u_o}^{i^{\prime }jk}+\underset{j^{\prime }=1}{\overset{ny}{\Sigma }}{D_{ny}}^{{\mathrm{jj}}^{\prime }}{v_o}^{{\mathrm{ij}}^{\prime }k}+\underset{k^{\prime }=1}{\overset{nz}{\Sigma }}{D_{nz}}^{{\mathrm{kk}}^{\prime }}{w_o}^{{\mathrm{ijk}}^{\prime }}$

其中，D_nx^mn,D_ny^mn,D_nz^mn分别是D_nx,D_ny,D_nz的m行n列元素。

原始速度场散度的概率密度分布曲线如图c所示。可以看到未经过修正的速度场存在不可忽略的散度残差。本发明的任务是在对原始速度场尽可能少修正的前提下，消除这些残差。

步骤三：

对矩阵D_nxD_nx^T,D_nyD_nx^T,D_nzD_nx^T进行特征值分解，如公式(11)所示。得到相应的特征值Λ_nx,Λ_ny,Λ_nz与特征向量Φ_nx,Φ_ny,Φ_nz。Φ_nx,Φ_ny,Φ_nz分别为nx×nx，ny×ny，nz×nz的标准正交阵，这有利于后面的计算。Λ_nx,Λ_ny,Λ_nz是对角阵，且其对角元均是不小于0的实数。由于D_nxD_nx^T,D_nyD_nx^T,D_nzD_nx^T维度为nx×nx，ny×ny，nz×nz，远小于速度场的维度nx×ny×nz，该步骤计算较快。

步骤四：

按照步骤104的公式(11)-(14)依次计算Γ^ijk，μ^ijk与修正速度场u_c^ijk,v_c^ijk,w_c^ijk。Γ^ijk的结果中有一个0元素。为了避免公式(12)对0元素取倒数的运算，这里将零元素改写为1后进行公式(12)的计算。这一操作不会改变修正速度场满足连续方程的性质，而且能够保证对速度场修正量最小。该步骤得到的修正后的速度场u_c^ijk在中间截面上的分布如图b所示。

本发明通过对速度场做最小修正，得到散度处处为零(<10^-14)，即离散形式下满足流体力学的连续方程的速度场。在MATLAB环境下，利用本发明的方法，完成对测试速度场的修正仅需0.1s。经过本发明实施后，速度场噪声有所减小，对原始速度场的修正量较小，如图2a,b。但散度残差的概率密度曲线有明显变化，由图2c变为图2d，基本上消除了原始速度场散度的残差。当测试速度场扩大到的1000×1000×160数据点，完成修正的时间约为60s，计算机占用内存小于20Gb。而这样大小的速度场在传统的最小二乘法修正方法下是很难实现的，因为计算机内存消耗过大。

实施例三、

本发明实施例还提供了一种基于流体力学连续方程的速度场快速修正装置，如图3所示，该装置可位于计算机内部，包括连续方程离散形式模块301、速度场散度计算模块302、离散求导算子分解模块303、速度场重构模块304。

连续方程离散形式模块301通过对速度场网格坐标进行分析，得到连续方程的离散形式，将散度算子通过一维求导算子表示，至少得到一维求导算子矩阵的形式。

速度场散度计算模块302通过连续方程的离散形式，计算原始速度场散度的残差，原始速度场散度的残差包括所有测量结点上速度场散度的取值。

离散求导算子分解模块303通过对一维求导算子进行特征值分解，至少得到该算子的特征值与特征向量，为后面散度算子求逆提供方便。

速度场重构模块304根据原始速度场散度的残差，利用算子的特征值与特征向量重构出满足连续方程的速度场。该模块通过计算机的并行计算实现，大大提高计算的效率。

本发明实施例中，所述各计算模块均可通过计算机中的中央处理器(Central Processing Unit，CPU)、数字信号处理器(Digital Signal Processor，DSP)或可编程逻辑阵列(Field－Programmable Gate Array，FPGA)实现。

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

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

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

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

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