专利名称：求解拉格朗日形式的欧拉方程的数值方法
技术领域：
本发明涉及一种模拟亚音速无粘流的流动和求解一类反问题的数值方法。这类反问题是指在亚音速流场中的空气动力学物体的几何形状的设计问题。本发明属于计算流体力学(CFD-Computational Fluid Dynamics)令页域。
2.
背景技术：
2. 1先前的工作在计算流体力学的计算中，大多数是求解欧拉形式的流体控制方程。这意味着在笛卡尔坐标下的欧拉平面中，计算网格必须根据物体的限定事先生成。网格构成网格单元。由于流体穿过网格单元的交界面，所以存在对流项的通量。正是这个通量在数值解中引起数值耗散，因为数值耗散直接与对对流项的数值逼近所引起的误差有关。从上世纪以来，CFD研究者致力于开发更精确、高效率的数值方法来降低这个数值耗散。迎风差分方法在求解流体的流动过程中取得显著成效，因为它合理表达了对流项的特征。典型代表有， Godimov方法[1]，它在网格单元交界面求解黎曼问题，给出非常精确的解；FVS方法[2] (通量分裂法)在网格单元交界面运用特征关系，使求解过程更加快捷。但是，作为以欧拉形式中不可避免的对流项的数值耗散，仍然存在于这些基于特征的数值方法中。另一方面，流体的拉格朗日描述强调流体颗粒在不同位置的运动和特点。对于流体的控制方程，例如，拉格朗日形式的欧拉方程，其中，必须有一个方向代表流线，数值上可以用流函数表示。另一个方向是流体颗粒运动的距离。这个坐标系统构成了拉格朗日平面， 在这个平面上，计算网格点理论上就是流体颗粒，网格线总是简单的直角网格。特别是稳定的流动中，流体迹线和流线是重合的，没有流体颗粒穿过流线，穿过网格单元交界面的对流项不存在。所以在拉格朗日平面上求解方程，数值耗散被降低至最小。应用拉格朗日形式的方程在求解一类反问题时体现了优势。在空气动力学中，一类典型的反问题是通过给定固体壁面上的压力分布，然后设计固体壁面形状以符合压力分布的要求。如果用基于欧拉平面的方法求解这类问题，例如adjoint法[3](共轭方程法)， 流程是首先估计这个未被确定的几何形状；然后在其周围生成网格；在对流场求解，下一步，是重要的、费时的，即求解共轭方程，改进几何形状。这个过程反复进行，直到找到目标为止。通常，这个过程持续较长时间。拉格朗日形式十分适合应用在这类几何形状不确定的问题中，因为在拉格朗日平面，不确定的几何边界，也就是固体壁面，也是由流线表示的， 而且，无论几何形状怎样变化，由于沿着一条流线的流函数是常数，所以流函数表示的流线在拉格朗日平面是直线。在拉格朗日平面求解这类几何形状的设计问题可以达到最优(最有效)的过程。尽管拉格朗日形式表现了如此卓越的优势，以前只是应用于超音速流动中W]。当欧拉方程在拉格朗日平面上被求解，即方程在空间上向下游方向发展，不需要考虑任何下游对上游的影响，这一切完美地符合超音速流动的特点。以前，拉格朗日形式也成功地应用于二维超音速流动中的固体壁面的几何形状的设计中[5]。到目前为止，应用流函数作为坐标进行求解几何形状设计的反问题的例子仅限于有势流和线性化的可压缩流[6]。2. 2目的和优势在工业界有明显的需求去应用拉格朗日形式的优势。例如在产品设计的初级阶段，需要对产品进行最小数值耗散的数值模拟研究和快速的几何形状设计。如机翼、喷管的流场计算、外形的设计等。亚音速也是工业界最常见的流动状态。用严格的拉格朗日概念求解亚音速流动会遇到障碍，因为物体的存在会对上游的流体颗粒产生影响。成功应用拉格朗日形式的数值方法的关键是如何正确使上游的流体颗粒感受到这种影响波动。本发明的目的在于(1)提供一个拉格朗日形式的欧拉方程，它在一个坐标方向上没有对流项，从而最大程度降低数值耗散；( 提供精确求解拉格朗日形式的欧拉方程的数值方法；(3)提供一个在亚音速流场中进行求解几何形状设计的反问题的最优方法。
3.
发明内容
3. 1 二维拉格朗日形式的欧拉方程本发明首先提供一个从欧拉平面的欧拉变量(t，χ, y)到拉格朗日平面上的拉格朗日变量(τ，λ, ξ)的坐标变换，其中变量τ为格拉朗日时间，和时间项t有相同的量纲，被引入作为时间历遍项，另外两个独立的变量分别是流函数ξ (量纲[πι2·。])和拉格朗日距离λ (量纲[m])。坐标ξ和流体颗粒的流线重合，λ被定义为不同流体颗粒在流向上的位置。本发明开始于描述二维、无粘、可压缩流体流动的欧拉平面上的欧拉方程
权利要求
1. 一种模拟无粘、亚音速流动的数值方法，其特征是，包括以下步骤(1)将二维欧拉平面上的欧拉方程变换为拉格朗日平面上的拉格朗日形式的欧拉方程(2)录入物体几何形状的表面坐标；(3)在所述物体周围建立计算网格；(4)用带有混合迎风差分格式通量算子的、Strang换方向的数值方法求解拉格朗日形式的欧拉方程。
2.根据权利要求1所述的方法，其中所述的拉格朗日形式的欧拉方程方程包括两个方，拉格朗日物理守恒方程和几何守恒方程t
3.根据权利要求1所述的方法，其中所述的拉格朗日平面有一个由拉格朗日时间τ表示的时间方向和由流函数ξ和拉格朗日距离λ表示的两个空间方向。
4.根据权利要求1所述的方法，其中所述的二维欧拉平面上的欧拉方程的变换包括一 个转换矩阵，
5.根据权利要求2所述，其中拉格朗日物理守恒方程表示为
6.根据权利要求2所述，其中几何守恒方程表示为
7.根据权利要求1所述的方法，其中所述的计算网格是拉格朗日平面上以λ和ξ为两个方向的二维直角网格。
8.根据权利要求1所述的方法，其中所述的在拉格朗日平面上求解拉格朗日形式的欧拉方程方程是沿着拉格朗日时间τ方向进行守恒变量的历遍，以找到其稳定解。
9.根据权利要求1所述的方法，其中所述的用^rang换方向法及其混合的通量运算方案求解拉格朗日形式的欧拉方程在守恒变量&沿着拉格朗日时间τ方向进行历遍时的每一步更新过程中包括以下步骤(1)用FVS方法沿λ方向以时间步长Δ τ来计算网格单元交界的对流项通量；(2)用求解跨过流线的黎曼问题的Godimov方法沿ξ方向以时间步长Δτ的二分之一来计算网格单元交界的对流项通量；(3)再次用FVS方法沿λ方向以时间步长Δτ来计算网格单元交界的对流项通量。
10.根据权利要求9所述的方法，其中所述沿ξ方向的跨过流线的黎曼问题包括在计算单元交界面左侧变量或右侧变量是激波或者是膨胀波，在其中间是中间变量，又可被分为左中间变量和右中间变量。
11.根据权利要求9所述的方法，其中所述的求解跨过流线的黎曼问题包括以下步骤(1)沿着拉格朗日物理守恒方程的特征方程积分，将左侧变量或右侧变量与中间变量连接；(2)在中间变量中恢复流动速度的幅值；(3)求解组合函数以找到中间变量的流动角；(4)在中间变量中求解速度分量。
12.根据权利要求11所述的方法，其中所述的在中间变量中恢复流动速度的幅值，是利用通过跨过激波的Rankine-Hugoniot关系和跨过膨胀波的焓不变关系来获得的。
13.根据权利要求11所述的方法，其中所述的组合函数f(u，ν)表示为
14.一种求解固体壁面几何形状设计的反问题的数值方法，其特征是，包括以下步骤(1)初始化流场变量参数和固体壁面角度；(2)录入规定的分布在物体壁面上的压力值；(3)求解反黎曼问题以找到物体壁面边界的镜像变量；(4)计算固体壁面角度；(5)检测固体壁面角度的收敛性，如果收敛则计算结束；(6)更新流场变量参数；(7)重复步骤(3)。
15.根据权利要求14所述的方法，其中所述的求解反黎曼问题包括求解一个已知压力分布的组合函数，该组合函数表示为
16.根据权利要求14所述的方法，其中所述的求解反黎曼问题包括在固体壁面上用 Newton-Raphson历遍的方法求解权利要求15中所述的已知压力分布的组合函数，以找到流动角度θρ
全文摘要
这个发明是关于一种数值方法，它提供并求解一种新的拉格朗日形式的二维欧拉方程来模拟亚音速流动和求解固体壁面几何形状设计的反问题。该发明提供了一个推导拉格朗日平面上的欧拉方程的变换方式，从而简化计算网格、最大程度降低了对流项的数值耗散。该发明构造了一个混合的通量运算方案，其中，二维黎曼问题在拉格朗日平面得以求解。使用本发明的方法可以同时得到流场物理量的解和固体壁面几何形状的设计的解。
文档编号G06F17/50GK102203782SQ201080001554
公开日2011年9月28日 申请日期2010年9月9日 优先权日2010年9月9日
发明者路明 申请人:天津空中代码工程应用软件开发有限公司