CFX处理弹丸在超音速下的烧蚀问题
2018/12/16 22:49:32
弹丸在高速飞行的过程中，由于稠密的大气，弹丸的表面将经受剧烈的气动加热，当弹丸的温度达到材料熔点时，弹丸局部就会开始熔化，熔化的部分将会立即被除去，这就是所谓的气动烧蚀现象....为了简化计算，本次算例采用的金属为铝，计算时间为1.6秒。
问题描述：
弹丸在高速飞行的过程中，由于稠密的大气，弹丸的表面将经受剧烈的气动加热，当弹丸的温度达到材料熔点时，弹丸局部就会开始熔化，熔化的部分将会立即被除去，这就是所谓的气动烧蚀现象。实际的烧蚀不单单是高温熔化的过程，它还包括金属的氧化反应。金属氧化物的熔点一般为金属纯净物熔点的一半低一点。一般弹丸采用高熔点的钨作为材料。为了简化计算，本次算例采用的金属为铝，计算时间为1.6秒。
涉及到的问题：
超音速问题，超音速的求解与一般流动的求解方法略有不同，超音速问题尤其对求解时间步与网格质量的要求相当高，非常容易出现求解发散的问题。其二就是动网格问题，要使用动网格来模拟材料的烧蚀。其三是网格的重构，当壁面发生较大变形时，网格质量不能维持计算，需要进行网格重构。CFX中并没有自带的网格重构功能，必须依靠第三方软件完成。
这里先开一个头，我会慢慢来完成这个帖子，这个帖子主要包括两部分。第一部分就是超音速问题的处理方法，第二部分是CFX调用user fortran实现动网格。
案例介绍
首先把模型发上来，计算域的大小和形状都是靠经验选取的。外流场的形状为半抛物线形，这样生成的网格方向与流场的流动方向能够高度保持一致，对于提高数值计算的精度有很大帮助。
对于边界条件的设置，CFX处理超音速问题与FLUENT略有不同，CFX进口可以选择速度入口、压力入口，可以定义流速、总压或者静压，总温或者静温，出口边界只需设置为超音速即可，没有多余的选项。
这里简单介绍一下超音速问题，其实超音速问题可以分为两种，风洞内的超音速问题与大气环境下的超音速问题。它们之间的区别在于，被测物体是否真实的发生了运动。风洞环境下的入口边界中我们一般会指定总压，总温，以及当地的音速，当地的音速只与当地的静温和材料有关，有兴趣的可以查看ansys帮助文件。而对于外流场问题，我们能给定的是静压，静温以及入口速度（已知的，与被测物体的相对速度）。静压与静温表示的是环境温度，外流场条件下会比较容易给定。风洞条件下的静温与静压是ma的函数，而ma又与温度等有关，所以风洞条件下会直接给定总压与总温。
介绍网格的生成方法
下面介绍网格的生成方法，首先要解除一个误区，网格的稀密程度只与计算结果的精度有关系，对于计算的收敛特性是没有绝对关系的。一个case的收敛好坏，与边界的设置有关，还有一个就是网格的质量。这里所说的质量，并不是简单的正交角和雅阁比矩阵的质量。它包括：边界层的布置是否合理、网格由疏到密过渡的控制以及网格的均匀程度。这些都会影响到计算的收敛特性。下面看网格的生成，要求要有边界层，网格有稀到密有很均匀的过渡。
对于超音速问题，边界层网格的划分对计算的影响非常大，所以这一部分的处理要花费一些精力。
二维网格准备
本次算例准备的是二维网格，我是简单进行了一次拉伸，生成两层节点。实际理想的二维计算应该是取计算域的一半，把对称轴取出来，然后旋转3度或者5度，这样的网格是与FLUENT中的旋转轴对称相一致的。但是我拉伸成平面的二维是有好处的，这是为后续的动网格做准备。如果选择轴对称的二维模型，在对称轴上，每个单元都是三角形的，而弹丸的其它部位（轴之外）则都是四边形的，这样USER CEL的循环位置既有四边形又有三角形，不容易能控制循环位置，而采取拉伸操作的话，弹丸的壁面全部为四边形。
CFXPRE中的设置，本次测试的工况为8MA，环境为海平面10m处的环境，环境压力101325pa，温度300K。气体选择理想气体，湍流模型选用SST模型（SST模型的兼容性比较好，壁面函数的选择与否也是自动的）。
进口边界条件：如下图所示
出口的设置直接设置为超音速，没有其它的附加选项。另外设置好对称面。
因为流体计算域与固体计算域都是在ICEM中一起划分的，她们的交接面是共节点的，一旦设置好两个计算域，交界面会自动生成。
对于共轭传热的流固交界面，不论网格是否共节点，网格的匹配类型都要使用GGI格式。首先建立的是一个稳态的，没有动网格的case，我们首先要把流场计算准确，否则后续的工作都是徒劳的。
边界条件设置好之后，尽管是稳态的计算，这里一定要给定一个流场初始化，分别指定速度、压力、温度。（对于收敛有极大的帮助）
初始化完毕之后，要进行时间步的控制。虽然是稳态的算法，但是CFX使用的是伪瞬态算法，这里也存在时间步的概念。与瞬态不同的是，对于共轭传热的问题，稳态算法可以对流体域和固体域指定不同的时间步长。一般的共轭传热问题，固体域的时间步至少是流体域时间步的100倍以上，甚至更大。这里我们只需要把流体域时间步长的选项设置为‘保守的’。或者把时间步系数设置为0.1或更小。
我们要反复的观察残差的收敛曲线，如果残差的收敛曲线呈现上扬或者持续不下降的趋势，说明时间步还是过大了，这样反复的调整。
看一个比较好的残差收敛曲线
计算开始一般都是波动的，但是一段时间之后，残差曲线能够稳定的下降。这就说明我们的网格、边界条件、时间步设置的比较合理。这样的计算结果才有可信度。
看一下没有考虑烧蚀的计算结果，
压力场：
温度场：
温度场近景：显示一下温度范围
可以看到，8Ma条件下，弹丸的头部最高温度可以达到3300以上。
流场的计算调试好之后，下面就要考虑烧蚀问题了，下面我会简单介绍一下烧蚀的基本原理，以及实现烧蚀在CFX中的遇到的难点。
首先介绍弹丸烧蚀的原理，在流固共轭传热的交界面上，热流量会从流体计算域流向固体域，这些热量会产生两个效果：1，以热传导的方式（固体域不存在对流）传入固体域的内部，导致内部温度升高。2，产生气动烧蚀，气动烧蚀也是一种吸收能量的方式。
在这里有几个物理概念：q，热流量单位W/m^2；F,熔化热单位KJ/kg；K,热导率单位W/m×K。还有一个表征弹丸表面切向与来流方向的夹角，beta。
假设Ds为弹丸在水平方向上烧蚀掉的厚度，在deta t时间内，如果控制体是无限薄的，则固体域内热量的传导方向为弹丸表面的法向（无切向分量），这样只须在这个方向上列热量守恒方程即可。
这个方程中，Q是流体域提供的热流量，在CFX中对应的变量为 Heat Flux；q（in）表示流入固体域的热流量，因为固体域不涉及对流项，所以q（in）可以表示为温度沿弹丸法向（指向固体域）的梯度与热传导率的乘积（温度的梯度和单元法向量是CEL功能实现不了的，这里就要用到USER FORTRAN）；另外，Ds为水平的后退距离，为了保证网格的质量，我们要求弹丸沿表面法向的后退速度，这个很简单，沿法向的后退速度为 Ds/deta t×sin（beta）。
最终形式为:
每个单元节点的位移（displacement）需要我们知道热流量Q，对应固体域在此节点处的温度梯度，以及单元的法向量。
难点来了
热量的守恒方程包含了两个计算域：流体域与固体域。我们知道，USER CEL是基于Nloc的循环，Nloc是不能同时在两个计算域上循环的。
例如，如果要求流体域交界面上的节点位移，不单要有此节点的热流量，还要知道对应此节点的固体域的温度梯度。这个就是利用CFX计算这类问题的难点。
我介绍一下我的求解思路：（在用到这些变量之前，把他们写为关于节点坐标的函数，调用的时候只要有对应节点的坐标就可以访问两个计算域的量）
首先利用junction box routine在MMS中创建几个数据单元（makdat命令），单元的维度为交界面上节点的个数。还有几个单元数据，用来存储稍后得到的多项式的系数。调用的位置可以指定为user input
然后继续创建一个junction box routine，要分别访问几个变量（getvar命令），包括：流体域内的热流量、固体域内的两个温度梯度、还有一个就是对应变量的节点坐标。
问题描述：
弹丸在高速飞行的过程中，由于稠密的大气，弹丸的表面将经受剧烈的气动加热，当弹丸的温度达到材料熔点时，弹丸局部就会开始熔化，熔化的部分将会立即被除去，这就是所谓的气动烧蚀现象。实际的烧蚀不单单是高温熔化的过程，它还包括金属的氧化反应。金属氧化物的熔点一般为金属纯净物熔点的一半低一点。一般弹丸采用高熔点的钨作为材料。为了简化计算，本次算例采用的金属为铝，计算时间为1.6秒。
涉及到的问题：
超音速问题，超音速的求解与一般流动的求解方法略有不同，超音速问题尤其对求解时间步与网格质量的要求相当高，非常容易出现求解发散的问题。其二就是动网格问题，要使用动网格来模拟材料的烧蚀。其三是网格的重构，当壁面发生较大变形时，网格质量不能维持计算，需要进行网格重构。CFX中并没有自带的网格重构功能，必须依靠第三方软件完成。
这里先开一个头，我会慢慢来完成这个帖子，这个帖子主要包括两部分。第一部分就是超音速问题的处理方法，第二部分是CFX调用user fortran实现动网格。
案例介绍
首先把模型发上来，计算域的大小和形状都是靠经验选取的。外流场的形状为半抛物线形，这样生成的网格方向与流场的流动方向能够高度保持一致，对于提高数值计算的精度有很大帮助。
对于边界条件的设置，CFX处理超音速问题与FLUENT略有不同，CFX进口可以选择速度入口、压力入口，可以定义流速、总压或者静压，总温或者静温，出口边界只需设置为超音速即可，没有多余的选项。
这里简单介绍一下超音速问题，其实超音速问题可以分为两种，风洞内的超音速问题与大气环境下的超音速问题。它们之间的区别在于，被测物体是否真实的发生了运动。风洞环境下的入口边界中我们一般会指定总压，总温，以及当地的音速，当地的音速只与当地的静温和材料有关，有兴趣的可以查看ansys帮助文件。而对于外流场问题，我们能给定的是静压，静温以及入口速度（已知的，与被测物体的相对速度）。静压与静温表示的是环境温度，外流场条件下会比较容易给定。风洞条件下的静温与静压是ma的函数，而ma又与温度等有关，所以风洞条件下会直接给定总压与总温。
介绍网格的生成方法
下面介绍网格的生成方法，首先要解除一个误区，网格的稀密程度只与计算结果的精度有关系，对于计算的收敛特性是没有绝对关系的。一个case的收敛好坏，与边界的设置有关，还有一个就是网格的质量。这里所说的质量，并不是简单的正交角和雅阁比矩阵的质量。它包括：边界层的布置是否合理、网格由疏到密过渡的控制以及网格的均匀程度。这些都会影响到计算的收敛特性。下面看网格的生成，要求要有边界层，网格有稀到密有很均匀的过渡。
边界层的控制
网格由密到疏的均匀过渡。
对于超音速问题，边界层网格的划分对计算的影响非常大，所以这一部分的处理要花费一些精力。
二维网格准备
本次算例准备的是二维网格，我是简单进行了一次拉伸，生成两层节点。实际理想的二维计算应该是取计算域的一半，把对称轴取出来，然后旋转3度或者5度，这样的网格是与FLUENT中的旋转轴对称相一致的。但是我拉伸成平面的二维是有好处的，这是为后续的动网格做准备。如果选择轴对称的二维模型，在对称轴上，每个单元都是三角形的，而弹丸的其它部位（轴之外）则都是四边形的，这样USER CEL的循环位置既有四边形又有三角形，不容易能控制循环位置，而采取拉伸操作的话，弹丸的壁面全部为四边形。
CFXPRE中的设置，本次测试的工况为8MA，环境为海平面10m处的环境，环境压力101325pa，温度300K。气体选择理想气体，湍流模型选用SST模型（SST模型的兼容性比较好，壁面函数的选择与否也是自动的）。
进口边界条件：如下图所示
出口的设置直接设置为超音速，没有其它的附加选项。另外设置好对称面。
因为流体计算域与固体计算域都是在ICEM中一起划分的，她们的交接面是共节点的，一旦设置好两个计算域，交界面会自动生成。
对于共轭传热的流固交界面，不论网格是否共节点，网格的匹配类型都要使用GGI格式。首先建立的是一个稳态的，没有动网格的case，我们首先要把流场计算准确，否则后续的工作都是徒劳的。
边界条件设置好之后，尽管是稳态的计算，这里一定要给定一个流场初始化，分别指定速度、压力、温度。（对于收敛有极大的帮助）
初始化完毕之后，要进行时间步的控制。虽然是稳态的算法，但是CFX使用的是伪瞬态算法，这里也存在时间步的概念。与瞬态不同的是，对于共轭传热的问题，稳态算法可以对流体域和固体域指定不同的时间步长。一般的共轭传热问题，固体域的时间步至少是流体域时间步的100倍以上，甚至更大。这里我们只需要把流体域时间步长的选项设置为‘保守的’。或者把时间步系数设置为0.1或更小。
我们要反复的观察残差的收敛曲线，如果残差的收敛曲线呈现上扬或者持续不下降的趋势，说明时间步还是过大了，这样反复的调整。
看一个比较好的残差收敛曲线
计算开始一般都是波动的，但是一段时间之后，残差曲线能够稳定的下降。这就说明我们的网格、边界条件、时间步设置的比较合理。这样的计算结果才有可信度。
看一下没有考虑烧蚀的计算结果，
压力场：
温度场：
温度场近景：显示一下温度范围
可以看到，8Ma条件下，弹丸的头部最高温度可以达到3300以上。
流场的计算调试好之后，下面就要考虑烧蚀问题了，下面我会简单介绍一下烧蚀的基本原理，以及实现烧蚀在CFX中的遇到的难点。
首先介绍弹丸烧蚀的原理，在流固共轭传热的交界面上，热流量会从流体计算域流向固体域，这些热量会产生两个效果：1，以热传导的方式（固体域不存在对流）传入固体域的内部，导致内部温度升高。2，产生气动烧蚀，气动烧蚀也是一种吸收能量的方式。
在这里有几个物理概念：q，热流量单位W/m^2；F,熔化热单位KJ/kg；K,热导率单位W/m×K。还有一个表征弹丸表面切向与来流方向的夹角，beta。
假设Ds为弹丸在水平方向上烧蚀掉的厚度，在deta t时间内，如果控制体是无限薄的，则固体域内热量的传导方向为弹丸表面的法向（无切向分量），这样只须在这个方向上列热量守恒方程即可。
Q×deta t= q（in）×deta t+ Ds×sin（beta）×density×F
这个方程中，Q是流体域提供的热流量，在CFX中对应的变量为 Heat Flux；q（in）表示流入固体域的热流量，因为固体域不涉及对流项，所以q（in）可以表示为温度沿弹丸法向（指向固体域）的梯度与热传导率的乘积（温度的梯度和单元法向量是CEL功能实现不了的，这里就要用到USER FORTRAN）；另外，Ds为水平的后退距离，为了保证网格的质量，我们要求弹丸沿表面法向的后退速度，这个很简单，沿法向的后退速度为 Ds/deta t×sin（beta）。
最终形式为:
每个单元节点的位移（displacement）需要我们知道热流量Q，对应固体域在此节点处的温度梯度，以及单元的法向量。
难点来了
热量的守恒方程包含了两个计算域：流体域与固体域。我们知道，USER CEL是基于Nloc的循环，Nloc是不能同时在两个计算域上循环的。
例如，如果要求流体域交界面上的节点位移，不单要有此节点的热流量，还要知道对应此节点的固体域的温度梯度。这个就是利用CFX计算这类问题的难点。
我介绍一下我的求解思路：（在用到这些变量之前，把他们写为关于节点坐标的函数，调用的时候只要有对应节点的坐标就可以访问两个计算域的量）
首先利用junction box routine在MMS中创建几个数据单元（makdat命令），单元的维度为交界面上节点的个数。还有几个单元数据，用来存储稍后得到的多项式的系数。调用的位置可以指定为user input
然后继续创建一个junction box routine，要分别访问几个变量（getvar命令），包括：流体域内的热流量、固体域内的两个温度梯度、还有一个就是对应变量的节点坐标。