CN108920872A

CN108920872A - 针对dsmc方法的bcp粒子定位实现方法及系统

Info

Publication number: CN108920872A
Application number: CN201810836325.1A
Authority: CN
Inventors: 张斌; 刘洪�; 余彬
Original assignee: Shanghai Jiaotong University
Current assignee: Shanghai Jiaotong University
Priority date: 2018-07-26
Filing date: 2018-07-26
Publication date: 2018-11-30

Abstract

一种针对DSMC方法的BCP粒子定位实现方法及系统，通过在结构化网格上加上一层背景笛卡尔网格并分割结构化网格，然后根据粒子更新坐标和笛卡尔网格几何尺寸预定位到其所在的背景笛卡尔网格，然后通过网格转换关系精确定位得到粒子所在的结构化网格，本发明兼备结构化网格模拟复杂外形的能力以及笛卡尔网格快速定位的优点，在付出尽可能少的代价条件下(增加合适的背景笛卡尔网格)，使DSMC方法在模拟粒子迁移时，快速准确构建粒子与结构化网格的映射关系，实现针对非笛卡尔网格体系下DSMC计算过程中粒子快速定位。

Description

针对DSMC方法的BCP粒子定位实现方法及系统

技术领域

本发明涉及一种航空航天领域的技术，具体涉及一种针对直接模拟蒙特卡洛(Direct Simulation Monte Carlo-DSMC)方法的背景笛卡尔网格定位技术(BackgroundCartesian grids Positioning method-BCP)粒子定位实现方法及系统。

背景技术

在DSMC计算中，尤其是在模拟大努森数、大马赫数的来流时，粒子迁移是主要耗时部分之一。而计算采用网格的形式直接决定着模拟粒子迁移算法的计算效率。目前常用的有笛卡尔网格、结构化网格和非结构化网格三种。三种网格都有各自的优缺点。其中笛卡尔网格和结构化网格具有互补性的优点。

发明内容

本发明针对现有技术存在的上述不足，提出一种针对DSMC方法的将笛卡尔网格和结构化网格相结合的BCP粒子定位实现方法及系统，通过在计算前构建背景笛卡尔网格和结构化网格的几何关系，具备了结构化网格模拟复杂外形的能力以及笛卡尔网格快速定位的优点。本发明在付出尽可能少的代价条件下(增加合适的背景笛卡尔网格)，充分利用笛卡尔快速定位以及结构化网格模拟复杂外形的优点，使DSMC方法在模拟粒子迁移时，较于传统方法，在计算精度不变的前提下，能够更加快速、高效的实现针对非笛卡尔网格体系下的粒子定位。

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

本发明涉及一种针对DSMC方法的BCP粒子定位实现方法，通过在结构化网格上加上一层背景笛卡尔网格并分割结构化网格，然后根据粒子更新坐标和笛卡尔网格几何尺寸预定位到其所在的背景笛卡尔网格，然后通过网格转换关系精确定位得到粒子所在的结构化网格。

本发明涉及一种实现上述方法的系统，包括：结构化网格模块、背景网格划分模块、预定位模块以及精确定位模块，其中：结构化网格模块与流场信息输入相连并传输结构化网格几何信息，背景网格划分模块与结构化网格模块相连并传输两者之间的几何结构关系信息，预定位模块与背景网格划分模块相连并传输粒子与背景笛卡尔网格映射关系信息，精确定位模块与结构化网格模块相连并传输粒子与结构化网格映射关系信息。

技术效果

与现有技术相比，本发明在保证计算精度不变的前提下，能够大幅提高粒子迁移过程的计算耗时，大大提高DSMC方法的计算效率。

附图说明

图1为BCP方法定位示意图；

图2为BCP方法与传统RC定位方法流程示意图；

图2中：P为粒子，C为笛卡尔网格，S为结构网格，I为与粒子运动轨迹相交的网格面， Sa为与共用网格面I的相邻网格，椭圆形代表映射，矩形代表算法，菱形代表判断；

图3为P2L判断准则示意图；

图3中：P_i为网格节点i的坐标(x_i，y_i)，X(t)为t时刻粒子位置坐标(x(t)，y(t))；

图4为实施例1计算域示意图；

图5为不同网格计算结果温度分布云图；

图6为沿x轴线压力、温度分布线图(左)和壁面CP、CQ分布线图(右)；

图7为RM1与内存随BGR的变化曲线；

图8为计算耗时与内存随BGR的变化曲线；

图9为实施例2计算域示意图；

图10为计算结果温度分布云图对比示意图。

具体实施方式

实施例1

如图1所示，本实施例在结构化网格上加上一层背景笛卡尔网格，其背景直角网格的疏密程度，即笛卡尔网格的相对密度值为背景直角网格的平均尺寸和结构化网格平均尺寸的比值，具体为：其中：a、b表示笛卡尔网格的长和宽，V_i代表第i个结构化网格的面积；BGR值越小，表示背景笛卡尔网格越密。

本发明中背景直角网格的疏密程度BGR优选范围为：0.1-0.8

所述的背景笛卡尔网格比结构化网格尺寸小(BGR<1)，所以一个结构化网格将被至少两个直角网格分割。易知直角坐标网格被结构化网格分成两类，即：第一型网格，其背景网格完全处于任一结构化网格内部，如图1右图中的阴影部分，以及第二型网格，其背景网格在结构化网格边界处(即与多个结构化网格相交)，如图1右图中白色部分。

如图2所示，所述的BCP方法具体包括以下步骤：

i)计算某一粒子运动时间步长Δt后的新坐标为：x(t+Δt)＝x(t)+Δt×u，y(t+Δt)＝ y(t)+Δt×v；

ii)根据粒子的新坐标和笛卡尔网格几何尺寸，瞬间定位其所在的背景笛卡尔网格 (i，j):其中：表示：取数A的整数部分；

iii)根据背景笛卡尔网格与结构化网格关系，分两种情况进行相应判断；当该背景网格属于第一型网格，则通过直角坐标网格和结构化网格的对应关系确定该粒子所在的贴体结构化网格。当该背景网格属于第二型网格，则该粒子只需要与和该直角网格相交的结构化网格进行左侧粒子(Particle to the left，P2L)方法判断，进而找出粒子所在的结构化网格。

完成一个粒子的定位后，继续对下一个粒子进行步骤(i)(ii)(iii)，直至完成所有粒子的定位。图1为粒子从A运动到B的粒子BCP定位示意图。

所述的P2L方法是指：逆时针沿着组成网格单元的各个面，即连续相邻两个网格点的直线进行判断，当该粒子在所有这些面的左边，则可以确定该粒子在此网格单元内；如图3所示， P_iP_i+1×P_iX(t+Δt)的z向分量Ω_i＝(x_I+1-x_i)(y(t+Δt)-y_i)-(x(t+Δt)-x_i)(y_i+1-y_i)，其中：Ω_i>0表示粒子在网格面X(t+Δt)左边，记该网格面为P2L；反之Ω_i<0表示点在网格面X(t+Δt)右边，记该网格面为P2R。

如图1所示，对于网格单元所有面都进行z向分量Ω_i计算，当所有面z向分量Ω_i计算结果都为正数，则粒子就在这个网格里，当出现一个面使得z向分量Ω_i成负数，则停止在这个网格单元的计算，转移到共用这个面的网格单元内进行计算，直到有一个网格中的所有z向分量Ω_i都为正数。

具体地，本实施例采用经典的氩气二维圆柱绕流实施例，检验BCP方法的准确性和高效性，并检验本发明方法中核心参数BGR的判断准则。如图4所示，计算域圆柱直径L＝0.308m。来流的马赫数、密度、和静温分别为Ma＝10、ρ＝2.818E-5kg/m³、T₀＝300K。由努森数定义式(λ为平均分子自由程，L为宏观尺度)可计算出标准实施例来流努森数 Kn＝0.01。氩气的变径硬球参数设定为：温度指数α＝0.734，参考直径D＝3.595E-10m以及参考温度T＝1000K。标准实施例计算域见图6左。通过网格无关性验证研究，为确保网格的平均尺寸为自由流平均自由程的1/3，最终选取网格数为24万(沿圆柱表面和沿驻点线的网格分别为600和400个)。边界条件设置：X轴为对称边界条件、壁面为等温(TW＝500K)漫反射边界条件、其它边界为自由来流条件。最小网格中的粒子数为5，在初始化的过程中整个计算域中的粒子数约为500万。计算时间步长为5E-6s。采样次数为5000。

本实施例中采用四种背景网格数的网格：Mesh A:BGR＝0.8(1137×544)，Mesh B:BGR＝0.4(2274×1087)，Mesh C:BGR＝0.2(4547×2173)、Mesh D:BGR＝0.1(9094×4346)。传统DSMC算法即无背景网格(BGR＝0)计算。

图5为计算所得圆柱绕流流场温度分布云图，图中揭示的技术信息为：BCP定位方法与传统定位计算精度一致，背景网格相对密度不影响计算精度。

图6为沿x轴线压力、温度分布线图(左图)和壁面CP、CQ分布线图(右图)，图中揭示的技术信息为：BCP定位方法与传统定位计算精度一致，背景网格相对密度不影响计算精度。

图7为计算之前绘制的第一型网格面积占整个计算域面积的比值(RM1)与内存关系图 (BGR＝0.1-0.8)，图中揭示的技术信息为：可通过RM1和MEMO随BGR变化关系选出合适的 BGR值，对于本算例BGR取0.2附近比较合适。

图8为计算耗时与内存增加关系图(BGR＝0.1-0.8)，图中揭示的技术信息为：对于本算例， BGR取0.2附近比较合适。

表1本发明方法和传统方法计算圆柱扰流计算耗时对比结果

注：上表中MEMO代表增加的内存，DSMC传统实施例消耗内存为587MB.

本实施例证明了在计算精度不变的情况下，BCP方法具有高于传统定位方法的计算效率。还证明了背景笛卡尔网格相对密度控制参数BGR的判断准则的有效性。

实施例2

本实施例采用真空羽流以检验BCP方法的准确性和高效性。

用DSMC计算调姿发动机的喷管扩张段羽流场，计算域如图9所示，喷管内网格划分为沿壁面300，吼道200，真空区域网格划分为x轴300，y轴400。该情形下，喷管内收敛段为亚音流，扩张段内为超音流，流动在喷管喉部附近达到音速。由于喷管外的环境为真空环境，因而喷管羽流为自由膨胀流动，没有激波的出现。考虑到轴对称性，只计算轴线以上的流场参数。

由于喷管内外流动均为较低密度流动，分析时作如下基本假设:①流场中分子的碰撞均为二体碰撞；②不考虑分子内自由度激发和化学非平衡效应；③假设喷管壁面为恒温壁面(见表 1)；④气体流动为定常流动；⑤将分子视为硬球分子，即分子间作用力忽略不计。

其它计算参数如下表

喉部直径D	1.0×10^-3m
		面积比A_e/A_t(出口/喉部)	100
扩张角θ	20°
		气体	N₂
总温T	1000K
		总压P^*	9200Pa
喉部克鲁曾数K_n	5.76×10^-3
		喷管壁温T_W	1000K

图10为计算结果温度分布云图对比，图中揭示的技术信息为：BCP定位方法与传统定位计算精度一致。

表2本发明方法和传统方法计算真空羽流耗时对比结果

本实施例证明了BCP方法在计算精度不变的情况下，能够有效提高DSMC的计算效率。

上述具体实施可由本领域技术人员在不背离本发明原理和宗旨的前提下以不同的方式对其进行局部调整，本发明的保护范围以权利要求书为准且不由上述具体实施所限，在其范围内的各个实现方案均受本发明之约束。

Claims

1.一种针对DSMC方法的BCP粒子定位实现方法，其特征在于，通过在结构化网格上加上一层背景笛卡尔网格并分割结构化网格，然后根据粒子更新坐标和笛卡尔网格几何尺寸预定位到其所在的背景笛卡尔网格，然后通过网格转换关系精确定位得到粒子所在的结构化网格。

2.根据权利要求1所述的方法，其特征是，所述的背景笛卡尔网格，其背景直角网格的疏密程度，即笛卡尔网格的相对密度值为背景直角网格的平均尺寸和结构化网格平均尺寸的比值，具体为：其中：a、b表示笛卡尔网格的长和宽，V_i代表第i个结构化网格的面积。

3.根据权利要求1或2所述的方法，其特征是，所述的背景笛卡尔网格密度存在一合适范围：使内存增加不是很大，而计算机时大幅减少。其相对密度值范围为：0.1-0.8。

4.根据权利要求1所述的方法，其特征是，所述的预定位是指：计算某一粒子运动时间步长Δt后的新坐标为：x(t+Δt)＝x(t)+Δt×u，y(t+Δt)＝y(t)+Δt×v；根据粒子的新坐标和笛卡尔网格几何尺寸，瞬间定位其所在的背景笛卡尔网格其中：表示：取A的整数部分。

5.根据权利要求1所述的方法，其特征是，所述的精确定位是指：根据背景笛卡尔网格与贴体结构化网格关系，分两种情况进行相应判断；当该背景网格属于第一型网格，则通过直角坐标网格和结构化网格的对应关系确定该粒子所在的贴体结构化网格。当该背景网格属于第二型网格，则该粒子只需要与和该直角网格相交的结构化网格进行左侧粒子判断，进而实现粒子所在的结构化网格的精确定位。

6.根据权利要求1所述的方法，其特征是，所述的左侧粒子判断是指：逆时针沿着组成网格单元的各个面，即连续相邻两个网格点的直线进行判断，当该粒子在所有这些面的左边，则可以确定该粒子在此网格单元内；如图3所示，P_iP_i+1×P_iX(t+Δt)的z向分量Ω_i＝(x_I+1-x_i)(y(t+Δt)-y_i)-(x(t+Δt)-x_i)(y_i+1-y_i)，其中：Ω_i>0表示粒子在网格面X(t+Δt)左边，记该网格面为P2L；反之Ω_i<0表示点在网格面X(t+Δt)右边，记该网格面为P2R；

对于网格单元所有面都进行z向分量Ω_i计算，当所有面z向分量Ω_i计算结果都为正数，则粒子就在这个网格里，当出现一个面使得z向分量Ω_i成负数，则停止在这个网格单元的计算，转移到共用这个面的网格单元内进行计算，直到有一个网格中的所有z向分量Ω_i都为正数。

7.根据权利要求5或6所述的方法，其特征是，所述的第一型网格，其背景网格完全处于任一结构化网格内部；第二型网格，其背景网格在结构化网格边界处，即与多个结构化网格相交。

8.一种实现上述任一权利要求所述方法的系统，其特征在于，包括：结构化网格模块、背景网格划分模块、预定位模块以及精确定位模块，其中：结构化网格模块与流场信息输入相连并传输结构化网格几何信息，背景网格划分模块与结构化网格模块相连并传输两者之间的几何结构关系信息，预定位模块与背景网格划分模块相连并传输粒子与背景笛卡尔网格映射关系信息，精确定位模块与结构化网格模块相连并传输粒子与结构化网格映射关系信息。