一种改进的kgf-sph方法
【技术领域】
[0001] 本发明涉及一种改进的KGF-SPH方法，属于计算流体力学技术领域。
【背景技术】
[0002] SPH (Smoothed Particle Hydrodynamics)是一种拉格朗日型无网格粒子方法，由 Lucy，Gingold与Monaghan于1977年提出并用于求解三维开放空间的天体物理学问题，现 在已经被广泛用于多相流、重力流、热传导等流体力学问题的数值模拟中。在SPH方法中， 用一系列粒子来代替连续介质，这些粒子可以在计算域内自由移动，因此SPH方法可以很 好的处理传统有网格方法难以解决的问题，如具有移动边界、自由表面或边界极大变形的 流动问题。虽然SPH方法已经在流体力学和固体力学领域有了非常广泛的应用，但是SPH 方法在对流问题的应用仍较少。
[0003] 对流现象无处不在，如日常生活中的热水沸腾，自然界中的海洋环流，工程应用 中的核反应堆冷却等，因此研究对流问题具有非常的意义。但是传统SPH方法存在粒子 不一致性、应力不稳定性等缺点，无法很好的模拟对流问题。在SPH方法的基础上，Huang 等人在 2015 年提出了 一种 KGF-SPH 方法（kernel gradient free-smoothed particle hydrodynamics)。该方法在核近似和粒子近似过程中只用到了核函数本身，不包含核函数 的导数，因此无论核函数是否可导都可用于该方法，是一种无核梯度的SPH方法。KGF-SPH 方法降低了无网格方法对于核函数的要求，具有较高的精度。但是KGF-SPH方法中二阶导 数通过对一阶导数再求导来求解，这样虽然提高了计算精度，但是存在以下缺点：
[0004] (1)计算量较大。因为二阶导数的求解需要在对一阶导数求导，增加了一倍的计算 量。
[0005] (2)不利于处理边界条件。因为在求解边界附近流体粒子的物理量时，需要先求解 出边界上粒子物理量的导数值，然而为了精确的求出边界上粒子物理量的导数值，还需要 在边界外布置更多的虚粒子。
【发明内容】
[0006] 本发明的目的是为了提高数值模拟的精度和稳定性，减小计算量，提出了一种改 进的KGF-SPH方法。
[0007] 本发明是通过以下技术方案实现的：
[0008] -种改进的KGF-SPH方法，具体实现步骤如下：
[0009] 步骤1、在所研究问题的计算域中布置粒子；
[0010] 步骤2、根据所研究问题的初始条件和边界条件对粒子的物理属性进行初始化；
[0011] 步骤3、利用下述空间离散格式对所研究问题的控制方程进行近似，得到函数的时 间导数#的粒子近似式，其近似原理如下： dt
[0012] 对任意函数f (r])在点r =巧处进行泰勒展开，并且只保留到一阶导数，可得：
[0013] f (r j) = f (r；) +Vf (r；) · Cr jTi) +〇 (h2) (I)
[0014] 忽略式⑴中的高阶项o (h2)，在式⑴两边同时乘以核函数W(R，h)和（rfrj W (R, h)并在计算域上进行积分得：
[0015] / f(rj)ffdV = f (r;) f WdV+V f (r;) · f (r^ri) WdV (2)
[0016] f f (Tj) (T-Ti)WdV = f (r J f (rj-r^ffdV+V f (r;) · f (r ^ri) 2WdV (3)
[0017] 联立式⑵和式（3)，并求解方程组，可得：
[0019] 上式对应的粒子近似式为：
[0021] 式中，t= f(r ;)，fj= f(r .)，（ ▽ :〇; = ▽ f(r J，N表示在i粒子的支持域内的 j粒子（即i粒子的最近相邻粒子）总量，P ,分别表示j粒子的质量和密度。
[0022] 在一维空间下，式（5)的具体表达式为：
[0024] 其中，Xjl= Xj-X1Jlix为粒子i在X方向上的一阶导数。修正矩阵厶立一维空间 下的表达式为：
[0026] 在一维空间下，对函数f (Xj)在点Xj= X 1进行二阶泰勒展开并略去高阶项可得：
[0028] 式中，fliXX为粒子i在X方向上的二阶导数。式⑶两边同时乘以核函数W并移 项可得：
[0030] 式（9)两边同时积分可得：
[0032] 由核函数的对称特性可知，上式右边第一项为0,因此式（10)可化简为：
[0034] 对（11)式进行移项可以推导出一维空间下拉普拉斯算子的离散格式为：
[0036] 式（12)对应的粒子近似式为：
[0038] 在二维空间下，式（5)的具体表达式为：
[0040] 其中，Xji= X fXi，Yji= yAjP f。分别为粒子i在X和y方向上的一阶导 数。修正矩阵Ai在二维空间下的表达式为：
[0042] 在二维空间下，对函数f(Xj，yj)在点（Xl，yi)进行二阶泰勒展开并略去高阶项可 得：
[0044] 式中，fiiXy为粒子i的混合导数，f iiX)^P f iiyy分别为粒子i在X和y方向上的二阶 导数。式（16)两边同时乘以核函数W并移项可得：
[0046] 式（17)两边同时积分得：
[0048] 由核函数的对称特性可知，式（18)右边第一项和第二项为0,因此上式可化简为：
[0050] 当粒子分布均匀且规则时将
在极坐标下进行求解可得：
[0055] 式（19)移项并利用式（22)对其进行化简，可以推导出二维空间下拉普拉斯算子 的离散格式为：
[0057] 式（23)对应的粒子近似式为：
[0059] 为了减弱粒子分布不均匀的影响，采用下式代替式（24):
[0061] 同理可推导出三维空间下任意函数及其一阶导数和拉普拉斯算子的近似方案：
[0064]其中，Xji= X .j-Xi，yji= y .jii，Zji= Z fZi，f iiZ分别为粒子 i 在 x、y 和 z方向上的一阶导数，fiiXX、4"和f iiZZ分别为粒子i在x、y和z方向上的二阶导数，x k、yk 和zk分别为粒子k在x、y和z方向上的坐标，k = i或j。修正矩阵A ""在三维空间下的表 达式为：
[0066] 式（6)、式（7)和式（13)给出了一维空间下任意函数及其一阶导数和拉普拉斯算 子的近似方案，式（14)、式（15)和式（25)给出了二维空间下任意函数及其一阶导数和拉普 拉斯算子的近似方案，式（26)、式（27)和式（28)给出了三维空间下任意函数及其一阶导数 和拉普拉斯算子的近似方案。该方案具有较高的数值精度和稳定性，便于边界条件的处理。 [0067] 步骤4、根据所研究问题的控制方程进行时间迭代，迭代一定步数后，得到模拟结 果，其中迭代公式为：
[0069] 式中，dt为迭代的时间步长。
[0070] 有益效果
[0071] 本发明对比已有技术具有如下优点：
[0072] (1)改进的KGF-SPH方法计算量较小；
[0073] (2)改进的KGF-SPH方法提高了数值模拟的精度和稳定性；
[0074] (3)改进的KGF-SPH方法便于边界条件的处理。
【附图说明】
[0075] 图1为本发明实施例1初始粒子分布示意图；
[0076] 图2为本发明实施例1基于（a) KGF-SPH方法、（b)改进的KGF-SPH方法以及（c) 解析解得到的一维热传导问题温度随时间变化的分布图；
[0077] 图3为本发明实施例1基于KGF-SPH方法和改进的KGF-SPH方法模拟一维热传导 问题得到的t = 0. 2s时的温度分布曲线与解析解的对比；
[0078] 图4为本发明实施例2初始粒子分布示意图；
[0079] 图5为本发明实施例2基于（a) SPH方法、（b)改进的KGF-SPH方法和（c) FVM方 法模拟封闭方腔自然对流得到的温度分布云图对比；
[0080] 图6为本发明实施例2基于（a) SPH方法、（b)改进的KGF-SPH方法和（C) FVM方 法模拟封闭方腔自然对流得到的水平方向速度分布云图对比；
[0081] 图7为本发明实施例2基于（a) SPH方法、（b)改