用拉格朗日形式的欧拉方程求解一类反问题的数值方法

摘要

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

说明书

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{\partial f}{\partial t}+\frac{\partial F}{\partial x}+\frac{\partial G}{\partial y}=0,---\left(1\right)$

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

$f=\left(\begin{array}{c}\rho \\ \mathrm{\rho u}\\ \mathrm{\rho v}\\ \mathrm{\rho E}\end{array}\right),$$F=\left(\begin{array}{c}\mathrm{\rho u}\\ {\mathrm{\rho u}}^2+p\\ \mathrm{\rho uv}\\ (\mathrm{\rho E}+p)u\end{array}\right),$$G=\left(\begin{array}{c}\mathrm{\rho v}\\ \mathrm{\rho uv}\\ {\mathrm{\rho v}}^2+p\\ (\mathrm{\rho 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}{\gamma -1}\frac{p}{\rho }$and$H=E+\frac{p}{\rho }=\frac{u^2+v^2}{2}+\frac{\gamma }{\gamma -1}\frac{p}{\rho }.---\left(3\right)$

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

且它的绝对值J成为

其中h是拉格朗日行为参数；θ是流动方向角；被称作拉格朗日几何状态变量， 单位量纲是[m²s·kg^-1]，被定义为

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

$\frac{\partial f_L}{\partial \tau }+\frac{\partial F_L}{\partial \lambda }+\frac{\partial G_L}{\partial \xi }=0,---\left(7\right)$

$\frac{\partial f_g}{\partial \tau }+\frac{\partial F_g}{\partial \lambda }+\frac{\partial G_g}{\partial \xi }=0,---\left(8\right)$

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

$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 \theta \\ p\cos \theta \\ 0\end{array}\right),---\left(9\right)$

$F_g=\left(\begin{array}{c}0\\ 0\end{array}\right),$$G_g=\left(\begin{array}{c}-\mathrm{hu}\\ -\mathrm{hv}\end{array}\right),---\left(10\right)$

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

$\theta ={\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{\partial F_L}{\partial Q},$$C_L=\frac{\partial G_L}{\partial Q},---\left(14\right)$

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

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

$K_L^{\left(1\right)}=\left(\begin{array}{c}1\\ 0\\ 0\\ 0\end{array}\right),$(16)

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

$\left(\begin{array}{c}{\mu }_{L_{\mathrm{1,2}}}=0,\\ {\mu }_{L_{\mathrm{3,4}}}=\pm 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 \theta /\rho \\ \sin \theta /\rho \\ 0\end{array}\right),$$K_L^{\left(3\right)}=\left(\begin{array}{c}1\\ -\sin \theta \left(\frac{a}{\rho }\right)\\ \cos \theta \left(\frac{a}{\rho }\right)\\ a^2\end{array}\right),$$K_L^{\left(4\right)}=\left(\begin{array}{c}1\\ \sin \theta \left(\frac{a}{\rho }\right)\\ -\cos \theta \left(\frac{a}{\rho }\right)\\ a^2\end{array}\right),---\left(18\right)$

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

沿着特征方程可被写为

$\mathrm{dp}\pm \mathrm{\rho 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{\partial f}{\partial \tau }+\frac{\partial F}{\partial \lambda }=0.---\left(23\right)$

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

$f_{i,j}^{n+1}=f_{i,j}^n-\frac{\mathrm{\Delta \tau }}{\mathrm{\Delta \lambda }}[L_i^{+}{f}_{i,j}^n-L_{i-1}^{+}{f}_{i-1,j}^n+L_{i+1}^{-}{f}_{i+1,j}^n-L_i^{-}{f}_{i,j}^n],---\left(24\right)$

其中n是时间步长的序号；i，j是网格单元的序号；Δτ是时间步长；Δλ是空间步长。

上式中分裂的矩阵表示为

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

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

${\Lambda }^{+}=\left(\begin{array}{cccc}{\mu }_{L_1}& & & \\ & {\mu }_{L_2}& & \\ & & {\mu }_{L_3}& \\ & & & 0\end{array}\right)$and${\Lambda }^{-}=\left(\begin{array}{cccc}0& & & \\ & 0& & \\ & & 0& \\ & & & {\mu }_{L_4}\end{array}\right);---\left(26\right)$

以及矩阵

和它的逆矩阵

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

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

$\frac{\partial f}{\partial \tau }+\frac{\partial G\left(f\right)}{\partial \xi }=0,$$f=\left(\begin{array}{cc}{f}_L,& \xi <0,\\ f_R,& \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}{cc}{p}_{*}+C_L·f(u_{*},v_{*})=p_L+C_L·f(u_L,v_L),& \left(31\right)\\ p_{*}-C_R·f(u_{*},v_{*})=p_R-C_R·(u_R,v_R),& \left(32\right)\\ {\rho }_{*L}+f(u_{*},v_{*})({\rho }_L/a_L)={\rho }_L+f(u_L,v_L)({\rho }_L/a_L),& \left(33\right)\\ {\rho }_{*R}-f(u_{*},v_{*})({\rho }_R/a_R)={\rho }_R-f(u_R,v_R)({\rho }_R/a_R),& \left(34\right)\end{array}\right)$

其中C_L＝ρ_Lα_L，C_R＝ρ_Rα_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_{\Delta }(V,\theta )=\frac{V}{2}(\mathrm{tg\theta }+\cos \theta \ln \left((1+\sin \theta \right))+\frac{\cos \theta }{2}V\ln V,---\left(36\right)$

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

$\left(\begin{array}{cc}{p}_{*}=\frac{C_R{p}_L+C_L{p}_R+C_L{C}_R[f(u_L,v_L)-f(u_R,v_R)]}{C_L+C_R},& \left(37\right)\\ f(u_{*},v_{*})=\frac{C_L·f(u_L,v_L)+C_R·f(u_R,v_R)+[p_L-p_R]}{C_L+C_R},& \left(38\right)\\ {\rho }_{*L}={\rho }_L+(f(u_L,v_L)-f(u_{*},v_{*}))({\rho }_L/a_L),& \left(39\right)\\ {\rho }_{*R}={\rho }_R+(f(u_{*},v_{*})-f(u_R,v_R))({\rho }_R/a_R).& \left(40\right)\end{array}\right)$

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

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

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

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

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

tgθ_*＝ε，       (43)

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

$u_{*}=V_{*}\cos {\theta }_{*}=\frac{V_{*}}{\sqrt{1+{ϵ}^2}}$和$v_{*}=V_{*}\sin {\theta }_{*}=\frac{ϵV_{*}}{\sqrt{1+{ϵ}^2}},---\left(44\right)$

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

$F\left(ϵ\right)=\frac{V_{*}}{2}(ϵ+\frac{\ln V_{*}+\ln (\sqrt{1+{ϵ}^2}+ϵ)-\ln \sqrt{1+{ϵ}^2}}{\sqrt{1+{ϵ}^2}})-B=0.---\left(45\right)$

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

$F^{\prime }\left(ϵ\right)=\frac{V_{*}}{2}\left(\frac{(2+{ϵ}^2)\sqrt{1+{ϵ}^2}-ϵ-ϵ(\ln V_{*}+\ln (\sqrt{1+{ϵ}^2}+ϵ)-\ln \sqrt{1+{ϵ}^2})}{{(1+{ϵ}^2)}^{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关系可以用来建立 激波和对应的中间变量的关系。它们是，

μ_s(ρ_LJ_LE_L-ρ_*LJ_*LE_*L)＝0，      (49)

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

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

和

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

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

跨过膨胀波的焓不变关系

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

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

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

$V_{*L}=\sqrt{{2H}_L-\frac{2\gamma }{\gamma -1}\frac{p_{*}}{{\rho }_{*L}}}---\left(53\right)$

或$V_{*R}=\sqrt{{2H}_R-\frac{2\gamma }{\gamma -1}\frac{p_{*}}{{\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)

${\rho }_{*L}={\rho }_{*R}={\rho }_R+\frac{1}{2}(f(V_R,{\theta }_L)-f(V_R,{\theta }_R))({\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{{\rho }_R+a_R}{2}(f(V_R,{\theta }_L)-f(V_R,{\theta }_R)),---\left(70\right)$

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

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

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

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

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

tgθ_L＝ε，        (73)

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

$u_L=V_R\cos {\theta }_L=\frac{V_R}{\sqrt{1+{ϵ}^2}}$和$v_L=V_R\sin {\theta }_L=\frac{ϵV_R}{\sqrt{1+{ϵ}^2}},---\left(74\right)$

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

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

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

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

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

${ϵ}^{n+1}={ϵ}^n-\frac{F\left({ϵ}^n\right)}{F^{\prime }\left({ϵ}^n\right)}.---\left(77\right)$

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

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

4.附图说明

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

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

①${\left[Q\right]}_{i,j}^n,$②${\left[Q\right]}_{i\pm 1/2,j}^n,$③${\left[F\right]}_{i\pm 1/2,j}^n,$

④$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],$⑤${\left[Q\right]}_{i,j}^{n+1/3},$⑥$[Q{]}_{i,j\pm 1/2}^{n+1/3},$⑦${\left[G\right]}_{i\pm 1/2,j}^{n+1/3},$

⑧$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],$⑨${\left[Q\right]}_{i,j}^{n+2/3},$⑩${\left[Q\right]}_{i\pm 1/2,j}^{n+2/3},$${\left[F\right]}_{i\pm 1/2,j}^n,$

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

图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}-a{x}^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}·\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)找到新时刻的原始变量从守恒变量解码后，原始 变量为

$\rho =\frac{f_1}{J},$$u=\frac{f_2}{f_1},$$v=\frac{f_3}{f_1},$$p=\frac{(\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}\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}\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{th}}{e}^{-\beta x^2},-1\le x\le 1,---\left(84\right)$

其中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)计算固体壁面角度。固体壁面角度计算式根据镜像条件，所以有

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

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

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

如果δθ_w小于收敛标准，流场物理量和固体壁面的角度均达到最终的解；否则，继 续下一步。

(6)更新流场变量。这一步的任务是根据更新的固体壁面角度计算流场的物理参数从 更新到这一步包括上例中的步骤(1)～(7)。完成这一步后，获得流场新的物 理参数变量(ρ_R，p_R，V_R，θ_R，)。

(7)重复步骤(2)。

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