搜索為什么我用FLUENT算的題總是發(fā)散？？求大神??！【轉(zhuǎn)發(fā)】
2017-03-30 by:CAE仿真在線 來源:互聯(lián)網(wǎng)
太長不看版
對于很多可壓縮流動問題,如果入口邊界和出口邊界的壓力差別太大,導(dǎo)致計算發(fā)散,可以先算出一個入口邊界和出口邊界壓力差別小一些的流場,然后以此作為計算原問題的迭代初值,在此基礎(chǔ)上算出原問題的解答。
作者簡介
葉漢玉,本科畢業(yè)于北京航空航天大學(xué)飛行器動力工程專業(yè)?，F(xiàn)為北京航空航天大學(xué)在讀博士,研究方向為流動穩(wěn)定性,在Physics of fluids、Journal of fluid mechanics等流體力學(xué)著名刊物上發(fā)表文章數(shù)篇。
對于很多可壓縮流動問題,入口邊界和出口邊界的壓力差別很大。例如在計算拉伐爾噴管的流動時,噴管入口的絕對壓力常常比出口絕對壓力大一到兩個數(shù)量級。但是,這樣給數(shù)值計算帶來了困難。由于流場中絕對壓力的變化非常劇烈,計算前的均勻的初始壓力場往往與實際差別很大,這使得計算的時候容易發(fā)散。
我們來看一個例子。我們來計算一個拉伐爾噴管內(nèi)部的流動以及其排氣羽流。這類流動常見的例子是噴氣發(fā)動機(航空發(fā)動機、火箭發(fā)動機)噴管的流動。拉伐爾噴管是一種先收縮后擴張的噴管,氣流在收縮段從亞聲速加速到聲速,壓力下降,在擴張段壓力進一步下降,流動加速到超聲速(圖1)。在噴管出口,由于氣流的壓力與外界環(huán)境壓力往往不相等,因此常常會呈現(xiàn)由激波、膨脹波組成的復(fù)雜的波系結(jié)構(gòu),而且這種波系會重復(fù)若干次,從外觀上能很明顯地看出來(圖2),直到遠處才因為粘性耗散而逐漸消失。
圖1 拉伐爾噴管內(nèi)部的流動及其排氣羽流(過膨脹狀態(tài))
圖2 AIM-120空空導(dǎo)彈的排氣羽流
(圖片來源:Wikipedia)
我們所計算的噴管的形狀如圖3(a)所示,圖中的單位為mm。噴管的出口面積為0.008m2,喉部面積為0.002m2。計算域如圖3(b)所示,在噴管出口的下游放置一個半徑約等于10倍噴管出口半徑,長度約等于40倍噴管出口半徑的圓柱形區(qū)域,以便容納羽流流場。計算的時候使用二維軸對稱模型(Axisymmetric)。
(a)噴管局部
(b)計算域
圖3 計算域和噴管幾何形狀
噴管入口總壓為1000kPa,總溫為500K,環(huán)境壓力為10kPa;工質(zhì)為空氣。通過拉伐爾噴管的一維模型理論分析[1]可以知道,這時噴管處于欠膨脹工作狀態(tài),在噴管出口截面上超聲速氣流的壓力大于外界反壓(環(huán)境壓力),氣流會在出口處產(chǎn)生膨脹波。
在FLUENT 14.5中計算這個問題。計算所使用的網(wǎng)格(圖4)為分塊結(jié)構(gòu)化網(wǎng)格,用ICEM CFD生成,網(wǎng)格數(shù)量約為2萬。邊界(1)為噴管入口,使用pressure-inlet邊界條件,總壓(Gauge Total Pressure)設(shè)為1000kPa,初始靜壓(Supersonic/Initial Gauge Pressure)設(shè)為984840Pa,總溫(Total Temperature)設(shè)為500K,湍流條件按照湍流強度(Turbulent Intensity)等于5%、水力直徑(Hydraulic Diameter)等于0.1m設(shè)定。初始靜壓僅是為了在計算前初始化流場使用。通過一維模型理論分析可以知道噴管的流量大約為3.61kg/s,與之對應(yīng)的噴管入口流速約為66m/s,因此可以推算出噴管入口靜壓的近似值。邊界(2)為對稱軸,使用axis邊界條件。邊界(3)、(4)為壁面,使用wall邊界條件。邊界(5)、(6)為出口,使用pressure-outlet邊界條件。湍流模型使用k-ω SST。使用基于密度(Density-Based)的求解器,求解算法選取默認的隱式算法。工質(zhì)的狀態(tài)方程用完全氣體模型(ideal-gas),工作壓力(Operating pressure)設(shè)為0。
圖4 網(wǎng)格
首先我們按照一般的方法來嘗試一下,即把邊界(5)、(6)的壓力直接設(shè)為10kPa(回流的湍流條件設(shè)定為與入口邊界相同的數(shù)值)。我們用入口邊界的變量的數(shù)值初始化整個流場(在“SolutionInitialization”頁面上,在“Compute from”組合框中選擇入口邊界,如圖5所示)。從初始化的流場開始迭代的時候,我們使用默認的空間離散格式(連續(xù)方程、動量方程和能量方程(三者合稱Flow)使用二階迎風(fēng)格式,湍動能(Turbulent Kinetic Energy)和比耗散率(Specific Dissipation Rate)使用一階迎風(fēng)格式)。求解的Courant數(shù)保持默認值(即5)。
圖5 用入口邊界的變量的數(shù)值初始化整個流場
很不幸的是,計算發(fā)散了(圖6)。根據(jù)FLUENT的用戶手冊[2],剛開始計算的時候,為了使得計算穩(wěn)定可以調(diào)小Courant數(shù)的值,并使用一階迎風(fēng)格式。因此我們嘗試將Courant數(shù)減小到1,并將Flow的離散格式改成一階迎風(fēng)格式,但是仍然無濟于事。計算發(fā)散的原因是這個流動問題中,壓力的變化范圍太大而且非常豐富。在噴管入口,流動馬赫數(shù)很低,壓力接近于流動的總壓(1000kPa);在噴管喉部,流動馬赫數(shù)等于1,壓力降低到約為總壓的0.53倍(臨界壓力比[1]);在噴管擴張段,流動加速到超聲速,壓力進一步下降;在噴管出口截面,按照一維模型理論分析可知壓力降低到約30kPa,流動馬赫數(shù)約為3;然后,氣流在出口處通過膨脹波繼續(xù)降壓,最終達到與環(huán)境壓力(10kPa)一致。對于這樣復(fù)雜的壓力變化,從初始的均勻的壓力場開始迭代顯然是過于困難了。
圖6 計算發(fā)散
為了克服這種困難,我們可以使用逐次降低出口邊界壓力的方法。在使用一階迎風(fēng)格式并將Courant數(shù)減小到1的條件下,我們先把出口邊界(5)、(6)的壓力設(shè)為100kPa,計算收斂后,再將出口邊界壓力改為10kPa,然后再次計算收斂。這實際上就是用出口壓力100kPa的計算結(jié)果作為出口壓力10kPa計算時的迭代初值。這樣分開兩次計算,每次計算時迭代初值與計算結(jié)果的差別都比較小,因此計算就不容易發(fā)散了。圖7是出口邊界壓力設(shè)為100kPa時的馬赫數(shù)云圖。圖8是將壓力改成10kPa后的結(jié)果。
圖7 出口邊界壓力為100kPa。一階迎風(fēng)格式。流動在噴管出口通過斜激波增壓到環(huán)境壓力。
圖8 出口邊界壓力為10kPa。一階迎風(fēng)格式。流動在噴管出口通過膨脹波減壓到環(huán)境壓力。
最終,我們得到了出口邊界壓力為10kPa的收斂的解答。最后我們可以將空間離散格式改成二階迎風(fēng)格式,算出最終的結(jié)果。從馬赫數(shù)云圖(圖9)可以清楚地看出噴管出口下游重復(fù)的波系結(jié)構(gòu)。當(dāng)然,這里的重點是說明避免計算發(fā)散的技巧,因此采用了較粗的網(wǎng)格,也沒有進一步做網(wǎng)格無關(guān)性驗證。
圖9 最終的計算結(jié)果的馬赫數(shù)云圖
在其它可壓縮流動或者復(fù)雜流動(如帶有空化的流動)的模擬中,如果遇到計算發(fā)散,也可以嘗試本文的技巧。而且,有時邊界上的壓力分成兩步來設(shè)置可能還嫌少,可能要分成好幾步,逐次地降低(或者升高)壓力。
其它FLUENT版本也可以參考本文。
作者非常感謝北京航空航天大學(xué)宇航學(xué)院童曉艷老師。作者正是在和她的討論中了解到本文所述的避免計算發(fā)散的技巧。
轉(zhuǎn)自:流體那些事兒
相關(guān)標簽搜索：為什么我用FLUENT算的題總是發(fā)散？？求大神??！【轉(zhuǎn)發(fā)】 Fluent培訓(xùn) Fluent流體培訓(xùn) Fluent軟件培訓(xùn) fluent技術(shù)教程 fluent在線視頻教程 fluent資料下載 fluent分析理論 fluent化學(xué)反應(yīng) fluent軟件下載 UDF編程代做 Fluent、CFX流體分析 HFSS電磁分析
