CN102203782A - 求解拉格朗日形式的欧拉方程的数值方法 - Google Patents

Abstract

这个发明是关于一种数值方法，它提供并求解一种新的拉格朗日形式的二维欧拉方程来模拟亚音速流动和求解固体壁面几何形状设计的反问题。该发明提供了一个推导拉格朗日平面上的欧拉方程的变换方式，从而简化计算网格、最大程度降低了对流项的数值耗散。该发明构造了一个混合的通量运算方案，其中，二维黎曼问题在拉格朗日平面得以求解。使用本发明的方法可以同时得到流场物理量的解和固体壁面几何形状的设计的解。

Description

求解拉格朗日形式的欧拉方程的数值方法

1.技术领域

本发明涉及一种模拟亚音速无粘流的流动和求解一类反问题的数值方法。这类反问题是指在亚音速流场中的空气动力学物体的几何形状的设计问题。本发明属于计算流体力学(CFD-Computational Fluid Dynamics)领域。

2.背景技术

2.1先前的工作

在计算流体力学的计算中，大多数是求解欧拉形式的流体控制方程。这意味着在笛卡尔坐标下的欧拉平面中，计算网格必须根据物体的限定事先生成。网格构成网格单元。由于流体穿过网格单元的交界面，所以存在对流项的通量。正是这个通量在数值解中引起数值耗散，因为数值耗散直接与对对流项的数值逼近所引起的误差有关。从上世纪以来，CFD研究者致力于开发更精确、高效率的数值方法来降低这个数值耗散。迎风差分方法在求解流体的流动过程中取得显著成效，因为它合理表达了对流项的特征。典型代表有，Godunov方法[1]，它在网格单元交界面求解黎曼问题，给出非常精确的解；FVS方法[2](通量分裂法)在网格单元交界面运用特征关系，使求解过程更加快捷。但是，作为以欧拉形式中不可避免的对流项的数值耗散，仍然存在于这些基于特征的数值方法中。

另一方面，流体的拉格朗日描述强调流体颗粒在不同位置的运动和特点。对于流体的控制方程，例如，拉格朗日形式的欧拉方程，其中，必须有一个方向代表流线，数值上可以用流函数表示。另一个方向是流体颗粒运动的距离。这个坐标系统构成了拉格朗日平面，在这个平面上，计算网格点理论上就是流体颗粒，网格线总是简单的直角网格。特别是稳定的流动中，流体迹线和流线是重合的，没有流体颗粒穿过流线，穿过网格单元交界面的对流项不存在。所以在拉格朗日平面上求解方程，数值耗散被降低至最小。

应用拉格朗日形式的方程在求解一类反问题时体现了优势。在空气动力学中，一类典型的反问题是通过给定固体壁面上的压力分布，然后设计固体壁面形状以符合压力分布的要求。如果用基于欧拉平面的方法求解这类问题，例如adjoint法[3](共轭方程法)，流程是首先估计这个未被确定的几何形状；然后在其周围生成网格；在对流场求解，下一步，是重要的、费时的，即求解共轭方程，改进几何形状。这个过程反复进行，直到找到目标为止。通常，这个过程持续较长时间。拉格朗日形式十分适合应用在这类几何形状不确定 的问题中，因为在拉格朗日平面，不确定的几何边界，也就是固体壁面，也是由流线表示的，而且，无论几何形状怎样变化，由于沿着一条流线的流函数是常数，所以流函数表示的流线在拉格朗日平面是直线。在拉格朗日平面求解这类几何形状的设计问题可以达到最优(最有效)的过程。

尽管拉格朗日形式表现了如此卓越的优势，以前只是应用于超音速流动中[4]。当欧拉方程在拉格朗日平面上被求解，即方程在空间上向下游方向发展，不需要考虑任何下游对上游的影响，这一切完美地符合超音速流动的特点。以前，拉格朗日形式也成功地应用于二维超音速流动中的固体壁面的几何形状的设计中[5]。到目前为止，应用流函数作为坐标进行求解几何形状设计的反问题的例子仅限于有势流和线性化的可压缩流[6]。

2.2目的和优势

在工业界有明显的需求去应用拉格朗日形式的优势。例如在产品设计的初级阶段，需要对产品进行最小数值耗散的数值模拟研究和快速的几何形状设计。如机翼、喷管的流场计算、外形的设计等。亚音速也是工业界最常见的流动状态。用严格的拉格朗日概念求解亚音速流动会遇到障碍，因为物体的存在会对上游的流体颗粒产生影响。成功应用拉格朗日形式的数值方法的关键是如何正确使上游的流体颗粒感受到这种影响波动。本发明的目的在于(1)提供一个拉格朗日形式的欧拉方程，它在一个坐标方向上没有对流项，从而最大程度降低数值耗散；(2)提供精确求解拉格朗日形式的欧拉方程的数值方法；(3)提供一个在亚音速流场中进行求解几何形状设计的反问题的最优方法。

3.发明内容

3.1二维拉格朗日形式的欧拉方程

本发明首先提供一个从欧拉平面的欧拉变量(t，x，y)到拉格朗日平面上的拉格朗日变量(τ，λ，ξ)的坐标变换，其中变量τ为格拉朗日时间，和时间项t有相同的量纲，被引入作为时间历遍项，另外两个独立的变量分别是流函数ξ(量纲[m²·s^-1])和拉格朗日距离λ(量纲[m])。坐标ξ和流体颗粒的流线重合，λ被定义为不同流体颗粒在流向上的位置。本发明开始于描述二维、无粘、可压缩流体流动的欧拉平面上的欧拉方程

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

其中f是守恒变量矢量；F，G是对流项通量在两个笛卡尔方向上的分量。具体地，

$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在两个笛卡尔方向上的分量、流体密度和压力。总能E和总焓分别是H

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

上式中的γ是气体比热比。从坐标(t，x，y)到(τ，λ，ξ)的坐标变换的雅各比矩阵可写为

且它的绝对值J成为

其中h是拉格朗日行为参数；θ是流动方向角； 

被称作拉格朗日几何状态变量，

单位量纲是[m²s·kg^-1]，被定义为

最终，二维不稳定欧拉方程(1)被转换成拉格朗日平面上的拉格朗日形式，它包括两组偏微分方程

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

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

其中f_L，F_L和G_L与(2)中的f，F，G称谓相同，是在拉格朗日平面。具体地，

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

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

0＜h≤1。                                                            (13)

方程(7)被称作拉格朗日物理守恒方程；(8)是几何守恒方程。在(9)中，可以发现矢量G_L的第一和第四个元素是零，这意味着在拉格朗日平面上，沿着ξ方向，对于连续方程和能量方程没有通量穿过各单元的交界面，所以数值耗散被减小。此外，拉格朗日形式的欧拉方程还有其他特性。它们是，

1.方程(7)满足齐次特性，这意味着FVS方法可被用来求解通量项。

2.方程(7)不满足旋转不变性，这意味着在不同的坐标方向上应该用不同的数值方法求解对流项通量，以保证各自的特征不被模糊掉。

3.偏微分方程(7)和(8)都是强双曲线型的，这意味着对于初值问题，方程可以用时间历遍的方法来求解。

3.2拉格朗日形式的欧拉方程的特征和特征方程

对于拉格朗日物理守恒方程(7)，可以找到其特征(特征值和特征向量)。按照雅各比矩阵的定义，

$B_L=\frac{{\∂F}_L}{\∂Q},C_L=\frac{{\∂G}_L}{\∂Q},---\left(14\right)$

其中Q＝[ρ，u，v，p]^T。雅各比矩阵B_L的四个特征值是，

及其对应的、沿着λ方向的四个特征向量是，

雅各比矩阵C_L的四个特征值是，

$\left.\begin{array}{c}{\mathrm{\μ}}_{L_{\mathrm{1,2}}}=0,\\ {\mathrm{\μ}}_{L_{\mathrm{3,4}}}=\±a/J;\end{array}\right\}---\left(17\right)$

及其对应的、沿着ξ方向的四个特征向量是，

$K_L^{\left(1\right)}=\left[\begin{array}{c}1\\ 0\\ 0\\ 0\end{array}\right],K_L^{\left(2\right)}=\left[\begin{array}{c}0\\ \cos \mathrm{\θ}/\mathrm{\ρ}\\ \sin \mathrm{\θ}/\mathrm{\ρ}\\ 0\end{array}\right],K_L^{\left(3\right)}=\left[\begin{array}{c}1\\ -\sin \mathrm{\θ}\left(\frac{a}{\mathrm{\ρ}}\right)\\ \cos \mathrm{\θ}\left(\frac{a}{\mathrm{\ρ}}\right)\\ a^2\end{array}\right],K_L^{\left(4\right)}=\left[\begin{array}{c}1\\ \sin \mathrm{\θ}\left(\frac{a}{\mathrm{\ρ}}\right)\\ -\cos \mathrm{\θ}\left(\frac{a}{\mathrm{\ρ}}\right)\\ a^2\end{array}\right],---\left(18\right)$

其中a是声速。沿着 拉格朗日物理守恒方程(7)的特征方程可被写为

沿着 

特征方程可被写为

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

根据坐标变换的雅各比矩阵(4)，在拉格朗日平面沿着λ方向的黎曼不变量为

沿着ξ方向的为

方程(15)～(22)被用来构造求解控制方程(7)和(8)的数值方法。

3.3数值方法的构造

为在拉格朗日平面为求解方程(7)和(8)，本发明采用有限容积法，变量存储在网格单元的中心。根据方程的齐次、旋转、双曲线特性，同时考虑到计算的效率和精度，本发明产生一种新的求解方法，带有混合迎风差分格式通量算予的、Strang换方向的数值方法，意味着，在计算网格交界面的对流项通量时，在λ方向用FVS方法，在ξ方向用求解黎曼问题的Godunov方法，体现了Godunov方法和FVS方法的各自优势。Strang换方向是二阶精度时间离散。从现在起，为方便起见，方程和变量中省去拉格朗日含义的下标L。

3.3.1沿着λ方向的FVS法

首先，将方程(7)写成仅沿着λ方向的形式

$\frac{\∂f}{\∂\mathrm{\τ}}+\frac{\∂F}{\∂\mathrm{\λ}}=0.---\left(23\right)$

通过应用FVS(通量矢量分裂)法，可以获得各网格单元(i，j)中的更新的守恒变量

$f_{i,j}^{n+1}=f_{i,j}^n-\frac{\mathrm{\Δ\τ}}{\mathrm{\Δ\λ}}[L_i^{+}{f}_{i,j}^n-L_{i-1}^{+}{f}_{i-1,j}^n+L_{i+1}^{-}{f}_{i+1,j}^m-L_i^{-}{f}_{i,j}^n],---\left(24\right)$

其中n是时间步长的序号；i，j是网格单元的序号；Δτ是时间步长；Δλ是空间步长。上式中分裂的矩阵表示为

$L^{+}=K_{\mathrm{Lcv}}{\mathrm{\Λ}}^{+}{K}_{\mathrm{Lcv}}^{-1}$ 和 $L^{-}=K_{\mathrm{Lcv}}{\mathrm{\Λ}}^{-}{K}_{\mathrm{Lcv}}^{-1},---\left(25\right)$

其中，分裂的特征值矩阵为

${\mathrm{\Λ}}^{+}=\left[\begin{array}{cccc}{\mathrm{\μ}}_{L_1}& & & \\ & {\mathrm{\μ}}_{L_2}& & \\ & & {\mathrm{\μ}}_{L_3}& \\ & & & 0\end{array}\right]\mathrm{and}{\mathrm{\Λ}}^{-}=\left[\begin{array}{cccc}0& & & \\ & 0& & \\ & & 0& \\ & & & {\mathrm{\μ}}_{L_4}\end{array}\right];---\left(26\right)$

以及矩阵

和它的逆矩阵

3.3.2沿着ξ方向求解黎曼问题的Godunov法

对于Godunov方法，核心内容是求解黎曼问题。本发明沿着ξ方向求解黎曼问题，被称作跨过流线的黎曼问题。图1表示了跨过流线的黎曼问题的定义。将方程(7)写成仅沿着ξ方向的形式，

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

其中f_L和f_R分别是网格交界面左侧和右侧变量中的守恒变量。从现在开始，下标L、R和*分别代表黎曼问题中的左、右和中间变量。如图1所示，左侧变量或右侧变量是激波或者是膨胀波，在其中间是中间变量，又可被接触波分为左中间变量和右中间变量。时间更新的解可以从方程(29)的离散方程获得，

$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(30\right)$

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

$\left\{\begin{array}{c}{p}_{*}+C_L\&CenterDot;f(u_{*},v_{*})=p_L+C_L\&CenterDot;f(u_L,v_L),---\left(31\right)\\ p_{*}-C_R\&CenterDot;f(u_{*},v_{*})=p_R-C_R\&CenterDot;f(u_R,v_R),---\left(32\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(33\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(34\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(35\right)$

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

$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(36\right)$

其中 

u＝Vcosθandv＝V sinθ。方程(31)～(34)的解表达了黎曼问题中的中间变量的值

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

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

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

$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(42\right)$

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

tgθ_*＝ε，                                (43)

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

$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(44\right)$

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

$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(45\right)$

下面的工作是通过给定的速度V_*求解方程F(ε)＝0以发现ε，再求出θ_*。方程(45)是可微分的，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})}{{\left({1+\mathrm{\&epsiv;}}^2\right)}^{3/2}}\right).---\left(46\right)$

数值试验表明，对于无量纲速度V_*≤5(亚音速区间)而且ε∈[-1，1]的范围内，F′(ε)＞0，这意味着函数F(ε)在此范围内是单调的。同时，F(-ε)·F(ε)＜0，所以Newton-Raphson历遍方法可用来找到方程(45)在给定V_*时的根ε。为找到V_*的值，需要用p_*，ρ_*L，ρ_*R的值去还原黎曼问题中的各种非线性波(接触波、膨胀波、激波)。本发明中考虑了下面的关系：

跨过激波的Rankine-Hugoniot关系

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

$\left\{\begin{array}{c}{\mathrm{\&mu;}}_s({\mathrm{\&rho;}}_L{J}_L-{\mathrm{\&rho;}}_{*L}{J}_{*L})=0,---\left(46\right)\\ {\mathrm{\&mu;}}_s({\mathrm{\&rho;}}_L{J}_L{u}_L-{\mathrm{\&rho;}}_{*L}{J}_{*L}{u}_{*L})=p_{*}\sin {\mathrm{\&theta;}}_{*}-p_L\sin {\mathrm{\&theta;}}_L,---\left(47\right)\\ {\mathrm{\&mu;}}_s({\mathrm{\&rho;}}_L{J}_L{v}_L-{\mathrm{\&rho;}}_{*L}{J}_{*L}{v}_{*L})=p_L\cos {\mathrm{\&theta;}}_L-p_{*}\cos {\mathrm{\&theta;}}_{*},---\left(48\right)\\ {\mathrm{\&mu;}}_s({\mathrm{\&rho;}}_L{J}_L{E}_L-{\mathrm{\&rho;}}_{*L}{J}_{*L}{E}_{*L})=0,---\left(49\right)\end{array}\right.$

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

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

和

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

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

跨过膨胀波的焓不变关系

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

$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(52\right)$

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

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

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

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

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

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

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

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

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

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

在获得θ_*之后，按照公式(9)和由公式(44)求得的u_*L，v_*L或u_*R，v_*R，在加上p_*，便可求得公式(30)中的G_L。同时，几何守恒方程也可被更新至新的时间，

和欧拉形式的方法相比，尽管这个步骤是额外的，但这个步骤简单，不占用过多的计算时间。

3.3.3边界条件

处理边界条件需要围绕在计算域周围再加上两层网格，被称为虚拟网格，它可以使空间离散算法扩展至边界。

固体壁面边界条件

构造的数值方法中，需要沿着ξ方向在物体壁面求解黎曼问题。物体壁面上的黎曼问题包括左侧或右侧变量和针对物体壁面的镜像变量。在物体壁面边界条件中，下标L表示物理变量ρ，p，u，v，θ在虚拟网格中，下标R表示这些物理变量在壁面网格中。θ_w是物体壁面的角度。根据(37)～(40)给出的黎曼问题的解，在物体壁面边界上有

p_*＝p_R，                                                (58)

${\mathrm{\&rho;}}_{*L}={\mathrm{\&rho;}}_{*R}={\mathrm{\&rho;}}_R+\frac{1}{2}(f(V_R,{\mathrm{\&theta;}}_L)-f(V_R,{\mathrm{\&theta;}}_R))({\mathrm{\&rho;}}_R/a_R).---\left(59\right)$

中间变量中的压力即可被认为是物体壁面上的压力p_w，因而

p_w＝p_*。                                                (60)

如果物体壁面上的压力已经被给定，这个给定压力可以直接被用做黎曼问题中的中间变量中的压力，所以

p_*＝p_w.                                                 (61)

出口、入口边界条件

第一个关系对应的是正向、进入的特征线。以下标b表示的边界上的值可以通过式(21)、(22)中的黎曼不变量获得，

其中，u_∞，a_∞是来流的速度和声速，下标i表示计算域内部的值。第二个关系是通过内部的值的插值获得的数值逼近。因而，速度u_b可以通过式(62)和(63)获得，于是有

相似的方法可以处理亚音速出口边界条件。但需要在出口边界规定外界压力值p_∞，有

拉格朗日几何状态变量的边界条件

在物体壁面边界上，

对于入口、出口边界，规定流动角度θ＝0，则

图2给出了从 

更新到 

的一个步骤的数值方法的流程图，这个方法在空间和时间上都是二阶精度，其中Q＝[ρ，u，v，p]^T是流场的物理变量。

3.4求解一类反问题的方法

使用拉格朗日形式时，固体壁面由流函数代表的流线表示，所以在拉格朗日平面，由流线代表的网格线是直线。本发明的另一个目的是探索拉格朗日形式的数值方法在求解一类反问题的优势。这类反问题被定义为几何形状的设计问题。

第2部分引入了拉格朗日平面上沿着ξ方向的黎曼问题的解和固体壁面边界条件。从(37)和(57)，可以获得

$p_w=p_R+\frac{{\mathrm{\&rho;}}_R{a}_R}{2}(f(V_R,{\mathrm{\&theta;}}_L)-f(V_R,{\mathrm{\&theta;}}_R)),---\left(70\right)$

其中p_R，ρ_R，a_R，V_R，θ_R是边界上已知的量，f是公式(35)和(36)给出的组合函数。从(70)推导出

$f(V_R,{\mathrm{\&theta;}}_L)=f(V_R,{\mathrm{\&theta;}}_R)+\frac{2(p_w-p_R)}{{\mathrm{\&rho;}}_R{a}_R}.---\left(71\right)$

该式在规定壁面压力分布p_w后可以求解出虚拟网格中流动角度θ_L的大小。根据镜像边界条件，固体避免的角度可以从下式求得

${\mathrm{\&theta;}}_w=\frac{1}{2}({\mathrm{\&theta;}}_L+{\mathrm{\&theta;}}_R).---\left(72\right)$

几何形状的设计问题意味着要根据给定的固体壁面压力分布，找到固体壁面的角度。这个过程需要求解一个反黎曼问题。一个反黎曼问题的求解器按照下面的步骤建立，首先

定义

tgθ_L＝ε，                            (73)

然后，虚拟网格中的速度分量按下式求得

$u_L=V_R\cos {\mathrm{\&theta;}}_L=\frac{V_R}{\sqrt{1+{\mathrm{\&epsiv;}}^2}}$ 和 $v_L=V_R\sin {\mathrm{\&theta;}}_L=\frac{\mathrm{\&epsiv;}{V}_R}{\sqrt{1+{\mathrm{\&epsiv;}}^2}},---\left(74\right)$

然后被带入三角形式给出的组合函数(71)中，并被进一步简化为ε的函数，

$F\left(\mathrm{\&epsiv;}\right)=\frac{V_R}{2}(\mathrm{\&epsiv;}+\frac{\ln V_R+\ln (\sqrt{1+{\mathrm{\&epsiv;}}^2}+\mathrm{\&epsiv;})-\ln \sqrt{1+{\mathrm{\&epsiv;}}^2}}{\sqrt{1+{\mathrm{\&epsiv;}}^2}})-f(V_R,{\mathrm{\&theta;}}_R)-\frac{2(p_w-p_R)}{{\mathrm{\&rho;}}_R{a}_R}.---\left(75\right)$

下面的任务是求解规定物理参数ρ_R，p_R，V_R，θ_R和p_w条件下的方程F(ε)＝0，解得ε后再求得θ_L。方程(75)是可微分的，F针对ε的一阶导数是

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

相同于求解方程(45)的方法，更新的ε值表示为

${\mathrm{\&epsiv;}}^{n+1}={\mathrm{\&epsiv;}}^n-\frac{F\left({\mathrm{\&epsiv;}}^n\right)}{F^{\&prime;}\left({\mathrm{\&epsiv;}}^n\right)}.---\left(77\right)$

虚拟网格中的流动角度可由θ_L＝tg^-1(εⁿ⁺¹)求得。固体壁面角度θ_w按式(72)求得。

图3给出了求解几何外形的设计的反问题的算法流程图。这个方法起始于假设的流场物理量的参数值和固体壁面角度。按固体壁面上按照的指定的压力分布，求解反黎曼问题，在虚拟网格单元中得出流动角度，然后壁面形状可获得并被带回到流场求解器中求解更新流场物理量。这个过程持续进行，直到历遍的固体壁面的角度达到收敛为止。很明显，收敛标准是壁面角度而不是任何一个物理变量，这意味着，壁面角度也变成了一个流场变量。壁面边界的变化并不改变方程的特征，方程仍然是强双曲线型，仍然是一个初值问题。与应用在欧拉平面上的方法(如adjoint法)相比，现在的方法既不需要在更新的几何形状周围生成网格，也不需要求解adjoint方程，它同时获得收敛的流场物理变量的解和几何形状的解。针对这类问题，总体的CPU时间降低一个量级。

4.附图说明

图1拉格朗日平面上沿ξ方向跨过流线的黎曼问题

图2使用二阶精度方法，带有混合迎风差分格式通量算子的、Strang换方向的数值方法进行一步流场更新的流程图，其中的标号代表以下变量或公式

① ② 

③ 

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

⑥ 

⑦ 

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

⑩ 

图3几何形状设计的反问题的计算流程图

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

(a)拉格朗日平面上的网格

(b)欧拉平面上的网格

(c)拉格朗日方法解得的流线

(d)欧拉方法解得的流线

(e)拉格朗日方法解得的固体壁面上压力分布和精解的比较

(f)拉格朗日方法构造的网格

图5根据给定压力分布计算的固体壁面的几何形状

(a)欧拉平面上的网格

(b)固体壁面上压力分布

(c)用本发明的方法计算的固体壁面形状

5.具体实施方式

5.1扩张喷管中的无粘亚音速流动

按照本发明介绍的方法，下面给出一个完整的用拉格朗日形式的欧拉方程模拟无粘、亚音速内流。这个例子是关于亚音速流通过一个二维抛物线形、长度为L、入口高度为H_in的扩张型喷管。其几何尺寸是由两段抛物线定义的，

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

其中H_in＝L/3，a＝H_in/2，b＝H_in/L。无粘流在喷管入口的马赫数是M_in＝0.5。流动是纯亚音速。拉格朗日形式的欧拉方程(7)和(8)将用第二节介绍的方法求解，其中，带有混合迎风差分格式通量算子在空间上是二阶精度，采用带有minmod通量限制器的MUSCL插值。Strang换方向在时间上是二阶精度。计算参数如下：CFL＝0.48；拉格朗日行为参数h＝0.98；总共60×20个计算网格单元在拉格朗日平面。计算结果将和基于欧拉平面上的有限体积法的JST方法[7]得到的结果及其精解做比较。从现在起，定义本发明给出的解为拉格朗日解，由欧拉方法在欧拉平面上的到的解为欧拉解。参考图2，具体地讲，计算采用以下步骤，

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

均为初始值，其中Q⁰取无穷远出的值， 上标0表示流动问题在初始时刻，在x-y平面和λ-ξ平面，t＝0(τ＝0)。进而可求得守恒变量f⁰。然后在λ-ξ平面生成以均匀的空间步长(Δλ，Δξ)为单元的直角网格。对应的x-y平面上的网格可按照下式生成

dx＝cosθdλ和dy＝sinθdλ                   (79)

(2)计算时间步长。按照变量在时间n(n＝0，1，2，...)，CFL＜0.5，

$\mathrm{\&Delta;\&tau;}=\mathrm{CFL}\&CenterDot;\max (\frac{\mathrm{\&Delta;\&lambda;}}{|(1-h)\sqrt{u^2+v^2}+(a/J)\sqrt{U^2+V^2}|},\frac{\mathrm{\&Delta;\&xi;}}{|a/J|}),---\left(80\right)$

其中， 

θ^n-1从上一个时刻获得，所以有

(3)沿着λ-方向用FVS方法求解方程(23)。通过已知的 

和 

值，守恒变量从 

更新至 

对于二阶精度的数值方法，使用带有TVD通量限制器的MUSCL插值。

(4)沿着ξ-方向用Godunov方法求解方程(29)。用n(n＝0，1，2，...)时刻的已知的f^n+1/2和Q^n+1/2的值，用Godunov方法更新守恒变量。当地的黎曼问题求解器给出网格单元交界面上的原始变量[ρ，u，v，p]_i，j±1/2值。对于二阶精度的数值方法，使用带有TVD通量限制器的MUSCL插值。

(5)更新拉格朗日几何变量。因为网格交界面上的速度 

已知，更新的拉格朗日几何变量可以由方程(57)获得。精度等级由速度 

决定，尽管公示(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(82\right)$

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

$\left.\begin{array}{c}{x}_{i,j}^{n+1}=x_{i,j}^n+\frac{1}{2}{\mathrm{\&Delta;\&tau;}}_{i,j}^n({\mathrm{hu}}_{i,j}^n+{\mathrm{hu}}_{i,j}^{n+1})\\ y_{i,j}^{n+1}=y_{i,j}^n+\frac{1}{2}{\mathrm{\&Delta;\&tau;}}_{i,j}^n({\mathrm{hv}}_{i,j}^n+{\mathrm{hv}}_{i,j}^{n+1})\end{array}\right\}---\left(83\right)$

每一个网格点式代表一个由计算网格单元表示的流体颗粒。本发明的原理讲述的是，在拉格朗日平面构造的网格线其实就是流线本身。

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

图4中，拉格朗日的解起始于图4(a)中表示的由拉格朗日距离λ和流函数ξ构成的直角网格。图4(c)中，拉格朗日方法获得的流线与图4(d)中欧拉方法获得的流线吻合良好，图4(e)表示固体壁面上的压力分布与精解完全一致。图4(f)中由拉格朗日生成的流线与图4(b)略有不同，沿着y方向的网格线不垂直于流线方向，以保证流线的长度。因该指出的是，图4(f)中表示的最终的网格是由图4(a)中的网格进化而来，尽管这个最终的网格并不需要，但是这个进化过程可以清楚说明拉格朗日方法的原理。

5.2几何形状设计的反问题

下面的例子是介绍用求解二维拉格朗日形式的欧拉方程去解决固体壁面几何形状的设计问题，是一类反问题。例子中是一个收缩-扩张喷管。例子中先规定几何形状，求得固体壁面压力分布，然后用着个压力分布作为输入，用本发明的方法求得固体壁面的几何形状，结果应该与最初的几何形状一致。喷管固体壁面外形由下面的函数给出

$y=H_{\mathrm{yh}}{e}^{-\mathrm{\&beta;}{x}^2},$ -1≤x≤1，                            (84)

其中H_th取0.1，表明收缩喉口高度和入口高度之比为10％。β取6.5以控制喉口长度约占喷管长度的三分之一。入口马赫数为0.5，是纯亚音速流动。喷管的几何形状和欧拉平面上的计算网格如图5(a)所示。拉格朗日平面上的网格与上例中图中4(a)所示一致。拉格朗日形式的欧拉方程(7)和(8)将用第二节介绍的方法求解，其中，带有混合迎风差分格式通量算子在空间上是二阶精度，采用带有minmod通量限制器的MUSCL插值。Strang换方向在时间上是二阶精度。计算参数如下：CFL＝0.48；拉格朗日行为参数h＝0.98；两组 网格数目60×20和90×30个计算网格单元在拉格朗日平面。首先按照上例中给出的计算流程完成计算，获得固体壁面压力分布，如图5(b)所示，并作为求解几何形状设计的反问题的压力分布的输入。同时，adjoint方法将同时应用这个例子以检验本发明的方法的计算效率。具体步骤是：

(1)初始化流场和固体壁面角度。流场初始变量的设置于上例相同。固体壁面初始角度假设θ_w＝0°。因为流场中的流动角度设为0，所以虚拟网格单元中的流动角度θ_L＝0°。

(2)输入规定的固体壁面的压力分布p_w。

(3)在固体壁面上求解基于给定压力分布p_w的反黎曼问题，公式按照第二节(70)～(77)。

(4)计算固体壁面角度。固体壁面角度计算式根据镜像条件，所以有

${\mathrm{\&theta;}}_w=\frac{1}{2}({\mathrm{\&theta;}}_R+{\mathrm{\&theta;}}_L)---\left(84\right)$

(5)检验固体壁面角度的收敛性。固体壁面角度θ_w的变化量是两个时刻n+1和n的差值，

${\mathrm{\&delta;\&theta;}}_w=({\mathrm{\&theta;}}_w^{n+1}-{\mathrm{\&theta;}}_w^n).---\left(85\right)$

如果δθ_w小于收敛标准，流场物理量和固体壁面的角度 

均达到最终的解；否则，继续下一步。

(6)更新流场变量。这一步的任务是根据更新的固体壁面角度 

计算流场的物理参数从 

更新到 

这一步包括上例中的步骤(1)～(7)。完成这一步后，获得流场新的物理参数变量(ρ_R，p_R，V_R，θ_R，)。

(7)重复步骤(2)。

图5(c)表示了计算的固体几何形状和原始的形状的比较。从结果可以估计计算精度，计算的固体几何形状和原始的形状吻合良好，最大误差在之内0.5％。最后，在计算效率方面，本发明的方法使用的计算的时间是adjoint方法的十分之一。

结论

本发明提供了一个模拟无粘亚音速流动的新技术。它包括一个坐标变换，将二维欧拉平面上的欧拉形式的欧拉方程转换到拉格朗日平面上。在拉格朗日平面上不仅可以使用简单的直角网格，而却在采用基于特征的数值方法求解方程时产生最小的数值耗散。使用这个新技术的一个优势在于求解一类几何形状设计的反问题。在固体壁面上求解反黎曼问题可以降低问题的复杂性，将计算时间降低一个量级。它同时给出流场的物理参数的解和固体壁面的几何形状的解，从而使这类问题达到了最高效率。

Reference

(1)Godunov，S.K.：A difference scheme for numerical computation discontinuous solution of hydrodynamic equations.Translated US Publ.Res.service，JPRS 7226，1969.

(2)Van Leer，B.：Flux vector splitting for the 1990′s.Invited Lecture for the CFD Symposium on Aero propulsion，Cleveland，Ohio，1990.

(3)Jameson，A.：Aerodynamic design via control theory.Journal of Scientific Computing，3：233-260，1988.

(4)Loh，C.Y.and Liu，M.S.：A new Lagrangian method for three-dimensional steady supersonic flows.Journal of Computational Physics，113，pp.224-248，1994.

(5)Mateescu，D.： Analysis of aerodynamic problem with geometrically unspecified boundaries using an enhanced Lagrang

(6)Stanitz，D.John：A review of certain inverse methods for the design of duct with 2-or3-dimensional potential flow.Appl.Mech.Rev.Vol.41，No.6，June 1988.

(7)Jameson，A.，Schmidt，W.，Turkel，E.：Numerical solutions of the Euler equations by finite volume methods using Runge-Kutta time-stepping schemes.AIAA-paper，81-1259.，1981.

Claims (16)

1.一种模拟无粘、亚音速流动的数值方法，其特征是，包括以下步骤：

(1)将二维欧拉平面上的欧拉方程变换为拉格朗日平面上的拉格朗日形式的欧拉方程；

(2)录入物体几何形状的表面坐标；

(3)在所述物体周围建立计算网格；

(4)用带有混合迎风差分格式通量算子的、Strang换方向的数值方法求解拉格朗日形式的欧拉方程。

2.根据权利要求1所述的方法，其中所述的拉格朗日形式的欧拉方程方程包括两个方程，拉格朗日物理守恒方程和几何守恒方程。

3.根据权利要求1所述的方法，其中所述的拉格朗日平面有一个由拉格朗日时间τ表示的时间方向和由流函数ξ和拉格朗日距离λ表示的两个空间方向。

4.根据权利要求1所述的方法，其中所述的二维欧拉平面上的欧拉方程的变换包括一个转换矩阵，

其中，h被称作拉格朗日行为参数；θ是流动方向角；u，v是流动速度V的笛卡尔坐标分量；

是拉格朗日几何状态变量。

5.根据权利要求2所述，其中拉格朗日物理守恒方程表示为

其中，f_L是守恒变量矢量；F_L，G_L分别是拉格朗日平面上的λ和ξ方向上的对流项通量；

$f_L=\left[\begin{array}{c}\mathrm{\&rho;J}\\ \mathrm{\&rho;Ju}\\ \mathrm{\&rho;Jv}\\ \mathrm{\&rho;JE}\end{array}\right],$

$G_L=\left[\begin{array}{c}0\\ -p\sin \mathrm{\&theta;}\\ p\cos \mathrm{\&theta;}\\ 0\end{array}\right],$ $\mathrm{\&theta;}={\mathrm{tg}}^{-1}\left(\frac{v}{u}\right),$

，0＜h≤1，其中，变量t，ρ，p，E和H分别是时间、流体密度、压力、总能、总焓。

6.根据权利要求2所述，其中几何守恒方程表示为

其中

$F_g=\left[\begin{array}{c}0\\ 0\end{array}\right],$ $G_g=\left[\begin{array}{c}-\mathrm{hu}\\ -\mathrm{hv}\end{array}\right].$

7.根据权利要求1所述的方法，其中所述的计算网格是拉格朗日平面上以λ和ξ为两个方向的二维直角网格。

8.根据权利要求1所述的方法，其中所述的在拉格朗日平面上求解拉格朗日形式的欧拉方程方程是沿着拉格朗日时间τ方向进行守恒变量的历遍，以找到其稳定解。

9.根据权利要求1所述的方法，其中所述的用Strang换方向法及其混合的通量运算方案求解拉格朗日形式的欧拉方程在守恒变量f_L沿着拉格朗日时间τ方向进行历遍时的每一步更新过程中包括以下步骤：

(1)用FVS方法沿λ方向以时间步长Δτ来计算网格单元交界的对流项通量；

(2)用求解跨过流线的黎曼问题的Godunov方法沿ξ方向以时间步长Δτ的二分之一来计算网格单元交界的对流项通量；

(3)再次用FVS方法沿λ方向以时间步长Δτ来计算网格单元交界的对流项通量。

10.根据权利要求9所述的方法，其中所述沿ξ方向的跨过流线的黎曼问题包括在计算单元交界面左侧变量或右侧变量是激波或者是膨胀波，在其中间是中间变量，又可被分为左中间变量和右中间变量。

11.根据权利要求9所述的方法，其中所述的求解跨过流线的黎曼问题包括以下步骤：

(1)沿着拉格朗日物理守恒方程的特征方程积分，将左侧变量或右侧变量与中间变量连接；

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

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

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

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

13.根据权利要求11所述的方法，其中所述的组合函数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}).$

14.一种求解固体壁面几何形状设计的反问题的数值方法，其特征是，包括以下步骤：

(1)初始化流场变量参数和固体壁面角度；

(2)录入规定的分布在物体壁面上的压力值；

(3)求解反黎曼问题以找到物体壁面边界的镜像变量；

(4)计算固体壁面角度；

(5)检测固体壁面角度的收敛性，如果收敛则计算结束；

(6)更新流场变量参数；

(7)重复步骤(3)。

15.根据权利要求14所述的方法，其中所述的求解反黎曼问题包括求解一个已知压力分布的组合函数，该组合函数表示为

其中，p_R，ρ_R，a_R，V_R，θ_R是在固体壁面边界上的已知的流动参数；f是权力要求13中所述的组合函数；p_w是规定的固体壁面的压力分布；

θ_L是虚拟网格单元中的流动角度。

16.根据权利要求14所述的方法，其中所述的求解反黎曼问题包括在固体壁面上用Newton-Raphson历遍的方法求解权利要求15中所述的已知压力分布的组合函数，以找到流动角度θ_L。

Images