FLUENT不收斂案例及解決方法：壓力出口導(dǎo)致的不收斂【轉(zhuǎn)發(fā)】

2017-06-25 by:CAE仿真在線來源:互聯(lián)網(wǎng)

今天分享一個FLUENT的不收斂案例及其解決方法。計算的對象是一個新型的渦扇發(fā)動機(jī)加力燃燒室(圖1)。在這種新型加力燃燒室中,火焰穩(wěn)定器被整合到整流支板上,因此整流支板和整流錐都需要冷卻。在整流支板和整流錐上開了很多小孔,冷卻氣從這些孔滲出,形成冷卻氣膜。

FLUENT不收斂案例及解決方法：壓力出口導(dǎo)致的不收斂【轉(zhuǎn)發(fā)】fluent結(jié)果圖片1

圖1 加力燃燒室

這個算例模擬的是實驗的工況。實驗中沒有在加力燃燒室內(nèi)燃燒,而只是在“主流入口”處引入高溫氣體,在“冷卻氣入口”處引入冷卻氣,以檢驗氣膜冷卻的效能。

整流支板共有15塊,為了減小計算量,只計算其中的一塊(360°/15=24°)。主流入口和冷卻氣入口都采用“mass-flow-inlet”條件,其中主流入口的流量是0.8kg/s,總溫是1241.3K,冷卻氣入口的流量是0.024kg/s,總溫是490.3K。出口采用“pressure-outlet”條件,反壓是101325Pa(絕對壓力)。由于形狀比較復(fù)雜,特別是其中有很多小孔,所以采用非結(jié)構(gòu)網(wǎng)格,四面體單元。流體的狀態(tài)方程采用理想氣體(ideal-gas)模型,湍流模型采用Realizable k-ε模型。采用基于壓力的求解器。

采用定常算法計算不收斂(圖2;這里我們使用FLUENT默認(rèn)的收斂條件,即能量方程的殘差降低到1e-6以下,其余方程降低到1e-3以下)?？紤]到可能是分離流誘發(fā)的非定常效應(yīng)導(dǎo)致不收斂(見公眾號先前的文章“為何我這個流動總是算不收斂?我要砸電腦!”),我們嘗試使用非定常算法。但是不幸的是非定常算法仍然不能在每個時間步內(nèi)收斂。非定常計算的典型殘差曲線如圖3所示。

FLUENT不收斂案例及解決方法：壓力出口導(dǎo)致的不收斂【轉(zhuǎn)發(fā)】fluent結(jié)果圖片2

圖2 定常計算的殘差曲線

FLUENT不收斂案例及解決方法：壓力出口導(dǎo)致的不收斂【轉(zhuǎn)發(fā)】fluent分析案例圖片3

圖3 非定常計算仍然不收斂。此圖是時間步長設(shè)為3×10-6秒時的結(jié)果。圖中顯示了5個時間步的殘差曲線。

為了弄清不收斂的原因,我們用MATLAB編寫兩個小程序。第一個程序用來產(chǎn)生一個命令文件j1.jou,讓FLUENT迭代30次,并把每次迭代后的流場都保存到文件里面:

clear fid=fopen('j1.jou','wt'); % 打開命令文件 n=30; fprintf(fid,'solve update-physical-time\n'); % 下一個時間步 for i=1:n fprintf(fid,'solve iterate 1\n'); % 迭代一次 fprintf(fid,'file interpolate write-data "d:\\case1-%d" yes yes \n',i); % 將計算結(jié)果保存到文件 end fclose(fid); % 關(guān)閉命令文件

在FLUENT中打開非定常計算的case和data,在菜單欄選擇[File]->[Read]->[Journal…],選取命令文件j1.jou,FLUENT便會更新到下一個時間步并迭代30次,并在D盤根目錄下產(chǎn)生30個計算結(jié)果文件(圖4)。

FLUENT不收斂案例及解決方法：壓力出口導(dǎo)致的不收斂【轉(zhuǎn)發(fā)】fluent分析案例圖片4

圖4 計算結(jié)果文件

第二個程序分析這些計算結(jié)果,找出壓力波動最劇烈的點:

clear n=30; m=375045; % 網(wǎng)格數(shù)量 n_points=15; fid1=fopen('j2.jou','wt'); % 打開命令文件 L=1; for i=27:n filename=sprintf('d:\\case1-%d.ip',i); fid2=fopen(filename); % 打開計算結(jié)果文件 % ASCII碼中,40代表左圓括號 % FLUENT的計算結(jié)果文件中,每一組數(shù)據(jù)都是以左圓括號開頭的 % 因此,可以通過查找左圓括號的方法找到每一組數(shù)據(jù)的起點 % 需要了解更多關(guān)于FLUENT計算結(jié)果文件格式的內(nèi)容,請 % 閱讀User's Guide中“Format of the Interpolation File”這一節(jié) % until函數(shù)需要自己編寫,見注1 until(fid2,40); arrx=fread(fid2,[m,1],'double'); % 讀取x坐標(biāo) until(fid2,40); arry=fread(fid2,[m,1],'double'); % 讀取y坐標(biāo) until(fid2,40); arrz=fread(fid2,[m,1],'double'); % 讀取z坐標(biāo) until(fid2,40); arrp=fread(fid2,[m,1],'double'); % 讀取壓力場 fclose(fid2); if i>27 % 只關(guān)注最后三次迭代 % 比較相鄰兩次迭代的結(jié)果,找出壓力波動最大的15個點 [sorted,ix]=sort(abs(arrp_old-arrp),'descend'); x_maxchange=arrx(ix(1:n_points)); y_maxchange=arry(ix(1:n_points)); z_maxchange=arrz(ix(1:n_points)); for j=1:n_points % 讓FLUENT標(biāo)出這些點 fprintf(fid1,'surface point-surf point_%d %e %e %e\n',L,x_maxchange(j),y_maxchange(j),z_maxchange(j)); L=L+1; end end arrp_old=arrp; end fclose(fid1); % 關(guān)閉命令文件

運行這個程序后將生成一個FLUENT命令文件j2.jou,在FLUENT中執(zhí)行它,便將最后三次迭代中壓力波動最劇烈的一些點標(biāo)記了出來(每次迭代標(biāo)記15個,共45個點)。

在FLUENT菜單欄選擇[Display]->[Mesh…],將這些點顯示出來,可以發(fā)現(xiàn),壓力波動最劇烈的點都位于出口截面上(圖5)。因此推測可能是出口邊界設(shè)置不當(dāng)導(dǎo)致不收斂。

FLUENT不收斂案例及解決方法：壓力出口導(dǎo)致的不收斂【轉(zhuǎn)發(fā)】fluent培訓(xùn)課程圖片5

圖5 壓力波動最劇烈的點

嘗試對出口的邊界條件進(jìn)行修改,發(fā)現(xiàn)當(dāng)使用“無反射”選項(Non-Reflecting Boundary,圖6)的時候,不收斂的問題就得以解決(圖7)。

FLUENT不收斂案例及解決方法：壓力出口導(dǎo)致的不收斂【轉(zhuǎn)發(fā)】fluent培訓(xùn)課程圖片6

圖6 “無反射”選項。

由于參考壓力設(shè)為101325Pa,所以表壓(Gauge Pressure)是0。

FLUENT不收斂案例及解決方法：壓力出口導(dǎo)致的不收斂【轉(zhuǎn)發(fā)】fluent仿真分析圖片7

圖7 非定常計算的殘差曲線。時間步長設(shè)為3×10-6秒。出口邊界使用“無反射”選項。圖中顯示了約25個時間步的殘差曲線?？梢钥闯?在每一個時間步內(nèi),殘差都能降低到默認(rèn)的收斂標(biāo)準(zhǔn)以下。

究其原因,在這個算例中出口邊界已經(jīng)達(dá)到了壅塞狀態(tài),這可以從馬赫數(shù)云圖上看出來(圖8)。從圖中可以看出,出口附近有馬赫數(shù)超過1的局部超音速區(qū)域。馬赫波在出口邊界的反射導(dǎo)致了出口截面的參數(shù)不斷振蕩,不能收斂。這種反射是邊界條件的數(shù)學(xué)處理造成的——因為我們強(qiáng)制地讓出口截面的壓力等于指定的數(shù)值,而這是不符合物理事實的。“無反射”的具體處理方法涉及特征線理論,這里不予以敘述,有興趣的讀者可以看計算流體力學(xué)原理方面的資料(例如[1];以及FLUENT的User's Guide中的“General Non-Reflecting Boundary Conditions”這一節(jié))。

FLUENT不收斂案例及解決方法：壓力出口導(dǎo)致的不收斂【轉(zhuǎn)發(fā)】fluent仿真分析圖片8