CN102890751A - 求解二维黎曼问题模拟亚音速无粘流的数值方法 - Google Patents

Abstract

本发明属于计算流体力学领域，具体是一种求解二维黎曼问题的数值方法，应用于求解欧拉方程，进行亚音速、无粘流流动的数值模拟。本发明将欧拉平面上的欧拉方程变换到流函数平面，在二维直角网格中，求解跨国流线的黎曼问题和沿着流线的黎曼问题，减小了数值解的误差。

Description

求解二维黎曼问题模拟亚音速无粘流的数值方法

1.技术领域

本发明属于计算流体力学领域，具体是一种求解二维黎曼问题的数值方法，应用于亚音速、无粘流的流动的数值模拟。

2.背景技术

黎曼问题(Riemann Problem)最初在空气动力学领域提出，原来是指一维激波管中的可压缩流体的一种特殊的流动现象。具体是指可压缩流体，如空气，在一维管道中被压缩直至出现激波。激波的出现表明流体的流动出现不连续。跨过激波，流体的速度、密度、压力、熵都出现跳跃变化。1858年，黎曼研究了这种不连续流动现象的特点，提出并解决了一维欧拉方程的一种最简单的初值问题，被后人称为黎曼问题。

图1(a)是黎曼问题的物理模型。图中表示了一维管道(1)中间，有一个隔膜(2)，以此为分界线，管道(1)被分为左、右两个状态。图中的ρ，u，p分别表示流体的密度、压力和速度，下标L和R分别表示左右状态。流体在左右两个状态的参数不同，此时隔膜(2)是流体流动中的不连续。当突然打开隔膜(2)后，在管道(1)中的流体流动会出现以不连续处为中心的、随时间发展变化的波动。德国数学家波恩哈德·黎曼利用该物理模型构造出了描述无粘、可压缩流动的一维欧拉方程在流动不连续处的解析解。如图1(b)，黎曼问题的定义所示。图中x-方向是管道(1)内流体流动的方向，纵坐标t是时间轴。从时间0时刻开始，即隔膜(2)被打开的时刻，在当地形成以激波或者是膨胀波形式存在的左侧变量Q_L或右侧变量Q_R，在它们中间是中间变量，以下标*表示，中间变量又可被接触波分为左中间变量Q_*L和右中间变量Q_*R。黎曼给出了此四类解析解，它们分别由左、右膨胀波和左、右激波组装而成，同时给出了的判别解的条件。这一工作具有极大的超前性和创造性，为非线性双曲型守恒定律的数学理论奠定了基石，开创了计算流体力学(CFD，Computational Fluid Dynamics)领域中一类基于特征的求解流体流动控制方程的数值方法之先河。

计算流体力学的核心任务在于为求解流体流动的控制方程，如欧拉方程，而构造数值方法。其中最重要的是控制方程中的对流项的空间离散，也是最的困难的工作。因为在数值计算中，大多是在笛卡尔坐标下进行的。计算网格必须根据物体的形状事先生成。网格构成网格单元。由于有流体穿过网格单元的交界面，所以存在控制方程的对流项的通量。对对流项的数值逼近会引起数值解的误差。从上世纪以来，CFD研究者致力于开发更精确、高效率的数值方法来降低数值解的误差。迎风差分方法，其典型代表是Godunov方法，在求解流体的流动过程中取得显著成效，因为它合理表达了对流项的特征。

上世纪六十年代出现的Godunov方法的核心技术是在网格单元交界面求解黎曼问题，给出控制方程在网格单元交界面上精确的解。对于多维空间计算中，如二维空间，需要在两个笛卡尔坐标方向逐次求解黎曼问题。然而，对于更一般的情况，如图2(a)所示，网格单元ABCD和ABEF的交界面AB，与速度矢量V不垂直。在交界面处求解黎曼问题时，需要按照图2(b)将V在笛卡尔坐标系(x-y)下的速度分量u，v投影到与AB正交的坐标系(ξ-η)中。在ξ-方向以u₁、v₁之和作为黎曼问题的左变量或右变量。因为黎曼问题在空间上是一维问题，在另一个方向的u₂、v₂之和作为常数处理，跨过激波，这个速度不变。很明显，违反了激波的熵增原理，会造成数值解一定的误差。而激波愈强，引起的误差越大。

为了最大程度减小误差，本发明给出了一个在求解二维欧拉方程时使用的，直接求解二维黎曼问题的数值方法。

3.发明内容

已知的描述二维、无粘、可压缩流体流动的欧拉方程为

$\frac{\∂f}{\∂t}+\frac{\∂F}{\∂x}+\frac{\∂G}{\∂y}=0,---\left(1\right)$

其中，f是守恒变量矢量；F，G是对流项通量在笛卡尔坐标系下在x，y方向上的分量，且

$f=\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{\ρE}+p)u\end{array}\right],$ $G=\left[\begin{array}{c}\mathrm{\ρv}\\ \mathrm{\ρuv}\\ {\mathrm{\ρv}}^2+p\\ (\mathrm{\ρE}+p)v\end{array}\right].---\left(2\right)$

变量t，u，v，ρ和p分别是时间、流动速度V在两个笛卡尔x，y方向上的分量、流体密度和压力。总能E和总焓H分别是

$E=\frac{1}{2}(u^2+v^2)+\frac{1}{\mathrm{\γ}-1}\frac{p}{\mathrm{\ρ}},$ $H=E+\frac{p}{\mathrm{\ρ}}=\frac{u^2+v^2}{2}+\frac{\mathrm{\γ}}{\mathrm{\γ}-1}\frac{p}{\mathrm{\ρ}},---\left(3\right)$

上式中的γ是气体比热比。

欲求解上述欧拉方程(1)，必须获得计算网格单元交界面上的对流项通量F和G的值。本发明的目的是开法一种二维黎曼问题求解器，更精确地获得这个对流项通量的值。数值方法的构造

按照图3所示，在笛卡尔坐标系x-y平面下的一条流线上的某个质点A，它的速度矢量是V。首先，建立另一个坐标系，λ-ξ平面。该平面被称作流函数平面，该平面的横坐标λ，代表流体颗粒运行的距离，平行于速度矢量V，纵坐标是流函数ξ，垂直于V。该平面时间变量τ与t相同。从坐标(t，x，y)到(τ，λ，ξ)的坐标变换的雅各比矩阵可写为

且它的绝对值J成为

其中，θ是流动方向角；

被称作流函数几何状态变量，单位量纲是[m²s·kg^-1]，被定义为

最终，二维不稳定欧拉方程(1)被转换成流函数平面上的形式，

$\frac{{\∂f}_S}{\∂\mathrm{\τ}}+\frac{{\∂F}_S}{\∂\mathrm{\λ}}+\frac{{\∂G}_S}{\∂\mathrm{\ξ}}=0,---\left(7\right)$

其中f_S，F_S和G_S与(2)中的f，F，G称谓相同，只是在流函数平面。具体地，

加上(5)中J的定义，且

$\mathrm{\θ}={\mathrm{tg}}^{-1}\left(\frac{v}{u}\right),---\left(9\right)$

方程(7)被称作成流函数形式的二维不稳定欧拉方程。其中，F_S和G_S向量的第一、四项为0，意味着连续方程和能量方程没有对流项穿过网格交界面，比方程(1)更加简化。

按照公知的偏微分方程的特征分析方法，可以写出方程(7)的特征线方程。在λ方向，

dλ/dτ＝0；         (10)

在ξ方向，

dξ/dτ＝±a/J，dξ/dτ＝0，            (11)

其中，a为声速。沿着特征线(10)和(11)的特征方程分别为

$\mathrm{dp}\±\mathrm{\ρa}\frac{\sqrt{u^2+v^2}}{u}\mathrm{dv}=0.---\left(20\right)$

此外，定义(6)的数值逼近可以写成，

可以获得方程(7)的重要特性，它与方程(1)一样具有旋转不变性，即将除时间坐标以外的坐标方向上的变量互换，方程性质不变。因而特征方程(19)可以重新写做将(21)带入方程(19)

$\mathrm{dp}\±\mathrm{\ρa}\frac{\sqrt{u^2+v^2}}{v}\mathrm{du}=0.---\left(22\right)$

本发明中，沿着ξ方向是与流线垂直的方向。沿着ξ方向求解黎曼问题，被称作跨过流线的黎曼问题。将方程(1)写成仅沿着ξ方向的形式

$\frac{\&PartialD;f}{\&PartialD;\mathrm{\&tau;}}+\frac{\&PartialD;G\left(f\right)}{\&PartialD;\mathrm{\&xi;}}=0,$ $f=\left\{\begin{array}{cc}{f}_L,& \mathrm{\&xi;}<0,\\ f_R,& \mathrm{\&xi;}>0,\end{array}\right.---\left(23\right)$

其中，f_L和f_R分别是网格交界面左侧和右侧变量中的守恒变量。从现在开始，下标L、R、*分别代表黎曼问题中的左、右和中间变量。图4表示了跨过流线的二维黎曼问题的定义。区别于图1(b)中的一维黎曼问题的定义，二维黎曼问题中的左右变量和左右中间变量中的速度包含了两个笛卡尔坐标系下的分量u和v。这两个速度分量同时跨过用特征方程(11)表示的激波、膨胀波和接触波。同时考虑两个速度分量，避免了一维黎曼问题在处理多维流动时会出现的、如前所述的误差。

下面的任务是通过f找到原始变量p，ρ，u，v再求得网格交界面上的G。所述的原始变量也是黎曼问题中的中间变量。一个黎曼问题求解器按照以下流程导出。首先特征方程(20)积分，将左侧变量和右侧变量与中间变量连接，于是有

$\left\{\begin{array}{cc}{p}_{*}+C_L\&CenterDot;f(u_{*},v_{*})=p_L+C_L\&CenterDot;f(u_L,v_L),& \left(24\right)\\ p_{*}-C_R\&CenterDot;f(u_{*},v_{*})=p_R-C_R\&CenterDot;f(u_R,v_R),& \left(25\right)\\ {\mathrm{\&rho;}}_{*L}+f(u_{*},v_{*})({\mathrm{\&rho;}}_L/a_L)={\mathrm{\&rho;}}_L+f(u_L,v_L)({\mathrm{\&rho;}}_L/a_L),& \left(26\right)\\ {\mathrm{\&rho;}}_{*R}-f(u_{*},v_{*})({\mathrm{\&rho;}}_R/a_R)={\mathrm{\&rho;}}_R-f(u_R,v_R)({\mathrm{\&rho;}}_R/a_R),& \left(27\right)\end{array}\right.$

其中C_L＝ρ_La_L，C_R＝ρ_Ra_R，并且

$f(u,v)=\frac{1}{2}\frac{v\sqrt{u^2+v^2}}{u}+\frac{1}{2}u\ln (v+\sqrt{u^2+v^2}).---\left(28\right)$

f(u，v)被称作组合函数，它可以被认为是速度分量(u，v)在流函数平面上的组合，与速度有同样的量纲[m/s]。公式(28)有它的极坐标形式

$f_{\mathrm{\&Delta;}}(V,\mathrm{\&theta;})=\frac{V}{2}(\mathrm{tg\&theta;}+\cos \mathrm{\&theta;}\ln \left((1+\sin \mathrm{\&theta;})\right)+\frac{\cos \mathrm{\&theta;}}{2}V\ln V,---\left(29\right)$

其中

u＝Vcosθ，v＝Vsinθ。方程(24)-(27)的解表达了黎曼问题中的中间变量的值，

为发现速度分量(u_*，v_*)，需要求解方程(31)。首先可将(31)写为

F(u_*，v_*)＝f(u_*，v_*)-B＝0，                     (34)

其中，考虑了已知的左侧和右侧变量，于是有

$B=\frac{C_L\&CenterDot;f(u_L,v_L)+C_R\&CenterDot;f(u_R,v_R)+[p_L-p_R]}{C_L+C_R}.---\left(35\right)$

组合函数(28)可以被进一步改写。定义

tgθ_*＝ε，                   (36)

其中，θ_*是左侧或右侧中间变量中的流动方向角，于是有

$u_{*}=V_{*}\cos {\mathrm{\&theta;}}_{*}=\frac{V_{*}}{\sqrt{1+{\mathrm{\&epsiv;}}^2}};$ $v_{*}=V_{*}\sin {\mathrm{\&theta;}}_{*}=\frac{{\mathrm{\&epsiv;V}}_{*}}{\sqrt{1+{\mathrm{\&epsiv;}}^2}},---\left(37\right)$

其中，V_*是左侧或右侧中间变量中的速度幅值。进而，方程(34)可以最终写成ε的函数

$F\left(\mathrm{\&epsiv;}\right)=\frac{V_{*}}{2}(\mathrm{\&epsiv;}+\frac{\ln V_{*}+\ln (\sqrt{1+{\mathrm{\&epsiv;}}^2}+\mathrm{\&epsiv;})-\ln \sqrt{1+{\mathrm{\&epsiv;}}^2}}{\sqrt{1+{\mathrm{\&epsiv;}}^2}})-B=0.---\left(38\right)$

下面的工作是通过给定的速度V_*求解方程F(ε)＝0，以发现ε，再求出θ_*。方程(38)是可微分的，F的一阶导数是

$F^{\&prime;}\left(\mathrm{\&epsiv;}\right)=\frac{V_{*}}{2}\left(\frac{(2+{\mathrm{\&epsiv;}}^2)\sqrt{1+{\mathrm{\&epsiv;}}^2}-\mathrm{\&epsiv;}-\mathrm{\&epsiv;}(\ln V_{*}+\ln (\sqrt{1+{\mathrm{\&epsiv;}}^2}+\mathrm{\&epsiv;})-\ln \sqrt{1+{\mathrm{\&epsiv;}}^2})}{{(1+{\mathrm{\&epsiv;}}^2)}^{3/2}}\right).---\left(39\right)$

数值试验表明，对于无量纲速度V_*≤5(亚音速区间)而且ε∈[-1，1]的范围内，F′(ε)＞0，这意味着函数F(ε)在此范围内是单调的。同时，F(-ε)·F(ε)＜0，所以Newton-Raphson历遍方法可用来找到方程(45)在给定V_*时的根ε。

为找到V_*的值，需要用p_*，ρ_*L，ρ_*R的值去还原黎曼问题中的各种非线性波(接触波、膨胀波、激波)。本发明中考虑了下面的关系：

(1)跨过激波的Rankine-Hugoniot关系

如果激波出现在一侧(例如，左侧)，跨过激波，Rankine-Hugoniot关系可以用来建立激波和对应的中间变量的关系。它们是，

$\left\{\begin{array}{cc}{\mathrm{\&mu;}}_s({\mathrm{\&rho;}}_L{J}_L-{\mathrm{\&rho;}}_{*L}{J}_{*L})=0,& \left(40\right)\\ {\mathrm{\&mu;}}_s({\mathrm{\&rho;}}_L{J}_L{E}_L-{\mathrm{\&rho;}}_{*L}{J}_{*L}{E}_{*L})=0,& \left(41\right)\end{array}\right.$

其中，μ_s激波速度，J_*L和E_*L分别是中间变量中的J(坐标变换雅各比矩阵绝对值)和E(总能)。从(40)和(41)，可以获得

$V_{*L}=\sqrt{V_L+\frac{1}{\mathrm{\&gamma;}-1}(\frac{p_L}{{\mathrm{\&rho;}}_L}-\frac{p_{*}}{{\mathrm{\&rho;}}_{*L}})}---\left(42\right)$

和

$V_{*R}=\sqrt{V_R+\frac{1}{\mathrm{\&gamma;}-1}(\frac{p_R}{{\mathrm{\&rho;}}_R}-\frac{p_{*}}{{\mathrm{\&rho;}}_{*R}})}.---\left(43\right)$

速度绝对值V_*L或V_*R可以被带入到公式(38)以求解θ_*L或θ_*R。

(2)跨过膨胀波的焓不变关系

如果膨胀波出现在一侧(例如，左侧)，跨过膨胀波，焓不变关系可以用来建立膨胀波和对应的中间变量的关系。它们是，

$H_L=\frac{1}{2}{V}_L+\frac{\mathrm{\&gamma;}}{\mathrm{\&gamma;}-1}\frac{p_L}{{\mathrm{\&rho;}}_L}=\frac{1}{2}{V}_{*L}+\frac{\mathrm{\&gamma;}}{\mathrm{\&gamma;}-1}\frac{p_{*}}{{\mathrm{\&rho;}}_{*L}}=H_{*L},---\left(44\right)$

从这个公式可以求得速度幅值

$V_{*L}=\sqrt{2H_L-\frac{2\mathrm{\&gamma;}}{\mathrm{\&gamma;}-1}\frac{p_{*}}{{\mathrm{\&rho;}}_{*L}}}---\left(45\right)$

或 $V_{*R}=\sqrt{2H_R-\frac{2\mathrm{\&gamma;}}{\mathrm{\&gamma;}-1}\frac{p_{*}}{{\mathrm{\&rho;}}_{*R}}}.---\left(46\right)$

具体那一种非线性波出现，是要靠黎曼问题中的左侧、右侧和中间变量中的压力值判断的。例如，假设

p_min＝min(p_L，p_R)；                     (47)

p_max＝max(p_L，p_R)，                     (48)

及其p_*，构成了下面的波形选择条件：

(1)如果p_min＜p_*＜P_max，意味着在黎曼问题中一个激波出现在p_min一侧，一个膨胀波出现在P_max一侧；

(2)如果p_*≥p_max，意味着在黎曼问题中两个激波出现；

(3)如果p_*≤p_min，意味着在黎曼问题中两个膨胀波出现。

在获得θ_*之后，按照公式(36)求得的u_*L，v_*L或u_*R，v_*R，再加上p_*，便可求得公式(23)中的G_L。时间更新的解可以从方程(23)的离散方程获得，

$f_{i,j}^{n+1/2}=f_{i,j}^n-\frac{\mathrm{\&Delta;\&tau;}}{\mathrm{\&Delta;\&xi;}}[G\left(f_{i,j+1/2}^n\right)-G\left(f_{i,j-1/2}^n\right)],---\left(49\right)$

其中，i，j分别为λ-ξ方向上计算网格的标记，n是时间离散标记。同时，流函数几何状态变量更新至新的时间n+1/2，

同理，沿着λ方向是流线的方向。沿着λ方向也可求解黎曼问题，被称作沿着流线方向的黎曼问题。将方程(7)写成仅沿着λ方向的形式，

$\frac{\&PartialD;f}{\&PartialD;\mathrm{\&tau;}}+\frac{\&PartialD;F\left(f\right)}{\&PartialD;\mathrm{\&lambda;}}=0,$ $f=\left\{\begin{array}{cc}{f}_L,& \mathrm{\&lambda;}<0,\\ f_R,& \mathrm{\&lambda;}>0,\end{array}\right.---\left(51\right)$

其中的符号与公式(23)中的意义相同。图5表示了沿着流线的二维黎曼问题的定义。与图4中的定义相同，二维黎曼问题中的左右变量和左右中间变量中的速度包含了两个笛卡尔坐标系下的分量。这两个速度分量同时跨过用特征方程(10)表示的激波、膨胀波和接触波。根据方程(7)的旋转不变性，沿着流线方向的二维黎曼问题的求解过程与前述(24)-(50)相同，只是将所有公式中的u与v互换、sinθ与cosθ互换、tgθ与ctgθ互换即可。

时间更新的解可以从方程(51)的离散方程获得，

$f_{i,j}^{n+1}=f_{i,j}^{n+1/2}-\frac{\mathrm{\&Delta;\&tau;}}{\mathrm{\&Delta;\&lambda;}}[F\left(f_{i+1/2,j}^n\right)-F\left(F_{i-1/2,j}^n\right)].---\left(52\right)$

同时，流函数几何状态变量更新至新的时间n+I，

4.附图说明

图1(a)黎曼问题的物理模型

图2(b)黎曼问题的定义

图2(a)计算网格交界面以及其上的速度矢量

图2(b)黎曼问题处理多维空间问题时的速度分解

图3从笛卡尔坐标系到流函数平面坐标系的变换示意图

图4跨过流线的二维黎曼问题的定义

图5沿着流线的二维黎曼问题的定义

图6入口马赫数为0.5的抛物型型喷管的流场的拉格朗日方法的解和JST方法的解

(a)流函数平面上的网格

(b)欧拉平面上的网格

(c)流函数解的流线

(d)欧拉方法解的流线

(e)流函数方法解得的固体壁面上压力分布和精解的比较

(f)流函数方法构造的网格

5.具体实施方式

按照本发明介绍的方法，下面给出一个完整的，用流函数的欧拉方程模拟无粘、亚音速喷管内的流动例子。其中，用到本发明给出的求解二维黎曼问题的数值方法计算计算网格边界上的对流项通量值。这个例子中的喷管是一个二维的、壁面呈抛物线形、长度为L、入口高度为H_in的扩张型喷管。其几何尺寸是由两段抛物线定义的，

$H\left(x\right)=\left\{\begin{array}{cc}-a{x}^2,& 0\&le;x<L/2;\\ a{(x-L)}^2-b,& L/2\&le;x\&le;L,\end{array}\right.---\left(54\right)$

其中H_in＝L/3，a＝H_in/2，b＝H_in/L。无粘、可压缩流在喷管入口的马赫数是M_in＝0.5。流动是纯亚音速。流函数形式的欧拉方程(7)的求解的方法中，空间离散采用有限体积法，其中，按照Godunov方法获得的计算网格边界上的对流项通量算子在空间上是二阶精度，采用带有minmod通量限制器的MUSCL插值；时间离散采用Strang换方向的二阶精度。计算参数如下：CFL＝0.48；总共60x20个计算网格单元在流函数平面。计算结果将和基于公知的JST方法直接对方程(1)得到的结果、及其本实施例子的精解做比较。从现在起，定义本发明给出的解为流函数解，在欧拉平面上的JST方法到的解为欧拉解。具体地，计算采用以下步骤，

(1)初始化。流场变量Q⁰＝[ρ⁰，u⁰，v⁰，p⁰]^T及

均为初始值，其中Q⁰取无穷远出的值，上标0表示流动问题在初始时刻，在x-y平面和λ-ξ平面，

t＝0(τ＝0)。进而可求得守恒变量f⁰。然后在λ-ξ平面生成以均匀的空间步长

(Δλ，Δξ)为单元的直角网格。对应的x-y平面上的网格可按照下式生成

dx＝cosθdλ；dy＝sinθdλ                    (55)

(2)计算时间步长，CFL＜0.5，

其中，

其中，

θ^n-1从上一个时刻获得。

(3)沿着λ-方向用Godunov方法求解方程(51)。计算网格边界上的对流项通量值采用本

发明中提出的二维黎曼问题的求解方法。对于二阶精度的插值方法，使用带有TVD通量限制器的MUSCL插值。通过已知的

和

值，将守恒变量从

更新至

(4)沿着ξ-方向用Godunov方法求解方程(23)。计算网格边界上的对流项通量值采用本

发明中提出的二维黎曼问题的求解方法。对于二阶精度的插值方法，使用带有TVD通量限制器的MUSCL插值。用n时刻的已知的f^n+1/2和Q^n+1/2的值，将守恒变量从

更新至

(5)更新流函数几何变量。因为网格交界面上的速度已知，更新的拉格朗日几何变量可以由方程(57)获得。精度等级由速度

决定。

(6)找到新时刻的原始变量从守恒变量

解码后，原始变量为

$\mathrm{\&rho;}=\frac{f_1}{J};$ $u=\frac{f_2}{f_1};$ $v=\frac{f_3}{f_1};$ $p=\frac{(\mathrm{\&gamma;}-1)}{J}(f_4-\frac{f_2^2+f_3^2}{2f_1}).---\left(58\right)$

(6)构造网格。如果需要在每一次历遍迭代的最后一步构造网格。网格点由下式给出，

$x_{i,j}^{n+1}=x_{i,j}^n+\frac{1}{2}\mathrm{\&Delta;}{\mathrm{\&tau;}}_{i,j}^n(h{u}_{i,j}^n+h{u}_{i,j}^{n+1});$ $y_{i,j}^{n+1}=y_{i,j}^n+\frac{1}{2}\mathrm{\&Delta;}{\mathrm{\&tau;}}_{i,j}^n(h{v}_{i,j}^n+h{v}_{i,j}^{n+1}).---\left(59\right)$

(7)重复步骤(2)。在完成一次时间步长的更新，回到步骤(2)，重复步骤(3)-(7)，直到守恒变量或是原始变量收敛。

图6中，流函数的解起始于图6(a)中表示的由流体颗粒运行距离λ和流函数ξ构成的直角网格。图6(c)中，流函数方法获得的流线与图6(d)中JST方法获得的流线吻合良好，图6(e)表示固体壁面上的压力分布与精解完全一致。图6(f)中由流函数法生成的流线与图4(b)JST法生成的略有不同，沿着y方向的网格线不垂直于流线方向，以保证流线的长度。因该指出的是，图6(f)中表示的最终的网格是由图6(a)中的网格进化而来，尽管这个最终的网格并不需要，但是这个进化过程可以清楚说明流函数方法的原理。每一个网格点式代表一个由计算网格单元表示的流体颗粒。本发明的原理讲述的是，在流函数平面构造的网格线其实就是流线本身。

Claims (8)

1.一种通过求解二维黎曼问题来模拟亚音速、无粘流的数值方法，其特征是包括以下步骤：

(1)将二维欧拉平面上的欧拉方程经过以雅各比矩阵

为特征的坐标变换，变成为时间τ表示的时间方向、由流函数ξ和流体颗粒运行的距离λ表示的两个空间方向所决定的流函数平面上的流函数形式的欧拉方程，它的形式为：

其中，f_S是守恒变量矢量；F_S和G_S分别是流函数平面上的λ和ξ方向上的对流项通量，而且，

$G_S=\left[\begin{array}{c}0\\ -p\sin \mathrm{\&theta;}\\ p\cos \mathrm{\&theta;}\\ 0\\ -u\\ -v\end{array}\right],$ $\mathrm{\&theta;}={\mathrm{tg}}^{-1}\left(\frac{v}{u}\right),$ 以上公式中，ρ，p和E分别

是流体密度、压力、总能；u，v是流动速度的笛卡尔坐标分量；

是流函数几何

状态变量；

(2)建立计算网格；

(3)用求解计算网格单元边界上的二维黎曼问题的Godunov方法求解流函数形式的欧拉方程。

2.根据权利要求1所述的一种通过求解二维黎曼问题来模拟无粘、亚音速流动的数值方法，其中所述的计算网格是流函数平面上以λ和ξ为两个方向的二维直角网格。

3.根据权利要求1所述的一种通过求解二维黎曼问题来模拟无粘、亚音速流动的数值方法，其中所述的求解流函数形式的欧拉方程是沿着时间τ方向进行守恒变量f_S的历遍，以找到其稳定解。

4.根据权利要求1所述的一种通过求解二维黎曼问题来模拟无粘、亚音速流动的数值方法，其中所述的求解计算网格单元边界上的二维黎曼问题的Godunov方法需要求解跨过流线的黎曼问题和沿着流线的黎曼来计算网格单元交界的对流项通量。

5.根据权利要求4所述的跨过流线的黎曼问题和沿着流线的黎曼问题，其特征是：在计算单元交界面两侧形成以激波或者是膨胀波形式存在的左侧变量或右侧变量，在其中间是中间变量；中间变量又可被分为左中间变量和右中间变量。

6.根据权利要求4所述的求解跨过流线的黎曼问题和沿着流线的问题的方法包括以下步骤：

(1)沿着流函数形式的欧拉方程的特征方程积分，将左侧变量或右侧变量与中间变量连接，其中，左侧变量、右侧变量、中间变量由权利要求5给出；

(2)在中间变量中恢复流动速度的幅值；

(3)求解组合函数f(u，v)以找到中间变量的流动角；

(4)在中间变量中求解速度分量。

7.根据权利要求6所述的求解跨过流线的黎曼问题和沿着流线的黎曼问题的方法，其中所述的在中间变量中恢复流动速度的幅值，是利用通过跨过激波的Rankine-Hugoniot关系和跨过膨胀波的焓不变关系来获得的。

8.根据权利要求6所述的求解跨过流线的黎曼问题和沿着流线的黎曼问题的方法，其中所述的组合函数f(u，v)表示为

$f(u,v)=\frac{1}{2}\frac{v\sqrt{u^2+v^2}}{u}+\frac{1}{2}u\ln (v+\sqrt{u^2+v^2}).$

Images