一种适用于五次样条重构格式的自适应人工粘性控制方法

技术领域

本发明涉及计算流体力学中数值计算方法领域，特别涉及一种适用于五次样条重构格式的自适应人工粘性控制方法。

背景技术

在计算流体力学中，数值格式的耗散起着非常重要的作用。数值耗散越大，计算越稳定，但对流场的刻画能力越弱；数值格式的耗散越小，流场的刻画能力就越强，计算越不稳定。因此如何调节格式的数值耗散一直是人们研究的重点。传统调节耗散的方法是调节格式的总体耗散，这种调节耗散的方式并不能完全满足人们的需求。理想的数值格式应当满足以下条件：在流动光滑区域数值耗散尽量小以保证格式精度和分辨率，在流动非光滑区或者梯度较大的区域数值耗散要足够大以保证计算稳定。因此发展根据局部流场信息确定格式耗散大小的自适应耗散方法具有广泛的科研以及应用前景。

发明内容

针对现有技术中存在的问题，提供了一种适用于五次样条重构格式的自适应人工粘性方法；该方法将流场的局部流动特征同格式的色散特性结合起来，由两者共同决定局部流场需要的人工粘性系数的大小，实现人工粘性系数的自适应改变。首先根据波数识别方法，由表征局部流场特征的密度的高阶导数求出该处流场的局部等效无量纲波数，再将等效无量纲波数与格式的修正无量纲波数进行对比。如果等效无量纲波数在格式的完全解析范围内，则人工粘性系数为零或者取一个很小的值；如果等效无量纲波数不在格式的完全解析范围内，则人工粘性系数根据一定的规则逐渐增大，直至最大值。这个判定过程即为自适应人工粘性判据。这种自适应人工粘性方法改变了以往的全流场使用统一的人工粘性系数且人工粘性系数的大小根据经验确定的方式，使得人工粘性系数的取值更为合理。

本发明采用的技术方案如下：一种适用于五次样条重构格式的自适应人工粘性控制方法，包括以下步骤：

步骤1、读取流场数据，求解五次样条重构方程组得到单元界面处物理量的状态值及其一阶导数；

步骤2、根据单元界面处物理量的状态值及一阶导数，计算出单元界面处的二阶至五阶导数；

步骤3、根据波数识别方法，由单元界面处流场密度的高阶导数求出单元界面处流场的等效无量纲波数，代入自适应人工粘性判据确定该处流场的人工粘性系数；

步骤4、根据单元界面处物理量的状态值计算无粘数值通量，并根据人工粘性系数添加六阶人工粘性项，最后采用相对应的时间离散方法在时间上进行推进。

进一步的，在结构网格有限体积方法的基本框架下，不计质量力和源项，在直坐标系下欧拉方程为：

其中，

其中，

其中

将欧拉方程在三维结构网格的控制体单元

其中

每个面上

其中

在求解过程中根据单元的平均值重构出单元界面处的物理量的状态值，五次样条重构方程组为：

其中，

下标

进一步的，所述步骤2具体包括，根据单元界面

从而得到单元界面

进一步的，波数识别方法为：

其中

进一步的，自适应人工粘性判据为：

其中，

进一步的，所述步骤4具体包括：将单元界面处物理量的状态值带入控制方程计算无粘数值通量，无粘通量计算过程为：

其中，

其中，

其中，

其中

与现有技术相比，采用上述技术方案的有益效果为：本发明的方案中人工粘性系数不再是全场统一的值，而是由自适应人工粘性判据根据局部流场的特征来确定,自适应人工粘性判据局部瞬时的流场特征和格式的色散特性共同来决定，此种方法将局部瞬时流场特征与格式的解析能力结合起来，通过判断局部瞬时的流场特征是否在格式的可解析范围内，从而确定格式的人工粘性系数的取值，使格式的人工粘性控制更为合理。

附图说明

图1是本发明的适用于五次样条重构格式的自适应人工粘性方法流程图。

图2是本发明中的等效无量纲波数与修正无量纲波数对比图。

图3是本发明中自适应人工粘性系数随波数变化曲线图。

图4是本发明中与自适应人工粘性系数对应的数值格式耗散曲线。

图5是采用一种固定人工粘性的二维黎曼问题密度等值线图。

图6是采用自适应人工粘性的二维黎曼问题密度等值线图。

图7是采用另一固定人工粘性的二维黎曼问题密度等值线图。

具体实施方式

下面结合附图对本发明做进一步描述。

本发明提出了一种适用于五次样条重构格式的自适应人工粘性方法，包括下述理论依据：

（1）波数识别方法

考虑傅里叶模态

由式 (1)可得：

将式（2）进行无量纲化，为防止导数为零造成分母为零，等效无量纲波数可写出如下形式：

其中

（2）波数识别在五次样条重构格式中的验证

对于五次样条重构格式，采用傅里叶分析，单元界面

上式中

而根据傅里叶分析得到的五次样条格式的修正无量纲波数为：

其中实部为格式的色散特性，虚部为格式的耗散特性，

图2显示了由高阶导数求出的等效无量纲波数与五次样条格式修正无量纲波数对比图，其中纵坐标表示色散特性，由图中可以看出在修正无量纲波数的可识别范围内，等效无量纲波数与修正无量纲波数符合的很好，说明由波数识别方法得到的等效无量纲波数在一定的波数范围内可以用来代表格式的修正无量纲波数。

（3）自适应人工粘性判据

由于等效无量纲波数在一定波数范围内可以代表格式的修正无量纲波数，因此可 以根据图2进行自适应人工粘性判据的设定：根据修正的无量纲波数的特性设一阈值

图3为人工粘性系数随波数的变化曲线，可以看出随着波数的增大，人工粘性随之变化。图4为自适应人工粘性系数对应的耗散曲线，其中纵坐标表示格式的耗散特性，可以看出自适应耗散在低波数区域耗散尽量小，不引入耗散误差；在高波数区域尽量大以抑制色散误差带来的波动。

基于上述理论依据，具体方案如下：

如图1所示，一种适用于五次样条重构格式的自适应人工粘性控制方法，包括以下步骤：

步骤1、读取流场数据，求解五次样条重构方程组得到单元界面处物理量的状态值及其一阶导数；

步骤2、根据单元界面处物理量的状态值及一阶导数，计算出单元界面处的二阶至五阶导数；

步骤3、根据波数识别方法，由单元界面处流场密度的高阶导数求出单元界面处流场的等效无量纲波数，代入自适应人工粘性判据确定此处流场的人工粘性系数；

步骤4、根据单元界面处物理量的状态值计算无粘数值通量，并根据人工粘性系数添加六阶人工粘性项，最后采用相对应的时间离散方法在时间上进行推进。

具体的，

步骤（1），

在结构网格有限体积方法的基本框架下，不计质量力和源项，在直坐标系下欧拉方程为：

其中，

其中，

其中

将欧拉方程在三维结构网格的控制体单元

其中

每个面上

其中

在求解过程中根据单元的平均值重构出单元界面处的物理量的状态值，五次样条重构方程组为：

其中，

下标

步骤（2）根据单元界面

从而得到单元界面

步骤（3）：式(18)得到单元界面处物理量的二阶至五阶导数后，由式(5)可得到流场界面处的等效无量纲波数，此处求解过程中使用密度的高阶导数。根据自适应人工粘性判据式(7)，确定界面处的人工粘性系数。

步骤（4）：根据控制方程求解无粘数值通量

将单元界面处物理量的状态值带入控制方程计算无粘数值通量，无粘通量计算过程为：

其中，

其中，

其中，

其中

最后求出通量以后，采用相对应的时间离散方法在时间上进行推进。

在本实施例中，给出以下自适应人工粘性的数值验证，如图5、6、7为采用固定人工 粘性的二维黎曼问题的密度等值线图和使用自适应人工粘性的二维黎曼问题密度等值线 图，由图中可以看出：图5

本发明并不局限于前述的具体实施方式。本发明扩展到任何在本说明书中披露的新特征或任何新的组合，以及披露的任一新的方法或过程的步骤或任何新的组合。如果本领域技术人员，在不脱离本发明的精神所做的非实质性改变或改进，都应该属于本发明权利要求保护的范围。

本说明书中公开的所有特征，或公开的所有方法或过程中的步骤，除了互相排斥的特征和/或步骤以外，均可以以任何方式组合。

本说明书中公开的任一特征，除非特别叙述，均可被其他等效或具有类似目的的替代特征加以替换。即，除非特别叙述，每个特征只是一系列等效或类似特征中的一个例子而已。