CN112182909A - 一种用于工业cae方向的流动求解器建立方法 - Google Patents

Abstract

本发明提出了一种用于工业CAE方向的流动求解器建立方法，包括：采用以下核心求解器：LBGK基础求解器、自由表面模型、SMA亚格子湍流模型，对工业CAE在应用于压铸分析时，出现的大紊流超高雷诺数内部流动进行模拟计算，包括：构建LBGK基础求解器；在所述LBGK基础求解器的基础上，采用多松弛时间MRT模型对步骤S1计算的数据进行稳定性增强；采用亚格子SMA湍流模型，对MRT模型的处理MRT模型的松弛矩阵中主对角线元素进一步处理；采用自由表面模型对所述SMA湍流模型处理后的数据，通过单相流并耦合自由表面模型实现两相流流动过程两相界面演化做出准确的模拟计算；采用基于消息传递界面机制MPI并行模型优化。

Description

一种用于工业CAE方向的流动求解器建立方法

技术领域

本发明涉及数字化模拟仿真技术领域，特别涉及一种用于工业CAE方向的流动求解器建立方法。

背景技术

工业CAE领域的内流场流动求解器通常基于Navier-Stokes方程采用不同的离散方法进行求解，如有限元方法、有限体积法、有限差分法等。基于有限元法的传统流动求解器如Fluent等对前处理生成计算网格要求很高，大大限制了其应用领域。基于有限体积法和有限差分法的流动求解器在生成计算网格过程要求不高，但多受限于中等雷诺数流动，对高速，大紊流超高雷诺数流动求解存在较大限制。

发明内容

本发明的目的旨在至少解决所述技术缺陷之一。

为此，本发明的目的在于提出一种用于工业CAE方向的流动求解器建立方法。

为了实现上述目的，本发明的实施例提供一种用于工业CAE方向的流动求解器建立方法，采用以下核心求解器：LBGK基础求解器、自由表面模型、SMA亚格子湍流模型，对工业CAE在应用于压铸分析时，出现的大紊流超高雷诺数内部流动进行模拟计算，具体包括以下步骤：

步骤S1，构建LBGK基础求解器，利用所述LBGK基础求解器首先对空间进行离散，或读入计算网格模型数据，然后对速度和时间进行离散，计算碰撞项，更新密度分布函数、宏观场速度和密度；

步骤S2，在所述LBGK基础求解器的基础上，采用多松弛时间MRT模型对步骤S1计算的数据进行稳定性增强；

步骤S3，采用亚格子SMA湍流模型，对MRT模型的处理MRT模型的松弛矩阵中主对角线元素进一步处理，提高流动求解器在大雷诺数条件下计算稳定性；

步骤S4，采用自由表面模型对所述SMA湍流模型处理后的数据，通过单相流并耦合自由表面模型实现两相流流动过程两相界面演化做出准确的模拟计算；

步骤S5，采用基于消息传递界面机制MPI并行模型对上述步骤S1至步骤S4中的算法实现并行优化。

进一步，在所述步骤S1中，所述LBGK基础求解器对速度用D3Q19模型进行离散；采用显式欧拉方法对时间进行离散。

进一步，在所述步骤S4中，所述自由表面模型采用VOF模型。

进一步，在所述步骤S2中，所述MRT模型针对LBGK基础求解器计算的数据，将密度空间的松弛时间τ转换至矩空间中，然后对矩空间中多个松弛时间分量构成的松弛矩阵参数进行调整数学描述如下所述：

f_i(x+c_iΔt,t+Δt)-f_i(x,t)＝-M^-1SM[f_i(x,t)-f_i ^eq(x,t)]Δt

其中，M为变换矩阵，S为松弛矩阵；

S数学表达式为：S＝diag(0,ω_e,ω_ε,0,ω_q,0,ω_q,ω_v,ω_v,ω_π,ω_v,ω_v,ω_v,ω_m,ω_m,ω_m)。

进一步，在所述步骤S3中，所述SMA模型对于松弛矩阵中主对角线元素做进一步处理，过程如下：

其中

式中，C_s为格子声速，Δx为格子间距，τ₀是本征松弛时间，Q为流体力学Q函数。

根据本发明实施例的用于工业CAE方向的流动求解器建立方法，采用了以下技术：基于MRT和SMA模型的LBM大雷诺数流动求解方法；基于VOF模型的LBM自由表面流动求解方法；基于MPI的LBM高性能流动求解方法。通过采用基于MRT和SMA模型的LBM方法，可以对大雷诺数(Re-10⁶)紊流过程进行准确的模拟计算；同时本发明方案还耦合自由表面模型，对于大密度比(1000倍)的两相流流动过程两相界面演化做出准确的模拟计算。其次，本流动求解器通过基于消息传递界面机制的并行算法优化，显著提升流动求解器的计算效率。

本发明附加的方面和优点将在下面的描述中部分给出，部分将从下面的描述中变得明显，或通过本发明的实践了解到。

附图说明

本发明的上述和/或附加的方面和优点从结合下面附图对实施例的描述中将变得明显和容易理解，其中：

图1为根据本发明实施例的用于工业CAE方向的流动求解器建立方法的流程图；

具体实施方式

下面详细描述本发明的实施例，实施例的示例在附图中示出，其中自始至终相同或类似的标号表示相同或类似的元件或具有相同或类似功能的元件。下面通过参考附图描述的实施例是示例性的，旨在用于解释本发明，而不能理解为对本发明的限制。

本发明提出一种用于工业CAE方向的流动求解器建立方法，首先对本申请应用的格子玻尔兹曼算法(LatticeBoltzmannMethod,LBM)算法进行说明，LBM算法作为一种崭新的计算流体力学(CFD)算法，相较基于Naiver-Stokes方程的SIMPLE算法等具有更优异的并行性能。LBM算法在汽车、航空等大规模、大雷诺数空气动力学计算中已得到成功应用。

本发明提出的流动求解器建立方法，基于MPI并行模型，其核心求解器采用基于MRT的LBM模型、VOF自由表面模型、SMA亚格子湍流模型，可模拟雷诺数高达107的内部流动。通过对算法架构以及数据存储方式的多次优化，经测试在200核范围内，并行效率基本呈现线性增长，相较传统内流场算法和模拟计算程序，计算效率得到一个数量级提升。

如图1所示，本发明实施例的用于工业CAE方向的流动求解器建立方法，采用以下核心求解器：LBGK基础求解器、自由表面模型、SMA亚格子湍流模型，对工业CAE在应用于压铸分析时，出现的大紊流超高雷诺数内部流动进行模拟计算，具体包括以下步骤：

步骤S1，构建LBGK基础求解器。

LBM基础为Lattice-Bhatnagar-Gross-Krook(LBGK)方法。LBM基本模型如下所示：

f_i(x+c_iΔt,t+Δt)＝f_i(x,t)+Ω_i(x,t)

其中，i为第i离散速度方向，x为空间位置坐标，c_i为第i离散速度分量，t为当前时刻，Δt为时间步长。f_i(x+c_iΔt,t+Δt)为x+c_iΔt位置，t+Δt时刻密度分布函数，f_i(x,t)为x位置，t时刻密度分布函数。

Ω_i(x,t)为碰撞项。

经过LBGK近似后得到：

平衡分布函数定义为:

其中w_i为第i离散速度方向上密度权重，ρ是流体的密度，u是当前位置、时刻的流体速度，c_s是格子声速，τ为松弛时间。

宏观场密度ρ由下式给出：

宏观场速度u由下式给出：

LBGK基本计算过程如下步骤：

利用LBGK基础求解器首先对空间进行离散，或读入计算网格模型数据，然后对速度和时间进行离散，计算碰撞项，更新密度分布函数、宏观场速度和密度。

在本发明的实施例中，LBGK基础求解器对速度用D3Q19模型进行离散；采用显式欧拉方法对时间进行离散。

步骤S2，在LBGK基础求解器的基础上，采用多松弛时间MRT模型对步骤S1计算的数据进行稳定性增强。

LBGK模型在计算流体过程中稳定性很低，为提高计算过程稳定性，本发明方法在LBGK模型基础上采用多松弛时间(MRT)模型。LBGK为保证计算稳定，需使得松弛时间τ>0.5，这个限制条件非常严苛，大部分应用场景均为大雷诺数流动，无法满足此稳定性条件。

MRT模型针对LBGK基础求解器计算的数据，将密度空间的松弛时间τ转换至矩空间中，然后对矩空间中多个松弛时间分量构成的松弛矩阵参数进行调整数学描述如下：

f_i(x+c_iΔt,t+Δt)-f_i(x,t)＝-M^-1SM[f_i(x,t)-f_i ^eq(x,t)]Δt

其中，M为变换矩阵，S为松弛矩阵；

S数学表达式为：S＝diag(0,ω_e,ω_ε,0,ω_q,0,ω_q,ω_v,ω_v,ω_π,ω_v,ω_v,ω_v,ω_m,ω_m,ω_m)。

步骤S3，采用亚格子SMA湍流模型，对MRT模型的处理MRT模型的松弛矩阵中主对角线元素进一步处理，提高流动求解器在大雷诺数条件下计算稳定性。

在本步骤中，为进一步提高流动求解器在大雷诺数条件下计算稳定性，采用亚格子Smovofsky(SMA)紊流模型。

SMA模型对于松弛矩阵中主对角线元素做进一步处理，过程如下：

其中

式中，C_s为格子声速，Δx为格子间距，τ₀是本征松弛时间，Q为流体力学Q函数。

步骤S4，采用自由表面模型对SMA湍流模型处理后的数据，通过单相流并耦合自由表面模型实现两相流流动过程两相界面演化做出准确的模拟计算。

在本发明的实施例中，自由表面模型采用VOF模型。

具体的，在内流场计算流体力学中，自由表面是液体和气体之间的界面，VOF自由表面模型是气相和液相密度比大于1000时对于两相流的一种较好近似，可通过单相流并耦合自由表面模型实现大密度比两相流模型计算。本发明自由表面模型采用VOF模型，VOF模型存在三要素，即自由表面分布方式，自由表面追踪方法以及在自由表面上边界条件的添加。

步骤S5，采用基于消息传递界面机制MPI并行模型对上述步骤S1至步骤S4中的算法实现并行优化。

具体的，LBM方法碰撞项为局部项，大部分的计算时间集中在碰撞项，其余时间主要分布在密度函数值的更新。因此LBM算法具有天生的并行性。

在本发明中，采用基于消息传递机制MPI机制对算法做并行优化，使得同等计算规模下，计算效率提高越100倍。

根据本发明实施例的用于工业CAE方向的流动求解器建立方法，采用了以下技术：基于MRT和SMA模型的LBM大雷诺数流动求解方法；基于VOF模型的LBM自由表面流动求解方法；基于MPI的LBM高性能流动求解方法。通过采用基于MRT和SMA模型的LBM方法，可以对大雷诺数(Re-10⁶)紊流过程进行准确的模拟计算；同时本发明方案还耦合自由表面模型，对于大密度比(1000倍)的两相流流动过程两相界面演化做出准确的模拟计算。其次，本流动求解器通过基于消息传递界面机制的并行算法优化，显著提升流动求解器的计算效率。

在本说明书的描述中，参考术语“一个实施例”、“一些实施例”、“示例”、“具体示例”、或“一些示例”等的描述意指结合该实施例或示例描述的具体特征、结构、材料或者特点包含于本发明的至少一个实施例或示例中。在本说明书中，对上述术语的示意性表述不一定指的是相同的实施例或示例。而且，描述的具体特征、结构、材料或者特点可以在任何的一个或多个实施例或示例中以合适的方式结合。

尽管上面已经示出和描述了本发明的实施例，可以理解的是，上述实施例是示例性的，不能理解为对本发明的限制，本领域的普通技术人员在不脱离本发明的原理和宗旨的情况下在本发明的范围内可以对上述实施例进行变化、修改、替换和变型。本发明的范围由所附权利要求及其等同限定。

Claims (5)

1.一种用于工业CAE方向的流动求解器建立方法，其特征在于，采用以下核心求解器：LBGK基础求解器、自由表面模型、SMA亚格子湍流模型，对工业CAE在应用于压铸分析时，出现的大紊流超高雷诺数内部流动进行模拟计算，具体包括以下步骤：

步骤S1，构建LBGK基础求解器，利用所述LBGK基础求解器首先对空间进行离散，或读入计算网格模型数据，然后对速度和时间进行离散，计算碰撞项，更新密度分布函数、宏观场速度和密度；

步骤S2，在所述LBGK基础求解器的基础上，采用多松弛时间MRT模型对步骤S1计算的数据进行稳定性增强；

步骤S3，采用亚格子SMA湍流模型，对MRT模型的处理MRT模型的松弛矩阵中主对角线元素进一步处理，提高流动求解器在大雷诺数条件下计算稳定性；

步骤S4，采用自由表面模型对所述SMA湍流模型处理后的数据，通过单相流并耦合自由表面模型实现两相流流动过程两相界面演化做出准确的模拟计算；

步骤S5，采用基于消息传递界面机制MPI并行模型对上述步骤S1至步骤S4中的算法实现并行优化。

2.如权利要求1所述的用于工业CAE方向的流动求解器建立方法，其特征在于，在所述步骤S1中，所述LBGK基础求解器对速度用D3Q19模型进行离散；采用显式欧拉方法对时间进行离散。

3.如权利要求1所述的用于工业CAE方向的流动求解器建立方法，其特征在于，在所述步骤S4中，所述自由表面模型采用VOF模型。

4.如权利要求1所述的用于工业CAE方向的流动求解器建立方法，其特征在于，在所述步骤S2中，所述MRT模型针对LBGK基础求解器计算的数据，将密度空间的松弛时间τ转换至矩空间中，然后对矩空间中多个松弛时间分量构成的松弛矩阵参数进行调整数学描述如下所述：

f_i(x+c_iΔt,t+Δt)-f_i(x,t)＝-M^-1SM[f_i(x,t)-f_i ^eq(x,t)]Δt

其中，M为变换矩阵，S为松弛矩阵；

S数学表达式为：S＝diag(0,ω_e,ω_ε,0,ω_q,0,ω_q,ω_v,ω_v,ω_π,ω_v,ω_v,ω_v,ω_m,ω_m,ω_m)。

5.如权利要求1所述的用于工业CAE方向的流动求解器建立方法，其特征在于，在所述步骤S3中，所述SMA模型对于松弛矩阵中主对角线元素做进一步处理，过程如下：

其中

式中，C_s为格子声速，Δx为格子间距，τ₀是本征松弛时间，Q为流体力学Q函数。

Images