CN115994463A - 一种受限空间三维流场的数值模拟方法、存储介质及系统 - Google Patents

Abstract

本发明涉及一种受限空间三维流场的数值模拟方法，包括步骤：进行控制方程的推导及求解；进行参数无量纲化；建立数值模型，并确定初始条件及边界条件；利用数值模型进行方程组求解。本发明还提供一种存储介质及受限空间三维流场的数值模拟系统，采用一种受限空间三维流场的数值模拟方法、存储介质及系统可采用直接数值模拟模型，以有限差分方法为基础，基于涡量‑速度方法推导了三维非定常不可压缩Burgers模型方程的求解方法，得到受限空间三维流场的数值模拟程序算法。

Description

一种受限空间三维流场的数值模拟方法、存储介质及系统

技术领域

本发明属于核设施放射性气溶胶输运基础研究技术领域，具体涉及一种受限空间三维流场的数值模拟方法、存储介质及系统。

背景技术

在核设施的典型事故中，如泄漏事故、燃烧事故、爆炸事故等均可能涉及放射性气溶胶的输运。在核设施的设计、安全分析和事故应急中，放射性气溶胶在厂房内部输运至厂房外部之前，其中一个非常重要的环节是在受限空间内的输运，不同的事故情景、事故原因将导致受限空间内的建筑结构和流场状态参数随之发生变化，放射性气溶胶在受限空间内的输运过程也将受到严重影响，只有在准确得到流场状态参数之后才能准确计算放射性气溶胶的浓度分布。

目前，Navier-Stokes方程的求解思路有两种：一是求解一个含有速度和压力的动量方程和一个关于压力的泊松方程，即压力-速度方程；二是通过将方程两边取旋度，得到涡量输运方程，从动量方程中消去求解压力项带来的困难，即涡量-速度方程。涡量速度的形式的边界条件容易找到，并且解得的涡量能直接表示出流动的特性。

控制方程的数值求解方法有多种，其中以有限差分方法、有限体积方法和有限元方法应用最为广泛。有限差分方法思路简捷明了、编程容易和格式构造方便灵活、精度高；系数矩阵为稀疏矩阵，便于求解。有限体积方法对复杂物理模型的高度适应性；积分型守恒方程；控制方程无法做到维度解耦，构造高精度格式较为困难。有限元方法计算精度可提高，可求解复杂的物理模型；但其数值间断处存在虚假震荡，没有局部守恒，只能保证全局守恒。

发明内容

针对现有技术中存在的缺陷，本发明的目的是提供一种受限空间三维流场的数值模拟方法、存储介质及系统，采用直接数值模拟模型，得到受限空间三维流场的数值模拟程序算法。

为达到以上目的，本发明采用的技术方案是：一种受限空间三维流场的数值模拟方法，包括步骤：进行控制方程的推导及求解；进行参数无量纲化；建立数值模型，并确定初始条件及边界条件；利用数值模型进行方程组求解。

进一步，所述控制方程的推导及求解时，采用涡量-速度方法求解三维非定常不可压缩Burgers模型方程，以涡量和速度作为求解变量，通过超松弛迭代方法反复迭代求解涡量输运方程和速度泊松方程，得到涡量和速度的数值解。

进一步，所述进行参数无量纲化时，引入特征长度、特征速度及特征时间，将所述控制方程中坐标及速度表示为无量纲形式。

进一步，所述边界条件包括涡量的边界条件及速度的边界条件，所述速度的边界条件采用第一类边界，所述涡量的边界条件须采用第二类边界。

进一步，求解任意一组方程时，采用升级后的SOR进行求解，形式为：

其中，表示方程左端未知项；、分别表示方程的数量和未知数的数量；为松弛因子；为方程右端已知项；为方程左端系数矩阵的元素，i，j表示系数矩阵A的行和列，i_max＝j_max＝x_max×y_max×z_max，x_max、y_max、z_max分别表示笛卡尔坐标系下三维空间X、Y、Z方向网格数量的最大值。

进一步，所述升级后的SOR求解时，未知数的编号至多为7个。

进一步，所述数值模型包括1个主程序及7个子程序，所述子程序包括3个涡量子程序、3个速度子程序及1个超松弛迭代子程序。

进一步，所述主程序用于参数定义、网格划分、参数初始化、迭代循环及结果输出。

本发明还提供一种存储介质，其上存储有计算机程序，所述计算机程序被处理器执行时实现一种受限空间三维流场的数值模拟方法。

本发明还提供一种受限空间三维流场的数值模拟系统，包括：推导单元，用于进行控制方程的推导及求解；无量纲化单元，用于进行参数无量纲化；模型单元，用于建立数值模型，并确定初始条件及边界条件；求解单元，用于利用数值模型进行方程组求解。

本发明的效果在于：可针对核设施事故工况下放射性气溶胶在受限空间内的输运特点，采用直接数值模拟模型，以有限差分方法为基础，基于涡量-速度方法推导了三维非定常不可压缩Burgers模型方程的求解方法，得到受限空间三维流场的数值模拟程序算法。

附图说明

图1是本发明中一种受限空间三维流场的数值模拟方法的步骤流程图；

图2是数值模型中程序的结构示意图。

具体实施方式

下面结合附图和具体实施方式对本发明作进一步描述。

如图1-2所示，本发明提出了一种受限空间三维流场的数值模拟方法包括步骤：

S1，进行控制方程的推导及求解；

具体的，在进行后续的步骤前，首先需要对不可压缩流体控制方程进行推导，其中，不可压缩流体控制方程如下：

质量守恒方程，三维不可压缩质量守恒方程形式为：

其中，

为笛卡尔坐标系下的梯度算子；U(ux，uy，uz)是笛卡尔坐标系下流体的三维速度矢量，ux，uy，uz是t时刻流体在点(x,y,z)处的速度分量。

故连续性方程可表示为如下形式：

动量守恒方程：

Euler方程，其形式如下：

其中，t为时间，单位：s；ρ为流体密度，单位：kg/m3；U为流体速度矢量，单位：m/s；f为体积力矢量，单位：N；p为流体压力，单位：Pa。

Navier-Stokes方程(简称“N-S方程”)

其形式如下：

其中，Π为黏性应力张量。

Burgers模型方程

其形式如下：

其中，

为拉普拉斯算子，Re是雷诺数。

本实施例采用涡量-速度方法求解三维非定常不可压缩Burgers模型方程，以涡量和速度作为求解变量，只需通过超松弛迭代方法反复迭代求解涡量输运方程和速度泊松方程，就可得到相应的涡量和速度的数值解。

涡量-速度方法的详细推导过程如下：

其中，涡量的求解为：

对(1.5)式两边取旋度，可以得到笛卡儿坐标系下三维非定常不可压缩Burgers模型方程的涡量输运方程，形式如下：

其中，Vor(Vorx，Vory，Vorz)为涡量矢量，定义为：

将(0.7)式整理得：

对式(0.6)两边取旋度，并使用连续性方程(0.1)，可以得到直角坐标系下三维非定常不可压缩Burgers模型方程的涡量输运方程，形式如下：

涡量对时间的导数采用一阶差分进行离散，离散形式如下：

式中，dt表示时间步长，dt的取值可根据courant数计算：

式中，|U|为网格内的特征速度，h为速度方向的空间步长。

为使对流项能满足二阶精度格式要求，同时能满足其系数矩阵具有对角占优性，采用迎风格式对对流项进行离散，将离散后代数方程的截断误差保留至二阶导数项，并对截断误差的二阶导数项采用旧时刻的具有二阶精度格式中心差分进行修正，可得：

当ux(x,y,z)≥0，uy(x,y,z)≥0，uz(x,y,z)≥0时采用：

当ux(x,y,z)<0，uy(x,y,z)<0，uz(x,y,z)<0时采用：

其中：

对涡量的二阶导数采用中心差分进行离散，离散形式如下：

对速度的一阶导数同样采用中心差分进行离散，离散形式如下：

将(0.12)式至(0.27)式代入(0.9)式，并对含有

的项进行整理并合并同类项：

可得参数

的具体形式如下：

其它参数同理可得：

对含有

的项进行整理并合并同类项，可得参数

的具体形式如下：

其它参数同理可得：

将参数

代入(0.9)式并整理可得涡量在x方向的分量(Vorx)的迭代公式：

该二阶差分格式离散方程的系数矩阵具有对角占优性，可以进行迭代求解，在本实施例中，采用升级之后的超松弛(Successive Over Relaxation，SOR)迭代法进行求解。

同理对(0.10)、(0.11)式整理可得，涡量在y、z方向的分量的迭代公式。

速度的求解为：

对(0.6)式两边取旋度并根据连续性方程(0.1)式，可以得到直角坐标系下三维非定常不可压缩Burgers模型方程的速度泊松方程，形式如下：

对(0.44)式左边速度的二阶导数项采用中心差分进行离散，离散形式如下：

对(0.44)式右边涡量的一阶导数项同样采用具有二阶精度的中心差分格式进行离散，离散形式如下：

将(0.47)式至(0.50)式代入(0.44)式整理得：

速度在x方向的分量(ux)的迭代公式：

同理对(0.45)式整理可得：

速度在y方向的分量(uy)的迭代公式：

同理对(0.46)式整理可得：

速度在z方向的分量(uz)的迭代公式

S2，进行参数无量纲化；

具体的，引入特征长度(L)、特征速度(U)和特征时间(τ)，则控制方程中坐标(x，y，z)和速度表示为无量纲形式：

x＝x′/L   (0.54)

y＝y′/L   (0.55)

z＝z′/L   (0.56)

ux＝ux′/U   (0.57)

uy＝uy′/U   (0.58)

uz＝uz′/U   (0.59)

其中，x′、y′、z′、ux′、uy′、uz′和t分别表示实际的坐标、速度和时间。

S3，建立数值模型，并确定初始条件及边界条件；

具体的，速度的初始条件相同，其形式如下：

ux＝0   (0.61)

uy＝0   (0.62)

uz＝0   (0.63)

边界条件包括涡量的边界条件及速度的边界条件，速度的边界条件采用第一类边界(Dirichlet边界条件)，因为涡量是速度的偏导数的差值，故涡量的边界条件须采用第二类边界条件(Neumann边界条件)。

其中，涡量的边界条件

壁面边界条件如下：

当i＝0时

当i＝ni时

当j＝0时

当j＝nj时

当k＝0时

当k＝nk时

速度的边界条件

速度在壁面边界条件均相同，如下：

ux＝0   (0.82)

uy＝0   (0.83)

uz＝0   (0.84)

S4，利用数值模型进行方程组求解；

具体的，在给定了求解方程组所需的前置条件(边界条件及初始条件)之后，即可对方程组进行数值求解。采用升级之后的SOR进行求解，将(0.43)式的系数矩阵展开，形式如下：

其中，X表示方程左端未知项；m、n分别表示方程的数量和未知数的数量；ω为松弛因子；b_i为方程右端已知项；a_ij为方程左端系数矩阵A的元素，i，j表示系数矩阵A的行和列，i_max＝j_max＝x_max×y_max×z_max，x_max、y_max、z_max分别表示笛卡尔坐标系下三维空间X、Y、Z方向网格数量的最大值，因此传统SOR需确定系数矩阵A的每一行、列的编号。

在求解其中一组方程时，与任意网格节点(x，y，z)相关的节点至多只有(x-1，y，z)、(x+1，y，z)、(x,y-1，z)、(x,y+1，z)、(x,y,z-1)和(x,y,z+1)六个，所以在求解其中任一方程组时，未知数的编号至多只有7个，故将传统超松弛迭代公式升级为：

升级之后的逐次超松弛迭代法数据占用量为原来的7/(x_max×y_max×z_max)，同时计算速度也得到提升。

可以理解，未知数编号之多只有7个，可节省在三维与二维转化过程中所需用到的内存，进而节省空间。

需说明的是，数值模型分为1主程序和7子程序，子程序包括涡量子程序(3个)、速度子程序(3个)和超松弛迭代子程序(1个)。

主程序的主要功能是对模型进行网格划分、参数初始化以及调用涡量方程、速度泊松方程的数值计算结果，以进行下一时间步的迭代计算。

主程序的功能包括：

参数定义，包括模型尺寸参数(x_start,x_stop,y_start,y_stop,z_start,z_stop)、网格空间步长(h)、移动上壁面速度(u_upwall)、时间步长(dt)和雷诺数(Re)等参数。

网格划分，主要功能是对模型进行网格划分，需要依次对点、线、面进行网格划分。

参数初始化，适宜的参数初值可以极大的缩减计算时间，快速计算出满足线性方程组的最优解，主程序的参数初始化包括速度、涡量以及时间等。

迭代循环，迭代循环部分的主要功能是：设置动态涡量值和动态速度值；调用涡量、速度子程序；参数无量纲化。

结果输出，结果输出部分包括涡量、速度和时间等参数的输出。

涡量子程序与速度子程序：

涡量方程与速度泊松方程均采用SOR进行求解，二者程序结构基本相同，仅方程左右两端的系数不同，其主要功能都是在确定方程左右两端系数矩阵元素的具体形式之后，调用超松弛迭代子程序进行迭代求解，计算出涡量和速度的值并输出，故仅对涡量子程序进行说明，包括：

确定涡量方程左端系数矩阵元素

确定涡量方程左端系数矩阵元素部分的主要功能是确定涡量方程左端项系数对应的元素。

涡量方程右端项处理

涡量方程右端项处理主要功能是给出满足边界以及内部节点的涡量方程右端项的形式，包括内部节点以及次边界点右端项的形式。

涡量方程左端项处理

主要功能是给出满足边界以及内部节点的涡量方程左端项的形式，包括内部节点以及次边界点未知数系数的特殊处理。

涡量方程左端系数矩阵的形式

主要功能是将传统超松弛迭代系数矩阵进行升级，即方程组系数矩阵的每一行至多只有7个元素(未知数)。

调用超松弛迭代子程序

主要功能是调用超松弛迭代子程序求解涡量方程，并将计算结果输出。

超松弛迭代子程序的功能包括：

迭代参数设置

主要功能是设置迭代误差、迭代步数以及松弛因子等迭代控制参数。

迭代计算

将升级后的系数矩阵的元素与升级之前的系数矩阵的元素一一对应，采用SOR逐次对线性方程组进行求解。

结果输出

将迭代计算结果返回到主程序中。计算结果包括，涡量各分量结果Vorx、Bory、Vorz，速度分量结果ux、uy、uz以及时间time。

本发明还提供一种存储介质，其上存储有计算机程序，所述计算机程序被处理器执行时实现实施一种受限空间三维流场的数值模拟方法的步骤。

需要说明的是，本申请所示的存储介质可以是计算机可读信号介质或者存储介质或者是上述两者的任意组合。存储介质例如可以是——但不限于——电、磁、光、电磁、红外线、或半导体的系统、系统或器件，或者任意以上的组合。存储介质的更具体的例子可以包括但不限于：具有一个或多个导线的电连接、便携式计算机磁盘、硬盘、随机访问存储器(RAM)、只读存储器(ROM)、可擦式可编程只读存储器(EPROM或闪存)、光纤、便携式紧凑磁盘只读存储器(CD-ROM)、光存储器件、磁存储器件、或者上述的任意合适的组合。在本申请中，存储介质可以是任何包含或存储程序的有形介质，该程序可以被指令执行系统、系统或者器件使用或者与其结合使用。而在本申请中，存储介质可以包括在基带中或者作为载波一部分传播的数据信号，其中承载了计算机可读的程序代码。这种传播的数据信号可以采用多种形式，包括但不限于电磁信号、光信号或上述的任意合适的组合。存储介质还可以是存储介质以外的任何计算机可读介质，该计算机可读介质可以发送、传播或者传输用于由指令执行系统、系统或者器件使用或者与其结合使用的程序。计算机可读介质上包含的程序代码可以用任何适当的介质传输，包括但不限于：无线、电线、光缆、RF等等，或者上述的任意合适的组合。

本发明还提供一种受限空间三维流场的数值模拟系统，其包括：

推导单元，用于进行控制方程的推导及求解；

无量纲化单元，用于进行参数无量纲化；

模型单元，用于建立数值模型，并确定初始条件及边界条件；

求解单元，用于利用数值模型进行方程组求解。

通过上述实施例可以看出，本发明的有益效果在于：可针对核设施事故工况下放射性气溶胶在受限空间内的输运特点，采用直接数值模拟模型，以有限差分方法为基础，基于涡量-速度方法推导了三维非定常不可压缩Burgers模型方程的求解方法，得到受限空间三维流场的数值模拟程序算法。

本发明所述并不限于具体实施方式中所述的实施例，本领域技术人员根据本发明的技术方案得出其他的实施方式，同样属于本发明的技术创新范围。

Claims (10)

1.一种受限空间三维流场的数值模拟方法，其特征在于，包括：

进行控制方程的推导及求解；

进行参数无量纲化；

建立数值模型，并确定初始条件及边界条件；

利用数值模型进行方程组求解。

2.如权利要求1所述的一种受限空间三维流场的数值模拟方法，其特征在于：

所述控制方程的推导及求解时，采用涡量-速度方法求解三维非定常不可压缩Burgers模型方程，以涡量和速度作为求解变量，通过超松弛迭代方法反复迭代求解涡量输运方程和速度泊松方程，得到涡量和速度的数值解。

3.如权利要求1所述的一种受限空间三维流场的数值模拟方法，其特征在：

所述进行参数无量纲化时，引入特征长度、特征速度及特征时间，将所述控制方程中坐标及速度表示为无量纲形式。

4.如权利要求1所述的一种受限空间三维流场的数值模拟方法，其特征在于：

所述边界条件包括涡量的边界条件及速度的边界条件，所述速度的边界条件采用第一类边界，所述涡量的边界条件须采用第二类边界。

5.如权利要求1所述的一种受限空间三维流场的数值模拟方法，其特征在于：

求解任意一组方程时，采用升级后的SOR进行求解，形式为：

其中，X表示方程左端未知项；m、n分别表示方程的数量和未知数的数量；ω为松弛因子；b_i为方程右端已知项；a_ij为方程左端系数矩阵A的元素，i，j表示系数矩阵A的行和列，i_max＝j_max＝x_max×y_max×z_max，x_max、y_max、z_max分别表示笛卡尔坐标系下三维空间X、Y、Z方向网格数量的最大值。

6.如权利要求4所述的一种受限空间三维流场的数值模拟方法，其特征在于：

所述升级后的SOR求解时，未知数的编号至多为7个。

7.如权利要求1所述的一种受限空间三维流场的数值模拟方法，其特征在于：

所述数值模型包括1个主程序及7个子程序，所述子程序包括3个涡量子程序、3个速度子程序及1个超松弛迭代子程序。

8.如权利要求7所述的一种受限空间三维流场的数值模拟方法，其特征在于：

所述主程序用于参数定义、网格划分、参数初始化、迭代循环及结果输出。

9.一种存储介质，其上存储有计算机程序，其特征在于：所述计算机程序被处理器执行时实现权利要求1至8任一项所述的一种受限空间三维流场的数值模拟方法。

10.一种受限空间三维流场的数值模拟系统，其特征在于，包括：

推导单元，用于进行控制方程的推导及求解；

无量纲化单元，用于进行参数无量纲化；

模型单元，用于建立数值模型，并确定初始条件及边界条件；

求解单元，用于利用数值模型进行方程组求解。

Images