一種減縮CFD非定常分析計算量的方法

楊 鳴，趙世平，王漢平，張煜冰，王邵助

(1.北京理工大學 宇航學院，北京 100081;2.北京機電工程研究所 11室，北京 1000781)

非定常氣動力精確計算是飛行器動氣彈和氣動伺服動力學問題得到準確分析前提［1-8］，也是了解飛行器操縱性、機動性和穩定性等飛行品質的關鍵。

針對非定常氣動力計算，傳統上往往引入準定常假設，假定在任何非定常運動瞬間飛行器氣動特性都與同一飛行器以等速度進行平移或轉動時所顯現的氣動特性一致，這實際上忽略了運動和頻率的依賴關系［9］。準定常假設在工程上簡單易用，但用于顫振和動響應分析時精度不足，在理論上也是站不住腳的。對于非定常運動的正確處理方式，是借助CFD分析中動網格技術建立流場中的動邊界模型，讓物體在流場中以真實的形態“運動”，這種指導思想下建立的模型在理論上才是完備的，在工程實踐上也具有更高精度。要讓運動物體在計算域中“運動”起來，需要預先構造足夠大且覆蓋物體運動空間的網格計算區域。但是，以飛行器為代表的航行體運動速度快，達到每秒數百米乃至上千米的數量級，即使僅僅只分析飛行器數秒時間范圍內的運動流場，在其主要運動方向上也必須構造成百上千，甚至上萬米的運動區域，如果還要進一步考慮飛行器的三維幾何效應，那么其所需計算資源將是現有絕大多數計算機都無法提供的。減縮計算量將是非定常運動問題分析一個重要課題。

本文抓住飛行器運動流場特征，即在遠離飛行器的流場區域應為未擾動遠場區域，提出在計算時只夠造部分計算域，對飛行器已經經過遠離且恢復為遠場條件的區域實時刪除，同時又對飛行器將要經過的遠場區域實時擴展以保證計算進行。這種處理方式使得計算域始終之只限制在飛行器周圍，無需構建起整個計算域，從而大大減少了計算量且對計算精度基本沒有影響。

1 方法原理

1.1 域動分層技術

對六面體或者楔形網格區域，可以使用運動分層技術來生成或者刪除邊界網格。當網格邊界和鄰近一層網格有相對運動時，若緊鄰邊界網格高度增長到一定高度，就將其劃分為2個網格層;反之，當網格層高度降低到一定程度時，就將緊鄰邊界的2層網格合并為一層網格以杜絕負體積網格的出現。

網格層擴大準則:

式(1)中:Hmin為網格最小高度;h0為理想高度;θs為分割因子，在滿足該條件前提下，就對網格單元執行分割操作。

網格層縮減準則:

式中:θc為合并因子，在滿足這個條件時，就將該層網格與外面一層網格合并［10］。

1.2 非定常流場特征與計算方案

飛行器在自由空間中運動，其流場呈現出如下特征:在靠近動壁面區域，速度場、壓力場和溫度場變化劇烈，呈現高度非線性、非定常特點，是流場計算主要關心區域［11-13］。而在遠離飛行器特征長度以外的區域(一般可取為特征長度20～30倍)，物理場變化微小，屬于平穩的未擾動流場區域，其流動狀態簡單且可以直接給定，其并非計算分析關心區域［14-15］。由于這個特點，如果將飛行器所有運動區域都預先構造計算網格，實際上就等于引入了多余的無關計算。而未擾動區域網格數占整個計算域網格數百分比很大，往往將達到30%～50%，如果能夠在計算時忽略掉這部分網格，那將大大縮減計算規模，且對計算結果影響極小。

利用這個特點，本文將飛行器流場計算分割為若干個階段，在計算一開始只構造一部分計算網格，遠離運動物體區域一律劃分為結構化網格。當運動物體在計算域中運動一段時間以后，在物體運動反方向上利用前述域動分層技術縮減網格，而在沿著物體運動方向則相應擴展網格以保證計算得以繼續。通過這種處理方式，計算量被大大減少，這對于必須利用動網格技術的CFD計算具有重要意義。

2 算例

2.1 基本模型的選取

以二維NACA0006翼型的非定常流場計算作為示例將傳統方法和新方法進行對比。湍流模型使用在空氣動力學計算中常用的Spalart-Allmaras一方程模型，該模型將湍動黏度μt表示為湍動能k的函數，其輸運方程為

流場邊界為遠場條件，由于本算例假設飛行器是以一定速度在相對慣性系靜止流場中運動，故來流馬赫數Ma=0，壓力為50662 Pa，溫度216.6 T。其中動區域和靜區域設置如圖1，①、②、③、④、⑤、⑥分別為動區域和靜區域的邊界，⑤、⑥在計算過程中為靜止邊界。

圖1 動區域和靜區域示意圖

機翼運動規律為沿x負方向以10 m/s2加速度從300 m/s速度在2 s時間內加速運動到320 m/s，為了減少初始速度對計算結果的影響，前0.2 s以300 m/s速度勻速運動，因此實際運動時間為2.2 s。由運動公式可知，機翼實際位移為680 m。

2.2 傳統計算方案網格布置

計算域尺寸如圖2所示，網格為四邊形結構化網格，總的網格數為501688個。

圖2 計算域大小示意圖

2.3 新的計算方案網格布置

初始計算域大小如圖3所示，機翼每前進100 m，就對網格區域進行一次“減縮-擴展”操作。進行這一操作時，應當將動區域和靜區域設置為可變形域，而將①、②、③、④、⑤、⑥6組邊界設置為剛體運動沿x負方向運動100 m，這樣既使得了計算域規模基本不變，又防止了內部區域隨同“減縮-擴展”操作時一起運動，保證不因為該操作破壞機翼附近區域正常的流場形態。整個計算域網格數為281796個，較之傳統計算方案減少了43.83%。

圖3 計算域大小示意圖

3 計算結果

在2臺配置相同計算機上計算，2組模型時間步長設置一致，迭代情況也大體一致。傳統計算方案計算時間為156 h左右，新方案計算加上“減縮-擴展”操作所耗費時間，一共耗時110 h左右，減少了30%左右的計算時間。

圖4～圖6為2種計算方案在計算終了時刻各個物理場對比等值線圖?？梢钥闯觯?種計算方案差異很小。

圖4 壓力等值線圖對比

圖5 速度等值線圖對比

圖6 溫度等值線圖對比

圖7為新計算方案得到的阻力與傳統計算方案之間的相對誤差隨時間變化曲線，圖8則為升力的相對誤差隨時間變化曲線。

圖7 阻力相對誤差

圖8 升力相對誤差

由相對誤差圖可以看出，阻力計算誤差比升力誤差要小一個數量級，造成升力誤差相對較大的原因在于計算域縱向尺寸取的還不夠大，從而放大了對網格進行“壓縮-擴展”操作以后造成的誤差，對此可以采取適當放寬計算域縱向尺寸辦法來加以克服。

4 結束語

結合飛行器非定常運動流場的特點，采用新的計算方案可以減少計算量，如果再采取其他相關措施，例如網格控制等手段，計算量的減縮量是很可觀的。和傳統計算方式相比，新的計算方案并不會降低精度。新的計算方法操作簡單，具備工程可實施性和易用性，對非定常流動問題計算具有指導意義。

［1］陳桂彬，鄒叢青，楊超.氣動彈性設計基礎［M］.北京：北京航空航天大學出版社，2007：64.

［2］Friedmann P P.Hypersonic aerothermoelastic characteristics of a finned missle［J］.Journal of Space craft，1979，16(3)：187-190.

［3］Lighthill M J.Oscillating airfoils at high mach numbers［J］.Journal of the Aeronautical Sciences，1953，20(6)：402-406.

［4］McNamara J J，Thuruthimattam B J，Friedmann P P，et al.Hypersonic aerothmoelastic studies for reusable launch vehicles［R］.AIAA-2004-1590，2004.

［5］Lind R.Linear parameter-varying modeling and control of structural dynamics with aerothermoelastic effects［R］.AIAA-1999-1393，1999.

［6］楊超，許贇，謝長川.高超聲速飛行器氣動彈性力學研究綜述［J］.航空學報，2010，31(1)：1-11.

［7］陸洋，王浩文，高正.電控旋翼氣彈動力學建模研究［J］.航空動力學報，2006，21(6)：1021-1026.

［8］康浩，高正.艦面直升機旋翼瞬態氣彈響應分析［J］.航空動力學報，2000，15(1)：67-70.

［9］J.R.賴特，J.E.庫珀.飛進氣動彈性力學及載荷導論［M］.姚一龍，譯.上海：上海交通大學出版社，2010：146.

［10］吳光中，宋婷婷，張毅.FLUENT基礎入門與案例精通［M］.北京：電子工業出版社，2012：329-330.

［11］伍貽兆，田書玲，夏健.基于非結構動網格的非定常流數值模擬方法［J］.航空學報，2011，32(1)：15-26.

［12］許和勇，葉正寅，王剛，等.基于非結構運動對接網格的旋翼前飛流場數值模擬［J］.空氣動力學報，2007，25(3)：325-328.

［13］喬陽，劉勇，徐敏，等.基于多塊對接網格的多側噴流瞬態干擾特性研究［J］.固體火箭技術，2007，30(21)：552-554.

［14］宋琦，楊樹興，徐勇，等.滾轉狀態下卷弧翼火箭彈氣動特性的數值模擬［J］.固體火箭技術，2008，31(6)：98-101.

［15］倪同兵，招啟軍，趙國慶，等.應用多塊對接結構網格方法的直升機涵道尾槳氣動特性［J］.航空動力學報，2011，29(6)：688-695.