一种基于流体力学连续方程的速度场快速修正方法及装置的制造方法
【专利摘要】本发明公开了一种基于流体力学连续方程的速度场快速修正方法及装置，速度场快速修正装置，包括连续方程离散形式模块、速度场散度计算模块、求导算子分解模块、速度场重构模块，速度场快速修正方法通过对一维求导算子的矩阵的特征值分解，实现了对修正速度场的直接计算，得到的速度场在离散形式下处处满足连续方程，且是对原始速度场的最保守的修正，本发明能大大提高计算效率，降低了内存的消耗，使得修正方法具有更大的适用性。
【专利说明】
一种基于流体力学连续方程的速度场快速修正方法及装置
技术领域
[0001] 本发明涉及一种基于流体力学连续方程的速度场快速修正方法及装置，属于流体 力学速度测量技术领域。
【背景技术】
[0002] 图像粒子测速技术(Particle Image Velocimetry，PIV)是实验流体力学领域最 重要的定量测量手段之一。该技术通过分析示踪粒子在相机上的成像，获得测量区域内的 速度场。随着层析粒子图像测速（Tomographic Particle Image Velocimetry，TPIV)等测 量技术的发展，三维体空间内全部三个速度分量(Three-dimensional Three-component， 3D3C)的获取越来越方便。3D3C速度场包含所有的速度梯度分量，允许使用流体力学的基本 方程，对速度场测量精度进行检验。对于不可压缩流体，流动满足连续方程V_U = 0(设u，v， w是三个速度分量；
，表示速度场的散度），因此速度场的散度应该处处 为零。然而实验测量的结果由于存在误差，所以不能严格满足全场散度为零，有不容忽视的 散度残差存在。这些残差不仅会影响流场显示的结果，造成统计量的偏差，还会传播到由速 度场重构的压力场中。基于连续方程速度场修正的任务是消除原始速度场散度的残差，或 者说从原始测量速度场中重构出一个新的满足连续方程的速度场。前人的研究结果表明， 通过对速度场的修正可以减小测量误差1 〇 % -20 %。
[0003] 目前，对于速度场散度残差的修正目前已有诸多方法，包括求解泊松方程、无散基 函数重构和数学优化方法。一种基于非线性优化问题的最小二乘法修正方法是最简便的方 法，并被数值与实验测试证明能够有效地改善湍流统计量。这种最小二乘法修正方法以连 续方程为约束条件，目标函数为修正量的二范数的优化修正方法，旨在寻求满足连续方程 的最小修正。基于优化方法的修正方案，可以保证对原始速度场进行最小的修正，从而使修 正速度场忠于原始实验数据。然而，目前求解这种优化问题的方法既耗内存，效率又低，缺 乏实用性，对计算机硬件要求较高，很能难用于实际的速度场修正。其他的非优化的修正方 法不能保证是对原始速度场的最小修正，因此可能造成对速度场的过度修正，从而引入额 外的误差。
【发明内容】
[0004] 本发明的目的是为了解决上述优化问题，提出一种快速、精确的速度场修正方法 及装置，降低了对计算机硬件的要求，使速度场修正更容易实现。
[0005] 本发明的一种基于流体力学连续方程的速度场快速修正方法，通过对一维求导算 子的矩阵的特征值分解，实现了对修正速度场的直接计算，得到的速度场在离散形式下处 处满足连续方程，且是对原始速度场的最保守的修正，本发明能大大提高计算效率，降低了 内存的消耗，使得修正方法具有更大的适用性，具体包括以下几个步骤：
[0006] 步骤101:确定连续方程的离散形式；
[0007] 步骤102:计算原始速度场的散度；
[0008] 步骤103:对矩阵DnxDnxT，D nyDnxT，DnzDnx T进行特征值分解；
[0009]步骤104:完成对速度场的修正。
[0010] 本发明的一种基于流体力学连续方程的速度场快速修正装置，包括连续方程离散 形式模块、速度场散度计算模块、求导算子分解模块、速度场重构模块；
[0011] 连续方程离散形式模块通过对速度场网格坐标进行分析，得到连续方程的离散形 式，将散度算子通过一维求导算子表示，同时得到一维求导算子矩阵的形式。
[0012] 速度场散度计算模块通过连续方程的离散形式，计算原始速度场散度的残差。
[0013] 求导算子分解模块通过对一维求导算子进行特征值分解，得到该算子的特征值与 特征向量。
[0014] 速度场重构模块根据原始速度场散度的残差，利用算子的特征值与特征向量重构 出满足连续方程的速度场，即修正速度场。该模块可以通过计算机的并行计算实现，大大提 高了计算效率。
[0015] 本发明的优点在于：
[0016] (1)本发明方法能够有效的消除测量误差造成的速度场散度的残差，修正后的速 度场处处满足离散形式的连续性方程；
[0017] (2)本发明方法是对原始速度场的最小修正，修正速度场高度忠于原始实验数据；
[0018] (3)本发明方法通过对散度算子的降维分解，实现对速度场的快速修正；方法简 单、直接，容易程序实现；
[0019] (4)本发明方法对内存占用较低，允许对三维大空间内的速度场快速修正；
[0020] (5)本发明方法可以通过矩阵运算实现，便于更快速的并行计算；
[0021] (6)本发明提供的速度场修正方法及装置，能够实现对实验测量的三维速度场的 快速修正，修正速度场忠实于原始实验数据，又满足流体力学的控制方程之一:连续方程， 从而能更好地反应真实流场的特征。
【附图说明】
[0022]图1速度场修正流程图；
[0023]图2速度场修正效果图。
[0024]图3本发明实施例装置的结构示意图。
【具体实施方式】
[0025]下面将结合附图和实施例对本发明作进一步的详细说明。
[0026] 实施例一、
[0027]本发明的一种基于流体力学连续方程的速度场快速修正方法，流程如图1所示，具 体实施步骤如下：
[0028]步骤101:确定连续方程的离散形式；
[0029] 不可压缩流动的TPIV测量结果是分布在nx X ny X nz规则结构网格结点上的3D3C 瞬时速度场，△ X，A y, A z是x，y，z三个方向的网格间距。设u，v，w是三个速度分量，流体力 学连续方程的形式为
。该方程涉及速度场的沿着三个坐标方向的偏导数。 为了方便后面的运算，首先给出连续方程的离散形式。假设U在网格结点（i，j，k)上的取值 为U1#
-相应的取值为
。在网格内部结点采用中心差分格式，在边界上采用前向或者 后向差分格式，相应的偏导数计算公式得到
^2)
[0031] 可以简写为
(°>)
[0033]其中，E是求和符号，i'是遍历的指标索引，DJ1'是下列矩阵Dnx的i行i'列元素
..(4)
[0035]求导算子矩阵Dnx与速度场的偏导数计算有关，另外两项偏导数
：也可以类似 地写出：
(5) (6)
[0038]其中，j'，k'是遍历的指标索引：
分别为
在网格结点（i，j，k)上的 取值，DnyB'Dnxkk'分别是下面矩阵的j(k)行j'（k'）列元素 (7)
(§)
[0041] 于是连续方程
的离散形式为 nx nv m
[0042] %Dju ^ +%D/w^ =0 (9) / =I j -1 k -1
[0043]最后指出，Dnx，Dny，Dnz与速度场的偏导数计算有关，这里称为求导算子矩阵。采用 不同的离散格式，求导算子矩阵的具体形式也各异，例如采用更高阶精度的中心格式或者 迎风格式等。
[0044] 步骤102:计算原始速度场的散度。
[0045] 速度场散度通过公式
计算得到，S。1#是原始速度场散度▽ 在结点（i，j，k)处的数值，利用步骤101给出的Dn，Dn，Dnx k1^可以得到速度场散度； nx nv nz
[0046] ^ (l〇) i_=:l /-I ^-1
[0047]其中：11。，￥。，'\￥。分别是实验测量的原始速度场11。的三个速度分量，下标0表示原始 实验测量的结果。
[0048]速度场散度的计算方式同样可以采用不同的离散形式。速度场散度的大小反应实 验测量误差的大小。本发明的任务是消除散度的残差，因此首先需要对原始速度场散度的 残差进行计算。
[0049] 步骤103:对矩阵DnxDnxT，D nyDnxT，DnzDnx T进行特征值分解。
[0050] D"rD"rr =，n. A,，,,/ (11)
[0051 ]其中 Onx，Ony，Onz分另lj是DnxD nxT，DnyDnxT，D nzDnxT的特征向量矩阵，
[0052] 设相应的 Onxmn，Onymn，O"，分别表示 Onx，Ony，O"2的111行n列元素 。A nx，A ny，A nz 是特征值矩阵(对角阵）。Anx1，Any1, Aj分别表示Anx, Any, Anz的第1个对角元。
[0053] 通过对DnxDnxT，D nyDnxT，DnzDnx T三个矩阵进行分解，得到了正交的特征向量矩阵。利 用特征向量矩阵可以方便快速地对速度场进行修正。这三个矩阵维度分别为nxX nx，nyX ny，nz X nz，远远小于速度场的维度nx X ny X nz，因此可以很快地计算特征值与特征向量。 而且，由于三个矩阵均是实对称阵，它们的所有特征值均为实数，所有特征向量正交，这为 后面的计算提供了方便。
[0054]步骤104:依次按照如下公式完成对速度场的修正：
[0055] rijk= Anxi+Anyj+Anzk (12)
[0058] 其中，0"/'0"广，0"严与八1?1，八耶1，八112 1为求导算子的特征值与特征向量，由 步骤103得到。S〇n'k'是原始实验速度场散度，由步骤102得到。（i，j，k)(i'，j'，k'）（l，m，n) 均为三个方向的坐标索引，取值范围均为（1~nx，l~ny，l~nz)。r ljk与]iljk是计算的中间 变量，与速度场散度S。1#类似。另外，-般会出现一个零元素。为了方便后面的公式计 算，令r$所有零元素均改为1。实际上，r$零元素改成任意非零数字均不影响修正结果。 可以证明，对r $零元素修改后，按照上述公式计算的速度场是满足连续方程的最小修正 的结果。ufk， vjk，w 是修正速度场，即本发明的最终目标。
[0059]该步骤运用了步骤103矩阵分解的结果，将对速度场的修正过程通过明确的公式 计算来实现。这是本发明得以快速实现的重要原因。另一方面，上述计算公式均是对原始速 度场的线性运算，可以通过矩阵乘法运算来实现，这非常利于并行运算。
[0060] 本发明得到的修正速度场^~^。^、^~，消除了原始速度场散度的残差^寸测量 误差具有抑制作用，提供了一个更能反映真实流场结构的速度场，为后面对流动机理的研 究提供了更好的素材。
[0061] 实施例二、
[0062] 下面结合具体应用场景对本发明方法进行描述。
[0063]下面速度场修正的实施是针对如图2a所示的施加高斯噪声的瞬时速度场进行的。 也可以采用不同人工噪音源进行测试。
[0064]在直接数值模拟(Direct Numerical Simulation，DNS)的一个 100 X 100 X 16速度 场上施加高斯随机噪声。A x= A y= A z = l。高斯随机噪声的均值y为0,方差〇是速度方差 的0.1倍。高斯随机噪声的概率密度函数f(x)为
(15)
[0066] 施加噪声的速度场用以模拟TPIV得到的原始速度场U。，V。，W。。图2a展示了施加噪 声后速度场U。分量云图。下面对该速度场进行修正。
[0067] 步骤一：
[0068]根据连续方程的离散形式，生成三个求导算子矩阵Dnx，Dny，Dnz。采用中心差分格 式，边界上采用前向或后向差分格式。Dnx，Dny，Dnz的形式如下
[0071]步骤二：
[0072 ]根据步骤102的公式(10)计算速度场的散度S。1 jk。 nx nv nz-
[0073] Sf = J]DjuJjk ^Djv/k+ i ~\ j -1 k =1
[0074] 其中，Dnxmn，Dnymn，Dnz mr^别是 Dnx，Dny，Dnz 的 m行 n 列元素。
[0075] 原始速度场散度的概率密度分布曲线如图c所示。可以看到未经过修正的速度场 存在不可忽略的散度残差。本发明的任务是在对原始速度场尽可能少修正的前提下，消除 这些残差。
[0076] 步骤三：
[0077] 对矩阵DnxDnxT，DnyDnx T，DnzDnxT进行特征值分解，如公式（11)所示。得到相应的特征 值 Anx, Any, Anz与特征向量C>nx，C>ny，C>nz。C> nx，C>ny，C>nz分别为nxXnx，ny Xny，nz Xnz的 标准正交阵，这有利于后面的计算。Anx，Any，Anz是对角阵，且其对角元均是不小于0的实 数。由于DnxD nxT，DnyDnxT，D nzDnxT维度为nx X nx，ny X ny，nz X nz，远小于速度场的维度nx X ny X nz，该步骤计算较快。
[0078] 步骤四：
[0079] 按照步骤104的公式（11)_(14)依次计算修正速度场Uc^k，Vc^ k，Wc^k。 的结果中有一个0元素。为了避免公式（12)对0元素取倒数的运算，这里将零元素改写 为1后进行公式（12)的计算。这一操作不会改变修正速度场满足连续方程的性质，而且能够 保证对速度场修正量最小。该步骤得到的修正后的速度场ujk在中间截面上的分布如图b 所示。
[0080] 本发明通过对速度场做最小修正，得到散度处处为零(<1(T14)，即离散形式下满足 流体力学的连续方程的速度场。在MATLAB环境下，利用本发明的方法，完成对测试速度场的 修正仅需〇 . Is。经过本发明实施后，速度场噪声有所减小，对原始速度场的修正量较小，如 图2a，b。但散度残差的概率密度曲线有明显变化，由图2c变为图2d，基本上消除了原始速度 场散度的残差。当测试速度场扩大到的1000 X 1000 X 160数据点，完成修正的时间约为60s， 计算机占用内存小于20Gb。而这样大小的速度场在传统的最小二乘法修正方法下是很难实 现的，因为计算机内存消耗过大。
[0081] 实施例三、
[0082]本发明实施例还提供了一种基于流体力学连续方程的速度场快速修正装置，如图 3所示，该装置可位于计算机内部，包括连续方程离散形式模块301、速度场散度计算模块 302、离散求导算子分解模块303、速度场重构模块304。
[0083]连续方程离散形式模块301通过对速度场网格坐标进行分析，得到连续方程的离 散形式，将散度算子通过一维求导算子表示，至少得到一维求导算子矩阵的形式。
[0084]速度场散度计算模块302通过连续方程的离散形式，计算原始速度场散度的残差， 原始速度场散度的残差包括所有测量结点上速度场散度
：的取值。
[0085]离散求导算子分解模块303通过对一维求导算子进行特征值分解，至少得到该算 子的特征值与特征向量，为后面散度算子求逆提供方便。
[0086]速度场重构模块304根据原始速度场散度的残差，利用算子的特征值与特征向量 重构出满足连续方程的速度场。该模块通过计算机的并行计算实现，大大提高计算的效率。 [0087]本发明实施例中，所述各计算模块均可通过计算机中的中央处理器（Central Processing Unit，CPU)、数字信号处理器(Digital Signal Processor，DSP)或可编程逻辑 阵列（Field-Programmable Gate Array，FPGA)实现。
[0088] 本领域内的技术人员应明白，本发明的实施例可提供为方法、装置、或计算机程序 产品。因此，本发明可采用硬件实施例、软件实施例、或结合软件和硬件方面的实施例的形 式。而且，本发明可采用在一个或多个其中包含有计算机可用程序代码的计算机可用存储 介质(包括但不限于磁盘存储器和光学存储器等)上实施的计算机程序产品的形式。
[0089] 本发明是参照根据本发明实施例的方法、设备(装置）、和计算机程序产品的流程 图和/或方框图来描述的。应理解可由计算机程序指令实现流程图和/或方框图中的每一流 程和/或方框、以及流程图和/或方框图中的流程和/或方框的结合。可提供这些计算机程序 指令到通用计算机、专用计算机、嵌入式处理机或其他可编程数据处理设备的处理器以产 生一个机器，使得通过计算机或其他可编程数据处理设备的处理器执行的指令产生用于实 现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的装置。
[0090] 这些计算机程序指令也可存储在能引导计算机或其他可编程数据处理设备以特 定方式工作的计算机可读存储器中，使得存储在该计算机可读存储器中的指令产生包括指 令装置的制造品，该指令装置实现在流程图一个流程或多个流程和/或方框图一个方框或 多个方框中指定的功能。
[0091] 这些计算机程序指令也可装载到计算机或其他可编程数据处理设备上，使得在计 算机或其他可编程设备上执行一系列操作步骤以产生计算机实现的处理，从而在计算机或 其他可编程设备上执行的指令提供用于实现在流程图一个流程或多个流程和/或方框图一 个方框或多个方框中指定的功能的步骤。
[0092] 以上所述，仅为本发明的较佳实施例而已，并非用于限定本发明的保护范围。
【主权项】
1. 一种基于流体力学连续方程的速度场快速修正方法，具体包括以下几个步骤： 步骤101:确定连续方程的离散形式； 步骤102:计算原始速度场的散度； 步骤103:对矩阵DnxDnxT，DnyD nxT，DnzDnxT进行特征值分解； 步骤104:完成对速度场的修正。2. 根据权利要求1所述的一种基于流体力学连续方程的速度场快速修正方法，所述的 步骤101具体为： 设不可压缩流动的ΤΡΙV测量结果是分布在nx X ny X nz规则结构网格结点上的3D3C瞬 时速度场，Δχ, Ay, Δζ是x，y，z三个方向的网格间距， 设u，v，w是三个速度分量，流体力学连续方程为假设u在网格结点（i， 夂^上的取值为!!1',在网格内部结点采用中心差分格式，在边界上 采用前向或者后向差分格式，相应的偏导数计算公式为：其中，Σ是求和符号，i '是遍历的指标索引，D:1'是下列矩阵Dnx的i行i '列元素求导算子矩阵Dnx与速度场的偏导数计算有关，另外两项偏导数其中，j'，k'是遍历的指标索引在网格结点（i, j，k)上的取值， Dnyh'Dnxkk'分别是下面矩阵的j(k)行j'（k'）列元素3. 根据权利要求1所述的一种基于流体力学连续方程的速度场快速修正方法，所述的 步骤102具体为：速度场散度通过公式▽ ·:计算得到，Sjk是原始速度场散度V _ U(?在结点 (i，j，k)处的数值，根据'得到速度场散度；其中:U。，V。，W。分别是实验测量的原始速度场U。的三个速度分量，下标0表示原始实验测 量的结果。4. 根据权利要求1所述的一种基于流体力学连续方程的速度场快速修正方法，所述的 步骤103具体为： 对矩阵DnxDnxT，DnyDnx T，DnzDnxT进行特征值分解；其中：φηχ，Φηγ，Φη4ν别是DnxDnx T，DnyDnxT，DnzD nxT的特征向量矩阵， 设相应的Φηχ?，i>nymn，Φηζ·分别表示Φηχ, i>ny，Φηζ的m行η列元素；Ληχ, Ληγ, Ληζ是特 征值矩阵；Ληχ1, Λη/，Ληζ1分别表示Ληχ, Any, Ληζ的第1个对角元； 通过对DnxDnxT，DnyDnx T，DnzDnxT三个矩阵进行分解，得到了正交的特征向量矩阵；利用特 征向量矩阵对速度场进行修正。5. 根据权利要求1所述的一种基于流体力学连续方程的速度场快速修正方法，所述的 步骤104具体为： 完成对速度场的修正：其中，φηχ' φη/' φηζ?与Λμ1，Λη/，为求导算子的特征值与特征向量，s/ j'k'是 原始实验速度场散度，（i，j，k) (i '，j '，k '）（1，m，η)均为三个方向的坐标索引，取值范围均 为（1~ηχ，1~ny, 1~ηζ); Γ 是计算的中间变量，令Γ ^所有零元素均改为任意非 零数字； 最终得到修正速度场:u。1' v。1' wfk。6. -种基于流体力学连续方程的速度场快速修正装置，包括连续方程离散形式模块、 速度场散度计算模块、求导算子分解模块、速度场重构模块； 连续方程离散形式模块通过对速度场网格坐标进行分析，得到连续方程的离散形式， 将散度算子通过一维求导算子表示，同时得到一维求导算子矩阵的形式； 速度场散度计算模块通过连续方程的离散形式，计算原始速度场散度的残差； 求导算子分解模块通过对一维求导算子进行特征值分解，得到该算子的特征值与特征 向量； 速度场重构模块根据原始速度场散度的残差，利用算子的特征值与特征向量重构出满 足连续方程的速度场，即修正速度场。7. 根据权利要求6所示的一种基于流体力学连续方程的速度场快速修正装置，所述的 速度场散度计算模块通过连续方程的离散形式，计算原始速度场散度的残差，原始速度场 散度的残差包括所有测量结点上速度场散度^·:的取值。
【文档编号】G01P5/00GK105929193SQ201610236214
【公开日】2016年9月7日
【申请日】2016年4月15日
【发明人】高琪, 王成跃, 王晋军
【申请人】北京航空航天大学