CN103823916B - 一种基于多维黎曼解的任意拉格朗日欧拉方法 - Google Patents

Abstract

本发明公开了一种基于二维黎曼解的任意拉格朗日欧拉方法，用以求解二维可压缩流体力学方程组及以之为基础的相关问题。整个算法的流程包括：网格生成与物理量的初始分布；确定下一时刻网格的策略和方式；确定网格单元边界流量的二维黎曼解法器算法；网格单元边界流量的表达形式；获得下一时刻的计算结果的数值方法。新的发明提供了一种确定网格单元边界流量的新的二维黎曼解法器算法，能够克服传统的一维黎曼解法导致的数值激波不稳定、网格易扭曲等现象，改正了现有二维黎曼解方法过于复杂而难以实施的缺点。本发明是一种简洁、健壮和精确的数值算法组件，适合多介质大变形问题的数值模拟，适合现在流行的有限差分、有限体积和有限元方法的形式。

Description

一种基于多维黎曼解的任意拉格朗日欧拉方法

技术领域

本发明属于流体力学数值模拟技术领域，涉及一种基于多维黎曼解的任意拉格朗日欧拉方法，具体地说，涉及一种利用二维黎曼解模拟多介质可压缩流动的任意拉格朗日欧拉方法。

背景技术

多介质大变形问题的数值模拟，是流体力学问题数值模拟中具有挑战性的课题。它需要求解涉及多种介质相互作用、并发生剧烈运动的问题(如，介质撞击和作用的弹塑性问题，爆轰问题，高温高压的辐射流体力学问题)。在这类问题中，介质常常具有复杂的物态性质，状态量(密度、压强)会发生突然的改变(激波、爆轰波出现)，作用区域有可能会发生几百倍的变化。这种问题的计算方法包含了复杂的流程，如网格运动、介质界面的追踪和可压缩流体力学方程组的数值解法。

在工程中目前应用的主要方法有欧拉方法、拉格朗日方法和任意拉格朗日欧拉方法。欧拉方法使用固定网将，实施简便。主要困难一个是难以精确地计算区域变动巨大的问题，另一个是难以正确地描述不同介质处于同一个网格内的状态(特另是当介质的物态性质相差较大、或者存在辐射等多物理作用过程时)。在弹塑性问题和辐射流体力学问题中，人们较少使用这一方法，更多的是采用拉格朗日方法和任意拉格朗日欧拉方法。拉格朗日方法使用的网格随着介质的运动而运动，可以清晰地分辨介质界面，不需要计算介质在人为混合后的状态。但这种方法非常容易造成网格扭曲，导致计算中断。任意拉格朗日欧拉方法可以在计算过程中调节网格的运动方式，例如仅在介质界面处让网格随着介质运动，而在其它地方则调整网格，使之保持较好的几何品质。该方法在理论上同时兼具欧拉与拉格朗日方法的优点，是一种较理想的解决方案。但在实际的计算过程中，计算结果还是会出现严重的错误，如虚假的涡度误差和数值激波不稳定现象。这些错误扭曲了真实的物理图像，并常常导致计算过程的中断。计算出现错误的主要原因是流体力学方程组的解法器出现了问题。

流体力学方程组的解法器是一种利用当前时刻的物理量，获得下一时刻的物理量的算法实施组件。其关键之处是设计当流体通过网格边界时数值通量(流量)的算法。目前在工程应用中广泛使用的流体力学方程组解法器，都是沿网格边界法线方向，求解一维的黎曼问题，并以获得的结果作为流量的近似值使用的。这种方法对一般流体力学问题是适用的，但在多介质大变形等极端条件下，由于缺少多维信息，对多维效应引起的涡度误差和数值激波不稳定性等难点问题难以解决。依靠加密计算网格，减少时间步长，高精度技术等传统方法并不能改正这些缺陷。有时候加密网格还导致计算终止得更快。采用包含多维信息的多维黎曼解法器，是解决这类问题的可行方案。

遗憾的是，针对二维黎曼问题的解法器还没有在任何这类问题的应用软件中实施。这是由于二维黎曼问题的精确求解理论上还没有完全解决，即使使用近似解法器方法，实施过程也非常复杂。现有的研究成果，基本都停留在学术层面，大部分仅仅应用在四边形网格的欧拉方法中。法国的Despres和Maire等人，分别设计了一种可以在非结构网格(六边形或者三角形)上使用的二维近似黎曼解法器，但仅仅适用于拉格朗日方法。现有的任意拉格朗日欧拉方法在网格边界的通量计算中，还没有实施简单、适用于复杂网格体系、健壮有效的二维黎曼解法器。

发明内容

为了克服现有技术中的缺陷，本发明提供了一种利用二维黎曼解模拟多介质可压缩流动的任意拉格朗日欧拉方法，针对复杂的流体力学计算问题，研制的方法能适应复杂的网格体系，具备分辨多介质问题的能力，并取得良好的计算效果。与传统方法使用局部一维的黎曼解方法不同的是，本发明设计了一个两维的黎曼问题解法器，利用多维信息改正传统方法的不足。本发明是一种健壮、简洁的数值算法，可以明显地减少传统方法的误差。

技术方案

一种基于多维黎曼解的任意拉格朗日欧拉方法，包括以下步骤：

(一)网格生成与物理量的初始分布

所有的初始物理量(密度、速度、压强、能量)均定义在网格的中心。

(二)下一时刻网格的确定和生成

针对不同的网格策略(欧拉方法，自适应网格方法和拉格朗日方法)，求解特定的网格方程，并形成软件程序，主要包括椭圆型网格生成器，自适应网格生成器等。

(三)二维黎曼解法器组件

传统方法缺陷：工程应用中的一维黎曼解算法，不具备二维黎曼解的性质。文献中已有的二维黎曼解算法需要利用前一个时刻的物理量，求解复杂的波系相互作用情况。由于理论不完备，算法的实施非常困难。

新的发明想法：如果二维黎曼解法器仅仅提供经过网格边界的流量算法，而不计算复杂波系的相互作用，就可以大大减少问题的难度。

新的发明遇到的困难：如何从前一个时刻网格中心的物理量，获得直到下一个时刻的网格节点和网格边界的物理量。如果采用直接的插值等数值技术，由于物理量存在间断，将会完全失效。

新的发明主要内容包括：节点运动速度的确定方法和网格边界通量的计算方法。

1)节点运动速度的确定方法

在某一网格V_c的所有边界上(看图1(a))，求解沿着边界qq⁺法方向的局部一维黎曼问题，获得网格边界上的运动速度。当围绕所要求解的节点q(看图1(b))的所有边界上的边运动速度全部求出时，利用这些边速度的法方向分量，采用最小二乘方法求出网格节点的运动速度其具体规则是节点的运动速度在网格边的投影与这些边速度的法方向分量的差，在最小二乘意义下最小。计算过程使用的权函数，需要特殊的设计。

2)网格边界通量的计算方法

利用节点处的二维流体运动速度将此速度在网格边界qq⁺处的投影，更新为相应边界的新的流动速度，再以两边的状态和此流动速度，设计沿qq⁺法向的新的一维黎曼解。此时由于网格节点处的运动速度，进入了边界通量的计算中，传统一维黎曼解中的间断跳跃条件不再成立。一种单边的新的跳跃条件被构造出来，一个边界上具有了两个流量(看图2(b))，分别对应不同的网格单元。这种方法能保证网格的运动方式与边界通量完全相容。沿边界的切向动量，也由于引进了节点速度的平均作用，导致计算更加稳定。

(四)边界通量表达形式

给定了两层网格节点的移动速度后，根据黎曼解的结构，设计出边界通量的表达形式。

(五)获得下一时刻的计算结果

如果计算时间已经到了需要的时刻，则计算终止，否则回到第(二)步。

本发明的有益效果：

本发明适合复杂的网格，特别是无结构网格及结构和非结构混合网格；清晰地分辨介质界面；明显地改善现有流体力学算法的健壮性，提高计算精度。

附图说明

图1是网格及其节点的几何示意图，其中图1(a)是网格单元V_c及相邻单元，图1(b)是网格节点及相邻单元；

图2是构造边界流量的算法示意图，其中图2(a)是传统算法，图2(b)是新提出的边界流量算法；

图3是新发明的局部一维黎曼解的波系分解图，其中图3(a)是传统算法，图3(b)是新提出的边界通量；

图4是激波衍射问题密度轮廓图，其中图4(a)是采用传统的HLLC一维近似黎曼解算法，图4(b)是采用新的二维黎曼解算法；

图5是初始网格一致倾斜时(Saltzmann网格)，其中图5(a)是采用传统的HLLC一维近似黎曼解算法，图5(b)是采用新的二维黎曼解算法；

图6是初始网格为三角形网格的双马赫问题，其中图6(a)是采用传统的HLLC一维近似黎曼解算法，图6(b)是采用新的二维黎曼解算法；

图7是初始网格为三角形网格的双马赫问题，其中图7(a)是为密度轮廓，图7(b)是为网格图；

图8是初始网格为三角形网格的多介质激波管问题，图8(a)为初始网格，图8(b)为计算过程中的拉格朗日网格，图8(c)为计算过程中的自适应移动网格；图8(d)是沿着横向的密度分布；图8(e)是沿着横向的压强分布。

具体实施方式

下面结合附图和具体实施方式对本发明的技术方案作进一步详细地说明。

使用一种基于多维黎曼解的任意拉格朗日欧拉方法，求解如下的流体力学方程组

$\frac{\∂U}{\∂t}+\frac{\∂F\left(U\right)}{\∂x}+\frac{\∂G\left(U\right)}{\∂y}=H$

其中H表示多物理过程的源项，U，F，G分别是

$U=\left[\begin{array}{c}\mathrm{\ρ}\\ \mathrm{\ρ}u\\ \mathrm{\ρ}v\\ \mathrm{\ρ}E\end{array}\right],F=\left[\begin{array}{c}\mathrm{\ρ}u\\ \mathrm{\ρ}{u}^2+p\\ \mathrm{\ρ}uv\\ \mathrm{\ρ}Eu+pu\end{array}\right],G=\left[\begin{array}{c}\mathrm{\ρ}v\\ \mathrm{\ρ}uv\\ \mathrm{\ρ}{v}^2+p\\ \mathrm{\ρ}Ev+pv\end{array}\right],$

这里ρ，u＝(u，v)，p，E分别是流体运动中的密度、速度、压强和总能量。

整个数值模拟过程包括以下步骤：

(一)网格生成与物理量的初始分布

在图1(a)显示的网格单元_Vc上，物理量的分布(密度压强速度内能和总能量)都定义在网格的中心，而网格节点位置q，和网格的几何信息(边长，面积，网格边的切方向和法方向)都已经确定。这里符号n表明时刻tⁿ，此时n＝0。

(二)网格节点处运动速度的确定

在图1(b)显示的节点邻域图上，环绕节点q的网格V_c，内的物理量，如密度、压强，速度等都给定为常数。

先计算沿着边界qq^-的法方向的局部一维黎曼解，其法向速度是

${\mathrm{\ν}}^{*}=\frac{p_L-p_R-{\left(\mathrm{\ρ}u\right)}_L(S_L-u_L)+{\left(\mathrm{\ρ}u\right)}_R(S_R-u_R)}{{\left(\mathrm{\ρ}u\right)}_L-{\left(\mathrm{\ρ}u\right)}_R+S_R{\mathrm{\ρ}}_R-S_L{\mathrm{\ρ}}_L},$

这里L，R分别指的是网格V_c和分别是网格内部的速度在网格边界的投影，S_L，S_R是相应的符号速度。

网格节点q处的流体运动速度应满足如下的最小二乘问题

$Func\left(u_q^{*}\right)=\underset{k\∈K\left(q\right)}{\mathrm{\Σ}}{L}_k{\mathrm{\α}}_k{(u_q^{*}\·N_k-{\mathrm{\ν}}_k^{*})}^2,$

这里K(q)是围绕着节点q的所有边的总数，L_k是第k条边的边长，a_k是相应的权函数，为了满足相应的性质，如质量、动量和能量守恒等性质而需要特定的设计。

(三)下一时刻网格的确定和生成

针对特定的网格策略：欧拉方法，自适应网格方法和拉格朗日方法，选取不同的方法。

如果采用欧拉方法，下一时刻的网格位置不变，节点坐标

如果采用自适应网格方法，需要求解特定的椭圆类型网格方程

$\▿\·(W\▿x)=0,$

这里是梯度算子，W是与密度或者压强相关权函数。这种算法形成的tⁿ⁺¹时刻的网格节点位置能在物理量变化大(梯度大)的地方自动汇集。现在已经形成软件；

如果采用拉格朗日方法，下一时刻的网格位置是这里由第二步计算给出。

至此，网格节点q的移动速度(不是其附近流体的运动速度)为：

(四)构造边界流量

如图2所示，构造边界流量的算法。与传统算法图2(a)仅仅利用边界两侧物理量不同的是，新提出的边界流量算法图2(b)利用了节点处的二维流体运动速度将此速度在网格边界qq⁺处的投影，作为相应边界更新后的流动速度。再以两边的状态和此流动速度，构造沿qq⁺法向的新的一维黎曼解。

由于网格节点处的运动速度，进入了边界通量的计算中，接触处的压强不再相同。看图3的一维波系分解图。

$p_L^{*}=p_L+{\mathrm{\ρ}}_L(S_L-u_L)(S_{*}-u_L),$

$p_R^{*}=p_R+{\mathrm{\ρ}}_R(S_R-u_R)(S_{*}-u_R).$

图3(a)是传统的方法，图3(b)是新发明的方法。新发明最主要的贡献之一是在每个波系k中设计了两个通量F_L，k，F_R，k，k＝1，2，3，使之满足单边的间断条件，

$F_{L,2}-F_{L,1}=S_L(U_L^{*}-U_L)=F_{R,2}-F_{R,1},$

$F_{L,3}-F_{L,2}=S_{*}(U_R^{*}-U_L^{*})=F_{R,3}-F_{R,2},$

$F_{L,4}-F_{L,3}=S_R(U_R-U_R^{*})=F_{R,4}-F_{R,3}.$

这两个通量分别对应两个网格V_c，V_f在其相邻边界上的通量，

$F_f^c(w_q^{*},U_L,U_R)=F_{L,k},F_c^f(w_q^{*},U_L,U_R)=F_{R,k},$

其中，UL，，UR分别是图3四个波系中的守恒物理量向量。

这种方法能保证网格的运动方式与边界通量完全相容(传统的方法没有这种性质)。此外，沿边界的切向动量，也引进了节点速度的作用，因此增加了更多的数值粘性，导致计算更加稳定。

(五)边界通量表达形式

在给定了两层网格节点的移动速度后，根据黎曼解的结构，设计边界通量的表达形式。

(六)获得下一时刻的计算结果

此时可以获得密度压强速度内能和总能量如果计算时间已经到了需要的时刻，则计算终止，否则将第n+1时刻的计算结果作为计算初值，回到第二步。

传统方法与新方法计算结果的比对

激波衍射(后台阶流动问题)：参照图4激波衍射问题密度轮廓图。激波自左向右传播，经过一个台阶，形成激波衍射形式。此例网格在计算过程中固定不动。左采用传统的HLLC一维近似黎曼解方法，右采用新的二维黎曼解方法。传统方法的激波会由于不稳定受到破坏。

活塞推动冷气体问题：如图5所示，初始网格一致倾斜时(Saltzmann网格)，由活塞驱动的激波自左向右传播的密度图。此例网格在计算过程中接近拉格朗日运动。(a)采用传统的HLLC一维近似黎曼解方法，(b)采用新的二维黎曼解方法。传统方法的伪涡度误差非常严重，而新的方法，保持了较好的图像。

双马赫问题：参照图6初始网格为三角形网格的双马赫问题。此例网格在计算过程中固定不动，(a)采用传统的HLLC一维近似黎曼解方法，(b)采用新的二维黎曼解方法。传统方法存在数值激波的不稳定性。

如图7初始网格为三角形网格的双马赫问题。此例网格在计算过程采用自适应网格运动方式。(a)为密度轮廓，(b)为网格图。较少的网格达到与密集网格同样的效果。

多介质强稀疏波问题(拉格朗日方法与任意拉格朗日欧拉方法的比较)：如图8初始网格为三角形网格的多介质激波管问题(Abgrall提出)。此例网格在计算过程采用自适应网格运动方式。(a)为初始网格，(b)为计算过程中的拉格朗日网格，(c)为计算过程中的自适应移动网格。(d)和(e)在y＝0.02处密度与压强计算结果的横切面。由于此问题包含了非常强的稀疏波，拉格朗日方法会出现震荡，而自适应网格的新方法表现良好。

以上所述，仅为本发明较佳的具体实施方式，本发明的保护范围不限于此，任何熟悉本技术领域的技术人员在本发明披露的技术范围内，可显而易见地得到的技术方案的简单变化或等效替换均落入本发明的保护范围内。

Claims (1)

1.一种基于多维黎曼解的任意拉格朗日欧拉方法，其特征在于，包括以下步骤：

1)网格生成与物理量的初始分布；

2)下一时刻网格的确定和生成：需要求解与求解特定的网格方程，并形成软件程序，主要包括椭圆型网格生成器，自适应网格生成器；

3)黎曼解法器组件：

构造一个全新的两维近似黎曼解法，该解法器包括节点速度的确定方法和网格边界通量的计算方法，在任意网格下，网格节点的速度算法：已知网格单元中心的物理量，需要求解2×2的代数方程组来获得节点处的流动速度，点处的流动速度与网格边界通过一些特定方法获得的法向速度差在最小二乘意义下最小；

构造边界流量的方法：利用了节点处的二维流体速度，以此速度在网格边界qq⁺处的投影，作为相应边界的流动速度，再以两边的状态和此流动速度，构造沿qq⁺法向的一维黎曼解，由于网格节点处的运动速度，进入了边界通量的计算中，这种方法能保证网格的运动方式与边界通量完全相容 ，沿边界的切向动量，也引进了节点速度的作用；

4)边界通量表达形式：给定了两层网格节点的移动速度后，根据黎曼解的结构，写出边界通量的表达形式；

5)获得下一时刻的计算结果：如果计算时间已经到了需要的时刻，则计算终止，否则回到第2)步。