获取不可压缩流动的流场的数值模拟方法与流程

本公开涉及一种获取不可压缩流动的流场的数值模拟方法。

背景技术

计算流体力学在工程实际中的应用是针对navier–stokes方程组进行离散化求解以模拟各种流动现象。根据密度在流动中的变化程度，一般将流动分为可压缩流动和不可压缩流动。对于液体流动，如广泛的水力学问题，以及马赫数小于0.3的气体流动问题，一般都认为是不可压缩流动。

求解不可压缩流动的主要困难在于如何处理“速度-压力”的耦合问题，为了求解不可压缩流动问题，研究人员发展了多种数值方法，有非原始变量方法和原始变量方法。在非原始变量方法中，涡量-流函数方法应用较多，但该方法只适用于二维不可压缩流动，由于三维流动的流函数无法定义，因而该方法不能扩展到三维实际流动。原始变量方法包括人工压缩性方法、压力poisson方程法以及压力修正法等，其中simple系列的压力修正法应用最为广泛。为了解决流场中“速度-压力”的耦合问题，一般需要使用交错网格，因为对流速的控制容积做积分导出该流速的离散方程时，该流速两侧的节点压力自然地进入离散方程中的压力梯度项，进而可以得到压力修正和速度修正的关系式，再通过连续方程对压力控制容积做积分就得到了压力修正方程。保证了速度与压力之间的偶合关系，实质上保证了连续方程与动量方程的耦合求解。但随着计算流体力学、计算传热学的发展，所计算的问题越发复杂化，特别是在计算三维问题时，这时交错网格的复杂与不便越来越突出。于是，研究人员又发展了在同一套网格上吸收交错网格的成功经验来完成流场与温度场的计算，即对控制体界面通量进行动量插值，界面通量的动量插值方法就是将在节点上离散的动量方程插值到界面上，进而得到界面速度的表达式，此表达式中就含有两相邻节点的压力信息，从而可以得到压力修正和速度修正的关系，进一步得到压力修正方程。采用动量插值法重构界面速度来引入压力修正的过程可以理解为界面动量方程的重构过程。采用动量插值法后，simple系列压力修正方法能够较好的解决压力的非物理振荡，但这种方法对网格质量要求高，而且把非线性动量方程线性化并迭代求解，线化后的动量方程收敛性不好，编程比较复杂，计算量大。因此，当前需要发展收敛性好、编程复杂性小且简单易用的不可压缩流动的数值计算方法。

技术实现要素：

为了解决至少一个上述技术问题，本公开提出了一种获取不可压缩流动的流场的数值模拟方法，包括如下步骤：

s1：读取模拟对象的计算网格信息、流体工质的物性参数和流动初边值条件信息，计算网格参数和时间步长参数；

s2：初始化速度场和压力场；

s3：使用m步显式runge-kutta法中的前m-1步求解不可压缩流体的动量方程，获得模拟对象的虚拟速度场，其中m表示使用的步数且m≥2；

s4：基于上述的虚拟速度场以及在步骤s1中读取的参数和信息，由连续方程推导压力方程，确定压力方程的系数，通过该压力方程确定下一步的压力；

s5：进行m步显式runge-kutta法的最后一步，利用在步骤s4中获得的压力修正速度场，获得下一步的速度；以及

s6：使用收敛准则判断数值求解是否收敛，若收敛，则输出流场数据信息；若未收敛，则重复步骤s3至s5，直至数值求解收敛，输出流场数据信息。

优选地，步骤s3包括使用节点中心型有限体积法离散动量方程。

动量方程具有如下形式：

其中，i,j为网格节点编号，wi,j为速度项，ci,j为对流通量项，di,j为扩散通量项，pi,j为压力项，ωi,j为第i,j个控制体的面积。

可选地，对流通量项通过中心差分格式、迎风格式、hlpa格式或quick格式进行离散；

可选地，扩散通量项和压力项通过中心差分格式进行离散。

在对流通量项采用中心差分格式离散的情况下，扩散通量项中包括自适应人工粘性项。

优选地，m步显式runge-kutta法满足如下公式：

其中，αk表示各步的系数，δt表示每个控制体单元对应的时间步长，m表示使用的步数且m≥2，n表示时间节点编号。

可选地，采用下式计算时间步长：

δt＝λω²

其中，λ表示比例系数。

优选地，压力方程为具有如下五点格式的方程：

appp＝aepe+awpw+anpn+asps+b

其中，pp、pe、pw、pn和ps分别表示控制体中p、e、w、n和s处的压力，ap、ae、aw、an和as分别表示pp、pe、pw、pn和ps的压力系数，b为一常数。

在定常流动的流场的数值模拟中，可以使用当地时间步长法、残值光顺法或多重网格方法加速数值求解的收敛。

在非定常流动的流场的数值模拟中，可以使用双时间步长方法加速数值求解的收敛。

附图说明

附图示出了本公开的示例性实施方式，并与其说明一起用于解释本公开的原理，其中包括了这些附图以提供对本公开的进一步理解，并且附图包括在本说明书中并构成本说明书的一部分。

图1是根据本公开的一个实施方式的获取不可压缩流动的流场的数值模拟方法的流程图。

图2是根据本公开的一个实施方式的节点中心型有限体积法的控制体示意图。

图3是根据本公开的一个实施方式的不可压缩平板湍流边界层流动计算网格示意图。

图4是根据本公开的一个实施方式的不可压缩平板湍流边界层流动速度剖面示意图。

图5是根据本公开的一个实施方式的不可压缩平板湍流边界层内层平均速度分布与理论值的比较。

图6是根据本公开的一个实施方式的不可压缩平板表面摩擦系数分布。

图7是根据本公开的一个实施方式的naca4412翼型绕流计算网格示意图。

图8是根据本公开的一个实施方式的naca4412翼型绕流流线计算图。

图9是根据本公开的一个实施方式的naca4412翼型绕流压力等值线计算分布图。

图10是根据本公开的一个实施方式的naca4412翼型绕流压力系数计算分布图。

图11是根据本公开的一个实施方式的naca4412翼型绕流计算速度剖面图(x/c＝0.675)。

图12是根据本公开的一个实施方式的naca4412翼型绕流计算速度剖面图(x/c＝0.786)。

图13是根据本公开的一个实施方式的naca4412翼型绕流计算速度剖面图(x/c＝0.897)。

具体实施方式

下面结合附图和实施方式对本公开作进一步的详细说明。可以理解的是，此处所描述的具体实施方式仅用于解释相关内容，而非对本公开的限定。另外还需要说明的是，为了便于描述，附图中仅示出了与本公开相关的部分。

需要说明的是，在不冲突的情况下，本公开中的实施方式及实施方式中的特征可以相互组合。下面将参考附图并结合实施方式来详细说明本公开。

图1是根据本公开的一个实施方式的获取不可压缩流动的流场的数值模拟方法的流程图，下面参照该图以二维不可压缩流动为例，详细说明本公开的数值模拟方法。

1、确定描述不可压缩流动的流场的控制方程：

以下以二维不可压缩navier–stokes动量方程组为例来进一步说明本公开的获取不可压缩流动的流场的数值模拟方法，该动量方程组的守恒形式如下式1所示：

其中，q为守恒变量矢量，f、g为无粘通量项，r、s为粘性通量项，▽p为压力梯度项，式1中q、f、g、r、s和▽p的具体含义如下：

其中，ρ为流体密度，p为流体压力，u为x方向的流体速度，v为y方向的流体速度，μ表示动力学粘性系数。

在用上式2的方程组计算层流流动时，动力学粘性系数μ＝μl，其中μl表示层流粘性系数，其可由如下的sutherland公式计算获得：

其式3中，t表示输入温度(k)，t0表示参考温度(k)，ts表示材料的sutherland常数，μ0表示参考温度t0下的参考粘性系数，μ表示输入温度t下的粘性系数。

在用上式2的方程组计算湍流流动时，动力学粘性系数μ＝μl+μt，其中μt表示湍流粘性系数。一般需要采用湍流模型得出湍流粘性系数，例如使用b-l、k-ε、k-ω等湍流模型。

2、对上述动量方程组进行离散：

在本公开的一个优选实施方式中，使用节点中心型有限体积法对上述动量方程组进行离散。在节点中心型有限体积法中，网格节点即为积分单元的中心，在网格节点周围构筑辅助积分单元作为控制体，做一些构筑辅助积分单元的计算，该方法易于给定物面边界条件。图2为节点中心型有限体积法的控制体示意图。在控制体上对式1进行积分，并应用格林公式可以得到下式：

其中，w为速度项，ω为控制体的面积，γ为相应控制体的周线。

上式4可以被半离散化为：

其中，i,j为网格节点编号，wi,j为速度项，ci,j为对流通量项，di,j为扩散通量项，pi,j为压力项，ωi,j为第i,j个控制体的面积。

针对不可压缩流动，可以采用很多格式来离散对流通量项和扩散通量项。ci,j离散后可以写成下式6的形式：

其中，下标k表示控制体的边界代号(1-e，2-n，3-w，4-s)，如图2所示。

ci,j的具体离散方法有多种，可以使用中心差分格式、迎风格式、quick格式、hlpa格式等。同样di,j离散后可以写成下式7的形式：

di,j项和pi,j项可以都采用中心差分格式计算，值得注意的是,在离散扩散项di,j之前，要先求出扩散项中的速度偏导数在控制体界面上的值。

在本公开的一个可选实施方式中，如果对流项采用中心格式计算，为了保证计算稳定，则di,j项中还应该包括自适应人工粘性项。

3、建立压力方程：

由于不可压缩流动中没有独立求解压力的方程，而且压力有可能产生很强的非物理振荡。由动量方程求解出的速度是满足动量方程，但是并不一定能够满足连续方程，压力的影响恰恰是要求速度既满足动量方程，又要满足连续方程，所以压力的影响是可以通过连续方程得出的。

在本公开的一个实施方式中，将通过jameson改进的五步显式runge-kutta法中的前4步求解动量方程，但不限于采用该五步显式runge-kutta法，可根据需要选择采用其他多步显式runge-kutta法，求解得到的速度作为虚拟速度，将该虚拟速度带入到连续方程中，得到关于压力的方程，再迭代求解该压力方程，求出新的压力。

具体地，首先，在连续方程的基础上直接推导求解压力的方程，具体推导过程如下：

由式5略去下标可写成：

其中，为虚拟速度，hⁿ＝dⁿ-cⁿ，n为时间节点编号。

控制体界面上速度的关系式可以表示为下式10：

其中，是用五步显式runge-kutta法经过前4步求解得到的相邻节点上的虚拟速度的平均值。

把连续方程在控制体上积分可以得到下式11：

应用高斯公式进行变换，最终可以得到下式12：

其中，ak为控制体边界面积，为界面单位法向矢量。

把上式10代入到上式12中整理，可以得到一个关于节点压力的五点格式的方程，这样我们就得到了求解压力的方程，即式13：

appp＝aepe+awpw+anpn+asps+b式13

其中，pp、pe、pw、pn和ps分别表示控制体中p、e、w、n和s处的压力，ap、ae、aw、an和as分别表示pp、pe、pw、pn和ps的压力系数，b为一常数。

通过使用五步显式runge-kutta法的前4步求得的虚拟速度以及网格参数等条件，可以确定压力方程系数。通过求解该压力方程，可以得到在虚拟速度场上的新的流体压力。

在基于获得的虚拟速度求得新的流体压力后，再通过jameson的五步显式runge-kutta法中的最后一步的计算，用上述新的流体压力修正流体速度，得到新的速度场。

下面详细介绍该算法的数值求解过程。对于每一个控制体单元，方程式5相当于一个常微分方程。在本公开的一个优选的实施方式中，采用jameson改进的五步显式runge-kutta法求解该常微分方程，则上式5可以写为下式14和式15：

其中，各步的系数分别为α5＝1。

在保证计算稳定的前提下，为加快收敛速度，可以采取当地时间步长法，隐式残值光顺和多重网格法。

在本公开的一种优选的实施方式中，采用下式16计算当地时间步长：

δt＝λω²式16

其中，δt表示每个控制体单元对应的时间步长，ω表示每个控制体单元的面积，λ表示比例系数。在保证计算稳定的前提下，时间步长的取值尽可能的大，采用当地时间步长的方法可以使收敛速度大幅提高。

五步显式runge-kutta法可以在同一计算框架下实现定常流动和非定常流动的数值求解。在定常流动的数值求解中，本公开的数值方法可以使用当地时间步长法、残值光顺法或多重网格方法加速收敛。在非定常流动的数值求解中，本公开的数值方法可以使用双时间步长方法获得非定常不可压缩流场结果。

下面给出采用本公开的方法获得的不可压缩流动的流场的数值模拟结果的两个实施例：

在第一个实施例中，模拟对象为不可压缩平板湍流边界层流动。在湍流边界层中，随着离壁面的距离的变化，粘性切应力和湍流附加切应力各自对流动的影响也发生变化。离开壁面后，粘性切应力的影响逐渐减小，而湍流附加切应力的影响开始不断增大，而后逐渐减小。

首先读入计算网格信息，计算网格为81x51的结构化网格，如图3所示为不可压缩平板湍流边界层流动计算网格示意图。读入无量纲化的粘性系数，即reynolds＝2.2×10⁶；读入设置边界条件选项，壁面采用无滑移条件，出口及上边界速度，压力均外推；对流通量的差分格式采用三阶精度的quick格式；扩散通量采用二阶精度的中心格式；湍流模型使用b-l模型。

其次，进行流场的初始化。

之后，使用上文详细描述的jameson的五步显式runge-kutta法进行计算。

图4为计算收敛后得到的速度剖面示意图。图5是使用对数率整理的平板湍流边界层内层平均速度分布与理论值的比较。其中选取x＝0.3和0.8的两个截面。

图6反映了表面摩擦系数的计算情况，它直接关系到流动的阻力，是工程上最关心的内容，其中τw为壁面剪切应力，ue为边界层外缘速度，一般也可以用来流速度表示。层流部分与湍流部分参考的公式分别是：层流部分湍流部分cf＝0.455/[ln(0.06rex)]²，其中rex是沿x方向的当地雷诺数。该实施例有效说明了本公开提出的不可压缩流动的流场的数值模拟方法的有效性。

在第二个具体实施例中，模拟对象为不可压缩低速大攻角naca4412翼型绕流。

首先，计算程序读入计算网格：349×61的c型结构化网格，naca4412翼型附近的计算网格如图7所示；无量纲翼型弦长取为c＝1.0；计算域的外边界各方向均取20倍翼型弦长；翼型壁面y+的值在0.05～3.0区间；程序读入翼型流动攻角为13.87°，雷诺数为1.5×10⁶；计算域外边界用无反射远场边界条件，对于翼型物面边界给定无滑移条件；湍流模型使用k-ω模式，湍流模式方程如下式17和式18所示：

其中k为湍流脉动动能，ω为比耗散率，表示湍流的时间尺度，或者频率尺度，m∞为来流的马赫数，re为来流雷诺数，α，β为方程的系数。

如图8所示，为计算的流线图。图9是计算的压力等值线分布示意图。图10是naca4412翼型绕流计算压力系数分布示意图。图11～图13分别是在x/c＝0.675、0.786、0.897三个位置上计算的速度分布情况，并与实验值进行了比较，结果显示本公开的数值方法获得了很好的效果。

需要说明的是，采用上述多步显式runge-kutta法不仅可以用于求解上述一维和二维流动问题，还可以应用于求解三维流动问题，且比传统求解方法更简便。

本公开的获取不可压缩流动的流场的数值模拟方法，是一种简单、高效的求解不可压缩流动的数值计算方法，该方法能够有效解决压力的非物理震荡问题，并直接求解非线性动量方程获得好的收敛特性。此外，该数值方法还可以方便地求解非定常流动。

本领域的技术人员应当理解，上述实施方式仅仅是为了清楚地说明本公开，而并非是对本公开的范围进行限定。对于所属领域的技术人员而言，在上述公开的基础上还可以做出其它变化或变型，并且这些变化或变型仍处于本公开的范围内。