專利名稱：求解拉格朗日形式的歐拉方程的數(shù)值方法
技術領域：
本發(fā)明涉及一種模擬亞音速無粘流的流動和求解一類反問題的數(shù)值方法。這類反問題是指在亞音速流場中的空氣動力學物體的幾何形狀的設計問題。本發(fā)明屬于計算流體力學(CFD-Computational Fluid Dynamics)令頁域。
2.
背景技術：
2. 1先前的工作在計算流體力學的計算中，大多數(shù)是求解歐拉形式的流體控制方程。這意味著在笛卡爾坐標下的歐拉平面中，計算網(wǎng)格必須根據(jù)物體的限定事先生成。網(wǎng)格構成網(wǎng)格單元。由于流體穿過網(wǎng)格單元的交界面，所以存在對流項的通量。正是這個通量在數(shù)值解中引起數(shù)值耗散，因為數(shù)值耗散直接與對對流項的數(shù)值逼近所引起的誤差有關。從上世紀以來，CFD研究者致力于開發(fā)更精確、高效率的數(shù)值方法來降低這個數(shù)值耗散。迎風差分方法在求解流體的流動過程中取得顯著成效，因為它合理表達了對流項的特征。典型代表有， Godimov方法[1]，它在網(wǎng)格單元交界面求解黎曼問題，給出非常精確的解；FVS方法[2] (通量分裂法)在網(wǎng)格單元交界面運用特征關系，使求解過程更加快捷。但是，作為以歐拉形式中不可避免的對流項的數(shù)值耗散，仍然存在于這些基于特征的數(shù)值方法中。另一方面，流體的拉格朗日描述強調流體顆粒在不同位置的運動和特點。對于流體的控制方程，例如，拉格朗日形式的歐拉方程，其中，必須有一個方向代表流線，數(shù)值上可以用流函數(shù)表示。另一個方向是流體顆粒運動的距離。這個坐標系統(tǒng)構成了拉格朗日平面， 在這個平面上，計算網(wǎng)格點理論上就是流體顆粒，網(wǎng)格線總是簡單的直角網(wǎng)格。特別是穩(wěn)定的流動中，流體跡線和流線是重合的，沒有流體顆粒穿過流線，穿過網(wǎng)格單元交界面的對流項不存在。所以在拉格朗日平面上求解方程，數(shù)值耗散被降低至最小。應用拉格朗日形式的方程在求解一類反問題時體現(xiàn)了優(yōu)勢。在空氣動力學中，一類典型的反問題是通過給定固體壁面上的壓力分布，然后設計固體壁面形狀以符合壓力分布的要求。如果用基于歐拉平面的方法求解這類問題，例如adjoint法[3](共軛方程法)， 流程是首先估計這個未被確定的幾何形狀；然后在其周圍生成網(wǎng)格；在對流場求解，下一步，是重要的、費時的，即求解共軛方程，改進幾何形狀。這個過程反復進行，直到找到目標為止。通常，這個過程持續(xù)較長時間。拉格朗日形式十分適合應用在這類幾何形狀不確定的問題中，因為在拉格朗日平面，不確定的幾何邊界，也就是固體壁面，也是由流線表示的， 而且，無論幾何形狀怎樣變化，由于沿著一條流線的流函數(shù)是常數(shù)，所以流函數(shù)表示的流線在拉格朗日平面是直線。在拉格朗日平面求解這類幾何形狀的設計問題可以達到最優(yōu)(最有效)的過程。盡管拉格朗日形式表現(xiàn)了如此卓越的優(yōu)勢，以前只是應用于超音速流動中W]。當歐拉方程在拉格朗日平面上被求解，即方程在空間上向下游方向發(fā)展，不需要考慮任何下游對上游的影響，這一切完美地符合超音速流動的特點。以前，拉格朗日形式也成功地應用于二維超音速流動中的固體壁面的幾何形狀的設計中[5]。到目前為止，應用流函數(shù)作為坐標進行求解幾何形狀設計的反問題的例子僅限于有勢流和線性化的可壓縮流[6]。2. 2目的和優(yōu)勢在工業(yè)界有明顯的需求去應用拉格朗日形式的優(yōu)勢。例如在產品設計的初級階段，需要對產品進行最小數(shù)值耗散的數(shù)值模擬研究和快速的幾何形狀設計。如機翼、噴管的流場計算、外形的設計等。亞音速也是工業(yè)界最常見的流動狀態(tài)。用嚴格的拉格朗日概念求解亞音速流動會遇到障礙，因為物體的存在會對上游的流體顆粒產生影響。成功應用拉格朗日形式的數(shù)值方法的關鍵是如何正確使上游的流體顆粒感受到這種影響波動。本發(fā)明的目的在于(1)提供一個拉格朗日形式的歐拉方程，它在一個坐標方向上沒有對流項，從而最大程度降低數(shù)值耗散；( 提供精確求解拉格朗日形式的歐拉方程的數(shù)值方法；(3)提供一個在亞音速流場中進行求解幾何形狀設計的反問題的最優(yōu)方法。
3.

發(fā)明內容
3. 1 二維拉格朗日形式的歐拉方程本發(fā)明首先提供一個從歐拉平面的歐拉變量(t，χ, y)到拉格朗日平面上的拉格朗日變量(τ，λ, ξ)的坐標變換，其中變量τ為格拉朗日時間，和時間項t有相同的量綱，被引入作為時間歷遍項，另外兩個獨立的變量分別是流函數(shù)ξ (量綱[πι2·。])和拉格朗日距離λ (量綱[m])。坐標ξ和流體顆粒的流線重合，λ被定義為不同流體顆粒在流向上的位置。本發(fā)明開始于描述二維、無粘、可壓縮流體流動的歐拉平面上的歐拉方程
權利要求
1. 一種模擬無粘、亞音速流動的數(shù)值方法，其特征是，包括以下步驟(1)將二維歐拉平面上的歐拉方程變換為拉格朗日平面上的拉格朗日形式的歐拉方程(2)錄入物體幾何形狀的表面坐標；(3)在所述物體周圍建立計算網(wǎng)格；(4)用帶有混合迎風差分格式通量算子的、Strang換方向的數(shù)值方法求解拉格朗日形式的歐拉方程。
2.根據(jù)權利要求1所述的方法，其中所述的拉格朗日形式的歐拉方程方程包括兩個方，拉格朗日物理守恒方程和幾何守恒方程t
3.根據(jù)權利要求1所述的方法，其中所述的拉格朗日平面有一個由拉格朗日時間τ表示的時間方向和由流函數(shù)ξ和拉格朗日距離λ表示的兩個空間方向。
4.根據(jù)權利要求1所述的方法，其中所述的二維歐拉平面上的歐拉方程的變換包括一 個轉換矩陣，
5.根據(jù)權利要求2所述，其中拉格朗日物理守恒方程表示為
6.根據(jù)權利要求2所述，其中幾何守恒方程表示為
7.根據(jù)權利要求1所述的方法，其中所述的計算網(wǎng)格是拉格朗日平面上以λ和ξ為兩個方向的二維直角網(wǎng)格。
8.根據(jù)權利要求1所述的方法，其中所述的在拉格朗日平面上求解拉格朗日形式的歐拉方程方程是沿著拉格朗日時間τ方向進行守恒變量的歷遍，以找到其穩(wěn)定解。
9.根據(jù)權利要求1所述的方法，其中所述的用^rang換方向法及其混合的通量運算方案求解拉格朗日形式的歐拉方程在守恒變量&沿著拉格朗日時間τ方向進行歷遍時的每一步更新過程中包括以下步驟(1)用FVS方法沿λ方向以時間步長Δ τ來計算網(wǎng)格單元交界的對流項通量；(2)用求解跨過流線的黎曼問題的Godimov方法沿ξ方向以時間步長Δτ的二分之一來計算網(wǎng)格單元交界的對流項通量；(3)再次用FVS方法沿λ方向以時間步長Δτ來計算網(wǎng)格單元交界的對流項通量。
10.根據(jù)權利要求9所述的方法，其中所述沿ξ方向的跨過流線的黎曼問題包括在計算單元交界面左側變量或右側變量是激波或者是膨脹波，在其中間是中間變量，又可被分為左中間變量和右中間變量。
11.根據(jù)權利要求9所述的方法，其中所述的求解跨過流線的黎曼問題包括以下步驟(1)沿著拉格朗日物理守恒方程的特征方程積分，將左側變量或右側變量與中間變量連接；(2)在中間變量中恢復流動速度的幅值；(3)求解組合函數(shù)以找到中間變量的流動角；(4)在中間變量中求解速度分量。
12.根據(jù)權利要求11所述的方法，其中所述的在中間變量中恢復流動速度的幅值，是利用通過跨過激波的Rankine-Hugoniot關系和跨過膨脹波的焓不變關系來獲得的。
13.根據(jù)權利要求11所述的方法，其中所述的組合函數(shù)f(u，ν)表示為
14.一種求解固體壁面幾何形狀設計的反問題的數(shù)值方法，其特征是，包括以下步驟(1)初始化流場變量參數(shù)和固體壁面角度；(2)錄入規(guī)定的分布在物體壁面上的壓力值；(3)求解反黎曼問題以找到物體壁面邊界的鏡像變量；(4)計算固體壁面角度；(5)檢測固體壁面角度的收斂性，如果收斂則計算結束；(6)更新流場變量參數(shù)；(7)重復步驟(3)。
15.根據(jù)權利要求14所述的方法，其中所述的求解反黎曼問題包括求解一個已知壓力分布的組合函數(shù)，該組合函數(shù)表示為
16.根據(jù)權利要求14所述的方法，其中所述的求解反黎曼問題包括在固體壁面上用 Newton-Raphson歷遍的方法求解權利要求15中所述的已知壓力分布的組合函數(shù)，以找到流動角度θρ
全文摘要
這個發(fā)明是關于一種數(shù)值方法，它提供并求解一種新的拉格朗日形式的二維歐拉方程來模擬亞音速流動和求解固體壁面幾何形狀設計的反問題。該發(fā)明提供了一個推導拉格朗日平面上的歐拉方程的變換方式，從而簡化計算網(wǎng)格、最大程度降低了對流項的數(shù)值耗散。該發(fā)明構造了一個混合的通量運算方案，其中，二維黎曼問題在拉格朗日平面得以求解。使用本發(fā)明的方法可以同時得到流場物理量的解和固體壁面幾何形狀的設計的解。
文檔編號G06F17/50GK102203782SQ201080001554
公開日2011年9月28日 申請日期2010年9月9日 優(yōu)先權日2010年9月9日
發(fā)明者路明 申請人:天津空中代碼工程應用軟件開發(fā)有限公司