FLUENT不收斂案例及解決方法

化佳為貴

2018年5月31日 13:34

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

圖1　加力燃燒室

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

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

采用定常算法計算不收斂（圖2；這里我們使用FLUENT默認的收斂條件，即能量方程的殘差降低到1e-6以下，其余方程降低到1e-3以下）。考慮到可能是分離流誘發的非定常效應導致不收斂（見公眾號先前的文章“為何我這個流動總是算不收斂？我要砸電腦！”），我們嘗試使用非定常算法。但是不幸的是非定常算法仍然不能在每個時間步內收斂。非定常計算的典型殘差曲線如圖3所示。

圖2　定常計算的殘差曲線

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

為了弄清不收斂的原因，我們用MATLAB編寫兩個小程序。第一個程序用來產生一個命令文件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); % 將計算結果保存到文件
end
fclose(fid); % 關閉命令文件

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

圖4　計算結果文件

第二個程序分析這些計算結果，找出壓力波動最劇烈的點：

clear
n=30;
m=375045; % 網格數量
n_points=15;
fid1=fopen('j2.jou','wt'); % 打開命令文件
L=1;
for i=27:n
    filename=sprintf('d:\\case1-%d.ip',i);
    fid2=fopen(filename); % 打開計算結果文件
    % ASCII碼中，40代表左圓括號
    % FLUENT的計算結果文件中，每一組數據都是以左圓括號開頭的
    % 因此，可以通過查找左圓括號的方法找到每一組數據的起點
    % 需要了解更多關于FLUENT計算結果文件格式的內容，請
    % 閱讀User's Guide中“Format of the Interpolation File”這一節
    % until函數需要自己編寫，見注1
    until(fid2,40); 
    arrx=fread(fid2,[m,1],'double'); % 讀取x坐標
    until(fid2,40);
    arry=fread(fid2,[m,1],'double'); % 讀取y坐標
    until(fid2,40);
    arrz=fread(fid2,[m,1],'double'); % 讀取z坐標
    until(fid2,40);
    arrp=fread(fid2,[m,1],'double'); % 讀取壓力場
    fclose(fid2);    
    if i>27 % 只關注最后三次迭代
        % 比較相鄰兩次迭代的結果，找出壓力波動最大的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標出這些點
            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); % 關閉命令文件

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

在FLUENT菜單欄選擇[Display]->[Mesh…]，將這些點顯示出來，可以發現，壓力波動最劇烈的點都位于出口截面上（圖5）。因此推測可能是出口邊界設置不當導致不收斂。

圖5　壓力波動最劇烈的點

嘗試對出口的邊界條件進行修改，發現當使用“無反射”選項（Non-Reflecting Boundary，圖6）的時候，不收斂的問題就得以解決（圖7）。

圖6　“無反射”選項。由于參考壓力設為101325Pa，所以表壓（Gauge Pressure）是0。

圖7　非定常計算的殘差曲線。

時間步長設為3×10^-6秒。出口邊界使用“無反射”選項。圖中顯示了約25個時間步的殘差曲線。可以看出，在每一個時間步內，殘差都能降低到默認的收斂標準以下。

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