今天分享一个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”这一节）。
图8 马赫数云图
如果计算域的出口没有局部超音速区域，就不必将出口边界设置为无反射的了。计算实践表明，这时不将出口设成无反射的也是可以收敛的。
来源：流体那些事儿（ID：happyfluid）转载分享
微信
QQ
微博