概述
我是CFD里可压流程序开发出身，一直研究超音速湍流、激波相关问题，博士阶段对不可压流算法接触较少。 工作之后才开始学习了解不可压流算法、程序，其中最知名的就是OpenFOAM。 OpenFOAM基于有限体积法，cell-centered格式。 与OpenFOAM相似，SU2也主要采用有限体积法，但是vertex-centered格式。 两者有限体积法的基本原理是一致的，但是具体实现上不同，也不好直接对比。
本文对比SU2中基于可压流的求解器与OpenFOAM中基于SIMPLE算法的simpleFoam求解器在马赫数$0.15$、雷诺数$6 \times 10^{6}$、攻角$10^{\circ}$的湍流NACA0012算例中的收敛性。 这个NACA0012算例是比较简单的算例，虽然高雷诺数湍流对边界层网格要求较高，但是我们此处不是为了算准翼型的$C_l, C_d$，而是为了检验非线性系统的收敛性。 攻角$10^{\circ}$会造成一定程度的流动分离，使得稳态收敛有一定程度的小困难。 但是这个程度的攻角一般认为收敛性问题不大。 但是simpleFoam使用的SIMPLE算法这种解耦法是否还能收敛？ 从本文的测试来看，simpleFoam这种不可压流解耦法求解器在NACA0012攻角$10^{\circ}$高雷诺数湍流工况下无法收敛，而SU2可压流耦合法可以顺利收敛。
算例描述
算例采用NASA Turbulence Modeling Resource网站上经典的NACA0012湍流测试算例的工况设定：马赫数$0.15$、雷诺数$6 \times 10^{6}$、攻角$10^{\circ}$。
OpenFOAM需要在初场、边界条件中指定有量纲的速度、压力。 当然不可压流下压力绝对值没意义，范围才有意义。 SU2可以使用有量纲的工况，也可以进行无量纲化。 虽然理论上用户可以自己对工况做无量纲化，但是暂时直接测试有量纲工况下的OpenFOAM和SU2。 使用如下脚本计算当前工况的各项有量纲参数，以便OpenFOAM中设定。
import numpy as np
T_infty = 300.0
GAMMA = 1.4
R = 287.015
Ref_Length = 1.0
Ma_infty = 0.15
Re = 6E6
AoA = 10.0
sound_speed = np.sqrt(GAMMA*R*T_infty)
U_infty = Ma_infty * sound_speed
Ux_infty = U_infty * np.cos(AoA/180.0*np.pi)
Uy_infty = U_infty * np.sin(AoA/180.0*np.pi)
nu = U_infty*Ref_Length/Re
print("Flow condition is as follows,")
print("Ma_infty = %.2f"%(Ma_infty))
print("Re_infty = %.1e"%(Re))
print("AoA = %.1f"%(AoA))
print("Ux_infty = %23.15e"%(Ux_infty))
print("Uy_infty = %23.15e"%(Uy_infty))
print("nu = %23.15e"%(nu))
工况参数如下
# python3 calcFlowCondition.py
Flow condition is as follows,
Ma_infty = 0.15
Re_infty = 6.0e+06
AoA = 10.0
Ux_infty = 5.128846016841165e+01
Uy_infty = 9.043539326682394e+00
nu = 8.679944556274538e-06
网格采用simpleFoam自带的tutorials/incompressible/simpleFoam/airFoil2D算例中的网格文件。
该网格文件并未针对高雷诺数湍流做边界层网格加密，如下图所示。
尽管当前网格无法实现$C_l,C_d$的精确计算，但是当前测试目标是收敛性。
为了使SU2使用上述OpenFOAM自带网格计算，使用jingchangshi/ConvertOpenFOAMToSU2脚本转换。
需要提及的是，OpenFOAM这个自带的网格中，翼型头尾两端分别位于x=-17.5472, x=17.5处，所以在网格格式转换前先平移缩放网格，使得翼型弦长为1。（一般弦长都会定为1，为了方便与其他人对比，故平移缩放为惯例的弦长。）
使用OpenFOAM自带的工具在airfoil2D算例根目录下执行如下命令，
transformPoints -translate "(17.5472 0 0)"
# 0.02853294985048734 = 1.0/(17.5+17.5472)
transformPoints -scale 0.02853294985048734
然后使用上述格式转换工具将OpenFOAM自带网格转换为SU2格式的网格文件。
值得一提的是，我尝试过OpenFOAM自带的网格格式转换工具plot3dToFoa和gmshToFoam，均存在各种问题，而且难以简便地解决，所以使用上述工具。
接下来就是按照OpenFOAM和SU2的流程运行该算例。 使用SA湍流模型。
结果对比
设定最大运行5000个非线性迭代步，收敛残差设定为$10^{-12}$。 SU2的密度残差达到$10^{-12}$后停止，simpleFoam运行5000个非线性迭代步后仍然无法达到指定收敛残差终止。 OpenFOAM和SU2的残差收敛历史如下
图中横轴是非线性迭代步数，左边纵轴是残差，取log10，右边纵轴是每个非线性迭代步时的$C_l$与终止计算时的$C_l$之间的差值，对差值的绝对值+$10^{-13}$取log10。
从图中可以看出，simpleFoam在600个非线性迭代步左右收敛滞止。
随后直至5000步，压力、速度和湍流量的残差稳定振荡。
同时$C_l$有规律振荡。
此处，我怀疑OpenFOAM在$C_l$的计算中引入了FP32浮点数，尽管全局计算是双精度的。
当然，这里也可能是$C_l$的有规律振荡是残差稳定振荡的表现。
我没有深入OpenFOAM这一块代码进行验证。
很多人、尤其是国内工业仿真工程师，将收敛滞止当作收敛处理。
从当前测试的NACA0012高雷诺数湍流算例可以看出，$C_l$仍然存在不超过$10^{-7}$的振荡，无法给出准确值。
而有量纲模式的SU2计算中能量、动量、湍流量的残差稳定收敛。 同时，$C_l$逐渐逼近最终值。
对比终止计算后的X方向速度云图如下
SU2的结果与OpenFOAM存在一定区别。
OpenFOAM的X方向速度在出口边界处存在振荡。
攻角$10^{\circ}$下，流动斜向上，但振荡区域存在于翼型斜向下区域。。
可能当前OpenFOAM自带算例设定的freestreamVelocity类型有问题。
另外，顺便对比了SU2下有量纲计算与无量纲计算的结果。 残差收敛历史对比如下
从图中可以看出，有量纲模式下，各变量残差幅值与变量本身幅值密切相关。 无量纲模式下的各项残差全面小于有量纲模式下的值。 无量纲模式下密度残差收敛至$10^{-12}$需2563步， 有量纲模式下需要4924步。 无量纲模式下2563步时，$C_l = 1.2349506076894$； 有量纲模式下2563步时，$C_l = 1.2349506164452$，与无量纲模式下2563步时结果相差$8.755800084969678 \times 10^{-9}$； 4924步时，$C_l = 1.234950591527$，无量纲模式下2563步时结果相差$-1.6162399862906796 \times 10^{-8}$。
无量纲模式下终止计算时的X方向速度云图与有量纲模式的结果肉眼观察无区别，这也是可以预见的。
源代码及结果文件
OpenFOAM和SU2的算例设定文件、网格文件及结果文件见 jingchangshi/CompareConvergenceOpenFOAMAndSU2。