CN104063554A - 一种通用的无网格数值结果后处理方法 - Google Patents

Abstract

本发明涉及一种通用的无网格数值结果后处理方法，属于计算力学技术领域。该方法包括以下步骤：1、将无网格方法得到的粒子物理信息储存起来；2、对所研究问题的计算域划分网格，将网格节点信息储存起来；3、利用粒子搜索方法，确定每个网格节点影响域中的无网格粒子；4、计算网格节点影响域内的核函数及其导数；5、计算系数矩阵的逆矩阵；6、计算网格节点上的物理量；7、将网格节点上的物理量导入绘图软件，展示结果。本发明方法能将数值模拟得到的物理信息差值到新的网格节点上，提高了无网格方法计算结果的可视性。

Description

一种通用的无网格数值结果后处理方法

技术领域

本发明涉及一种通用的无网格数值结果后处理方法，属于计算力学技术领域。

背景技术

流体力学是航空航天领域中一个十分重要的学科，如今的流体已不再局限于静态问题，更多的是关注动态问题。传统的CFD方法是有网格方法，这些方法通常基于动网技术或嵌套网格技术来解决网格变形问题，然而当所研究问题的边界变形过大时这两种方法可能会产生网格畸变，此时通常利用网格重构来消除畸变的网格，然而网格重构必然会降低精度和增大计算量。为此人们提出了无网格方法，无网格方法不需要网格因此不会遇到网格畸变的问题。但是现有的绘图软件都不能很好的展示无网格方法的数值结果，因此有必要对无网格方法的数值结果进行特殊的后处理。

2007年Price基于SPH(smooth particle hydrodynamics)的对称特性提出了一种无网格数值结果后处理方法。但是该方法存在以下缺点：

(1)通用性较差，在实际使用过程时需要对结果进行正则化；

(2)该方法精度偏低，只有零阶精度。

发明内容

本发明的目的是为了弥补上述缺陷，提高无网格方法数值结果的可视化特性，提出了一种通用的无网格数值结果后处理方法。

本发明的目的是通过下述技术方案实现的：

一种通用的无网格数值结果后处理方法，具体实现步骤如下：

步骤1、首先将无网格方法得到的粒子物理信息储存起来；

步骤2、对所研究问题的计算域划分网格，将网格节点信息储存起来；

划分网格时可以借助于现有的商业软件，比如Gambit，Icem等，利用商业软件划分网格时需要将所有相关的网格节点信息保存。对于简单的外形，也可以通过编写程序划分网格，利用该方法获得的网格只需要保存节点的坐标。

步骤3、利用粒子搜索方法，确定每个网格节点影响域中的无网格粒子；如图2所示；

网格节点的影响域半径r_a＝kh，其中k由核函数的类型决定，通常取为2或3，h是光滑长度，从无网格方法数值结果的粒子信息中获得。计算每个网格节点与所有无网格粒子之间的距离r＝|x_g-x_p|，x_g和x_p分别表示任意网格节点和无网格粒子的矢量坐标，当r<r_a时，则认为无网格粒子x_p在网格节点x_g的影响域内。

步骤4、计算网格节点影响域内的核函数及其导数；

核函数有很多种，这里仅给出三次样条核函数，公式(1)是三次样条核函数的表达式，公式(2)、(3)、(4)分别是三次样条核函数在x、y、z方向的导数表达式。公式中R是基于光滑长度h进行无量纲化后的粒子间距，R＝r/h，α_d是根据归一性条件确定的常数，α_d在一维、二维和三维中分别取为1/h，15/7h²和3/2h³。x_g、y_g和z_g分别表示网格节点在x、y、z方向的坐标，x_p、y_p和z_p分别表示无网格粒子在x、y、z方向的坐标。类似的，如果选择其它核函数，只需将对应的核函数及其在x、y、z方向的导数赋值给W、W_x、W_y、W_z即可。

$W=a_d\left\{\begin{array}{cc}\frac{2}{3}-R^2+\frac{1}{2}{R}^3& 0\&le;R<1\\ \frac{1}{6}{(2-R)}^3& 1\&le;R<2\\ 0& R\&GreaterEqual;2\end{array}\right.---\left(1\right)$

$W_x=\frac{x_g-x_p}{r}\frac{\&PartialD;W}{\&PartialD;r}---\left(2\right)$

$W_y=\frac{y_g-y_p}{r}\frac{\&PartialD;W}{\&PartialD;r}---\left(3\right)$

$W_z=\frac{z_g-z_p}{r}\frac{\&PartialD;W}{\&PartialD;r}---\left(4\right)$

步骤5、计算系数矩阵的逆矩阵，计算过程如下；

5.1求系数矩阵：

系数矩阵A用公式(5)表示，其中的元素值通过公式(6)求出，(6)式中N是网格节点影响域内无网格粒子的个数，通过步骤(3)计算得到；m_p和ρ_p分别表示无网格粒子的质量和密度，从无网格方法数值结果的粒子信息中获得；

$A=\left[\begin{array}{cccc}{a}_0& b_0& c_0& d_0\\ a_1& b_1& c_1& d_1\\ a_2& b_2& c_2& d_2\\ a_3& b_3& c_3& d_3\end{array}\right]---\left(5\right)$

$\begin{array}{cccc}{a}_0=\underset{p=1}{\overset{N}{\mathrm{\&Sigma;}}}W\frac{m_p}{{\mathrm{\&rho;}}_p}& b_0=\underset{p=1}{\overset{N}{\mathrm{\&Sigma;}}}(x_p-x_g)W\frac{m_p}{{\mathrm{\&rho;}}_p}& c_0=\underset{p=1}{\overset{N}{\mathrm{\&Sigma;}}}(y_p-y_g)W\frac{m_p}{{\mathrm{\&rho;}}_p}& d_0=\underset{p=1}{\overset{N}{\mathrm{\&Sigma;}}}(z_p-z_g)W\frac{m_p}{{\mathrm{\&rho;}}_p}\\ a_1=\underset{p=1}{\overset{N}{\mathrm{\&Sigma;}}}{W}_x\frac{m_p}{{\mathrm{\&rho;}}_p}& b_1=\underset{p=1}{\overset{N}{\mathrm{\&Sigma;}}}(x_p-x_g)W_x\frac{m_p}{{\mathrm{\&rho;}}_p}& c_1=\underset{p=1}{\overset{N}{\mathrm{\&Sigma;}}}(y_p-y_g)W_x\frac{m_p}{{\mathrm{\&rho;}}_p}& d_1=\underset{p=1}{\overset{N}{\mathrm{\&Sigma;}}}(z_p-z_g)W_x\frac{m_p}{{\mathrm{\&rho;}}_p}\\ a_2=\underset{p=1}{\overset{N}{\mathrm{\&Sigma;}}}{W}_y\frac{m_p}{{\mathrm{\&rho;}}_p}& b_2=\underset{p=1}{\overset{N}{\mathrm{\&Sigma;}}}(x_p-x_g)W_y\frac{m_p}{{\mathrm{\&rho;}}_p}& c_2=\underset{p=1}{\overset{N}{\mathrm{\&Sigma;}}}(y_p-y_g)W_y\frac{m_p}{{\mathrm{\&rho;}}_p}& d_2=\underset{p=1}{\overset{N}{\mathrm{\&Sigma;}}}(z_p-z_g)W_y\frac{m_p}{{\mathrm{\&rho;}}_p}\\ a_3=\underset{p=1}{\overset{N}{\mathrm{\&Sigma;}}}{W}_z\frac{m_p}{{\mathrm{\&rho;}}_p}& b_3=\underset{p=1}{\overset{N}{\mathrm{\&Sigma;}}}(x_p-x_g)W_z\frac{m_p}{{\mathrm{\&rho;}}_p}& c_3=\underset{p=1}{\overset{N}{\mathrm{\&Sigma;}}}(y_p-y_g)W_z\frac{m_p}{{\mathrm{\&rho;}}_p}& d_3=\underset{p=1}{\overset{N}{\mathrm{\&Sigma;}}}(z_p-z_g)W_z\frac{m_p}{{\mathrm{\&rho;}}_p}\end{array}---\left(6\right)$

5.2求系数矩阵的逆矩阵A^-1；

矩阵的逆可以直接求解，也可以利用迭代算法求解。直接求解的方法是先求出矩阵A对应的行列式值Γ，再求出矩阵A的伴随阵A^*，最后用伴随阵A^*除以行列式的值Γ即可求出A^-1。对于二维问题直接求解比较方便，但是对于三维问题直接求解比较复杂，建议采用迭代算法求解，常用的迭代方法有高斯消元法、LU分解法。无论采用那种方法都可以将逆矩阵A^-1表示为：

$A^{-1}=\left[\begin{array}{cccc}{A}_0& B_0& C_0& D_0\\ A_1& B_1& C_1& D_1\\ A_2& B_2& C_2& D_2\\ A_3& B_3& C_3& D_3\end{array}\right]---\left(7\right)$

步骤6、计算网格节点上的物理量；

利用公式(8)计算网格节点上的物理量，(8)式中f_g、f_p分别表示网格节点、无网格粒子的任意物理量，比如密度ρ、x方向速度u、y方向速度v、z方向速度w、压力p、温度T。

$f_g=A_0\underset{p=1}{\overset{N}{\mathrm{\&Sigma;}}}{f}_{p}W\frac{m_p}{{\mathrm{\&rho;}}_p}+B_0\underset{p=1}{\overset{N}{\mathrm{\&Sigma;}}}{f}_p{W}_x\frac{m_p}{{\mathrm{\&rho;}}_p}+C_0\underset{p=1}{\overset{N}{\mathrm{\&Sigma;}}}{f}_p{W}_y\frac{m_p}{{\mathrm{\&rho;}}_p}+D_0\underset{p=1}{\overset{N}{\mathrm{\&Sigma;}}}{f}_p{W}_z\frac{m_p}{{\mathrm{\&rho;}}_p}---\left(8\right)$

步骤7、将网格节点上的物理量导入绘图软件，展示结果。

自此，就完成了无网格数值结果的后处理。

有益效果

本发明具有如下优点：

(1)具有一阶精度；

下面从原理上来说明本方法是具有一阶精度的。对任意函数f(x,y,z)，如果在点(x_g,y_g,z_g)的某一邻域内具有n+1阶倒数，那么就可以在该点对函数f(x,y,z)进行n阶泰勒展开。为了表示方便，这里用f，f_g，f_g,x，f_g,y，f_g,z分别表示f(x,y,z)，f(x_g,y_g,z_g)，fx(x_g,y_g,z_g)，f_y(x_g,y_g,z_g)，fz(x_g,y_g,z_g)。对函数f(x,y,z)在点(x_g,y_g,z_g)进行泰勒展开，且只保留一阶项，余项用o(r)表示：

f＝f_g+(x-x_g)f_g,x+(y-y_g)f_g,y+(z-z_g)f_g,z+o(r) (9)

因为公式(9)保留了一阶项，因此该展开式具有一阶精度。在公式(9)两边同时乘以核函数及其三个方向的导数，得到下列方程组：

$\left\{\begin{array}{c}\mathrm{fW}=f_{g}W+f_{g,x}(x-x_g)W+f_{g,y}(y-y_g)W+f_{g,z}(z-z_g)W\\ {\mathrm{fW}}_x=f_g{W}_x+f_{g,x}(x-x_g)W_x+f_{g,y}(y-y_g)W_x+f_{g,z}(z-z_g)W_x\\ {\mathrm{fW}}_y=f_g{W}_y+f_{g,x}(x-x_g)W_y+f_{g,y}(y-y_g)W_y+f_{g,z}(z-y_g)W_y\\ {\mathrm{fW}}_z=f_g{W}_z+f_{g,x}(x-x_g)W_y+f_{g,y}(y-y_g)W_z+f_{g,z}(z-z_g)W_z\end{array}\right.---\left(10\right)$

方程组(10)对(x,y,z)进行体积分可得：

$\left\{\begin{array}{c}\&Integral;\mathrm{fWdV}={\&Integral;f}_g\mathrm{WdV}+f_{g,x}\&Integral;(x-x_g)\mathrm{WdV}+f_{g,y}\&Integral;(y-y_g)\mathrm{WdV}+f_{g,z}\&Integral;(z-z_g)\mathrm{WdV}\\ \&Integral;{\mathrm{fW}}_x\mathrm{dV}={\&Integral;f}_g{W}_x\mathrm{dV}+f_{g,x}\&Integral;(x-x_g)W_x\mathrm{dV}+f_{g,y}\&Integral;(y-y_g)W_x\mathrm{dV}+f_{g,z}\&Integral;(z-z_g)W_x\mathrm{dV}\\ \&Integral;{\mathrm{fW}}_y\mathrm{dV}=\&Integral;f_g{W}_y\mathrm{dV}+f_{g,x}\&Integral;(x-x_g)W_y\mathrm{dV}+f_{g,y}\&Integral;(y-y_g)W_y\mathrm{dV}+f_{g,z}\&Integral;(z-y_g)W_y\mathrm{dV}\\ \&Integral;{\mathrm{fW}}_z\mathrm{dV}=\&Integral;f_g{W}_z\mathrm{dV}+f_{g,x}\&Integral;(x-x_g)W_y\mathrm{dV}+f_{g,y}\&Integral;(y-y_g)W_z\mathrm{dV}+f_{g,z}\&Integral;(z-z_g)W_z\mathrm{dV}\end{array}\right.---\left(11\right)$

用求和符号代替方程组(11)的积分号，被积函数dV＝m_p/ρ_p，这样可以将方程组(11)写成矩阵形式的表达式：

$\left[\begin{array}{c}\underset{p=1}{\overset{N}{\mathrm{\&Sigma;}}}{f}_{p}W\frac{m_p}{{\mathrm{\&rho;}}_p}\\ \underset{p=1}{\overset{N}{\mathrm{\&Sigma;}}}{f}_p{W}_x\frac{m_p}{{\mathrm{\&rho;}}_p}\\ \underset{p=1}{\overset{N}{\mathrm{\&Sigma;}}}{f}_p{W}_y\frac{m_p}{{\mathrm{\&rho;}}_p}\\ \underset{p=1}{\overset{N}{\mathrm{\&Sigma;}}}{f}_p{W}_z\frac{m_p}{{\mathrm{\&rho;}}_p}\end{array}\right]=\left[\begin{array}{cccc}\underset{p=1}{\overset{N}{\mathrm{\&Sigma;}}}W\frac{m_p}{{\mathrm{\&rho;}}_p}& \underset{p=1}{\overset{N}{\mathrm{\&Sigma;}}}(x_p-x_g)W\frac{m_p}{{\mathrm{\&rho;}}_p}& \underset{p=1}{\overset{N}{\mathrm{\&Sigma;}}}(y_p-y_g)W\frac{m_p}{{\mathrm{\&rho;}}_p}& \underset{p=1}{\overset{N}{\mathrm{\&Sigma;}}}(z_p-z_g)W\frac{m_p}{{\mathrm{\&rho;}}_p}\\ \underset{p=1}{\overset{N}{\mathrm{\&Sigma;}}}{W}_x\frac{m_p}{{\mathrm{\&rho;}}_p}& \underset{p=1}{\overset{N}{\mathrm{\&Sigma;}}}(x_p-x_g)W_x\frac{m_p}{{\mathrm{\&rho;}}_p}& \underset{p=1}{\overset{N}{\mathrm{\&Sigma;}}}(y_p-y_g)W_x\frac{m_p}{{\mathrm{\&rho;}}_p}& \underset{p=1}{\overset{N}{\mathrm{\&Sigma;}}}(z_p-z_g)W_x\frac{m_p}{{\mathrm{\&rho;}}_p}\\ \underset{p=1}{\overset{N}{\mathrm{\&Sigma;}}}{W}_y\frac{m_p}{{\mathrm{\&rho;}}_p}& \underset{p=1}{\overset{N}{\mathrm{\&Sigma;}}}(x_p-x_g)W_y\frac{m_p}{{\mathrm{\&rho;}}_p}& \underset{p=1}{\overset{N}{\mathrm{\&Sigma;}}}(y_p-y_g)W_y\frac{m_p}{{\mathrm{\&rho;}}_p}& \underset{p=1}{\overset{N}{\mathrm{\&Sigma;}}}(z_p-z_g)W_y\frac{m_p}{{\mathrm{\&rho;}}_p}\\ \underset{p=1}{\overset{N}{\mathrm{\&Sigma;}}}{W}_z\frac{m_p}{{\mathrm{\&rho;}}_p}& \underset{p=1}{\overset{N}{\mathrm{\&Sigma;}}}(x_p-x_g)W_z\frac{m_p}{{\mathrm{\&rho;}}_p}& \underset{p=1}{\overset{N}{\mathrm{\&Sigma;}}}(y_p-y_g)W_z\frac{m_p}{{\mathrm{\&rho;}}_p}& \underset{p=1}{\overset{N}{\mathrm{\&Sigma;}}}(z_p-z_g)W_z\frac{m_p}{{\mathrm{\&rho;}}_p}\end{array}\right]\left[\begin{array}{c}{f}_g\\ f_{g,x}\\ f_{g,y}\\ f_{g,z}\end{array}\right]---\left(12\right)$

用A来表示公式(12)右面的系数矩阵，用A^-1表示其逆矩阵，对公式(12)进行移向后表示为：

$\left[\begin{array}{c}{f}_g\\ f_{g,x}\\ f_{g,y}\\ f_{g,z}\end{array}\right]=A^{-1}\left[\begin{array}{c}\underset{p=1}{\overset{N}{\mathrm{\&Sigma;}}}{f}_{p}W\frac{m_p}{{\mathrm{\&rho;}}_p}\\ \underset{p=1}{\overset{N}{\mathrm{\&Sigma;}}}{f}_p{W}_x\frac{m_p}{{\mathrm{\&rho;}}_p}\\ \underset{p=1}{\overset{N}{\mathrm{\&Sigma;}}}{f}_p{W}_y\frac{m_p}{{\mathrm{\&rho;}}_p}\\ \underset{p=1}{\overset{N}{\mathrm{\&Sigma;}}}{f}_p{W}_z\frac{m_p}{{\mathrm{\&rho;}}_p}\end{array}\right]---\left(13\right)$

其中A^-1的元素表示为：

$A^{-1}=\left[\begin{array}{cccc}{A}_0& B_0& C_0& D_0\\ A_1& B_1& C_1& D_1\\ A_2& B_2& C_2& D_2\\ A_3& B_3& C_3& D_3\end{array}\right]---\left(14\right)$

利用(13)式和(14)式我们最终可以求出f_g的表达式为：

$f_g=A_0\underset{p=1}{\overset{N}{\mathrm{\&Sigma;}}}{f}_{p}W\frac{m_p}{{\mathrm{\&rho;}}_p}+B_0\underset{p=1}{\overset{N}{\mathrm{\&Sigma;}}}{f}_p{W}_x\frac{m_p}{{\mathrm{\&rho;}}_p}+C_0\underset{p=1}{\overset{N}{\mathrm{\&Sigma;}}}{f}_p{W}_y\frac{m_p}{{\mathrm{\&rho;}}_p}+D_0\underset{p=1}{\overset{N}{\mathrm{\&Sigma;}}}{f}_p{W}_z\frac{m_p}{{\mathrm{\&rho;}}_p}---\left(15\right)$

证明到此结束。在推导过程中只有公式(9)是近似表达式，其它公式的推导过程没有采用任何近似，因此公式(15)与公式(9)拥有相同的精度，根据公式(9)的精度可知本方法具有一阶精度。

(2)可以对任意无网格方法的数值结果进行后处理；

(3)可以对任意复杂外形的数值结果进行后处理；

(4)显著的提高了计算结果的可视化特性。

附图说明

图1是一种通用的无网格数值结果后处理方法流程图；

图2是网格节点x_g的影响域；

图3是本发明实施例1由无网格方法得到的腔内粒子分布图；

图4是本发明实施例1划分的计算域网格；

图5是本发明实施例1后处理之前的压力云图；

图6是本发明实施例1后处理之前的x方向速度云图；

图7是本发明实施例1后处理之前的y方向速度云图；

图8是本发明实施例1后处理之后的压力云图；

图9是本发明实施例1后处理之后的x方向速度云图；

图10是本发明实施例1后处理之后的y方向速度云图；

图11是本发明实施例2由无网格方法得到的腔内粒子分布图；

图12是本发明实施例2划分的计算域网格；

图13是本发明实施例2后处理之前的温度云图；

图14是本发明实施例2后处理之前的x方向速度云图；

图15是本发明实施例2后处理之前的y方向速度云图；

图16是本发明实施例2后处理之后的温度云图；

图17是本发明实施例2后处理之后的x方向速度云图；

图18是本发明实施例2后处理之后的y方向速度云图。

具体实施方式

下面结合附图和实施例对本发明作详细说明。

实施例1下面以二维腔内剪切流为例来说明本发明的技术方案的过程，如图1所示。

步骤1

以二维腔内剪切流为例，方腔的边长L＝0.001m，初始时刻粒子间距dx＝dy＝0.00025m，光滑长度h＝λdx，这里λ＝1.3，粒子的密度为标况下空气的密度ρ_p＝1.27kg/m³，粒子的质量m_p＝ρdxdy，上壁面粒子的速度u_p＝0.001m/s保持不变，左右壁面和下壁面的速度为0。通过任意无网格方法(比如SPH、ISPH、FPM、FVPM)对腔内剪切流进行数值模拟，这里采用SPH方法获得腔内及其边界上粒子的所有物理量，其中包括密度ρ_p、质量m_p、x方向速度u_p、y方向速度v_p、压力p_p。图3给出了用SPH方法获得的腔内粒子的分布图。

步骤2

该问题的外形较为简单，这里通过编写程序方法划分网格，保存网格节点的坐标。图4是计算域的网格。

步骤3

计算每个网格节点与所有无网格粒子之间的距离r，当r<r_a时，则认为无网格粒子x_p在网格节点x_g的影响域内，网格节点的影响域半径r_a＝kh，这里k＝2，h＝λdx。

步骤4

通过公式(1)计算核函数W的值，通过(2)、(3)计算核函数分别在x、y方向的导数值，式(1)中的α_d取值为15/7h²。

步骤5

对于二维问题，系数矩阵A变为：

$A=\left[\begin{array}{ccc}{a}_0& b_0& c_0\\ a_1& b_1& c_1\\ a_2& b_2& c_2\end{array}\right]$

利用公式(6)求出系数矩阵中所有的元素，对于二维问题可以利用直接法求出系数矩阵的逆矩阵A^-1为：

$A^{-1}=\frac{\left[\begin{array}{ccc}{b}_1{c}_2-b_2{c}_1& c_0{b}_2-b_0{c}_2& b_0{c}_1-b_1{c}_0\\ a_2{c}_1-a_1{c}_2& a_0{c}_2-a_2{c}_0& a_1{c}_0-a_0{c}_1\\ a_1{b}_2-a_2{b}_1& a_2{b}_0-a_0{b}_2& a_0{b}_1-a_1{b}_0\end{array}\right]\begin{array}{}\end{array}}{a_0(b_1{c}_2-b_2{c}_1)-b_0(a_1{c}_2-a_2{c}_1)+c_0(a_1{b}_2-a_2{b}_1)}$

步骤6

利用公式(8)计算网格节点上的密度ρ_g、x方向速度u_g、y方向速度v_g、压力p_g。

步骤7

将网格节点上的物理量导入绘图软件，展示结果。图5、6、7分别是后处理之前的压力云图、x方向速度云图、y方向速度云图；图8、9、10分别是利用本专利提出的后处理方法得到的压力云图、x方向速度云图、y方向速度云图。将处理之前与处理之后的云图比较可以发现在处理之前方腔角落的云图有锯齿状，但是在处理之后这些锯齿状消失，云图变的更加光滑和漂亮，而且结果更加符合物理实际情况。

实施例2下面以自然对流为例来说明本发明的技术方案。

步骤1

以自然对流为例，方腔的边长L＝0.1m，初始时刻粒子间距dx＝dy＝0.002m，光滑长度h＝λdx，这里λ＝1.3，粒子的密度为标况下空气的密度ρ_p＝1.27kg/m³，粒子的质量m_p＝ρdxdy，所有壁面粒子的速度为0，上下壁面采用绝热边界条件，左右壁面温度恒定不变，左壁面为高温壁面，右壁面为低温壁面，左右壁面温差ΔT＝10k。通过任意无网格方法(比如SPH、ISPH、FPM、FVPM)对自然对流进行数值模拟，这里采用SPH方法获得腔内及其边界上粒子的所有物理量，其中包括密度ρ_p、质量m_p、x方向速度u_p、y方向速度v_p、温度T_p。如图11给出了用SPH方法获得的腔内粒子的分布图。

步骤2

该问题的外形较为简单，这里通过编写程序方法划分网格，保存网格节点的坐标。如图12是计算域的网格。

步骤3

计算每个网格节点与所有无网格粒子之间的距离r，当r<r_a时，则认为无网格粒子x_p在网格节点x_g的影响域内。网格节点的影响域半径r_a＝kh，这里k＝2，h＝λdx。

步骤4

通过公式(1)计算核函数W的值，通过(2)、(3)计算核函数分别在x、y方向的导数值，式一中的α_d取为15/7h²。

步骤5

对于二维问题，系数矩阵A变为：

$A=\left[\begin{array}{ccc}{a}_0& b_0& c_0\\ a_1& b_1& c_1\\ a_2& b_2& c_2\end{array}\right]$

利用公式(6)求出系数矩阵中所有的元素，对于二维问题可以利用直接法求出系数矩阵的逆矩阵A^-1为：

$A^{-1}=\frac{\left[\begin{array}{ccc}{b}_1{c}_2-b_2{c}_1& c_0{b}_2-b_0{c}_2& b_0{c}_1-b_1{c}_0\\ a_2{c}_1-a_1{c}_2& a_0{c}_2-a_2{c}_0& a_1{c}_0-a_0{c}_1\\ a_1{b}_2-a_2{b}_1& a_2{b}_0-a_0{b}_2& a_0{b}_1-a_1{b}_0\end{array}\right]\begin{array}{}\end{array}}{a_0(b_1{c}_2-b_2{c}_1)-b_0(a_1{c}_2-a_2{c}_1)+c_0(a_1{b}_2-a_2{b}_1)}$

步骤6

利用公式(8)计算网格节点上的密度ρ_g、x方向速度u_g、y方向速度v_g、温度T_g。

步骤7

将网格节点上的物理量导入绘图软件，展示结果。图13、14、15分别是后处理之前的温度云图、x方向速度云图、y方向速度云图；图16、17、18分别是利用本专利提出的后处理方法得到的温度云图、x方向速度云图、y方向速度云图。将处理之前与处理之后的云图比较可以发现在处理之前，整个方腔的云图呈现出团状，云图连续性非常差，但是在处理之后方腔的云图变的光滑，这有助于对数值结果的认识和分析，由此可见本发明方法是有效的。

以上所述的具体描述，对发明的目的、技术方案和有益效果进行了进一步详细说明，所应理解的是，以上所述仅为本发明的具体实施例，用于解释本发明，并不用于限定本发明的保护范围，凡在本发明的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (1)

1.一种通用的无网格数值结果后处理方法，其特征在于，具体实现步骤如下：

步骤1、首先将无网格方法得到的粒子物理信息储存起来；

步骤2、对所研究问题的计算域划分网格，将网格节点信息储存起来；

步骤3、利用粒子搜索方法，确定每个网格节点影响域中的无网格粒子；

步骤4、计算网格节点影响域内的核函数及其导数；

步骤5、计算系数矩阵的逆矩阵，其中系数矩阵A用下式表示：

$A=\left[\begin{array}{cccc}{a}_0& b_0& c_0& d_0\\ a_1& b_1& c_1& d_1\\ a_2& b_2& c_2& d_2\\ a_3& b_3& c_3& d_3\end{array}\right]$

矩阵中各元素值通过下式计算：

$\begin{array}{cccc}{a}_0=\underset{p=1}{\overset{N}{\mathrm{\&Sigma;}}}W\frac{m_p}{{\mathrm{\&rho;}}_p}& b_0=\underset{p=1}{\overset{N}{\mathrm{\&Sigma;}}}(x_p-x_g)W\frac{m_p}{{\mathrm{\&rho;}}_p}& c_0=\underset{p=1}{\overset{N}{\mathrm{\&Sigma;}}}(y_p-y_g)W\frac{m_p}{{\mathrm{\&rho;}}_p}& d_0=\underset{p=1}{\overset{N}{\mathrm{\&Sigma;}}}(z_p-z_g)W\frac{m_p}{{\mathrm{\&rho;}}_p}\\ a_1=\underset{p=1}{\overset{N}{\mathrm{\&Sigma;}}}{W}_x\frac{m_p}{{\mathrm{\&rho;}}_p}& b_1=\underset{p=1}{\overset{N}{\mathrm{\&Sigma;}}}(x_p-x_g)W_x\frac{m_p}{{\mathrm{\&rho;}}_p}& c_1=\underset{p=1}{\overset{N}{\mathrm{\&Sigma;}}}(y_p-y_g)W_x\frac{m_p}{{\mathrm{\&rho;}}_p}& d_1=\underset{p=1}{\overset{N}{\mathrm{\&Sigma;}}}(z_p-z_g)W_x\frac{m_p}{{\mathrm{\&rho;}}_p}\\ a_2=\underset{p=1}{\overset{N}{\mathrm{\&Sigma;}}}{W}_y\frac{m_p}{{\mathrm{\&rho;}}_p}& b_2=\underset{p=1}{\overset{N}{\mathrm{\&Sigma;}}}(x_p-x_g)W_y\frac{m_p}{{\mathrm{\&rho;}}_p}& c_2=\underset{p=1}{\overset{N}{\mathrm{\&Sigma;}}}(y_p-y_g)W_y\frac{m_p}{{\mathrm{\&rho;}}_p}& d_2=\underset{p=1}{\overset{N}{\mathrm{\&Sigma;}}}(z_p-z_g)W_y\frac{m_p}{{\mathrm{\&rho;}}_p}\\ a_3=\underset{p=1}{\overset{N}{\mathrm{\&Sigma;}}}{W}_z\frac{m_p}{{\mathrm{\&rho;}}_p}& b_3=\underset{p=1}{\overset{N}{\mathrm{\&Sigma;}}}(x_p-x_g)W_z\frac{m_p}{{\mathrm{\&rho;}}_p}& c_3=\underset{p=1}{\overset{N}{\mathrm{\&Sigma;}}}(y_p-y_g)W_z\frac{m_p}{{\mathrm{\&rho;}}_p}& d_3=\underset{p=1}{\overset{N}{\mathrm{\&Sigma;}}}(z_p-z_g)W_z\frac{m_p}{{\mathrm{\&rho;}}_p}\end{array}$

其中，W、W_x、W_y、W_z分别是核函数及其在x、y、z方向的导数，N是网格节点影响域内无网格粒子的个数，m_p和ρ_p分别表示无网格粒子的质量和密度，x_g、y_g和z_g分别表示网格节点在x、y、z方向的坐标，x_p、y_p和z_p分别表示无网格粒子在x、y、z方向的坐标；A的逆矩阵A^-1可以表示成下式：

$A^{-1}=\left[\begin{array}{cccc}{A}_0& B_0& C_0& D_0\\ A_1& B_1& C_1& D_1\\ A_2& B_2& C_2& D_2\\ A_3& B_3& C_3& D_3\end{array}\right];$

步骤6、通过下式计算网格节点上的物理量；

$f_g=A_0\underset{p=1}{\overset{N}{\mathrm{\&Sigma;}}}{f}_{p}W\frac{m_p}{{\mathrm{\&rho;}}_p}+B_0\underset{p=1}{\overset{N}{\mathrm{\&Sigma;}}}{f}_p{W}_x\frac{m_p}{{\mathrm{\&rho;}}_p}+C_0\underset{p=1}{\overset{N}{\mathrm{\&Sigma;}}}{f}_p{W}_y\frac{m_p}{{\mathrm{\&rho;}}_p}+D_0\underset{p=1}{\overset{N}{\mathrm{\&Sigma;}}}{f}_p{W}_z\frac{m_p}{{\mathrm{\&rho;}}_p}$

其中，f_g、f_p分别表示网格节点、无网格粒子的任意物理量；

步骤7、将网格节点上的物理量导入绘图软件，展示结果。