一种改进的KGF-SPH方法

摘要

本发明涉及一种改进的FPM方法，属于计算力学技术领域，该方法包括以下步骤：首先在所研究问题的计算域中布置粒子，其次根据所研究问题的初始条件和边界条件对粒子的物理属性初始化，接下来利用改进的KGF-SPH方法的离散格式对所研究问题的控制方程进行近似，最后对所研究问题的控制方程进行迭代，得出数值模拟结果。对比KGF-SPH方法，本发明方法保持了KGF-SPH方法无核梯度的特性，具有较高的数值精度和稳定性，便于边界条件的处理。

说明书

技术领域

本发明涉及一种改进的KGF-SPH方法，属于计算流体力学技术领域。

背景技术

SPH(SmoothedParticleHydrodynamics)是一种拉格朗日型无网格粒子方法，由Lucy，Gingold与Monaghan于1977年提出并用于求解三维开放空间的天体物理学问题，现在已经被广泛用于多相流、重力流、热传导等流体力学问题的数值模拟中。在SPH方法中，用一系列粒子来代替连续介质，这些粒子可以在计算域内自由移动，因此SPH方法可以很好的处理传统有网格方法难以解决的问题，如具有移动边界、自由表面或边界极大变形的流动问题。虽然SPH方法已经在流体力学和固体力学领域有了非常广泛的应用，但是SPH方法在对流问题的应用仍较少。

对流现象无处不在，如日常生活中的热水沸腾，自然界中的海洋环流，工程应用中的核反应堆冷却等，因此研究对流问题具有非常的意义。但是传统SPH方法存在粒子不一致性、应力不稳定性等缺点，无法很好的模拟对流问题。在SPH方法的基础上，Huang等人在2015年提出了一种KGF-SPH方法(kernelgradientfree-smoothedparticlehydrodynamics)。该方法在核近似和粒子近似过程中只用到了核函数本身，不包含核函数的导数，因此无论核函数是否可导都可用于该方法，是一种无核梯度的SPH方法。KGF-SPH方法降低了无网格方法对于核函数的要求，具有较高的精度。但是KGF-SPH方法中二阶导数通过对一阶导数再求导来求解，这样虽然提高了计算精度，但是存在以下缺点：

(1)计算量较大。因为二阶导数的求解需要在对一阶导数求导，增加了一倍的计算量。

(2)不利于处理边界条件。因为在求解边界附近流体粒子的物理量时，需要先求解出边界上粒子物理量的导数值，然而为了精确的求出边界上粒子物理量的导数值，还需要在边界外布置更多的虚粒子。

发明内容

本发明的目的是为了提高数值模拟的精度和稳定性，减小计算量，提出了一种改进的KGF-SPH方法。

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

一种改进的KGF-SPH方法，具体实现步骤如下：

步骤1、在所研究问题的计算域中布置粒子；

步骤2、根据所研究问题的初始条件和边界条件对粒子的物理属性进行初始化；

步骤3、利用下述空间离散格式对所研究问题的控制方程进行近似，得到函数的时间导数的粒子近似式，其近似原理如下：

对任意函数f(r_j)在点r＝r_i处进行泰勒展开，并且只保留到一阶导数，可得：

f(r_j)＝f(r_i)+▽f(r_i)·(r_j-r_i)+o(h²)(1)

忽略式(1)中的高阶项o(h²)，在式(1)两边同时乘以核函数W(R,h)和(r_j-r_i)W(R,h)并在计算域上进行积分得：

∫f(r_j)WdV＝f(r_i)∫WdV+▽f(r_i)·∫(r_j-r_i)WdV(2)

∫f(r_j)(r_j-r_i)WdV＝f(r_i)∫(r_j-r_i)WdV+▽f(r_i)·∫(r_j-r_i)²WdV(3)

联立式(2)和式(3)，并求解方程组，可得：

$>\left(\begin{array}{c}f\left(r_i\right)\\ ▿f\left(r_i\right)\end{array}\right)=A^{*}\left(\begin{array}{c}\int f\left(r_j\right)WdV\\ \int f\left(r_j\right)(r_j-r_i)WdV\end{array}\right)---\left(4\right)$>

上式对应的粒子近似式为：

$>\left(\begin{array}{c}{f}_i\\ {\left(▿f\right)}_i\end{array}\right)=A^{*}\left(\begin{array}{c}\underset{j=1}{\overset{N}{\Sigma }}{f}_{j}W\frac{m_j}{{\rho }_j}\\ \underset{j=1}{\overset{N}{\Sigma }}{f}_j(r_j-r_i)W\frac{m_j}{{\rho }_j}\end{array}\right)---\left(5\right)$>

式中，f_i＝f(r_i)，f_j＝f(r_j)，(▽f)_i＝▽f(r_i)，N表示在i粒子的支持域内的j粒子(即i粒子的最近相邻粒子)总量，m_j和ρ_j分别表示j粒子的质量和密度。

在一维空间下，式(5)的具体表达式为：

$>\left(\begin{array}{c}{f}_i\\ f_{i,x}\end{array}\right)=A^{*}\left(\begin{array}{c}\underset{j=1}{\overset{N}{\Sigma }}{f}_{j}W\frac{m_j}{{\rho }_j}\\ \underset{j=1}{\overset{N}{\Sigma }}{f}_j{x}_{ji}W\frac{m_j}{{\rho }_j}\end{array}\right)---\left(6\right)$>

其中，x_ji＝x_j-x_i，f_i,x为粒子i在x方向上的一阶导数。修正矩阵A^*在一维空间下的表达式为：

$>A^{*}=\left(\begin{array}{cc}{A}_0& B_0\\ A_1& B_1\end{array}\right)={\left(\begin{array}{cc}\underset{j=1}{\overset{N}{\Sigma }}W\frac{m_j}{{\rho }_j}& \underset{j=1}{\overset{N}{\Sigma }}{x}_{ji}W\frac{m_j}{{\rho }_j}\\ \underset{j=1}{\overset{N}{\Sigma }}{x}_{ji}W\frac{m_j}{{\rho }_j}& \underset{j=1}{\overset{N}{\Sigma }}{x}_{ji}^{2}W\frac{m_j}{{\rho }_j}\end{array}\right)}^{-1}---\left(7\right)$>

在一维空间下，对函数f(x_j)在点x_j＝x_i进行二阶泰勒展开并略去高阶项可得：

$>f\left(x_j\right)=f_i+f_{i,x}{x}_{ji}+\frac{1}{2}{f}_{i,xx}{x}_{ji}^2---\left(8\right)$>

式中，f_i,xx为粒子i在x方向上的二阶导数。式(8)两边同时乘以核函数W并移项可得：

$>(f_j-f_i)W=(f_{i,x}{x}_{ji}+\frac{1}{2}{f}_{i,xx}{x}_{ji}^2)W---\left(9\right)$>

式(9)两边同时积分可得：

$>\int (f_j-f_i)WdV=\int f_{i,x}{x}_{ji}WdV+\int \frac{1}{2}{f}_{i,xx}{x}_{ji}^{2}WdV---\left(10\right)$>

由核函数的对称特性可知，上式右边第一项为0，因此式(10)可化简为：

$>\int (f_j-f_i)WdV=\int \frac{1}{2}{f}_{i,xx}{x}_{ji}^{2}WdV---\left(11\right)$>

对(11)式进行移项可以推导出一维空间下拉普拉斯算子的离散格式为：

$>f_{i,xx}=\frac{2\int (f_j-f_i)WdV}{\int x_{ji}^{2}WdV}---\left(12\right)$>

式(12)对应的粒子近似式为：

$>f_{i,xx}=\frac{2\underset{j=1}{\overset{N}{\Sigma }}(f_j-f_i)W\frac{m_j}{{\rho }_j}}{\underset{j=1}{\overset{N}{\Sigma }}{x}_{ji}^{2}W\frac{m_j}{{\rho }_j}}---\left(13\right)$>

在二维空间下，式(5)的具体表达式为：

$>\left(\begin{array}{c}{f}_i\\ f_{i,x}\\ f_{i,y}\end{array}\right)=A^{*}\left(\begin{array}{c}\underset{j=1}{\overset{N}{\Sigma }}{f}_{j}W\frac{m_j}{{\rho }_j}\\ \underset{j=1}{\overset{N}{\Sigma }}{f}_j{x}_{ji}W\frac{m_j}{{\rho }_j}\\ \underset{j=1}{\overset{N}{\Sigma }}{f}_j{y}_{ji}W\frac{m_j}{{\rho }_j}\end{array}\right)---\left(14\right)$>

其中，x_ji＝x_j-x_i，y_ji＝y_j-y_i，f_i,x和f_i,y分别为粒子i在x和y方向上的一阶导数。修正矩阵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)={\left(\begin{array}{ccc}\underset{j=1}{\overset{N}{\Sigma }}W\frac{m_j}{{\rho }_j}& \underset{j=1}{\overset{N}{\Sigma }}{x}_{ji}W\frac{m_j}{{\rho }_j}& \underset{j=1}{\overset{N}{\Sigma }}{y}_{ji}W\frac{m_j}{{\rho }_j}\\ \underset{j=1}{\overset{N}{\Sigma }}{x}_{ji}W\frac{m_j}{{\rho }_j}& \underset{j=1}{\overset{N}{\Sigma }}{x}_{ji}^{2}W\frac{m_j}{{\rho }_j}& \underset{j=1}{\overset{N}{\Sigma }}{y}_{ji}{x}_{ji}W\frac{m_j}{{\rho }_j}\\ \underset{j=1}{\overset{N}{\Sigma }}{y}_{ji}W\frac{m_j}{{\rho }_j}& \underset{j=1}{\overset{N}{\Sigma }}{x}_{ji}{y}_{ji}W\frac{m_j}{{\rho }_j}& \underset{j=1}{\overset{N}{\Sigma }}{y}_{ji}^{2}W\frac{m_j}{{\rho }_j}\end{array}\right)}^{-1}---\left(15\right)$>

在二维空间下，对函数f(x_j,y_j)在点(x_i,y_i)进行二阶泰勒展开并略去高阶项可得：

$>f(x_j,y_j)=f_i+f_{i,x}{x}_{ji}+f_{i,y}{y}_{ji}+\frac{1}{2}{f}_{i,xx}{x}_{ji}^2+f_{i,xy}{x}_{ji}{y}_{ji}+\frac{1}{2}{f}_{i,yy}{y}_{ji}^2---\left(16\right)$>

式中，f_i,xy为粒子i的混合导数，f_i,xx和f_i,yy分别为粒子i在x和y方向上的二阶导数。式(16)两边同时乘以核函数W并移项可得：

$>(f_j-f_i)W=(f_{i,x}{x}_{ji}+f_{i,y}{y}_{ji}+\frac{1}{2}{f}_{i,xx}{x}_{ji}^2+f_{i,xy}{x}_{ji}{y}_{ji}+\frac{1}{2}{f}_{i,yy}{y}_{ji}^2)W---\left(17\right)$>

式(17)两边同时积分得：

$>\begin{array}{c}\int (f_j-f_i)WdV=\int (f_{i,x}{x}_{ji}+f_{i,y}{y}_{ji})MdV+\int f_{i,xy}{x}_{ji}{y}_{ji}WdV\\ +\int (\frac{1}{2}{f}_{i,xx}{x}_{ji}^2+\frac{1}{2}{f}_{i,yy}{y}_{ji}^2)WdV\end{array}---\left(18\right)$>

由核函数的对称特性可知，式(18)右边第一项和第二项为0，因此上式可化简为：

$>\int (f_j-f_i)WdV=\int (\frac{1}{2}{f}_{i,xx}{x}_{ji}^2+\frac{1}{2}{f}_{i,yy}{y}_{ji}^2)WdV---\left(19\right)$>

当粒子分布均匀且规则时将在极坐标下进行求解可得：

$>\int x_{ji}^{2}WdV=\int {\left(cos\theta \right)}^{2}d\theta \int r_{ji}^{3}Wdr=\pi \int r_{ji}^{3}Wdr---\left(20\right)$>

$>\int y_{ji}^{2}WdV=\int {\left(sin\theta \right)}^{2}d\theta \int r_{ji}^{3}Wdr=\pi \int r_{ji}^{3}Wdr---\left(21\right)$>

故有：

$>\int x_{ji}^{2}WdV=\int y_{ji}^{2}WdV---\left(22\right)$>

式(19)移项并利用式(22)对其进行化简，可以推导出二维空间下拉普拉斯算子的离散格式为：

$>f_{i,xx}+f_{i,yy}=\frac{2\int (f_j-f_i)WdV}{\int x_{ji}^{2}WdV}---\left(23\right)$>

式(23)对应的粒子近似式为：

$>f_{i,xx}+f_{i,yy}=\frac{2\underset{j=1}{\overset{N}{\Sigma }}(f_j-f_i)W\frac{m_j}{{\rho }_j}}{\underset{j=1}{\overset{N}{\Sigma }}{x}_{ji}^{2}W\frac{m_j}{{\rho }_j}}---\left(24\right)$>

为了减弱粒子分布不均匀的影响，采用下式代替式(24)：

$>f_{i,xx}+f_{i,yy}=\frac{2\underset{j=1}{\overset{N}{\Sigma }}(f_j-f_i)W\frac{m_j}{{\rho }_j}}{\frac{1}{2}(\underset{j=1}{\overset{N}{\Sigma }}{x}_{ji}^{2}W\frac{m_j}{{\rho }_j}+\underset{j=1}{\overset{N}{\Sigma }}{y}_{ji}^{2}W\frac{m_j}{{\rho }_j})}---\left(25\right)$>

同理可推导出三维空间下任意函数及其一阶导数和拉普拉斯算子的近似方案：

$>\left(\begin{array}{c}{f}_i\\ f_{i,x}\\ f_{i,y}\\ f_{i,z}\end{array}\right)=A^{*}\left(\begin{array}{c}\underset{j=1}{\overset{N}{\Sigma }}{f}_{j}W\frac{m_j}{{\rho }_j}\\ \underset{j=1}{\overset{N}{\Sigma }}{f}_j{x}_{ji}W\frac{m_j}{{\rho }_j}\\ \underset{j=1}{\overset{N}{\Sigma }}{f}_j{y}_{ji}W\frac{m_j}{{\rho }_j}\\ \underset{j=1}{\overset{N}{\Sigma }}{f}_j{z}_{ji}W\frac{m_j}{{\rho }_j}\end{array}\right)---\left(26\right)$>

$>f_{i,xx}+f_{i,yy}+f_{i,zz}=\frac{2\underset{j=1}{\overset{N}{\Sigma }}(f_j-f_i)W\frac{m_j}{{\rho }_j}}{\frac{1}{3}(\underset{j=1}{\overset{N}{\Sigma }}{x}_{ji}^{2}W\frac{m_j}{{\rho }_j}+\underset{j=1}{\overset{N}{\Sigma }}{y}_{ji}^{2}W\frac{m_j}{{\rho }_j}+\underset{j=1}{\overset{N}{\Sigma }}{z}_{ji}^{2}W\frac{m_j}{{\rho }_j})}---\left(27\right)$>

其中，x_ji＝x_j-x_i，y_ji＝y_j-y_i，z_ji＝z_j-z_i，f_i,x、f_i,y和f_i,z分别为粒子i在x、y和z方向上的一阶导数，f_i,xx、f_i,yy和f_i,zz分别为粒子i在x、y和z方向上的二阶导数，x_k、y_k和z_k分别为粒子k在x、y和z方向上的坐标，k＝i或j。修正矩阵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)={\left(\begin{array}{cccc}\underset{j=1}{\overset{N}{\Sigma }}W\frac{m_j}{{\rho }_j}& \underset{j=1}{\overset{N}{\Sigma }}{x}_{ji}W\frac{m_j}{{\rho }_j}& \underset{j=1}{\overset{N}{\Sigma }}{y}_{ji}W\frac{m_j}{{\rho }_j}& \underset{j=1}{\overset{N}{\Sigma }}{z}_{ji}W\frac{m_j}{{\rho }_j}\\ \underset{j=1}{\overset{N}{\Sigma }}{x}_{ji}W\frac{m_j}{{\rho }_j}& \underset{j=1}{\overset{N}{\Sigma }}{x}_{ji}^{2}W\frac{m_j}{{\rho }_j}& \underset{j=1}{\overset{N}{\Sigma }}{y}_{ji}{x}_{ji}W\frac{m_j}{{\rho }_j}& \underset{j=1}{\overset{N}{\Sigma }}{z}_{ji}{x}_{ji}W\frac{m_j}{{\rho }_j}\\ \underset{j=1}{\overset{N}{\Sigma }}{y}_{ji}W\frac{m_j}{{\rho }_j}& \underset{j=1}{\overset{N}{\Sigma }}{x}_{ji}{y}_{ji}W\frac{m_j}{{\rho }_j}& \underset{j=1}{\overset{N}{\Sigma }}{y}_{ji}^{2}W\frac{m_j}{{\rho }_j}& \underset{j=1}{\overset{N}{\Sigma }}{z}_{ji}{y}_{ji}W\frac{m_j}{{\rho }_j}\\ \underset{j=1}{\overset{N}{\Sigma }}{z}_{ji}W\frac{m_j}{{\rho }_j}& \underset{j=1}{\overset{N}{\Sigma }}{x}_{ji}{z}_{ji}W\frac{m_j}{{\rho }_j}& \underset{j=1}{\overset{N}{\Sigma }}{y}_{ji}{z}_{ji}W\frac{m_j}{{\rho }_j}& \underset{j=1}{\overset{N}{\Sigma }}{z}_{ji}^{2}W\frac{m_j}{{\rho }_j}\end{array}\right)}^{-1}---\left(28\right)$>

式(6)、式(7)和式(13)给出了一维空间下任意函数及其一阶导数和拉普拉斯算子的近似方案，式(14)、式(15)和式(25)给出了二维空间下任意函数及其一阶导数和拉普拉斯算子的近似方案，式(26)、式(27)和式(28)给出了三维空间下任意函数及其一阶导数和拉普拉斯算子的近似方案。该方案具有较高的数值精度和稳定性，便于边界条件的处理。

步骤4、根据所研究问题的控制方程进行时间迭代，迭代一定步数后，得到模拟结果，其中迭代公式为：

$>\left(\begin{array}{c}t=t+dt\\ f_i(t+dt)=f_i\left(t\right)+dt·{\left(\frac{df\left(t\right)}{dt}\right)}_i\end{array}\right)---\left(29\right)$>

式中，dt为迭代的时间步长。

有益效果

本发明对比已有技术具有如下优点：

(1)改进的KGF-SPH方法计算量较小；

(2)改进的KGF-SPH方法提高了数值模拟的精度和稳定性；

(3)改进的KGF-SPH方法便于边界条件的处理。

附图说明

图1为本发明实施例1初始粒子分布示意图；

图2为本发明实施例1基于(a)KGF-SPH方法、(b)改进的KGF-SPH方法以及(c)解析解得到的一维热传导问题温度随时间变化的分布图；

图3为本发明实施例1基于KGF-SPH方法和改进的KGF-SPH方法模拟一维热传导问题得到的t＝0.2s时的温度分布曲线与解析解的对比；

图4为本发明实施例2初始粒子分布示意图；

图5为本发明实施例2基于(a)SPH方法、(b)改进的KGF-SPH方法和(c)FVM方法模拟封闭方腔自然对流得到的温度分布云图对比；

图6为本发明实施例2基于(a)SPH方法、(b)改进的KGF-SPH方法和(c)FVM方法模拟封闭方腔自然对流得到的水平方向速度分布云图对比；

图7为本发明实施例2基于(a)SPH方法、(b)改进的KGF-SPH方法和(c)FVM方法模拟封闭方腔自然对流得到的垂直方向速度分布云图对比。

具体实施方式

为了使本发明的目的、技术方案和优点更加清晰明白，下面结合实施例和附图，对本发明实施例做进一步的详细说明。

实施例1下面以一维热传导问题来说明本发明的技术方案。

步骤1、在所研究问题的计算域中布置粒子；

以一维热传导问题为例，计算域为从0到1的线段L，粒子分布方式如图1所示，粒子间距为Δ，图中实心圆点为实粒子，空心圆点为虚粒子。

步骤2、根据所研究问题的初始条件和边界条件对粒子的物理属性进行初始化；

粒子的密度ρ＝1，粒子的质量m＝ρΔ，T表示粒子的温度。初始条件和边界条件分别为：T(x,0)＝sinπx和T(0,t)＝T(1,t)＝0。

步骤3、利用改进的KGF-SPH方法的离散格式对一维热传导问题的控制方程进行近似；

本实施例中一维热传导问题的控制方程为：

$>\frac{dT}{dt}=T_{xx}---\left(30\right)$>

式中，T_xx为温度在x方向上的二阶导数。式(30)对应的解析解为：

$>T(x,t)=e^{-{\pi }^{2}t}sin\pi x---\left(31\right)$>

利用式(13)对一维热传导问题的控制方程式(30)进行离散可得：

$>{\left(\frac{dT}{dt}\right)}_i=T_{i,xx}=\frac{2\underset{j=1}{\overset{N}{\Sigma }}(T_j-T_i)W\frac{m_j}{{\rho }_j}}{\underset{j=1}{\overset{N}{\Sigma }}{x}_{ji}^{2}W\frac{m_j}{{\rho }_j}}---\left(32\right)$>

式中，$>W=W(R,h)={\alpha }_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& 2\le R\end{array}\right),$>h是定义核函数W的影响区域(即支持域)的光滑长度，h＝1.3Δ，R＝|r_j-r_i|/h，在一维空间下α_d＝1/h。

步骤4、根据所研究问题的控制方程进行时间迭代，时间步长dt＝2×10^-4s，迭代一定步数后，得到模拟结果。

图2给出了分别采用(a)KGF-SPH方法和(b)改进的KGF-SPH方法对一维热传导问题进行模拟所得到的温度随时间(0～1s)变化的分布图，为对比分析各数值方法的精度，图中也给出了(c)解析解的结果。从图2中可以看出，KGF-SPH方法所得结果有明显的震荡，而本发明提出的改进的KGF-SPH方法能得到比较光滑的结果，且与解析解基本一致。

图3给出了基于KGF-SPH方法和改进的KGF-SPH方法迭代1000次(t＝0.2s)后得到的温度分布曲线，并与t＝0.2s时刻的解析解进行对比。从图中可以看出，KGF-SPH方法所得结果在解析解上方震荡，而本发明提出的改进的KGF-SPH方法所得结果与解析解吻合较好，表明改进的KGF-SPH方法比KGF-SPH方法稳定性好，且具有更高的精度。

实施例2下面以封闭方腔自然对流问题来说明本发明的技术方案。

步骤1、在所研究问题的计算域中布置粒子；

以封闭方腔自然对流问题为例，计算域为边长为L的方腔，粒子分布方式如图4所示，图中实心圆点为实粒子，空心圆点为虚粒子。

步骤2、根据所研究问题的初始条件和边界条件对粒子的物理属性进行初始化；

方腔的边长为L＝0.1m，初始时刻粒子间距dx＝dy＝0.0005m，光滑长度h＝1.3dx，粒子的质量m＝ρdxdy，粒子的初始密度为标况下空气的密度ρ₀＝1.225kg/m³，粘性系数μ＝1.707×10^-5kg/(m·s)，热扩散系数α＝1.963×10^-5m²/s，热膨胀系数β_T＝4.095×10^-3K^-1，重力加速度g＝6.680m/s²。上、下壁面为绝热边界，左、右壁面为恒温壁面，左壁面的温度T_H＝283K，右壁面的温度T_C＝273K，初始时刻方腔内部的温度T₀＝(T_H+T_C)/2。本实施例中，压力通过人工状态方程p＝(ρ-ρ₀)c²进行求解，其中c为人工声速，取c＝0.8m/s。

步骤3、利用改进的KGF-SPH方法的离散格式对封闭方腔自然对流的控制方程进行近似；

封闭方腔自然对流的流体力学控制方程为：

$>\left(\begin{array}{c}\frac{d\rho }{dt}=-\rho (u_x+v_y)\\ \frac{du}{dt}=-\frac{1}{\rho }{p}_x+\frac{\mu }{\rho }(u_{xx}+u_{yy})\\ \frac{dv}{dt}=-\frac{1}{\rho }{p}_y+\frac{\mu }{\rho }(v_{xx}+v_{yy})+{\mathrm{g\beta }}_T(T-T_0)\\ \frac{dT}{dt}-\alpha (T_{xx}+T_{yy})\end{array}\right)---\left(33\right)$>

式中，ρ表示密度，u和v分别表示速度在x方向和y方向上的分量，p表示压力，g表示重力加速度，T表示温度；下标x和y分别表示变量在x和y方向上的一阶导数；下标xx和yy分别表示变量在x和y方向上的二阶导数。

利用式(14)、式(15)和式(25)对上述控制方程进行离散可得：

$>{\left(\frac{d\rho }{dt}\right)}_i=-{\rho }_i\left(\begin{array}{c}\underset{j=1}{\overset{N}{\Sigma }}{A}_1{u}_{ji}W\frac{m_j}{{\rho }_j}+\underset{j=1}{\overset{N}{\Sigma }}{B}_1{u}_{ji}{x}_{ji}W\frac{m_j}{{\rho }_j}+\underset{j=1}{\overset{N}{\Sigma }}{C}_1{u}_{ji}{v}_{ji}W\frac{m_j}{{\rho }_j}\\ +\underset{j=1}{\overset{N}{\Sigma }}{A}_2{v}_{ji}W\frac{m_j}{{\rho }_j}+\underset{j=1}{\overset{N}{\Sigma }}{B}_2{v}_{ji}{x}_{ji}W\frac{m_j}{{\rho }_j}+\underset{j=1}{\overset{N}{\Sigma }}{C}_2{u}_{ji}{y}_{ji}W\frac{m_j}{{\rho }_j}\end{array}\right)---\left(34\right)$>

$>\begin{array}{c}{\left(\frac{du}{dt}\right)}_i=-\underset{j=1}{\overset{N}{\Sigma }}{A}_1{\mathrm{Wm}}_j\frac{p_{ji}}{{\rho }_i{\rho }_j}-\underset{j=1}{\overset{N}{\Sigma }}{B}_1{x}_{ji}{\mathrm{Wm}}_j\frac{p_{ji}}{{\rho }_i{\rho }_j}\\ -\underset{j=1}{\overset{N}{\Sigma }}{C}_1{y}_{ji}{\mathrm{Wm}}_j\frac{p_{ji}}{{\rho }_i{\rho }_j}+\frac{\mu }{{\rho }_i}\frac{2\underset{j=1}{\overset{N}{\Sigma }}(u_j-u_i)W\frac{m_j}{{\rho }_j}}{\frac{1}{2}\left(\underset{j=1}{\overset{N}{\Sigma }}{x}_{ji}^2\frac{m_j}{{\rho }_j}+\underset{j=1}{\overset{N}{\Sigma }}{y}_{ji}^2\frac{m_j}{{\rho }_j}\right)}\end{array}---\left(35\right)$>

$>\begin{array}{c}{\left(\frac{du}{dt}\right)}_i=-\underset{j=1}{\overset{N}{\Sigma }}{A}_2{\mathrm{Wm}}_j\frac{p_{ji}}{{\rho }_i{\rho }_j}-\underset{j=1}{\overset{N}{\Sigma }}{B}_2{x}_{ji}{\mathrm{Wm}}_j\frac{p_{ji}}{{\rho }_i{\rho }_j}\\ -\underset{j=1}{\overset{N}{\Sigma }}{C}_2{y}_{ji}{\mathrm{Wm}}_j\frac{p_{ji}}{{\rho }_i{\rho }_j}+\frac{\mu }{{\rho }_i}\frac{2\underset{j=1}{\overset{N}{\Sigma }}(v_j-v_i)W\frac{m_j}{{\rho }_j}}{\frac{1}{2}\left(\underset{j=1}{\overset{N}{\Sigma }}{x}_{ji}^2\frac{m_j}{{\rho }_j}+\underset{j=1}{\overset{N}{\Sigma }}{y}_{ji}^2\frac{m_j}{{\rho }_j}\right)}\\ +{\mathrm{g\beta }}_T(T_i-T_0)\end{array}---\left(36\right)$>

$>{\left(\frac{dT}{dt}\right)}_i=\alpha (T_{i,xx}+T_{i,yy})=\alpha \frac{2\underset{j=1}{\overset{N}{\Sigma }}(T_j-T_i)W\frac{m_j}{{\rho }_j}}{\frac{1}{2}(\underset{j=1}{\overset{N}{\Sigma }}{x}_{ji}^{2}W\frac{m_j}{{\rho }_j}+\underset{j=1}{\overset{N}{\Sigma }}{y}_{ji}^{2}W\frac{m_j}{{\rho }_j})}---\left(37\right)$>

式中，u_ji＝u_j-u_i，v_ji＝v_j-v_i，p_ji＝p_j-p_i；A_n、B_n、C_n(n＝1、2)为式(15)中矩阵A*的元素。本实施例中用到的核函数与实施例1相同，在二维空间下α_d＝15/(7πh²)。

步骤4、根据所研究问题的控制方程进行时间迭代，迭代的时间步长为dt＝5×10^-6s，迭代一定步数后，得到模拟结果。

图5给出了分别采用(a)SPH方法以及(b)改进的KGF-SPH方法对瑞利数Ra＝10⁶条件下的封闭方腔自然对流问题进行模拟所得到的稳定时刻的温度云图，为对比分析各数值方法的精度，图中也给出了(c)FVM方法(有限体积法)的结果。从图5中可以看出，SPH方法和本发明提出的改进的KGF-SPH方法均与FVM方法所得结果吻合较好，都能很好的模拟出流场温度的分布规律。

图6和图7给出了分别采用(a)SPH方法、(b)改进的KGF-SPH方法与(c)FVM方法模拟得到的流场速度分布云图，其中，图6为水平方向速度u的模拟结果，图7为垂直方向速度v的模拟结果。从图6和图7中可以看出，SPH方法所得模拟结果在方腔中部和壁面附近不太光滑，有明显的震荡，但本发明提出的改进的KGF-SPH方法所得模拟结果非常光滑，与FVM方法数值模拟结果基本一致，证明了改进的KGF-SPH方法具有更高的精度和稳定性。

以上所述仅为本发明的具体实施例，并不用于限定本发明的保护范围，凡在本发明的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。