《用拉格朗日形式的欧拉方程求解一类反问题的数值方法.pdf》由会员分享，可在线阅读，更多相关《用拉格朗日形式的欧拉方程求解一类反问题的数值方法.pdf（21页珍藏版）》请在专利查询网上搜索。
1、(10)申请公布号 CN 102880588 A (43)申请公布日 2013.01.16 C N 1 0 2 8 8 0 5 8 8 A *CN102880588A* (21)申请号 201210366939.0 (22)申请日 2011.02.15 201080001554.3 2011.02.15 G06F 17/10(2006.01) (71)申请人天津空中代码工程应用软件开发有 限公司 地址 300384 天津市华苑产业区华天道2号 火炬大厦辅楼302室 (72)发明人路明 (54) 发明名称 用拉格朗日形式的欧拉方程求解一类反问题 的数值方法 (57) 摘要 本发明是关于一种数值方。
2、法，它提供并求解 一种新的拉格朗日形式的二维欧拉方程来求解固 体壁面几何形状设计的反问题。该发明提供了一 个推导拉格朗日平面上的欧拉方程的变换方式， 从而简化计算网格、最大程度降低了对流项的数 值耗散。使用本发明的方法可以同时得到流场物 理量的解和固体壁面几何形状的设计的解。 (62)分案原申请数据 (51)Int.Cl. 权利要求书1页 说明书15页 附图4页 (19)中华人民共和国国家知识产权局 (12)发明专利申请 权利要求书 1 页 说明书 15 页 附图 4 页 1/1页 2 1.一种用拉格朗日形式的欧拉方程求解一类反问题的数值方法，其特征是，包括以下 步骤： (1)初始化流场变量参。
3、数和固体壁面角度； (2)录入规定的分布在固体壁面上的压力值； (3)求解反黎曼问题以找到物体壁面边界的镜像变量； (4)计算固体壁面角度； (5)检测固体壁面角度的收敛性，如果收敛则计算结束； (6)更新流场变量参数； (7)重复步骤(2)。 2.根据权利要求1所述的一种用拉格朗日形式的欧拉方程求解一类反问题的数值方 法，其中所述的固体壁面上的反黎曼问题，其特征是：以固体壁面为对称，形成计算域一侧 的壁面网格单元中的变量及其镜像变量，其中一侧形成以激波或者是膨胀波形式存在的左 侧变量或右侧变量，在其中间是中间变量；中间变量又可被分为左中间变量和右中间变量。 3.根据权利要求1所述的一种用；格。
4、朗日形式的欧拉方程求解一类反问题的数值方 法，其中所述的找到物体壁面边界的镜像变量，包括以下步骤： (1)沿着流函数形式的欧拉方程的特征方程积分，将左侧变量或右侧变量与中间变量 连接，其中，左侧变量、右侧变量、中间变量由权利要求2给出； (2)在中间变量中恢复流动速度的幅值； (3)求解组合函数以找到中间变量的流动角，组合函数在已知压力分布的情况下，表示 为 其中，p R ， R ， R ，V R ， R 是在固体壁面边界上的已知 的流动参数；p w 是规定的固体壁面的压力分布； L 是镜像中的流动角度； (4)在中间变量中求解速度分量。 权 利 要 求 书CN 102880588 A 1/1。
5、5页 3 用拉格朗日形式的欧拉方程求解一类反问题的数值方法 1. 技术领域 0001 本发明涉及一种模拟亚音速无粘流的流动和求解一类反问题的数值方法。这类反 问题是指在亚音速流场中的空气动力学物体的几何形状的设计问题。本发明属于计算流体 力学(CFD-Computational Fluid Dynamics)领域。 2. 背景技术 0002 2.1 先前的工作 0003 在计算流体力学的计算中，大多数是求解欧拉形式的流体控制方程。这意味着在 笛卡尔坐标下的欧拉平面中，计算网格必须根据物体的限定事先生成。网格构成网格单 元。由于流体穿过网格单元的交界面，所以存在对流项的通量。正是这个通量在数值解。
6、中 引起数值耗散，因为数值耗散直接与对对流项的数值逼近所引起的误差有关。从上世纪以 来，CFD研究者致力于开发更精确、高效率的数值方法来降低这个数值耗散。迎风差分方法 在求解流体的流动过程中取得显著成效，因为它合理表达了对流项的特征。典型代表有， Godunov方法1，它在网格单元交界面求解黎曼问题，给出非常精确的解；FVS方法2 (通量分裂法)在网格单元交界面运用特征关系，使求解过程更加快捷。但是，作为以欧拉 形式中不可避免的对流项的数值耗散，仍然存在于这些基于特征的数值方法中。 0004 另一方面，流体的拉格朗日描述强调流体颗粒在不同位置的运动和特点。对于流 体的控制方程，例如，拉格朗日形。
7、式的欧拉方程，其中，必须有一个方向代表流线，数值上可 以用流函数表示。另一个方向是流体颗粒运动的距离。这个坐标系统构成了拉格朗日平面， 在这个平面上，计算网格点理论上就是流体颗粒，网格线总是简单的直角网格。特别是稳定 的流动中，流体迹线和流线是重合的，没有流体颗粒穿过流线，穿过网格单元交界面的对流 项不存在。所以在拉格朗日平面上求解方程，数值耗散被降低至最小。 0005 应用拉格朗日形式的方程在求解一类反问题时体现了优势。在空气动力学中，一 类典型的反问题是通过给定固体壁面上的压力分布，然后设计固体壁面形状以符合压力分 布的要求。如果用基于欧拉平面的方法求解这类问题，例如adjoint法3(共。
8、轭方程法)， 流程是首先估计这个未被确定的几何形状；然后在其周围生成网格；在对流场求解，下一 步，是重要的、费时的，即求解共轭方程，改进几何形状。这个过程反复进行，直到找到目标 为止。通常，这个过程持续较长时间。拉格朗日形式十分适合应用在这类几何形状不确定 的问题中，因为在拉格朗日平面，不确定的几何边界，也就是固体壁面，也是由流线表示的， 而且，无论几何形状怎样变化，由于沿着一条流线的流函数是常数，所以流函数表示的流线 在拉格朗日平面是直线。在拉格朗日平面求解这类几何形状的设计问题可以达到最优(最 有效)的过程。 0006 尽管拉格朗日形式表现了如此卓越的优势，以前只是应用于超音速流动中4。当。
9、 欧拉方程在拉格朗日平面上被求解，即方程在空间上向下游方向发展，不需要考虑任何下 游对上游的影响，这一切完美地符合超音速流动的特点。以前，拉格朗日形式也成功地应用 于二维超音速流动中的固体壁面的几何形状的设计中5。到目前为止，应用流函数作为坐 说 明 书CN 102880588 A 2/15页 4 标进行求解几何形状设计的反问题的例子仅限于有势流和线性化的可压缩流6。 0007 2.2 目的和优势 0008 在工业界有明显的需求去应用拉格朗日形式的优势。例如在产品设计的初级阶 段，需要对产品进行最小数值耗散的数值模拟研究和快速的几何形状设计。如机翼、喷管的 流场计算、外形的设计等。亚音速也是工。
10、业界最常见的流动状态。用严格的拉格朗日概念 求解亚音速流动会遇到障碍，因为物体的存在会对上游的流体颗粒产生影响。成功应用拉 格朗日形式的数值方法的关键是如何正确使上游的流体颗粒感受到这种影响波动。本发明 的目的在于(1)提供一个拉格朗日形式的欧拉方程，它在一个坐标方向上没有对流项，从 而最大程度降低数值耗散；(2)提供精确求解拉格朗日形式的欧拉方程的数值方法；(3)提 供一个在亚音速流场中进行求解几何形状设计的反问题的最优方法。 3. 发明内容 0009 3.1 二维拉格朗日形式的欧拉方程 0010 本发明首先提供一个从欧拉平面的欧拉变量(t，x，y)到拉格朗日平面上的拉格 朗日变量(，)的坐。
11、标变换，其中变量为格拉朗日时间，和时间项t有相同的量 纲，被引入作为时间历遍项，另外两个独立的变量分别是流函数(量纲m 2 s -1 )和拉格 朗日距离(量纲m)。坐标和流体颗粒的流线重合，被定义为不同流体颗粒在流 线上的位置。本发明开始于描述二维、无粘、可压缩流体流动的欧拉平面上的欧拉方程 0011 0012 其中f是守恒变量矢量；F，G是对流项通量在两个笛卡尔方向上的分量。具体地， 0013 0014 变量t，u，v，和p分别是时间、流动速度V在两个笛卡尔方向上的分量、流体密 度和压力。总能E和总焓分别是H 0015 and 0016 上式中的是气体比热比。从坐标(t，x，y)到(，)的坐。
12、标变换的雅各 比矩阵可写为 0017 说 明 书CN 102880588 A 3/15页 5 0018 且它的绝对值J成为 0019 0020 其中h是拉格朗日行为参数；是流动方向角；被称作拉格朗日几何状态变 量，单位量纲是m 2 skg -1 ，被定义为 0021 0022 最终，二维不稳定欧拉方程(1)被转换成拉格朗日平面上的拉格朗日形式，它包 括两组偏微分方程 0023 0024 0025 其中f L ，F L 和G L 与(2)中的f，F，G称谓相同，是在拉格朗日平面。具体地， 0026 0027 0028 加上(5)中J的定义，且 0029 0030 说 明 书CN 10288058。
13、8 A 4/15页 6 0031 0h1。 (13) 0032 方程(7)被称作拉格朗日物理守恒方程；(8)是几何守恒方程。在(9)中，可以发 现矢量G L 的第一和第四个元素是零，这意味着在拉格朗日平面上，沿着方向，对于连续 方程和能量方程没有通量穿过各单元的交界面，所以数值耗散被减小。此外，拉格朗日形式 的欧拉方程还有其他特性。它们是， 0033 1.方程(7)满足齐次特性，这意味着FVS方法可被用来求解通量项。 0034 2.方程(7)不满足旋转不变性，这意味着在不同的坐标方向上应该用不同的数值 方法求解对流项通量，以保证各自的特征不被模糊掉。 0035 3.偏微分方程(7)和(8)都是。
14、强双曲线型的，这意味着对于初值问题，方程可以用 时间历遍的方法来求解。 0036 3.2拉格朗日形式的欧拉方程的特征和特征方程 0037 对于拉格朗日物理守恒方程(7)，可以找到其特征(特征值和特征向量)。按照雅 各比矩阵的定义， 0038 0039 其中Q，u，v，p T 。雅各比矩阵BL的四个特征值是， 0040 0041 及其对应的、沿着方向的四个特征向量是， 0042 (16) 0043 雅各比矩阵C L 的四个特征值是， 0044 0045 及其对应的、沿着方向的四个特征向量是， 说 明 书CN 102880588 A 5/15页 7 0046 0047 其中a是声速。沿着拉格朗日物。
15、 理守恒方程(7)的特征方程可被写为 0048 0049 沿着特征方程可被写为 0050 0051 根据坐标变换的雅各比矩阵(4)，在拉格朗日平面沿着方向的黎曼不变量为 0052 0053 沿着方向的为 0054 0055 方程(15)(22)被用来构造求解控制方程(7)和(8)的数值方法。 0056 3.3 数值方法的构造 0057 为在拉格朗日平面为求解方程(7)和(8)，本发明采用有限容积法，变量存储在网 格单元的中心。根据方程的齐次、旋转、双曲线特性，同时考虑到计算的效率和精度，本发明 说 明 书CN 102880588 A 6/15页 8 产生一种新的求解方法，带有混合迎风差分格式通。
16、量算子的、Strang换方向的数值方法，意 味着，在计算网格交界面的对流项通量时，在方向用FVS方法，在方向用求解黎曼问 题的Godunov方法，体现了Godunov方法和FVS方法的各自优势。Strang换方向是二阶精 度时间离散。从现在起，为方便起见，方程和变量中省去拉格朗日含义的下标L。 0058 3.3.1 沿着方向的FVS法 0059 首先，将方程(7)写成仅沿着方向的形式 0060 0061 通过应用FVS(通量矢量分裂)法，可以获得各网格单元(i，j)中的更新的守恒变 量 0062 0063 其中n是时间步长的序号；i，j是网格单元的序号；是时间步长；是空间 步长。 0064 上。
17、式中分裂的矩阵表示为 0065 和 0066 其中，分裂的特征值矩阵为 0067 and 0068 以及矩阵 0069 0070 和它 的逆矩阵 0071 说 明 书CN 102880588 A 7/15页 9 0072 0073 3.3.2 沿着方向求解黎曼问题的Godunov法 0074 对于Godunov方法，核心内容是求解黎曼问题。本发明沿着方向求解黎曼问题， 被称作跨过流线的黎曼问题。图1表示了跨过流线的黎曼问题的定义。将方程(7)写成仅 沿着方向的形式， 0075 0076 其中f L 和f R 分别是网格交界面左侧和右侧变量中的守恒变量。从现在开始，下标 L、R和*分别代表黎曼问。
18、题中的左、右和中间变量。如图1所示，左侧变量或右侧变量是激 波或者是膨胀波，在其中间是中间变量，又可被接触波分为左中间变量和右中间变量。时间 更新的解可以从方程(29)的离散方程获得， 0077 0078 下面的任务是通过f找到原始变量p，u，v，再求得网格交界面上j1/2的G。 所述的原始变量也是黎曼问题中的中间变量。一个黎曼问题求解器按照以下流程导出。首 先沿着拉格朗日物理守恒方程的特征方程(19)和(20)积分，将左侧变量和右侧变量与中 间变量连接，于是有 0079 0080 其中C L L L ，C R R R ，并且 0081 0082 f(u，v)被称作组合函数，它可以被认为是速度。
19、分量(u，v)在拉格朗日平面上的组 合，与速度有同样的量纲m/s。公式(35)有它的极坐标形式 说 明 书CN 102880588 A 8/15页 10 0083 0084 其中uV cos and vV sin。方程(31)(34)的解表达 了黎曼问题中的中间变量的值 0085 0086 为发现速度分量(u * ，v * )，需要求解方程(38)。首先可将(38)写为 0087 F(u * ，v * )f(u * ，v * )-B0， (41) 0088 其中，考虑了已知的左侧和右侧变量，于是有 0089 0090 组合函数(35)可以被进一步改写。定义 0091 tg * ， (43) 0。
20、092 其中 * 是左侧或右侧中间变量中的流动方向角，于是有 0093 和 0094 其中V * 是左侧或右侧中间变量中的速度幅值。进而，方程(41)可以最终写成 的函数 0095 0096 下面的工作是通过给定的速度V * 求解方程F()0以发现，再求出 * 。方 程(45)是可微分的，F的一阶导数是 0097 0098 数值试验表明，对于无量纲速度V * 5(亚音速区间)而且-1，1的范围 内，F()0，这意味着函数F()在此范围内是单调的。同时，F(-)F()0， 所以Newton-Raphson历遍方法可用来找到方程(45)在给定V * 时的根。为找到V * 的值， 需要用p * ， 。
21、*L ， *R 的值去还原黎曼问题中的各种非线性波(接触波、膨胀波、激波)。本 发明中考虑了下面的关系： 0099 跨过激波的Rankine-Hugoniot关系 说 明 书CN 102880588 A 10 9/15页 11 0100 如果激波出现在一侧(例如，左侧)，跨过激波，Rankine-Hugoniot关系可以用来 建立激波和对应的中间变量的关系。它们是， 0101 0102 s ( L J L E L - *L J *L E *L )0， (49) 0103 其中 s 激波速度，J *L 和E *L 分别是中间变量中的J(坐标变换雅各比矩阵绝对值， 见公式(5)和E(总能)。从(4。
22、6)和(49)，可以获得 0104 0105 和 0106 0107 速度绝对值V *L 或V *R 可以被带入到公式(45)以求解 *L 或 *R 。 0108 跨过膨胀波的焓不变关系 0109 如果膨胀波出现在一侧(例如，左侧)，跨过膨胀波，焓不变关系可以用来建立膨 胀波和对应的中间变量的关系。它们是， 0110 0111 从这个公式可以求得速度幅值 0112 0113 或 0114 具体那一种非线性波出现，是要靠黎曼问题中的左侧、右侧和中间变量中的压力 值判断的，例如，假设 0115 P min min(p L ，p R )； (55) 0116 p max max(p L ，p R )。
23、， (56) 0117 及其p * ，构成了下面的波形选择条件： 0118 (1)如果p min p * P max ，意味着在黎曼问题中一个激波出现在p min 一侧，一个膨 胀波出现在P max 一侧； 0119 (2)如果p * p max ，意味着在黎曼问题中两个激波出现； 0120 (3)如果p * p min ，意味着在黎曼问题中两个膨胀波出现。 说 明 书CN 102880588 A 11 10/15页 12 0121 在获得 * 之后，按照公式(9)和由公式(44)求得的u *L ，v *L 或u *R ，v *R ，在加上p * ， 便可求得公式(30)中的G L 。同时，几。
24、何守恒方程也可被更新至新的时间， 0122 0123 和欧拉形式的方法相比，尽管这个步骤是额外的，但这个步骤简单，不占用过多的 计算时间。 0124 3.3.3 边界条件 0125 处理边界条件需要围绕在计算域周围再加上两层网格，被称为虚拟网格，它可以 使空间离散算法扩展至边界。 0126 固体壁面边界条件 0127 构造的数值方法中，需要沿着方向在物体壁面求解黎曼问题。物体壁面上的黎 曼问题包括左侧或右侧变量和针对物体壁面的镜像变量。在物体壁面边界条件中，下标L 表示物理变量，p，u，v，在虚拟网格中，下标R表示这些物理变量在壁面网格中。 w 是物体壁面的角度。根据(37)(40)给出的黎曼。
25、问题的解，在物体壁面边界上有 0128 p * p R ， (58) 0129 0130 中间变量中的压力即可被认为是物体壁面上的压力p w ，因而 0131 p w p * 。 (60) 0132 如果物体壁面上的压力已经被给定，这个给定压力可以直接被用做黎曼问题中的 中间变量中的压力，所以 0133 p * p w (61) 0134 出口、入口边界条件 0135 第一个关系对应的是正向、进入的特征线。以下标b表示的边界上的值可以通过 式(21)、(22)中的黎曼不变量获得， 0136 0137 0138 其中，u ，a 是来流的速度和声速，下标i表示计算域内部的值。第二个关系是通 过内部。
26、的值的插值获得的数值逼近。因而，速度u b 可以通过式(62)和(63)获得，于是有 0139 0140 相似的方法可以处理亚音速出口边界条件。但需要在出口边界规定外界压力值 说 明 书CN 102880588 A 12 11/15页 13 p ，有 0141 0142 拉格朗日几何状态变量的边界条件 0143 在物体壁面边界上， 0144 0145 0146 对于入口、出口边界，规定流动角度0，则 0147 0148 0149 图2给出了从更新到的一个步骤的数值方法的流程图，这个方法在空间 和时间上都是二阶精度，其中Q，u，v，p T 是流场的物理变量。 0150 3.4 求解一类反问题的方。
27、法 0151 使用拉格朗日形式时，固体壁面由流函数代表的流线表示，所以在拉格朗日平面， 由流线代表的网格线是直线。本发明的另一个目的是探索拉格朗日形式的数值方法在求解 一类反问题的优势。这类反问题被定义为几何形状的设计问题。 0152 第2部分引入了拉格朗日平面上沿着方向的黎曼问题的解和固体壁面边界条 件。从(37)和(57)，可以获得 0153 0154 其中p R ， R ，a R ，V R ， R 是边界上已知的量，f是公式(35)和(36)给出的组合函 数。从(70)推导出 0155 0156 该式在规定壁面压力分布p w 后可以求解出虚拟网格中流动角度 L 的大小。根据 镜像边界条件。
28、，固体避免的角度可以从下式求得 0157 0158 几何形状的设计问题意味着要根据给定的固体壁面压力分布，找到固体壁面的角 度。这个过程需要求解一个反黎曼问题。一个反黎曼问题的求解器按照下面的步骤建立， 首先定义 0159 tg L ， (73) 0160 然后，虚拟网格中的速度分量按下式求得 说 明 书CN 102880588 A 13 12/15页 14 0161 和 0162 然后被带入三角形式给出的组合函数(71)中，并被进一步简化为的函数， 0163 0164 下面的任务是求解规定物理参数 R ，p R ，V R ， R 和p w 条件下的方程F()0，解 得后再求得 L 。方程(7。
29、5)是可微分的，F针对的一阶导数是 0165 0166 相同于求解方程(45)的方法，更新的值表示为 0167 0168 虚拟网格中的流动角度可由 L tg -1 ( n+1 )求得。固体壁面角度 w 按式(72) 求得。 0169 图3给出了求解几何外形的设计的反问题的算法流程图。这个方法起始于假设的 流场物理量的参数值和固体壁面角度。按固体壁面上按照的指定的压力分布，求解反黎曼 问题，在虚拟网格单元中得出流动角度，然后壁面形状可获得并被带回到流场求解器中求 解更新流场物理量。这个过程持续进行，直到历遍的固体壁面的角度达到收敛为止。很明 显，收敛标准是壁面角度而不是任何一个物理变量，这意味着。
30、，壁面角度也变成了一个流场 变量。壁面边界的变化并不改变方程的特征，方程仍然是强双曲线型，仍然是一个初值问 题。与应用在欧拉平面上的方法(如adjoint法)相比，现在的方法既不需要在更新的几 何形状周围生成网格，也不需要求解adjoint方程，它同时获得收敛的流场物理变量的解 和几何形状的解。针对这类问题，总体的CPU时间降低一个量级。 4. 附图说明 0170 图1拉格朗日平面上沿方向跨过流线的黎曼问题 0171 图2使用二阶精度方法，带有混合迎风差分格式通量算子的、Strang换方向的数 值方法进行一步流场更新的流程图，其中的标号代表以下变量或公式 0172 0173 0174 0175。
31、 说 明 书CN 102880588 A 14 13/15页 15 0176 图3几何形状设计的反问题的计算流程图 0177 图4入口马赫数为0.5的抛物型型喷管的流场的拉格朗日方法的解和欧拉方法的 解 0178 (a)拉格朗日平面上的网格 0179 (b)欧拉平面上的网格 0180 (c)拉格朗日方法解得的流线 0181 (d)欧拉方法解得的流线 0182 (e)拉格朗日方法解得的固体壁面上压力分布和精解的比较 0183 (f)拉格朗日方法构造的网格 0184 图5根据给定压力分布计算的固体壁面的几何形状 0185 (a)欧拉平面上的网格 0186 (b)固体壁面上压力分布 0187 (c)。
32、用本发明的方法计算的固体壁面形状 5. 具体实施方式 0188 5.1 扩张喷管中的无粘亚音速流动 0189 按照本发明介绍的方法，下面给出一个完整的用拉格朗日形式的欧拉方程模拟无 粘、亚音速内流。这个例子是关于亚音速流通过一个二维抛物线形、长度为L、入口高度为 H in 的扩张型喷管。其几何尺寸是由两段抛物线定义的， 0190 0191 其中H in L/3，aH in /2，bH in /L。无粘流在喷管入口的马赫数是M in 0.5。流 动是纯亚音速。拉格朗日形式的欧拉方程(7)和(8)将用第二节介绍的方法求解，其中，带 有混合迎风差分格式通量算子在空间上是二阶精度，采用带有minmod。
33、通量限制器的MUSCL 插值。Strang换方向在时间上是二阶精度。计算参数如下：CFL0.48；拉格朗日行为参 数h0.98；总共6020个计算网格单元在拉格朗日平面。计算结果将和基于欧拉平面 上的有限体积法的JST方法7得到的结果及其精解做比较。从现在起，定义本发明给出 的解为拉格朗日解，由欧拉方法在欧拉平面上的到的解为欧拉解。参考图2，具体地讲，计算 采用以下步骤， 0192 (1)初始化。流场变量Q 0 0 ，u 0 ，v 0 ，p 0 T 及均为初始值，其中Q 0 取无 穷远出的值，上标0表示流动问题在初始时刻，在x-y平面和-平面，t 0(0)。进而可求得守恒变量f 0 。然后在-。
34、平面生成以均匀的空间步长(， )为单元的直角网格。对应的x-y平面上的网格可按照下式生成 0193 dxcosd和dysind (79) 0194 (2)计算时间步长。按照变量在时间n(n0，1，2，.)，CFL0.5， 0195 说 明 书CN 102880588 A 15 14/15页 16 0196 其中， n-1 从上一个时刻获得，所以有 0197 0198 (3)沿着-方向用FVS方法求解方程(23)。通过已知的和 值，守恒变量从更新至对于二阶精度的数值方法，使用带有TVD通量限制器的 MUSCL插值。 0199 (4)沿着-方向用Godunov方法求解方程(29)。用n(n0，1，。
35、2，.)时刻的 已知的f n+1/2 和Q n+1/2 的值，用Godunov方法更新守恒变量。当地的黎曼问题求解器给出网格 单元交界面上的原始变量，u，v，p i，j1/2 值。对于二阶精度的数值方法，使用带有TVD 通量限制器的MUSCL插值。 0200 (5)更新拉格朗日几何变量。因为网格交界面上的速度已知，更新的拉 格朗日几何变量可以由方程(57)获得。精度等级由速度决定，尽管公示(57)表 示的是一阶精度。 0201 (6)找到新时刻的原始变量从守恒变量解码后， 原始变量为 0202 0203 (6)构造网格。如果需要在每一次历遍迭代的最后一步构造网格。网格点由下式 给出， 0204。
36、 0205 每一个网格点式代表一个由计算网格单元表示的流体颗粒。本发明的原理讲述的 是，在拉格朗日平面构造的网格线其实就是流线本身。 0206 (7)重复步骤(2)。在完成一次时间步长的更新，回到步骤(2)，重复步骤(3)-(7)， 直到守恒变量或是原始变量收敛。 0207 图4中，拉格朗日的解起始于图4(a)中表示的由拉格朗日距离和流函数构 成的直角网格。图4(c)中，拉格朗日方法获得的流线与图4(d)中欧拉方法获得的流线吻合 良好，图4(e)表示固体壁面上的压力分布与精解完全一致。图4(f)中由拉格朗日生成的 流线与图4(b)略有不同，沿着y方向的网格线不垂直于流线方向，以保证流线的长度。。
37、因 该指出的是，图4(f)中表示的最终的网格是由图4(a)中的网格进化而来，尽管这个最终的 网格并不需要，但是这个进化过程可以清楚说明拉格朗目方法的原理。 0208 5.2 几何形状设计的反问题 0209 下面的例子是介绍用求解二维拉格朗日形式的欧拉方程去解决固体壁面几何形 说 明 书CN 102880588 A 16 15/15页 17 状的设计问题，是一类反问题。例子中是一个收缩-扩张喷管。例子中先规定几何形状，求 得固体壁面压力分布，然后用着个压力分布作为输入，用本发明的方法求得固体壁面的几 何形状，结果应该与最初的几何形状一致。喷管固体壁面外形由下面的函数给出 0210 0211 其中。
38、H th 取0.1，表明收缩喉口高度和入口高度之比为10。取6.5以控制喉 口长度约占喷管长度的三分之一。入口马赫数为0.5，是纯亚音速流动。喷管的几何形状 和欧拉平面上的计算网格如图5(a)所示。拉格朗日平面上的网格与上例中图中4(a)所示 一致。拉格朗日形式的欧拉方程(7)和(8)将用第二节介绍的方法求解，其中，带有混合 迎风差分格式通量算子在空间上是二阶精度，采用带有minmod通量限制器的MUSCL插值。 Strang换方向在时间上是二阶精度。计算参数如下：CFL0.48；拉格朗日行为参数h 0.98；两组网格数目6020和9030个计算网格单元在拉格朗日平面。首先按照上例中 给出的计。
39、算流程完成计算，获得固体壁面压力分布，如图5(b)所示，并作为求解几何形状 设计的反问题的压力分布的输入。同时，adjoint方法将同时应用这个例子以检验本发明 的方法的计算效率。具体步骤是： 0212 (1)初始化流场和固体壁面角度。流场初始变量的设置于上例相同。固体壁面初 始角度假设 w 0。因为流场中的流动角度设为0，所以虚拟网格单元中的流动角度 L 0。 0213 (2)输入规定的固体壁面的压力分布p w 。 0214 (3)在固体壁面上求解基于给定压力分布p w 的反黎曼问题，公式按照第二节 (70)(77)。 0215 (4)计算固体壁面角度。固体壁面角度计算式根据镜像条件，所以有。
40、 0216 0217 (5)检验固体壁面角度的收敛性。固体壁面角度 w 的变化量是两个时刻n+1和 n的差值， 0218 0219 如果 w 小于收敛标准，流场物理量和固体壁面的角度均达到最终的解；否 则，继续下一步。 0220 (6)更新流场变量。这一步的任务是根据更新的固体壁面角度计算流场的物 理参数从更新到这一步包括上例中的步骤(1)(7)。完成这一步后，获得流场 新的物理参数变量( R ，p R ，V R ， R ，)。 0221 (7)重复步骤(2)。 0222 图5(c)表示了计算的固体几何形状和原始的形状的比较。从结果可以估计计算 精度，计算的固体几何形状和原始的形状吻合良好，最大误差在之内0.5。最后，在计算效 率方面，本发明的方法使用的计算的时间是adjoint方法的十分之一。 说 明 书CN 102880588 A 17 1/4页 18 图1 图2 说 明 书 附 图CN 102880588 A 18 2/4页 19 图3 说 明 书 附 图CN 102880588 A 19 3/4页 20 图4 说 明 书 附 图CN 102880588 A 20 4/4页 21 图5 说 明 书 附 图CN 102880588 A 21 。