【毕设奉献】STAR-CCM+ 外流场计算通用教程
0.前言
本教程主要是面向非专业用户想要了解 全机空气动力学普遍特征（尤其是超音速，捕捉激波间断的需要，这几乎必须要体积网格，高阶面元法只能部分模拟且观察不了流场）所写，面向CFD零基础但对飞机最基本概念有体会（如KSP飞机）玩家。
对于想要了解真正的真实气动的玩家，这显然远优于KSP原版或者FAR的空气动力学的。
本文基于我的本科团队毕设副组长经验所写，大后掠二段三角翼超音速飞机。
需要准备的原材料：
Simcenter STAR-CCM+已安装软件一份，模型一个，不小于8GB内存和至少CPU还算可以的电脑一个，不吃显卡，笔记本也可以。
我使用的是Star-CCM+ Single Precision 2022.1 (17.02.007)
STAR-CCM+主要的优点在于：支持自动网格可以将全机网格划分的手工成本几乎减到没有，以及有自带的专家初始化和收敛加速器来加速收敛和防止发散。建模请使用CATIA V5，体积只有Solidworks的六分之一，但是性能（尤其是曲面性能）远远高于Solidworks。
为了方便理解，这里就不用飞行力学学术上的坐标系了(直观看略阴间，至少美式机体系是这样，y轴向下，升力是负的)，我这里的坐标系是：X从机头指向机尾，Z指向上，Y按照右手系确定，这样有个优点，就是迎角可以直接与无侧滑的流向矢量对应，升力就是正值。
转换矩阵：
D=[cos(a)*cos(b),sin(b),sin(a)*cos(b)];//同时作为流向向量
C=[-cos(a)*sin(b),cos(b),-sin(a)*sin(b)];
L=[-sin(a),0,cos(a)]
1.基础设定
打开STAR-CCM+，建立文件，选择"在本地主机上并行"，核心数使用(N-2)。
注1：N是你的CPU核数，留出两个供系统和别的应用维持其运行就行了，不然太卡；STAR-CCM+会检测你的CPU的物理核心数目，如果超出你的CPU物理核心数则会解除CPU关联性并告诉你不要分配超过物理核心的数目）。
注2：在加载sim文件时，如果TEMP文件夹为假硬盘（内存盘放TEMP，或者还原卡系统）会告诉你java.io.FileNotFoundException；但是这个无须担心，不影响（但不可以将TEMP指向装满或者不存在的内存盘）。运行sim文件会产生几十MB的临时文件。
先别急着导入模型：提前进行物理连续体的设置，对于多个模型可以通用，节省设置步骤数。
连续体-新建-物理连续体，出现“物理1”，为灰色，重命名为air。
选择模型，三维，定常，气体，耦合流体，实际气体，平衡空气（耦合流体有利于超音速以及复杂流动计算），之后物理连续体将会被点亮。之后选择湍流，k-omega湍流（适用于全湍流计算，转捩模型不必设置，似乎无法直接使用）。
设置最小压力限制从1000Pa改为50Pa甚至10Pa（因为有高空计算的问题，如30km处静压就只有1200Pa，限制1000没得算了），当然这个高度应该还不算是滑移流。然后（这一步我当时没做，虽说影响不大。还好才1%阻力），STAR-CCM+的动力粘度默认是常值，需要在物理模型-材料-Air-属性里面找到动力粘度，手动调成sutherland。
此步骤完成时sim文件大小约为数百KB。保存一次，比较方便使用。
注1：尽量不要将有粘性的模型改成无粘（反过来可以）。我实测将k-omega湍流中的湍流项目逐个去掉然后在viscous regime改成无粘性，但是TurbModel依然被载入，疑似无影响但并不能确定对性能到底有没有隐性损耗。
注2：隐式非稳态大概是每5步出一个很高的残差尖峰（谷底的趋势与正常类似，而且力收敛按步数也没变快），对于未分离的外流场计算并无好处，而且很慢。定常就是OK的。
注3：会有自动选择的模型（湍流模型选择的时候，有全y+壁面处理，SST(Mentor) k-omega模型，且会载入ke和kw两个模型，属正常现象，搜索SST k-omega查原理便知，是基于壁面距离进行插值组合的一个拟合函数）。
2.导入模型
此时导入模型：3D-CAD模型，新建，进入3D-CAD界面后开始操作：注意，CAD转几何部件操作要一次性完成。实际上STAR-CCM+支持直接导入CATpart，也支持导入由CATIA等建模软件导出的igs文件。
导入igs文件则CATpart里面未合并的几何体就会变成每个几何体一个面（如果只有零件几何体那就是一个faces；但是似乎有比较大的精度损失，上面可能会出现凹坑什么的(0.1m/0.001m尺寸控制能明显看出，使用CATpart疑似可以减轻此问题。
而导入CATpart生成的曲面比较多（比如双圆弧机翼草图前缘圆角+多截面实体+实体后缘切角+镜像，就会有GSMLoft两块+surface一块+镜像一整块；镜像在ccm里面是单块曲面，但是被镜像的特征组不是，STAR-CCM+的原则是一个操作一块面，因此需要手动按需组合形面，并重命名（如fuselage，mainwing等等）。当然这个组合控制肯定效果比igs要好就是了，就是操作麻烦。）
导入模型后，需要生成几何零件，默认设置即可，实际上影响似乎不大。
此步骤完成时sim文件大小约为几MB。
3.设置计算域与流体域
这里选择长方体计算域，优点是设置方便。
几何-部件-新建形状部件-体，然后输入X、Y、Z坐标范围即可生成。
设置计算域需要大概估计一下计算域的参数，一般来说亚音速各个方向都是数倍，Z和Y设置为一样或者Z向尺寸设置为Y的60%即可。超音速可以在飞机前面缩短一些倍率，高超音速则可以进一步减低侧面的倍率（相对于飞机尺寸）。总之，超音速下保证最后出来的远场激波大概压在后侧面棱边线附近即可，特殊的超音速计算域处理（如锥形计算域等）先不说明。
设置完毕后，对计算域和上面的部件进行减运算操作，生成的部件重命名为fluid。之后选中减运算生成的部件，将部件分配给区域。选项为“为每个部件表面生成一个边界”，按需设置：
对于壁面，如果有分解的需要（如机身，主翼，鸭翼/尾翼等形面），按照导入模型处理时的结果设置多块壁面就行。对于计算域边界，重命名为FarField，设置为(六面)自由流，不用拆，也不用设置流体出入口。计算条件（压力、温度）按照所需高度依照ISA设置就可以。
迎角设置：单迎角的设置很简单，因为一般来说对非战斗机的运输类载具感兴趣的就是对称飞行（尤其是几乎为对称定直飞行的巡航）。
D=[cos(a),0,sin(a)]，L=[-sin(a),0,cos(a)]，完整版见上面的转换矩阵。
我直接使用Dev-C++用printf输出了，然后直接复制到对应的壁面设置的Java宏命令里面。
4.自动网格
重点。
需要按照以下方法设置，这个我摸索出来还算好：
基于切割体网格单元生成器（体积网格）和棱柱层网格生成器（边界层网格）的网格控制为：（括号内为相对基数的百分数值）
基础尺寸0.1m；对于物面，目标表面尺寸0.1m(基数)，最小尺寸0.002m(2)
曲率/圆点数为8（保证对机头机尾尖点小圆角处理区域的物面网格反映外形精度，限制网格数量不至于太大，也包括后续的前缘圆角的大量加密网格），
最大网格单元尺寸限制为10m(10000.0)（可以看清激波面）；对于Block.Block Surface（计算域边界），目标与最小尺寸均为10m(10000.0)，以防止计算域边界被细分浪费网格数；
体积增长率在切断尺寸0.1m下使用预设“中等”（小单元）和“慢”（大单元）[注：为了防止生成负体积网格]，后网格优化器两个选项（优化）均启用。
关于棱柱层的注意事项：
棱柱层的总厚度设置要小心，同样防止生成负体积网格。对于有粘计算要按照等比数列的思想求棱柱层的总厚度。当棱柱层数>5的时候近似公式可以为【1/(1-(1/k))-1】*100%，如对于延伸率2.0就是总厚度100%，尖前缘8层2.0是可以的。对于无粘计算的设置，采用2层棱柱层来过渡的时候统一采用棱柱层厚度75%（相对值），棱柱层延伸2，即可保证壁面附近体积网格精度足够。
目前尚无法克服的最大问题：
STAR-CCM+不支持各向异性加密，大曲率的小圆角是很容易出问题的，对于前缘圆角化的双圆弧机翼问题尤其大。因为我们希望的是前缘从机头看过去在YZ方向只有Z方向加密，但是实际上自动网格是按照各向同性生成的，因此Y方向也会加密，导致前缘网格数量多到离谱，尚无好的解决办法。所以给出的建议是，单纯为了探索超音速气动，机翼前缘可以不用圆角化，直接就是尖前缘，比较方便快速计算。
备注：当然机头建议用相切支持的封闭曲面+缝合曲面功能圆角化，因为机头如为尖头很容易在生成的网格中出现小范围的物面精度问题(头部前几个网格会变成一根很扭曲的尖劈)
5.求解器设置
在Part Meshes建立后，求解器即可开始设置。
给出一个比较方便的国际标准大气(ISA)计算器，在digitaldutch.com就有。当然你也可以用KSP-RSS的大气曲线或者按研究需要设置的条件（如KSP-RSS里面各个星球的大气参数）。
关于停止准则：
stop file我不知道咋用的，删了。这计算不慢的，对于未分离的常规较小攻角（个位数角度）巡航，可以设置一个200步避免超量跑计算。然后手工监控。跨音速设置个330步左右（Mach 0.91-1.09）就可以。其他设置，参见7.后处理栏目。
6.计算
点那个“人在跑”的图标就可以了。
计算开始的时候会先进行网格序列法初始化，该系列步骤计算很粗糙，而且不会被包括在残差图里。目的主要是最开始lv1就灌注了与边界条件相同的流体，lv2-9给出一个宽松的粗结果方便收敛。实际上这么设置过，物理连续体的初始条件就可以不用调，很松的。
之后进入正式计算，那就看收敛形状。一般曲线会有光滑段和粗糙段的区别，光滑段就是从最开始的残差1逐渐掉到最低的残差水平，因为未收敛，所以基本上就是向着收敛进行的过程。之后粗糙段基本上就是水平的每一步都有明显抖动的锯齿形曲线，图线趋势几乎无长线浮动，也就是收敛了。
注1：关于计算时操作
计算时大部分数据实际上并未冻结，都是可以改的，如修改条件会直接跳过网格序列法在已有基础上算导致残差回到1附近（因为边界附近不符）。如果发现计算条件设置错了，或者计算结束后想要开始下一组，务必记得清除所有计算结果（勾选域）。计算的时候你做的改动会在这一步完成的时候被输入进去。
注2：Mach 1.0附近收敛速率明显减缓。建议Mach 1±0.10设置330步起步，其他200步足够。
具体总结的收敛形状，我这里以我做的超声速巡航飞机为例：(Star-CCM+ Single Precision (17.02.007)，其余不保证准确性，有完整跨声速段可供参考)：
组1为Continuity，Energy，X-momentum(Xp)；组2为Y-momentum(Yp)，Z-momentum(Zp)；组3为Sdr(Specific dissipation rate)，Tke(Turbulent Kinetic Energy)
1.4-2.0：在单缩颈并收敛后明显张开并水平，步数约200，log RSD组1约在组2上方0.6
1.1：收敛形状类似1.4和2.0，但是凸包个数很多，收敛略困难，可能会组1类似正弦震荡而不完全收敛；log RSD组1和组2均偏大，且组1约在组2上方1.0。
1.05：残差极难收敛，可能需要400步以上才会收敛，但是会收敛。log RSD组1和组2均偏大，且组1约在组2上方1.5。力可以收敛，且力收敛比较快（超音速力收敛都比较快）
1.0：残差不收敛(log RSD组1可达0.0，组2也有-1.3左右)，但是力能收敛，大概残差稳定在上述水平时力即收敛，为每步浮动正负万分之一以内(<10N量级)。组1约在组2上方2.0。收敛约需270-330步。显然，因为超音速和亚音速气流几乎是本质区别，所以不收敛是正常的。
0.95：残差很难收敛，需要350步以上，但是也会收敛。log RSD组1和组2收敛后较为正常。力可以收敛，且力收敛比较慢（极高亚音速力收敛偏慢）。
0.9-0.7：残差较容易收敛，
0.7-0.3：高气压，缩颈后张开少，且Yp更接近组1紧挨下方，Zp在Yp下方离组1更远处。对于所有的速度范围，随着气压升高到接近1atm，最终组1和2几乎会变成一条线(>0.5atm时)。
所有速度组3一般在100步左右就要收敛，且明显收敛到组1和组2下方（组3也不属于耦合流体求解器）；超音速时收敛后极其稳定，亚音速&高气压下浮动大一点且无规律。
至于力的报告收敛曲线（见7.后处理），同样的收敛残差，超音速的力收敛的比亚音速更快，超音速在残差完全收敛前力就会完全收敛。
7.后处理
后处理大致可以分为数据采集和流场观测两类。
前者是基于报告的，STAR-CCM+里面报告直接使用力或者力矩报告就可以了。系数建议把ref.Area和DynPres(q)等扔进表格里面处理，就无需在生料中处理，以免弄错，也更直观。
新建报告，命名为Drag或者Lift，之后打开之后在里面的零部件里面选择你要研究的壁面部件（一般就选全部就可以，也就是整个飞机的所有形面）。然后对于力，需要设置方向（见前述“坐标系”段落；对于力矩，设置轴的指向和原点。报告可以用于生成监视器，而监视器之后就可以用于制作迭代过程中力随着迭代步数（在定常计算中）或物理时间（在隐式非稳态计算中）变化的状况，即我们关心的输出量的变化趋势（水平线就是收敛）。
后者是基于表面的，可以使用衍生部件，或者物面，进行二维面或曲面表面云图观测。
新建标量场景，然后准备建立衍生部件。建立衍生部件这里采用“等值面”然后选取Position-Laboratory-选择x或者y或者z（搜索pos就可以快速找到），之后建立的等值面会呈现深蓝色（因为默认是“新的表面显示器”）。之后找到标量场景自带的标量1，将衍生部件设置为标量1显示部件即可。
另一种是，直接将物面设置为标量1的显示部件，则会得到物面云图（压力等；只有无粘才能观测物面的速度，因为有粘计算中无滑移壁面处流体速度强制为0）。
8.结果判读
这里给出一些我大概的结果判读办法吧，学术上不好说，但从流体实际出发，是很好用的。
主要是流场观测的，用来理解飞机设计意义上的全机流场特性，主要在于实用而非理论价值。
(1)激波
基本上超音速流场里用恒定Y或者恒定Z的切片看到十分锐利(只有一层方格过渡的压力突增)的弓形或者直线边缘那肯定就是激波。包括在翼面与机身交汇的前缘处附近的机身压力云图上也可以看到一道>形状的锐利压力增加边线，这也是激波。包括Mach 0.9-0.95的时候你去看全机壁面压力云图看到一根明显由蓝绿变为黄橙的线那也是 上表面气流加速为超音速后通过激波在尾部回到亚音速产生的正激波线。
注1：超音速翼式(非乘波体)飞机的亚音速前缘(也就是Mach*sin(pi/2-Swept_Angle)<1，一般也就Mach 2.x及以下可以实现)
注2：学术上一般鉴定激波看马赫数突变，确实更准确，但是其实你这么看也是可以的，甚至更方便，因为空气动力学性能最终评判仅仅由全机表面（一般上和下）压力分布决定，这样对私下理解甚至更容易些，效应也更直接（当然压力分布确实会受到其他因素影响）。
(2)流动方向
这里其实不是很建议用矢量去做，因为矢量通常比较接近于x+方向(尤其流动未分离时)。所以其实可以用标量展示，不难理解，自己想象就看出来了。
去底下找到：工具-场函数，右键，新建场函数，将显示名称与实际函数名称都改为FlowAngleXZ或者FlowAngleYZ，之后输入以下公式即可：
纵向气流偏角Velocity angle XZ = atan2($$Velocity[2], $$Velocity[0])*57.295779513-【set_AoA】(角度/deg)(set_AoA是你在指定流向向量的时候的数值，这里需要手工输入)
已修正的横向气流内外流动角velocity angle YZ_Fix = atan2($$Velocity[1], $$Velocity[0])*57.295779513*$${Position}[1]/(abs($${Position}[1])+0.001) (角度/deg) (abs部分搓出符号函数，+0.001防止无穷大出现；57.3是弧度转角度)
注：学术上一般不用这个。但是这样表述对理解帮助是很大的，纵向气流偏角实际上可以让你很容易看出气流经过机翼后尾流的下洗程度，包括大后掠类三角翼构型前缘的内外流动效应、超声速双圆弧翼型的激波与膨胀波气流折转效应等，也能很容易分析出来。
(3)壁面Y+
检测网格生成效果的重要指标。关于Y+实际上是边界层无量纲化的一个距离参数，而u+是边界层无量纲化的速度参数。使用Wall Y+场函数可以查看，但壁面y+对于流场切片是无效的，只有壁面的表面直接在标量场景中设置为标量1的显示部件才可以显示。可以大概看看，估计数量级即可。
后记
终于能做全机全域气动了，我的确非常高兴。
以前确实基本上就是KSP里面玩一玩，飞行力学或者说启蒙感受可以，但是如果直接要于现实设计中使用，那还是肯定不行的。
至于FAR，我只能说，这空气动力根本就不真实，比如我捏了和CFD基本相同的上述飞机，但是在24km才能达到和CFD有粘17km相近的升阻比。波阻估计严重超标(CD0~205%，Mach 1.8)，明显亚音速性能严重偏高而超音速性能严重偏低。更何况KSP里面的飞机，基本上空天用途都是几乎全程在超音速飞行。
对于激波估测的不准确性以及积木式建造的不灵活性使得KSP在航空气动方面几乎无任何除了启蒙之外的意义，也包括FAR——可以不必认为FAR是真实气动了，性能不达标，这只是修正后的模块化气动，且FAR无翼型概念等。
另：
今天，特别写下此文，并祝贺本届本团队毕设圆满完成！
