一种适用于复杂流动的高效高精度数值模拟方法与流程

本发明公开了一种适用于复杂流动的高效高精度数值模拟方法，涉及计算流体力学领域。

背景技术：

在航空宇航领域，计算流体力学(cfd)已经成为了一种不可或缺的技术手段，为模拟真实流动现象和降低研究成本提供了有力的支持。相比于传统的风洞实验，cfd技术具有更低的计算成本和更高的飞行大气环境仿真度，从而以较小的计算代价获得较为精确的气动特性，为飞行器设计提供快速、准确的指导。随着飞行器功能需求和任务剖面的日益复杂多样，低耗散性、高分辨率特性以及高精度特性在流动模拟中的需求日益增加，人们对cfd也提出了更高的要求。

cfd最为关键的技术有三种：网格技术，数值离散方法和物理模型。这三个方面中任何一项的进步都能有效的推动cfd的发展以及在工程领域的应用。其中，通量格式和重构格式均属于空间离散格式，对cfd的求解精度和计算效率影响很大。在目前广泛应用于航空航天飞行器设计的较为成熟的cfd软件中，通量格式大多采用一维黎曼求解器，重构格式大多采用二阶格式。其中一维黎曼求解器在计算网格界面的数值通量时只考虑界面法向波系的传播而忽略横向波的影响，无法描述横向波传递至区域边界的流动特征，从而导致多维流动计算中波系结构分辨率降低以及计算允许的cfl数减小。而二阶重构格式虽说基本能够满足对于全机升阻力分析，气动外形设计优化方面的精度需求，但是对于湍流，分离等多尺度流动现象，仍难以给出满意的结果。

balsara通过推导界面角点处多维黎曼问题的求解公式，提出了一种基于角点框架模式的真正多维黎曼求解器。该求解器通过计算单元界面角点处的数值通量来描述横向波传递的流动特征，从而体现流动的多维效应，提高了波系分辨率以及计算允许的cfl数。该方法具备简单的封闭形式，计算效果良好且算法易于实现。但是，相比于一维黎曼求解器，多维黎曼求解器存在单步求解耗时长、效率低的问题。此外，目前针对多维黎曼求解器高阶格式的研究很少，由于角点处的数值通量求解需要用到角点四周的重构物理量值，传统的适用于结构网格的高阶重构方法无法直接对角点四周的物理量进行重构，从而成为多维求解器向高阶格式推广的主要困难。

技术实现要素：

为进一步提高针对复杂流动的数值模拟精度，本发明提供一种适用于复杂流动的高精度数值模拟方法，以期为更加精准的飞行器设计(例如高性能翼型设计)工作提供一定的技术支撑。本发明主要针对传统cfd求解器计算精度低，稳定cfl数小的弊端，采用多维黎曼求解器以及与之对应的高阶重构方案对流动控制方程进行离散并模拟；同时，针对多维黎曼求解器单步求解耗时长、效率低的问题，采用间断探测技术以提高重构效率、减少计算耗时，最终完成多维复杂流动的高效高精度数值模拟。

为实现上述目的，本发明采用的技术方案包括以下步骤：

步骤1：根据设计任务要求，构建设计对象模型，并建立计算网格，获得需要的网格单元信息；

步骤2：在步骤1建立的网格空间中构造空间离散形式的半离散控制方程：

当采用欧拉方程时，其微分形式的控制方程为

其中，

q是守恒形式的流场变量，f和g表示x和y方向的通量，ρ,u,v,p,^e分别表示流体密度、x方向速度，y方向速度，压强和能量；

对通量项进行空间离散得到：

其中i，j是单元节点编号；fi+1/2,j和gi,j+1/2分别为x方向和y方向的界面数值通量；

步骤3、采用hweno重构方案构造重构多项式

首先对控制方程分别沿x和y进行求导，求导后的控制方程转化为：

其中，fx(q，r)＝f′(q)r，gx(q，r)＝g′(q)r，fy(q，s)＝f′(q)s，gy(q，s)＝g′(q)s，r和s分别是变量q关于x和y的导数，

选择含有v0～v8在内共9个单元的模板，并将其分为8个子模板：

s1＝{v0，v1，v2，v8}，s2＝{v0，v2，v3，v4}，s3＝{v0，v4，v5，v6}，s4＝{v0，v6，v7，v8}

s5＝{v0，v1，v2，v3，v7，v8}，s6＝{v0，v1，v2，v3，v4，v5}，s7＝{v0，v3，v4，v5，v6，v7}，s8＝{v0，v1，v5，v6，v7，v8}

其中v0～v8分别指网格单元vi，j、vi+1，j+1、vi，j+1、vi1，j+1、vi-1，j、vi1，j-1、vi，j-1、vi+1，jj、vi+1，j；

在子模板s1，s2，s3，s4中，关于任意变量q的插值多项式pn需要满足如下约束条件：

其中，

n＝1，k＝0，1，2，8，kx＝8，ky＝2；n＝2，k＝0，2，3，4，kx＝4，ky＝2；

n＝3，k＝0，4，5，6，kx＝4，ky＝6；n＝4，k＝0，6，7，8，kx＝8，ky＝6.

在子模板s5，s6，s7，s8中，关于任意变量q的插值多项式pn需要满足如下约束条件：

其中，

n＝5，k＝0，1，2，3，7，8；n＝6，k＝0，1，2，3，4，5；

n＝7，k＝0，3，4，5，6，7；n＝8，k＝0，1，5，6，7，8.

在每个子模板中，具有三阶精度的插值多项式写为：

pn(x，y)＝a0+a1(x-x0)+a2(y-y0)+a3(x-x0)(y-y0)

+a4(x-x0)²+a5(y-y0)²，n＝1，2，3.4，5，6，7，8

将插值多项式代入约束条件，从而在每个子模板上得到一组关于多项式系数ak(k＝0，1，2，3，4，5)的线性代数方程组，求解该方程组得到各子模板中插值多项式pn的系数ak；

在得到各子模板上的插值多项式后，采用weno限制器的方法，通过光滑指示因子求得9个多项式的权重，加权组合成最终的重构多项式pi，j(x，y)：

光滑指示因子定义如下：

其中|α|＝α1+α2，根据光滑指示因子得到每个多项式的权重如下：

最后通过加权组合的方式得到单元vi,j上的空间重构多项式：

步骤4、根据步骤3得到的最终的重构多项式求解多维黎曼求解器所需的重构状态量；

对于单元vi,j，首先通过最终的重构多项式求得变量在包含界面中点和角点在内的8个插值点(xi+1,j，yi+1/2,j)，(xi,j，yi+1,j)，(xi-1,j，yi,j)，(xi,j，yi-1,j)，(xi+1,j，yi+1/2,j)，(xi-1/2,j，yi+1/2,j)，(xi-1/2,j，yi-1/2,j)，(xi+1/2,j，yi-1/2,j)处的重构值，继而求得界面i+1/2处中点以及上下角点处的多维黎曼求解器所需的状态量以及其中上标“r”和“l”分别表示界面中点两侧的重构变量值，上标“ru”,“lu”,“ld”和“rd”表示角点四周的重构变量值，下标“i+1/2，j”表示界面中点，下标“i+1/2，j+1/2”和“i+1/2，j-1/2”分别表示界面的上下角点；

步骤5、采用多维黎曼求解器进行界面通量求解：

在步骤2的半离散控制方程中，界面通量的具体求解公式如下：

其中，ω1＝1/6,ω2＝4/6,ω3＝1/6为权重系数，分别为x方向和y方向上的辛普森插值点，在x方向上，分别表示了界面i+1/2的上角点、中点和下角点；辛普森插值点处的数值通量为：

其中，通过经典的一维hlle格式求得：

其中，上表“m”表示与界面中点相关的物理量，下标“r”和“l”分别表示界面两侧的重构变量值，由步骤2中求解得到；和分别表示左右传播的最大波速，采用如下公式进行计算：

a是声速，上标“～”表示roe平均；

界面角点处的通量和则通过balsara的真正二维hlle格式求得；

步骤6、根据界面通量求解残差，并将半离散有限体积格式转化为时空全离散有限体积格式，全流场进行时间推进求解，得到最终的流场解。

进一步的，步骤3中采用间断探测器判断模板内是否含有间断，在含有间断的模板中采用weno方法，对各子模板构造得到的插值多项式进行加权平均，得到最终的重构多项式，而在不含间断的模板中，则只求解任一子模板上的插值多项式作为最终的重构多项式。

进一步的，所述间断探测器通过间断探测因子βk判断包含单元vk及其8个邻域单元在内的模板中是否存在间断；间断探测因子βk定义为

其中，(xk，yk)是单元vk的形心坐标，nneighbors是单元vk的相邻单元总数，在笛卡尔网格下，nneighbors＝8，p是重构多项式的阶数，qk和ql是未经限制的解的表达式，是单元均值；

利用以下判断式可以判断单元附近是否存在间断：

其中β0为设定的判断阈值。

进一步的，选取β0＝5。

进一步的，步骤5中，界面角点通量通过以下过程得到：

其中上表“c”表示与界面角点相关的物理量，下标“ru”,“lu”,“ld”和“rd”表示角点四周的重构变量值；其中的波速计算采用如下公式：

其中,表示状态qru处x方向的最大波速；表示状态qru处x方向的最小波速；表示(qlu，qru)之间roe平均状态沿x方向的最大波速；表示(qlu，qru)之间roe平均状态沿x方向的最小波速。

有益效果

本发明的有益效果是，提供一种适用于复杂流动的高精度数值模拟方法,通过采用二维空间模板插值的方式，完成高阶重构多项式的构造，解决了多维黎曼求解器中所需的重构变量无法由传统适用于结构化网格的高阶格式直接求解的弊端，提高波系结构的分辨率以及计算稳定cfl数；并优选通过采用间断探测技术，有效提高了程序的求解效率。本发明能够在解的光滑区域保持一致的时空高阶精度，基本无震荡地完成对流场间断的捕捉并保证流场解的多维特性保持良好。

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

附图说明

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

图1是本发明的实现流程图。

图2是hweno方法的模板示意图。

图3是多维黎曼求解器通量所需的插值点位置示意图。

图4是多维黎曼求解器通量求解的示意图。

图5是实施例中，rae2822翼型的计算网格示意图。

图6是实施例中，采用本发明方案求得的翼型压力系数等值线图。

图7是实施例中，采用本发明方案求得的翼型表面压力分布和风洞试验结果的对比。

具体实施方式

本发明目的是针对多维黎曼求解器单步求解耗时长、效率低的问题，采用间断探测技术以提高重构效率、减少计算耗时，最终完成对例如高性能翼型设计等应用需求下的多维复杂流动的高效高精度数值模拟。

主要包括以下步骤：

步骤1、根据设计任务要求，构建设计对象模型，如本实施例中rae2822翼型，模型，并建立复杂流场计算所需的网格，由给定网格得到需要的网格单元信息，如网格尺度，节点坐标等。

步骤2、构造空间离散形式的半离散控制方程。

在步骤1提供的网格空间上构造空间离散形式的半离散控制方程。以欧拉方程为例，其微分形式的控制方程如下所示：

其中，

q是守恒形式的流场变量，f和g表示x和y方向的通量，ρ,u,v,p,e分别表示流体密度、x方向速度，y方向速度，压强和能量。

对通量项进行空间离散可以得到：

其中i，j是单元节点编号；fi+1/2,j和gi,j+1/2分别为x方向和y方向的界面数值通量，由高阶重构格式和多维黎曼求解器求解得到，具体求解过程详见下述步骤。

步骤3、构造重构多项式

重构多项式的构造是高阶重构方案的关键步骤之一。通过构造空间上的二维插值多项式可以获得物理量在网格单元上的分布，从而求得所需位置处的重构变量值。本发明采用hweno重构方案进行空间二维插值多项式的构造。

为了获得步骤2中界面处的数值通量fi+1/2,j和gi,j+1/2，需要首先获得守恒流场变量在网格单元上的空间分布，进而用于下述步骤5,6中重构变量和界面通量的求解。本发明采用hweno重构方案进行空间二维插值多项式的构造。

hweno重构方案采用hermite插值多项式来构造单元vi,j上具备三阶精度的空间插值多项式。该方法既需要函数值，也需要函数值的导数，从而加强了单元之间的联系，使得格式更为紧凑和稳定。详细构造过程如下：

首先对控制方程分别沿x和y进行求导，求导后的控制方程转化为：

其中，fx(q,r)＝f'(q)r,gx(q,r)＝g'(q)r，fy(q,s)＝f'(q)s,gy(q,s)＝g'(q)s，r和s分别是变量q关于x和y的导数，

如图2所示，选择含有v0～v8在内共9个单元的模板，并将其分为8个子模板:

s1＝{v0，v1，v2，v8}，s2＝{v0，v2，v3，v4}，s3＝{v0，v4，v5，v6}，s4＝{v0，v6，v7，v8}

s5＝{v0，v1，v2，v3，v7，v8}，s6＝{v0，v1，v2，v3，v4，v5}，s7＝{v0，v3，v4，v5，v6，v7}，s8＝{v0，v1，v5，v6，v7，v8}

其中v0～v8分别指网格单元vi，j、vi+1，j+1、vi，j+1、vi-1，j+1、vi-1，j、vi-1，j-1、vi，j-1、vi+1，j-j、vi+1，j。

在子模板s1，s2，s3，s4中，关于任意变量q的插值多项式pn需要满足如下约束条件：

其中，

n＝1，k＝0，1，2，8，kx＝8，ky＝2；n＝2，k＝0，2，3，4，kx＝4，ky＝2；

n＝3，k＝0，4，5，6，kx＝4，ky＝6；n＝4，k＝0，6，7，8，kx＝8，ky＝6.

在子模板s5，s6，s7，s8中，关于任意变量q的插值多项式pn需要满足如下约束条件：

其中，

n＝5，k＝0，1，2，3，7，8；n＝6，k＝0，1，2，3，4，5；

n＝7，k＝0，3，4，5，6，7；n＝8，k＝0，1，5，6，7，8.

在每个子模板中，具有三阶精度的插值多项式可以写为：

将式(7)代入式(5)和(6)，每个子模板上都可得到一组关于多项式系数ak(k＝0，1，2，3，4，5)的线性代数方程组，求解该方程组便可得各子模板中插值多项式pn的系数ak。

在得到各子模板上的插值多项式后，类似tweno重构方案的方法，采用weno限制器的方法，通过光滑指示因子求得9个多项式的权重，加权组合成最终的重构多项式pi，j(x，y)。下面采用weno限制器，根据这五个多项式的光滑程度来赋予其不同的权重。多项式所在的区域越光滑，其所占权重越大，通过该方式能够抑制流场间断附近出现非物理震荡。光滑指示因子定义如下：

其中|α|＝α1+α2，根据光滑指示因子得到每个多项式的权重如下：

最后通过加权组合的方式得到单元vi,j上的空间重构多项式：

步骤4：间断探测器的构造。

步骤3中的weno格式通过自适应选取候选模板以降低间断所在区域对重构多项式的影响，从而避免流场中出现震荡。但在光滑流场区域，各候选模板的权重相近，继续采用weno格式会造成不必要的计算消耗。因此，本实施例通过引入间断探测器来进一步重构效率。该方法的核心思想是，通过间断探测技术判断模板内是否含有间断，在含有间断的模板中采用weno方法，对各子模板构造得到的插值多项式进行加权平均，得到最终的重构多项式p(x,y)，而在不含间断的模板中，只需求解任一子模板上的插值多项式(如p0(x,y))作为最终的重构多项式即可，无需再对其它子模板进行多项式构造，从而大幅降低计算耗时。

下面提出间断探测因子βk定义，通过该因子判断包含单元vk及其8个邻域单元在内的模板中是否存在间断：

其中，(xk，yk)是单元vk的形心坐标，nneighbors是单元vk的相邻单元总数，在笛卡尔网格下，nneighbors＝8，p是重构多项式的阶数，qk和ql是未经限制的解的表达式，是单元均值。

显然，利用以下判断式可以判断单元附近是否存在间断：

随着网格尺度h→0，在光滑区域βk→0，而在间断区域βk→∞，因此采用如下判据来监测间断位置：

在本实施例中，选取β0=5。

步骤5、根据重构多项式求解多维黎曼求解器所需的重构状态量；

下面利用步骤3中构造的空间重构多项式对多维黎曼求解器所需的输入重构变量值进行求解，包括如下两个步骤：

步骤5.1采用步骤3中的间断探测技术对流场间断位置进行探测，在光滑区域构造步骤2中的多项式p0作为最终的重构多项式在间断附近采用tweno或hweno方法构造出加权多项式pi,j作为最终的重构多项式

步骤5.2如图3所示，令m表示单元内的包含界面中点和角点在内的8个插值点编号，以单元vi,j为例，m＝1,2,…,8分别表示点(xi+1,j，yi+1/2,j)，(xi,j，yi+1,j)，(xi-1,j，yi,j)，(xi,j，yi-1,j)，(xi+1,j，yi+1/2,j)，(xi-1/2,j，yi+1/2,j)，(xi-1/2,j，yi-1/2,j)，(xi+1/2,j，yi-1/2,j)。通过可以求得变量在这些插值点处的重构值。

通过求解得到的插值点处的重构值，求得界面i+1/2处中点以及上下角点处的多维黎曼求解器所需的状态量以及其中上标“r”和“l”分别表示界面中点两侧的重构变量值，上标“ru”,“lu”,“ld”和“rd”表示角点四周的重构变量值，下标“i+1/2，j”表示界面中点，下标“i+1/2，j+1/2”和“i+1/2，j-1/2”分别表示界面的上下角点。

例如单元界面两侧的重构变量值为

单元角点c1四周的重构变量值为：

其中，表示单元vi,j内插值点m的坐标。

步骤6、采用多维黎曼求解器进行界面通量求解

图4给出了多维黎曼求解器计算通量所需要的初态，除了界面两侧的两个初态“r”和“l”以外，单元的角点处(点c1～c4)还需要四个初态“ru”,“lu”,“ld”和“rd”，这些初态值步骤5得到。下面给出采用多维黎曼求解器进行界面通量求解的详细步骤。

在半离散控制方程(3)中，界面通量的具体求解公式如下：

其中，ω1＝1/6,ω2＝4/6,ω3＝1/6为权重系数，分别为x方向和y方向上的辛普森插值点，在x方向上，分别表示了界面i+1/2的上角点、中点和下角点。辛普森插值点处的数值通量为：

其中，通过经典的一维hlle格式求得：

其中，上表“m”表示与界面中点相关的物理量，下标“r”和“l”分别表示界面两侧的重构变量值，由步骤2中求解得到。和分别表示左右传播的最大波速，采用如下公式进行计算：

a是声速，上标“～”表示roe平均。

界面角点处的通量则通过balsara的真正二维hlle格式求得，以的求解为例进行说明：

其中上表“c”表示与界面角点相关的物理量，下标“ru”,“lu”,“ld”和“rd”表示角点四周的重构变量值，传统适用于结构化网格的高阶重构方法无法求解这些位置的重构变量，但可由本发明步骤3中提出的重构方案求解得到。这里的波速计算采用如下公式：

其中,表示状态qru处x方向的最大波速；表示状态qru处x方向的最小波速；表示(qlu，qru)之间roe平均状态沿x方向的最大波速；表示(qlu，qru)之间roe平均状态沿x方向的最小波速。

同理可以求得下角点处通量然后通过式(16)加权得到x方向界面通量fi+1/2,j。y方向界面通量gi,j+1/2的求解同上述方法一致。

步骤7、根据界面通量求解残差，并将半离散有限体积格式转化为时空全离散有限体积格式，全流场进行时间推进求解，得到最终的流场解。

步骤6中求得界面通量fi+1/2,j和gi,j+1/2后，由步骤2的式(3)即可求解得到当前第n层时间步的残差l。对时间变量采用三阶runge-kutta离散公式将半离散有限体积格式转化为时空全离散有限体积格式：

l表示残差，在得到界面数值通量后，可由式(3)右端项求解得到，如果是ns方程，则额外计算粘性项的影响。上标“n”表示时间步。利用时空全离散有限体积格式求解下一时间步上的流场变量值。重复以上步骤，依次推进，直到得到全流场稳定的数值模拟结果。

本实施例中是以rae2822翼型跨音速绕流问题进行求解的。图6给出了采用本发明方案求得的压力系数等值线图，能精确捕捉到激波的位置以及激波的强弱。图7给出采用本发明方案求得的翼型表面压力分布和风洞试验结果的对比，可以看出计算结果和试验数据吻合良好。

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