為什么我用FLUENT算的題總是發(fā)散??求大神?。 巨D(zhuǎn)發(fā)】

2017-03-30  by:CAE仿真在線  來(lái)源:互聯(lián)網(wǎng)

太長(zhǎng)不看版

對(duì)于很多可壓縮流動(dòng)問(wèn)題,如果入口邊界和出口邊界的壓力差別太大,導(dǎo)致計(jì)算發(fā)散,可以先算出一個(gè)入口邊界和出口邊界壓力差別小一些的流場(chǎng),然后以此作為計(jì)算原問(wèn)題的迭代初值,在此基礎(chǔ)上算出原問(wèn)題的解答。


作者簡(jiǎn)介

葉漢玉,本科畢業(yè)于北京航空航天大學(xué)飛行器動(dòng)力工程專業(yè)。現(xiàn)為北京航空航天大學(xué)在讀博士,研究方向?yàn)榱鲃?dòng)穩(wěn)定性,在Physics of fluids、Journal of fluid mechanics等流體力學(xué)著名刊物上發(fā)表文章數(shù)篇。


對(duì)于很多可壓縮流動(dòng)問(wèn)題,入口邊界和出口邊界的壓力差別很大。例如在計(jì)算拉伐爾噴管的流動(dòng)時(shí),噴管入口的絕對(duì)壓力常常比出口絕對(duì)壓力大一到兩個(gè)數(shù)量級(jí)。但是,這樣給數(shù)值計(jì)算帶來(lái)了困難。由于流場(chǎng)中絕對(duì)壓力的變化非常劇烈,計(jì)算前的均勻的初始?jí)毫?chǎng)往往與實(shí)際差別很大,這使得計(jì)算的時(shí)候容易發(fā)散。

我們來(lái)看一個(gè)例子。我們來(lái)計(jì)算一個(gè)拉伐爾噴管內(nèi)部的流動(dòng)以及其排氣羽流。這類流動(dòng)常見(jiàn)的例子是噴氣發(fā)動(dòng)機(jī)(航空發(fā)動(dòng)機(jī)、火箭發(fā)動(dòng)機(jī))噴管的流動(dòng)。拉伐爾噴管是一種先收縮后擴(kuò)張的噴管,氣流在收縮段從亞聲速加速到聲速,壓力下降,在擴(kuò)張段壓力進(jìn)一步下降,流動(dòng)加速到超聲速(圖1)。在噴管出口,由于氣流的壓力與外界環(huán)境壓力往往不相等,因此常常會(huì)呈現(xiàn)由激波、膨脹波組成的復(fù)雜的波系結(jié)構(gòu),而且這種波系會(huì)重復(fù)若干次,從外觀上能很明顯地看出來(lái)(圖2),直到遠(yuǎn)處才因?yàn)檎承院纳⒍饾u消失。

為什么我用FLUENT算的題總是發(fā)散??求大神?。 巨D(zhuǎn)發(fā)】fluent分析圖片1

圖1 拉伐爾噴管內(nèi)部的流動(dòng)及其排氣羽流(過(guò)膨脹狀態(tài))

為什么我用FLUENT算的題總是發(fā)散??求大神??!【轉(zhuǎn)發(fā)】fluent圖片2

圖2 AIM-120空空導(dǎo)彈的排氣羽流

(圖片來(lái)源:Wikipedia)

我們所計(jì)算的噴管的形狀如圖3(a)所示,圖中的單位為mm。噴管的出口面積為0.008m2,喉部面積為0.002m2。計(jì)算域如圖3(b)所示,在噴管出口的下游放置一個(gè)半徑約等于10倍噴管出口半徑,長(zhǎng)度約等于40倍噴管出口半徑的圓柱形區(qū)域,以便容納羽流流場(chǎng)。計(jì)算的時(shí)候使用二維軸對(duì)稱模型(Axisymmetric)。

為什么我用FLUENT算的題總是發(fā)散??求大神!!【轉(zhuǎn)發(fā)】fluent圖片3

(a)噴管局部

為什么我用FLUENT算的題總是發(fā)散??求大神??!【轉(zhuǎn)發(fā)】fluent圖片4

(b)計(jì)算域

圖3 計(jì)算域和噴管幾何形狀

噴管入口總壓為1000kPa,總溫為500K,環(huán)境壓力為10kPa;工質(zhì)為空氣。通過(guò)拉伐爾噴管的一維模型理論分析[1]可以知道,這時(shí)噴管處于欠膨脹工作狀態(tài),在噴管出口截面上超聲速氣流的壓力大于外界反壓(環(huán)境壓力),氣流會(huì)在出口處產(chǎn)生膨脹波。

在FLUENT 14.5中計(jì)算這個(gè)問(wèn)題。計(jì)算所使用的網(wǎng)格(圖4)為分塊結(jié)構(gòu)化網(wǎng)格,用ICEM CFD生成,網(wǎng)格數(shù)量約為2萬(wàn)。邊界(1)為噴管入口,使用pressure-inlet邊界條件,總壓(Gauge Total Pressure)設(shè)為1000kPa,初始靜壓(Supersonic/Initial Gauge Pressure)設(shè)為984840Pa,總溫(Total Temperature)設(shè)為500K,湍流條件按照湍流強(qiáng)度(Turbulent Intensity)等于5%、水力直徑(Hydraulic Diameter)等于0.1m設(shè)定。初始靜壓僅是為了在計(jì)算前初始化流場(chǎng)使用。通過(guò)一維模型理論分析可以知道噴管的流量大約為3.61kg/s,與之對(duì)應(yīng)的噴管入口流速約為66m/s,因此可以推算出噴管入口靜壓的近似值。邊界(2)為對(duì)稱軸,使用axis邊界條件。邊界(3)、(4)為壁面,使用wall邊界條件。邊界(5)、(6)為出口,使用pressure-outlet邊界條件。湍流模型使用k-ω SST。使用基于密度(Density-Based)的求解器,求解算法選取默認(rèn)的隱式算法。工質(zhì)的狀態(tài)方程用完全氣體模型(ideal-gas),工作壓力(Operating pressure)設(shè)為0。

為什么我用FLUENT算的題總是發(fā)散??求大神!!【轉(zhuǎn)發(fā)】fluent圖片5

圖4 網(wǎng)格

首先我們按照一般的方法來(lái)嘗試一下,即把邊界(5)、(6)的壓力直接設(shè)為10kPa(回流的湍流條件設(shè)定為與入口邊界相同的數(shù)值)。我們用入口邊界的變量的數(shù)值初始化整個(gè)流場(chǎng)(在“SolutionInitialization”頁(yè)面上,在“Compute from”組合框中選擇入口邊界,如圖5所示)。從初始化的流場(chǎng)開(kāi)始迭代的時(shí)候,我們使用默認(rèn)的空間離散格式(連續(xù)方程、動(dòng)量方程和能量方程(三者合稱Flow)使用二階迎風(fēng)格式,湍動(dòng)能(Turbulent Kinetic Energy)和比耗散率(Specific Dissipation Rate)使用一階迎風(fēng)格式)。求解的Courant數(shù)保持默認(rèn)值(即5)。

為什么我用FLUENT算的題總是發(fā)散??求大神?。 巨D(zhuǎn)發(fā)】fluent培訓(xùn)的效果圖片6

圖5 用入口邊界的變量的數(shù)值初始化整個(gè)流場(chǎng)

很不幸的是,計(jì)算發(fā)散了(圖6)。根據(jù)FLUENT的用戶手冊(cè)[2],剛開(kāi)始計(jì)算的時(shí)候,為了使得計(jì)算穩(wěn)定可以調(diào)小Courant數(shù)的值,并使用一階迎風(fēng)格式。因此我們嘗試將Courant數(shù)減小到1,并將Flow的離散格式改成一階迎風(fēng)格式,但是仍然無(wú)濟(jì)于事。計(jì)算發(fā)散的原因是這個(gè)流動(dòng)問(wèn)題中,壓力的變化范圍太大而且非常豐富。在噴管入口,流動(dòng)馬赫數(shù)很低,壓力接近于流動(dòng)的總壓(1000kPa);在噴管喉部,流動(dòng)馬赫數(shù)等于1,壓力降低到約為總壓的0.53倍(臨界壓力比[1]);在噴管擴(kuò)張段,流動(dòng)加速到超聲速,壓力進(jìn)一步下降;在噴管出口截面,按照一維模型理論分析可知壓力降低到約30kPa,流動(dòng)馬赫數(shù)約為3;然后,氣流在出口處通過(guò)膨脹波繼續(xù)降壓,最終達(dá)到與環(huán)境壓力(10kPa)一致。對(duì)于這樣復(fù)雜的壓力變化,從初始的均勻的壓力場(chǎng)開(kāi)始迭代顯然是過(guò)于困難了。

為什么我用FLUENT算的題總是發(fā)散??求大神??!【轉(zhuǎn)發(fā)】fluent培訓(xùn)的效果圖片7

圖6 計(jì)算發(fā)散

為了克服這種困難,我們可以使用逐次降低出口邊界壓力的方法。在使用一階迎風(fēng)格式并將Courant數(shù)減小到1的條件下,我們先把出口邊界(5)、(6)的壓力設(shè)為100kPa,計(jì)算收斂后,再將出口邊界壓力改為10kPa,然后再次計(jì)算收斂。這實(shí)際上就是用出口壓力100kPa的計(jì)算結(jié)果作為出口壓力10kPa計(jì)算時(shí)的迭代初值。這樣分開(kāi)兩次計(jì)算,每次計(jì)算時(shí)迭代初值與計(jì)算結(jié)果的差別都比較小,因此計(jì)算就不容易發(fā)散了。圖7是出口邊界壓力設(shè)為100kPa時(shí)的馬赫數(shù)云圖。圖8是將壓力改成10kPa后的結(jié)果。

為什么我用FLUENT算的題總是發(fā)散??求大神??!【轉(zhuǎn)發(fā)】fluent培訓(xùn)的效果圖片8

圖7 出口邊界壓力為100kPa。一階迎風(fēng)格式。流動(dòng)在噴管出口通過(guò)斜激波增壓到環(huán)境壓力。

為什么我用FLUENT算的題總是發(fā)散??求大神??!【轉(zhuǎn)發(fā)】fluent培訓(xùn)的效果圖片9

圖8 出口邊界壓力為10kPa。一階迎風(fēng)格式。流動(dòng)在噴管出口通過(guò)膨脹波減壓到環(huán)境壓力。

最終,我們得到了出口邊界壓力為10kPa的收斂的解答。最后我們可以將空間離散格式改成二階迎風(fēng)格式,算出最終的結(jié)果。從馬赫數(shù)云圖(圖9)可以清楚地看出噴管出口下游重復(fù)的波系結(jié)構(gòu)。當(dāng)然,這里的重點(diǎn)是說(shuō)明避免計(jì)算發(fā)散的技巧,因此采用了較粗的網(wǎng)格,也沒(méi)有進(jìn)一步做網(wǎng)格無(wú)關(guān)性驗(yàn)證。

為什么我用FLUENT算的題總是發(fā)散??求大神??!【轉(zhuǎn)發(fā)】fluent分析案例圖片10

圖9 最終的計(jì)算結(jié)果的馬赫數(shù)云圖

在其它可壓縮流動(dòng)或者復(fù)雜流動(dòng)(如帶有空化的流動(dòng))的模擬中,如果遇到計(jì)算發(fā)散,也可以嘗試本文的技巧。而且,有時(shí)邊界上的壓力分成兩步來(lái)設(shè)置可能還嫌少,可能要分成好幾步,逐次地降低(或者升高)壓力。

其它FLUENT版本也可以參考本文。

作者非常感謝北京航空航天大學(xué)宇航學(xué)院童曉艷老師。作者正是在和她的討論中了解到本文所述的避免計(jì)算發(fā)散的技巧。

轉(zhuǎn)自:流體那些事兒


開(kāi)放分享:優(yōu)質(zhì)有限元技術(shù)文章,助你自學(xué)成才

相關(guān)標(biāo)簽搜索:為什么我用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電磁分析 

編輯
在線報(bào)名:
  • 客服在線請(qǐng)直接聯(lián)系我們的客服,您也可以通過(guò)下面的方式進(jìn)行在線報(bào)名,我們會(huì)及時(shí)給您回復(fù)電話,謝謝!
驗(yàn)證碼

全國(guó)服務(wù)熱線

1358-032-9919

廣州公司:
廣州市環(huán)市中路306號(hào)金鷹大廈3800
電話:13580329919
          135-8032-9919
培訓(xùn)QQ咨詢:點(diǎn)擊咨詢 點(diǎn)擊咨詢
項(xiàng)目QQ咨詢:點(diǎn)擊咨詢
email:kf@1cae.com