CN103793598A

CN103793598A - 一种改进的布置均匀粒子的方法

Info

Publication number: CN103793598A
Application number: CN201410019496.7A
Authority: CN
Inventors: 雷娟棉; 黄灿; 王锁柱; 何建东; 郑志伟
Original assignee: Beijing Institute of Technology BIT
Current assignee: Beijing Institute of Technology BIT
Priority date: 2014-01-15
Filing date: 2014-01-15
Publication date: 2014-05-14

Abstract

本发明属于计算力学技术领域，涉及一种改进的布置均匀粒子的方法。该方法包括1、布置流体粒子和边界粒子。2、对粒子的物理属性初始化。3、计算边界粒子对流体粒子的边界力、由于粒子不均匀性产生的速度变化率、粒子速度的修正值。4、求解粒子速度，然后进一步求解粒子的位置。5、重复步骤2、3、4，当所有的粒子速度u<10^-m的时候认为计算已经收敛，此时输出粒子的位置，m的取值范围为2～7。本发明方法能快速有效的布置均匀的粒子，可用于光滑粒子动力学计算前处理，提高其计算的精度。

Description

一种改进的布置均匀粒子的方法

技术领域

本发明属于计算力学技术领域，涉及一种改进的布置均匀粒子的方法。

背景技术

光滑粒子动力学(smooth particle hydrodynamic)是一种纯无网格方法，在众多的无网格方法中，光滑粒子动力学具有很多独特的优势，这些优势使得光滑粒子动力学被广泛的应用于固体和流体领域。本发明能够为光滑粒子动力学在计算之前获得一个均匀分布的粒子场，对提高光滑粒子动力学计算结果的精度和稳定性以及加快计算的收敛具有重要的意义。

1977年Gingold和Monaghan，及Lucy首次利用光滑粒子动力学解决天体物理问题。然而基于光滑粒子动力学方法的流动模拟目前仍主要是二维问题，而且通常是较为简单的边界，造成这一问题的众多原因之一就是光滑粒子动力学在处理复杂边界形状问题时难以得到一个初始分布均匀的粒子场。目前很少有人对光滑粒子动力学布置均匀粒子的方法进行研究。基于核函数的归一特性和对称特性，2012年Colagrossi等人提出了一种光滑粒子动力学粒子均匀化算法。然而Colagrossi等人提出的粒子均匀化算法存在着一些问题：

(1)固体边界周围较大范围内流体粒子的位置都会发生变动;

(2)粒子的均匀程度依赖光滑长度与初始粒子间距的比值;

(3)边界附近的流体粒子可能会穿透固体表面。

基于2012年Colagrossi等人提出的粒子均匀化算法，本发明通过速度修正以及添加边界力模型，构建了一种改进的布置均匀粒子的方法，利用该方法能够对任意边界形状的计算域快速地布置均匀粒子，并且有效的解决了流体粒子穿透固体表面的问题。

发明内容

本发明的目的是为了克服上述已有技术的缺陷，实现更快速、更均匀、更方便的分布均匀的粒子和避免前处理计算时流体粒子穿透固体壁面，提出了一种改进的布置均匀粒子的方法。

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

提出了一种改进的布置均匀粒子的方法，具体实现步骤如下：

步骤1、布置流体粒子和边界粒子，对粒子的物理属性初始化。

1.1、布置流体粒子和边界粒子

一般外形物体绕流计算的流场域可用一个封闭区域和物面之间的区域来表示。先在整个封闭区域内布上均匀的流体粒子；然后去除边界上和物面内的流体粒子；最后在外边界和固壁边界上布置均匀的边界粒子。

1.2、对粒子的物理属性初始化

给出粒子的质量m、密度ρ、光滑长度h(一般取为略大于初始粒子间距的值)，粒子的初始速度v=0。

步骤2、搜索粒子对，计算核函数。

计算所有粒子之间的距离r=|x_i-x_j|，x_i和x_j分别表示任意i粒子和j粒子的坐标。核函数有很多种，这里以三次样条核函数为例，通过公式(1)计算核函数W的值，R是基于光滑长度h进行无量纲化后的粒子间距，R=|x_i-x_j|/h，为了满足正则化条件，在一维、二维和三维中αd分别取为1/h，15/7h²和3/2h³。。

W_{ij} = a_{d} \{\begin{matrix} \frac{2}{3} - R^{2} + \frac{1}{2} R^{3} & 0 \leq R < 1 \\ \frac{1}{6} {(2 - R)}^{3} & 1 \leq R < 2 \\ 0 & R \leq 2 \end{matrix} - - - (1)

步骤3、计算边界粒子对流体粒子的边界力、由于粒子不均匀性产生的速度变化率、粒子速度的修正值。

3.1计算边界粒子对流体粒子的边界力

通过公式(2)求出边界粒子j对边界附近任意流体粒子i的边界力，其中fi_j ^B表示表示边界粒子j对边界附近流体粒子i的边界力；u_i和u_j分别表示i和j粒子的速度；为了避免边界力的震荡问题，选择(u_i-u_j)·n与-φ之间较小的值进行计算，φ的取值范围为0.1～1；n为边界的外法线方向；Hi_j=hi_j ^d，hi_j=(h_i+h_j)/2，d为求解问题的维数，对二维问题，d=2，对三维问题，d=3；V_max ²表示计算过程中流体粒子的最大速度。

3.2计算由于粒子不均匀性产生的速度变化率

利用粒子均匀化公式(3)求出任意粒子i由于粒子的不均匀度产生的速度变化率du_i/dt，其中β、α为可调参数，β的取值范围为0.1～0.01，α的取值范围为0.01～0.0001。

3.3计算粒子速度的修正值

利用公式(4)求出任意粒子的修正速度u′_i，其中ε为可调无量纲参数，ε的取值范围为0.01～0.3。

f_{ij}^{B} = \{\begin{matrix} - (V_{\max}^{2} \frac{\min ((u_{i} - u_{j}) \cdot n_{j}, - φ) W_{ij} H_{ij} n_{j}}{| r_{ij} \cdot n |}), (u_{i} - u_{j}) \cdot n_{j} < 0 \\ 0, (u_{i} - u_{j}) \cdot n_{j} &GreaterEqual; 0 \end{matrix} - - - (2)

\frac{{du}_{i}}{dt} = - β {&dtri; Γ}_{i} - α \frac{\sqrt{β}}{h} u - - - (3)

{u_{i}}^{'} = - ϵ \underset{j}{Σ} V_{j} u_{ij} W_{ij} - - - (4)

步骤4、求解粒子速度，然后进一步求解粒子的位置。

4.1求解粒子速度

先对步骤(3)求得的du_i/dt、fi_j ^B值通过公式(5)进行时间积分，求出下一个时间步的待修正的速度u_i ^*，再利用步骤(3)中求出的修正速度u′_i通过公式(6)对速度u_i ^*进行修正，求出下一个时间步的真实速度u_i。

4.2求解粒子的位置

对真实的速度u_i通过公式(7)进行时间积分求出下一个时间步的粒子位置x_i。依据CFL条件，时间步长Δt为通过公式(8)来确定，其中CFL=1。

u_{i}^{*} = (\frac{{du}_{i}}{dt} + f_{ij}^{B}) Δt - - - (5)

u_{i} = u_{i}^{*} + u_{i}^{'} - - - (6)

x_i＝u_iΔt (7)

Δt \leq CFL \frac{h}{\sqrt{β}} - - - (8)

步骤5、重复步骤2、3、4，当所有的粒子速度u<10^-m的时候认为计算已经收敛，此时输出粒子的位置，m的取值范围为2～7。

自此，就实现了一种改进的布置均匀粒子的方法。

有益效果

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

(1)可有效抑制前处理过程的数值震荡，加快计算的收敛，使流体粒子分布更加均匀。观察公式(4)可知修正速度u′_i是通过其邻近粒子的速度u_j加权平均得到的，因此当i粒子速度通过u′_i进行修正后能够减小该粒子速度与其周围粒子速度的差值，进而减小了由于速度波动造成的数值震荡问题。

(2)粒子的均匀程度并不依赖光滑长度与初始粒子间距的比值。当光滑长度与初始粒子间距的比值增大后，速度的波动相应的增大，但是在通过公式(4)得到的修正速度对粒子速度进行修正后就能够有效的抑制由于速度波动引起的数值震荡问题。

(3)有效的解决了流体粒子穿透固体边界的问题。当流体粒子靠近边界并将要穿过边界时，流体粒子i的速度u_i与边界的外法线方向相反，边界粒子j始终固定不动u_j=0，因此(u_i-u_j)·n<0，在这种条件下公式(2)将会计算出与边界外法线方向相同的边界力fi_j ^B，边界力会随着流体粒子i靠近边界粒子j而增大，这种边界力最终会改变流体粒子速度u_i的方向，进而阻止了流体粒子穿透边界。

附图说明

图1一种改进的布置均匀粒子的方法流程图；

图2粒子i的支持域；

图3为本发明实施例；其中，a为计算域及其边界示意图；b为在计算域内布置均匀粒子示意图；c为去除边界上和边界外的粒子示意图；d为在边界上布置均匀的粒子示意图；e为发明前的效果图；f为发明后的效果图。

具体实施方式

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

实施例1下面以矩形封闭域，二维翼型绕流场的粒子分布为例来说明本发明的技术方案。如图1所示。

步骤1

以二维翼型绕流场的粒子分布为例，对弦长b=0.1m的NACA0012翼型绕流场，取外边界为长和宽分别为0.2m和0.2m的矩形计算域，见图3a；首先在整个矩形域内及其边界上布置上均匀的流体粒子，粒子间距为dx=dy=0.005m，见图3b；然后去除矩形域边界上、翼型表面、翼型内部的流体粒子，见图3c；最后在矩形域边界上和上下翼面布置均匀的边界粒子，翼型边界上的粒子间距与初始流体粒子间距相接近（ds≈dx），矩形域边界上布置两排与初始流体粒子间距相接等边界粒子，见图3d。粒子的密度为标况下空气的密度ρ=1.27kg/m³，粒子的质量m=ρdxdy，光滑长度h=1.2dx，粒子的初始速度u=0。

步骤2

如图2所示，计算所有粒子之间的距离r=|x_i-x_j|，计算无量纲粒子间距R=|x_i-x_j|/h，当0<R<1时通过公式(1)中的1式计算核函数W，当1≤R<2时通过公式(1)中的2式计算核函数W，当R为其它值时W=0。

步骤3、先通过公式(2)求出边界粒子j对边界附近任意流体粒子i的边界力，其中φ=1。再利用公式(3)求出任意粒子i由于粒子的不均匀度产生的速度变化率du_i/dt，其中β=0.01，α的=0.001。最后利用公式(4)求出任意粒子的修正速度u′_i，其中ε=0.2。

步骤4、先对步骤(3)求得的du_i/dt、fi_j ^B值通过公式(5)进行时间积分，求出下一个时间步的待修正的速度u_i ^*，再利用步骤(3)中求出的修正速度u′_i通过公式(6)对速度u_i ^*进行修正，求出下一个时间步的真实速度u_i。对真实的速度u_i通过公式(7)进行时间积分求出下一个时间步的粒子位置x_i。依据CFL条件，时间步长Δt=0.005。

步骤5、重复步骤2、3、4，当所有的粒子速度u≤10^-3时认为计算已经收敛，此时输出粒子的位置。结果如图3e、f。图3e为采改进的前处理方法进行粒子均匀化的结果，由图可见翼型周围较大区域内粒子位置都发生了变化，而且有粒子穿透了翼型表面。图3f是采用本发明提出的前处理方法，得到的流体粒子分布结果，由图3e、f可见采用本文的前处理方法，不仅解决了粒子穿透翼型边界的问题，而且得到了更均匀分布的粒子场。

Claims

1.一种改进的布置均匀粒子的方法，其特征在于，具体实现步骤如下：

步骤1、布置流体粒子和边界粒子，对粒子的物理属性初始化

1.1、布置流体粒子和边界粒子

先在整个封闭区域内布上均匀的流体粒子；然后去除边界上和物面内的流体粒子；最后在外边界和固壁边界上布置均匀的边界粒子；

1.2、对粒子的物理属性初始化

给出粒子的质量m、密度ρ、光滑长度h，粒子的初始速度v=0；

步骤2、搜索粒子对，计算核函数

步骤3、计算边界粒子对流体粒子的边界力、由于粒子不均匀性产生的速度变化率、粒子速度的修正值

3.1计算边界粒子对流体粒子的边界力

通过公式

f_{ij}^{B} = \{\begin{matrix} - (V_{\max}^{2} \frac{\min ((u_{i} - u_{j}) \cdot n_{j}, - φ) W_{ij} H_{ij} n_{j}}{| r_{ij} \cdot n |}), (u_{i} - u_{j}) \cdot n_{j} < 0 \\ 0, (u_{i} - u_{j}) \cdot n_{j} &GreaterEqual; 0 \end{matrix}

求

出边界粒子j对边界附近任意流体粒子i的边界力，其中fi_j ^B表示表示边界粒子j对边界附近流体粒子i的边界力；u_i和u_j分别表示i和j粒子的速度；为了避免边界力的震荡问题，选择(u_i-u_j)·n与-φ之间较小的值进行计算；n为边界的外法线方向；Hi_j=hi_j ^d，hi_j=(h_i+h_j)/2，d为求解问题的维数，对二维问题，d=2，对三维问题，d=3；V_max ²表示计算过程中流体粒子的最大速度；

3.2计算由于粒子不均匀性产生的速度变化率

利用粒子均匀化公式