CN110032756B - 基于流函数分数坐标系变换的流动边界层数值分析方法 - Google Patents

Abstract

本发明公开了一种基于流函数分数坐标系变换的流动边界层数值分析方法，包括步骤：构建一个二维或三维流体区域作为仿真环境；设置所述仿真环境的初始参数，进行仿真环境的配置和初始化；对所要研究的流动现象进行抽象，建立物理模型，选择相应的控制方程；设置所述仿真环境的形状参数，根据所述形状参数贴体建立空间坐标系；定义所述仿真环境的流函数，根据边界层发生质量变动情况，将之标准化为流函数质量分数，将所述空间坐标系变换为空间‑流函数质量分数坐标系；定义所述仿真环境的求解器，进行控制方程的离散求解并输出至计算机进行可视化分析。本发明降低了离散化和理论模型失真造成的误差，简化了数值分析模型，降低了对计算资源的需求。

Description

基于流函数分数坐标系变换的流动边界层数值分析方法

技术领域

本发明涉及计算流体力学数值分析方法领域，尤其涉及一种基于流函数分数坐标系变换的流动边界层数值分析方法。

背景技术

计算流体力学(Computational Fluid Dynamics，CFD)是20世纪50年代以来，随着计算机的发展而产生的一个介于数学、流体力学和计算机之间的交叉学科，主要研究内容是通过计算机和数值方法来求解流体力学的控制方程，对流体力学问题进行模拟和分析。

CFD最基本的考虑是如何把连续流体在计算机上用离散的方式处理。一个方法是把空间区域离散化成小胞腔，以形成一个立体网格或者格点，然后应用合适的算法来解运动方程。这样的一个网格可以是不规则的或者是规则的，前者的特征是每个胞腔必须单独存储在内存中。如果是高度动态的并且在尺度上跨越很大的范围，网格本身应该可以动态随时间调整。

对于层流情况和对于所有相关的尺度都可以包含在格点中的湍流的情形，通过直接数值模拟求解纳维-斯托克斯方程是可能的。一般情况下，适合于问题的尺度的范围甚至大于今天的大型并行计算机可以建模的范围，需要引入湍流模型。

以二维不可压缩流体恒定流动边界层为例，需要求解的方程为：

连续性方程：

x方向动量方程：

一般标量输运方程：

针对上述方程在(x,y)坐标系进行离散。在使用中的一些离散化方法包括：有限体积法、有限元法(FEM)、有限差分方法等。将离散得到的矩阵(通常为稀疏)进行迭代求解。离散化完成后要进行的方程系统的基本求解由很多数值线性代数的为人熟知的许多算法完成。可以使用静态的迭代方法，譬如对称高斯-赛得尔(Gauss-Seidel)、渐进松弛(successive overrelaxation)或克雷洛夫子空间法(Krylov subspaces)。最后将计算结果输出为约定的可用于后处理可视化分析的数据格式。

传统的CFD方法存在以下缺点：

1.将数学模型在空间上进行离散化处理的时候，会引起误差。

2.理论模型没有完整包含实验中的各种复杂因素，在缺少实验数据修正的情况下会引起误差。

3.需要实现给定计算域，但是求解的问题本身实际上影响区域是在变动的，会引起误差。

4.需要求解完整的流体方程，使得计算资源需求十分巨大。

5.难以求解某些特定的工程问题，例如覆膜冷却等。

因此，本领域的技术人员致力于开发一种基于流函数分数坐标系变换的流动边界层数值分析方法，将物理问题在空间-流函数质量分数坐标系下进行离散，该坐标系是基于流通本身特性的天然坐标系，显著降低离散化造成的误差；在理论模型的基础上，加入了大量实验数据的修正，使得修正过的模型更加准确的反应求解问题的真实情况，显著降低因为理论模型失真造成的误差；通过采用空间-流函数质量分数坐标系，使得所求解问题的计算域在该坐标系下保持恒定不变，但是映射到真实世界视角下的空间中时，计算域会随着所求解问题的发展一起发展，从而显著降低误差；针对求解边界层占主导的区域内的问题，所求解的控制方程为该区域内特有的简化方程，合理地简化了问题，极大地降低了对计算资源的需求；设计通用的、开放的模型接口，有针对性地对特定问题进行订制开发求解。

发明内容

有鉴于现有技术的上述缺陷，本发明所要解决的技术问题是如何减少流动边界层数值分析的误差，简化数值分析模型，降低对计算资源的需求。

为实现上述目的，本发明提供了一种基于流函数分数坐标系变换的流动边界层数值分析方法，所述方法包括以下步骤：

步骤1、构建一个二维或三维流体区域作为仿真环境；

步骤2、设置所述仿真环境的初始参数，进行仿真环境的配置和初始化；

步骤3、对所要研究的流动现象进行抽象，建立物理模型，选择相应的控制方程；

步骤4、设置所述仿真环境的形状参数，根据所述形状参数贴体建立空间坐标系；

步骤5、定义所述仿真环境的流函数，根据边界层发生质量变动情况，将之标准化为流函数质量分数，将所述空间坐标系变换为空间-流函数质量分数坐标系；

步骤6、定义所述仿真环境的求解器，进行控制方程的离散求解；

步骤7、所述仿真环境根据计算域在沿边界表面前进方向的位置判断仿真是否结束，如未结束则会进行下一坐标的计算，同时在每一坐标的计算完成之后，进行当前坐标的部分后处理，并输出至计算机进行可视化分析。

进一步地，所述步骤2中所述仿真环境的初始参数包括流体的属性参数、边界条件、计算域。

进一步地，所述步骤2中设置的方式为手动输入或读取设置参数文件的任意一种。

进一步地，所述步骤2中，当所述仿真环境发出读取设置参数文件指令时，通过读取设置参数文件设置所述仿真环境的初始参数，当所述仿真环境发出写入设置参数文件指令时，将所述仿真环境的初始参数写入设置参数文件。

进一步地，所述步骤4中设置所述仿真环境的形状参数通过读取形状参数文件设置。

进一步地，所述步骤6中定义所述仿真环境的求解器包括针对边界层占主导的区域内的问题，所求解的控制方程根据边界层假设条件进行简化。

进一步地，所述步骤6中定义所述仿真环境的求解器包括在建立理论物理模型的基础上，加入实验数据的修正。

进一步地，所述步骤6中定义所述仿真环境的求解器包括根据速度梯度变化，加入流动方向的自适应网格加密算法。

进一步地，所述步骤6中定义所述仿真环境的求解器包括定义线性代数求解器和定义物理场求解器。

进一步地，所述步骤6中定义物理场求解器求解包括以下步骤：

步骤6.2.1、根据旧一列网格的数据，计算离散方程的系数，求解新一列网格的数据；

步骤6.2.2、将新一列网格的数据更新至中间迭代步的数据，使用旧一列网格的数据和中间迭代步的数据的平均值计算离散方程的系数，计算新的新一列网格的数据。

步骤6.2.3、比较新的新一列网格的数据和中间迭代步的数据的误差，

a)当误差大于给定阈值，且迭代步数小于限定步数，跳至步骤6.2.2；

b)当误差小于给定阈值，或迭代步数大于限定步数，则将旧一列网格的数据用新一列网格的数据更新，进入下一坐标的步骤6.2.1计算。

本发明的有益技术效果为：

1.本发明并非在真实世界视角下的空间进行离散化处理，而是通过做表转化，将物理问题在空间-流函数质量分数坐标系下进行离散，该坐标系是基于流通本身特性的天然坐标系，使用该坐标系，将会显著降低离散化造成的误差。

2.本发明在理论模型的基础上，加入了大量基于国外(NASA等)和国内科研院所提供的实验数据的修正，使得修正过的模型更加准确的反应求解问题的真实情况，显著降低因为理论模型失真造成的误差。

3.本发明通过采用空间-流函数质量分数坐标系，使得所求解问题的计算区域在该坐标系下保持恒定不变，但是映射到真实世界视角下的空间中时，计算区域会随着所求解问题的发展一起发展，从而显著降低误差。

4.本方法针对求解边界层占主导的区域内的问题，所求解的控制方程为该区域内特有的简化方程，合理地简化了问题，极大地降低了对计算资源的需求。

5.本方法设计了通用的，开放的模型接口，可以有针对性地对特定问题进行订制开发求解。

以下将结合附图对本发明的构思、具体结构及产生的技术效果作进一步说明，以充分地了解本发明的目的、特征和效果。

附图说明

图1是本发明的一个较佳实施例的外流分析坐标系描述示意图；

图2是本发明的一个较佳实施例的仿真软件系统操作流程图；

图3是本发明的一个较佳实施例的计算网格示意图；

图4是空间坐标系下的计算域；

图5是空间-流函数质量分数坐标系下的计算域；

图6是一般CFD方法计算域和边界层算法计算域对比。

具体实施方式

以下参考说明书附图介绍本发明的优选实施例，使其技术内容更加清楚和便于理解。本发明可以通过许多不同形式的实施例来得以体现，本发明的保护范围并非仅限于文中提到的实施例。

一、软件系统的基本框架：

首先建立本仿真软件系统(简称为AlphaBLAS)的起点和框架。在软件开始运行时，会首先尝试读取包含在执行命令中的第二个参数，(第一个参数为程序自身，第二个参数为一个路径)，该参数为由用户指定的设置文件的路径。如果该参数存在，则读取该路径下的设置参数文件，否则，读取本地目录下的设置参数文件。接着，AlphaBLAS根据获得的参数，进行仿真环境的配置和初始化。获得配置参数的途径有两条：读取代码中预先设置好的全局参数，或者调用AlphaBLAS的输入模块读取配置参数中用户设置的参数。

在运行的时候，AlphaBLAS程序会根据计算域在x方向的位置判断仿真是否结束，如未结束则会进行下一步的计算。同时，会在每一步的计算完成之后，进行当前步的部分后处理，并输出。

二、软件系统的数据输入输出模块：

AlphaBLAS为了让参数输入和配置更明确，设定了通用的抽象类型，由之派生出了各种子设置类型，这些设置类型的对象，共同作为输入模块的成员，集中管理和使用。在输出的功能实现上，AlphaBLAS可以通过读取内存，输出当前存储着的那一列网格的数据，也可以在硬盘上存储全部的数据，并将之转换成指定的格式，如tecplot等。

三、软件系统的数据结构：

AlphaBLAS中包含3级网格数据，即点网格，线网格和面网格，分别在系统中进行定义。其中，SJBL_线网格的类，包含了由SJBL_点网格组成的数据列，可以直接调用。出于内存使用效率的原因，没有将面网格数据存储于内存中，而是和输出模块_io协作，将数据存储在硬盘中。

四、软件系统的线性代数求解器：

AlphaBLAS定义了系统内部矩阵的数据格式、操作定义。AlphaBLAS求解的离散方程为3对角矩阵形式，故矩阵中包含3列系数(mA，mB，mC)和1列值(mX)。并定义了矩阵之间的加(+)和减(-)操作以及矩阵和系数之间的乘(*)操作。

五、软件系统的物理场求解器：

AlphaBLAS在运行时，在内存中始终同时存储着3条列网格的数据，分别是：

旧一列网格的数据，里面所有的信息已经计算完毕。

新一列网格的数据，里面的信息是最新的计算结果。

中间迭代步的数据，当新的数据计算完成后，将会被更新至其中。

物理场求解时的顺序为：

1)根据旧一列网格的数据，计算离散方程的系数，求解新一列网格的数据；

2)将新一列网格的数据更新至中间迭代步的数据，使用旧一列网格的数据和中间迭代步的数据的平均值计算离散方程的系数，计算新的新一列网格的数据。

3)比较新的新一列网格的数据和中间迭代步的数据的误差，

a)当误差大于给定阈值，且迭代步数小于限定步数，跳至步骤2；

b)当误差小于给定阈值，或迭代步数大于限定步数，则将旧一列网格的数据用新一列网格的数据更新，进入下一步的步骤1计算。

六、求解方程的简化：

在求解外流问题时，所关注的物理特性(壁面应力、热流)均主要取决于边界层的特性，而边界层以外的流场可以视为给定的外部边界条件。可以用边界层模型进行分析求解的问题有很多种，包括：平板流动，管道流动，叶片表面流动等。可以分析2D问题，也可以分析轴对称问题。

对燃气轮机叶片的外流场边界层简化假设为：

如图1所示，其中，x为沿着边界表面前进的距离，y为垂直于固体表面的距离，r为横截面的曲率半径，纵向的曲率半径认为变化可以忽略。使用雷诺应力平均方式(RANS)建立控制方程，在图1所示的坐标系中，基于边界层假设条件的控制方程形式为：

连续性方程：

动量方程：

压强梯度方程：

能量守恒方程：

式中I^*为滞止焓，I为静态焓：

控制方程中的S项为广义源项，μ_eff为有效粘度，由两部分组成，工质自身的粘度μ，和包含了湍流脉动效果的湍流粘度μ_t，依据选择的湍流模型有不同的表达，并可在此基础上定义湍流普朗特数Pr_t和有效普朗特数Pr_eff。

七、坐标系转换：

定义流场的固壁面边界1为I面，边界层上边界2为E面,如图1所示。定义流函数ψ，由于在具体的问题中，I面和E面都有可能会发生质量变动(m"_I和m"_E，如壁面蒸发，覆膜冷却，流动混合等)，因此流函数的取值范围可能会发生偏移，故将其标准化为流函数的质量分数ω，并将(x,y)坐标系的控制方程转写至(x,ω)坐标系上。

令ψ为流函数：

则有：

当不存在由固壁面边界1输入的质量和边界层上边界2质量交换时，随着流动沿x方向发展，边界层将会逐渐变厚，上边界所对应的y将逐渐变大，但是对应的流函数ψ始终保持不变，因此由(x,ψ)坐标系描述的流体控制区域始终保持恒定。当考虑到由固壁面边界1输入的质量(如冷气孔)和边界层上边界2质量交换(如混流)时，ψ发生变化，但是认为该变化为成比例变化的，使用流函数质量分数ω代替ψ：

将坐标系由(x,ψ)转换至(x,ω)

随着流动沿x方向发展，边界层上边界2对应的流函数质量分数ω始终保持不变，由(x,ω)坐标系描述的流体控制区域始终保持恒定。

如图1所示，(x,y)坐标系贴体建立，边界层上边界2为E面，曲率半径r_E，质量流出为m"_E，固壁面边界1为I面，曲率半径r_I，质量流入为m"_I，并考虑湍流效果：

令：

则有

动量方程：

能量方程：

以平板流为例，如图6所示，完整的流场区域3是一般CFD方法求解的区域，边界层区域4为边界层算法所计算的区域，该区域在(x,y)坐标系下逐渐拓展，如图4所示，但是在(x,ω)坐标系下，边界层上边界2(E面)为边界层处流线，ω为定值，故计算域在(x,ω)坐标系保持均匀分布，如图5所示。实际应用时，方程可以根据固体表面轮廓关系，求解问题是否轴对称等进行修改，对于传热，传质等过程采用边界层相似原理处理，湍流模型的修正，以及I面和E面的传质均可作为源项处理。

八、控制方程离散求解

输运方程的一般形式可写为：

其中A，B，C，D均为常数(D中包含源项)，该方程从左至右的四项，可以将其值认为是其网格节点控制容积5内的平均值，如图3所示：

将之各自展开获得差分形式：

该展开在流动方向上为一阶迎风格式，在垂直流动方向为二阶中心差分格式，也可以根据需要改成其它格式。将各差分格式带入控制方程，则控制方程可写成如下的形式并进行求解：

是上游x_u计算得到的参数。

如图2所示，本实施例的仿真软件系统操作流程如下：

首先，选择设置初始参数的方式，可以选择手动输入各项参数，包括流体的属性参数、求解器、用户自定义标量、固壁面边界、上表面开放边界、计算域等，也可以单击“读取设置参数文件”按钮，从设置参数文件中读入各项参数，根据需要对各项参数进行修改。然后，用户可以单击“写入设置参数文件”按钮将修改好的设置参数文件保存在本地路径以便下次直接调用，也可以完成本地路径设置参数文件的生成之后，进一步读取固壁面形状参数文件，单击“开始计算”按钮，开始仿真计算。最后，程序运行时的输出信息会显示在仿真软件系统界面下方。等待AlphaBLAS软件计算完成之后，本地路径文件夹内会生成三个结果文件：分别为seifcheck.cfg文件，提供用户做输入参数的检查；output.dat文件，用于tecplot绘图；facemesh.sjbl文件，提供用户最终的计算结果数据。

本实施例只求解2D边界层区域，对计算资源的需求极小，可以快速得到计算结果。即使是同样求解2D问题，由于CFD方法需要求解整个外部流动区域，STAN算法等效率仍然明显胜出。由于根据流线划分计算网格，使得计算资源的分配更加合理，可以在相关边界层区域内布置较密的计算节点，对于符合边界层方程的算例，精度较高。可以根据给定的参数，自动生成求解多种形状计算域(平板，翼型，管道流动，各种轴对称问题等)，避免划分CFD计算网格带来的一系列问题。在新求解器中采用的各种湍流、传热传质、壁面修正模型等，均采用了大量的经验模型参数，这些参数的选择和所求解的问题紧密相关，需要基于国内外大量同类型研究的数据库制定经验参数的选用方案，包含多种建立在理论研究和大量实践数据基础上的湍流，传质传热，表面粗糙度，覆膜模型和算法改进。在使用新求解器求解边界层方程时，由于外边界逆压梯度的影响，会造成计算难以收敛。为此，在算法中须根据速度梯度变化，进行相应的调整，如加入流动方向的自适应网格加密算法等功能，以改进算法的收敛性。本发明数据结构和算法架构比较简单，易于实现和改进。

以上详细描述了本发明的较佳具体实施例。应当理解，本领域的普通技术无需创造性劳动就可以根据本发明的构思作出诸多修改和变化。因此，凡本技术领域中技术人员依本发明的构思在现有技术的基础上通过逻辑分析、推理或者有限的实验可以得到的技术方案，皆应在由权利要求书所确定的保护范围内。

Claims (10)

1.一种基于流函数分数坐标系变换的流动边界层数值分析方法，其特征在于，所述方法包括以下步骤：

步骤1、构建一个二维或三维流体区域作为仿真环境；

步骤2、设置所述仿真环境的初始参数，进行仿真环境的配置和初始化；

步骤3、对所要研究的流动现象进行抽象，建立物理模型，选择相应的控制方程；

步骤4、设置所述仿真环境的形状参数，根据所述形状参数贴体建立空间坐标系；

步骤5、定义所述仿真环境的流函数，根据边界层发生质量变动情况，将之标准化为流函数质量分数，将所述空间坐标系变换为空间-流函数质量分数坐标系；

步骤6、定义所述仿真环境的求解器，进行控制方程的离散求解；

步骤7、所述仿真环境根据计算域在沿边界表面前进方向的位置判断仿真是否结束，如未结束则会进行下一坐标的计算，同时在每一坐标的计算完成之后，进行当前坐标的部分后处理，并输出至计算机进行可视化分析。

2.如权利要求1所述的基于流函数分数坐标系变换的流动边界层数值分析方法，其特征在于，所述步骤2中所述仿真环境的初始参数包括流体的属性参数、边界条件、计算域。

3.如权利要求1所述的基于流函数分数坐标系变换的流动边界层数值分析方法，其特征在于，所述步骤2中设置的方式为手动输入或读取设置参数文件的任意一种。

4.如权利要求1所述的基于流函数分数坐标系变换的流动边界层数值分析方法，其特征在于，所述步骤2中，当所述仿真环境发出读取设置参数文件指令时，通过读取设置参数文件设置所述仿真环境的初始参数，当所述仿真环境发出写入设置参数文件指令时，将所述仿真环境的初始参数写入设置参数文件。

5.如权利要求1所述的基于流函数分数坐标系变换的流动边界层数值分析方法，其特征在于，所述步骤4中设置所述仿真环境的形状参数通过读取形状参数文件设置。

6.如权利要求1所述的基于流函数分数坐标系变换的流动边界层数值分析方法，其特征在于，所述步骤6中定义所述仿真环境的求解器包括针对边界层占主导的区域内的问题，所求解的控制方程根据边界层假设条件进行简化。

7.如权利要求1所述的基于流函数分数坐标系变换的流动边界层数值分析方法，其特征在于，所述步骤6中定义所述仿真环境的求解器包括在建立理论物理模型的基础上，加入实验数据的修正。

8.如权利要求1所述的基于流函数分数坐标系变换的流动边界层数值分析方法，其特征在于，所述步骤6中定义所述仿真环境的求解器包括根据速度梯度变化，加入流动方向的自适应网格加密算法。

9.如权利要求1所述的基于流函数分数坐标系变换的流动边界层数值分析方法，其特征在于，所述步骤6中定义所述仿真环境的求解器包括定义线性代数求解器和定义物理场求解器。

10.如权利要求9所述的基于流函数分数坐标系变换的流动边界层数值分析方法，其特征在于，所述步骤6中定义物理场求解器求解包括以下步骤：

步骤6.2.1、根据旧一列网格的数据，计算离散方程的系数，求解新一列网格的数据；

步骤6.2.2、将新一列网格的数据更新至中间迭代步的数据，使用旧一列网格的数据和中间迭代步的数据的平均值计算离散方程的系数，计算新的新一列网格的数据；

步骤6.2.3、比较新的新一列网格的数据和中间迭代步的数据的误差，

a)当误差大于给定阈值，且迭代步数小于限定步数，跳至步骤6.2.2；

b)当误差小于给定阈值，或迭代步数大于限定步数，则将旧一列网格的数据用新一列网格的数据更新，进入下一坐标的步骤6.2.1计算。

Images