CFX处理弹丸在超音速下的烧蚀问题

coolwolf0100

问题描述：
      弹丸在高速飞行的过程中，由于稠密的大气，弹丸的表面将经受剧烈的气动加热，当弹丸的温度达到材料熔点时，弹丸局部就会开始熔化，熔化的部分将会立即被除去，这就是所谓的气动烧蚀现象。实际的烧蚀不单单是高温熔化的过程，它还包括金属的氧化反应。金属氧化物的熔点一般为金属纯净物熔点的一半低一点。一般弹丸采用高熔点的钨作为材料。为了简化计算，本次算例采用的金属为铝，计算时间为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命令），包括：流体域内的热流量、固体域内的两个温度梯度、还有一个就是对应变量的节点坐标。

游客，如果您要查看本帖隐藏内容请回复

    如何将这些变量和坐标关联起来呢，这里要用到IMSL库函数中的RCURV命令，我们把对应节点的要求变量（热流量，温度梯度）写作关于节点坐标的函数，RCURV可以指定其函数形式为多项式的形式，函数的返回值为多项式系数的数组。（IMSL的安装和使用可以在网上随便搜，很多）


上面的test中用到了两个附加变量GX，GY，这个是温度的梯度，代码如下


把这些量准备好之后，我们就可以定义网格运动的字程序了。在定义Function时，我们只需要把坐标传递过去，就可以根据坐标求得对应节点上两个计算域内的量。

您如果有更好的方法，不妨在这里讨论；

把最后的两个文件传上来，




说明一下，以上上传的附件都是程序调试初期的文件，使用的是计算域的一半，而且计算并不是0时刻开始计算的，是使用的稳态结果最为初始值的。很多地方会带来数值误差。

最后上传一个初期制作的动画



本帖已全部结束

40534427

不错 能否讲的详细点

faee0

请问如何考虑材料的去除？动网格是包含材料的去除过程还是仅仅只涉及到弹丸的运动？
打算用fortran来改变计算域边界形状？那么这个形状如何确定？

coolwolf0100

40534427 发表于 2012-10-30 12:24
不错 能否讲的详细点

多谢版主关注，我会慢慢把这个帖子完善的。

coolwolf0100

faee0 发表于 2012-10-30 12:59
请问如何考虑材料的去除？动网格是包含材料的去除过程还是仅仅只涉及到弹丸的运动？
打算用fortran来改变计 ...

材料的去除就是通过动网格来实现的，计算过程中弹丸是静止的。涉及到的气动烧蚀理论我会在后续的补充中全部给出的

水若无痕

coolwolf0100 发表于 2012-10-30 18:59
材料的去除就是通过动网格来实现的，计算过程中弹丸是静止的。涉及到的气动烧蚀理论我会在后续的补充中全 ...

你应该是把材料去除的经验公式或者与温度变量的关系定义在了网格运动里面吧？

我见过轴承磨损也是这样做的。

coolwolf0100

水若无痕 发表于 2012-10-31 20:53
你应该是把材料去除的经验公式或者与温度变量的关系定义在了网格运动里面吧？

我见过轴承磨损也是这样做 ...

对的，材料的去除会由热流量与温度梯度确定。

frank0079

不错，已经很久没有这样的帖子了，期待更新！

tianming1025

讲的非常详细了，很好！
Ma=4我计算过一个标模，1600万的四面体网格，cfx结果还可以！
Ma=8  cfx的计算结果是否可信，需要实验数据支撑

40534427

仔细看了楼主的 *.F 文件，由于没有注释，有些地方很难懂，如果能够加上注释就完美了。此外对于“
如何将这些变量和坐标关联起来呢，这里要用到IMSL库函数中的RCURV命令，我们把对应节点的要求变量（热流量，温度梯度）写作关于节点坐标的函数，RCURV可以指定其函数形式为多项式的形式，函数的返回值为多项式系数的数组。（IMSL的安装和使用可以在网上随便搜，很多）”
我在考虑是否有更简单的方法实现呢？

coolwolf0100

40534427 发表于 2012-11-3 10:58
仔细看了楼主的 *.F 文件，由于没有注释，有些地方很难懂，如果能够加上注释就完美了。此外对于“
如何将这 ...

我使用这个方法的原因是因为：user fortran对于节点循环的位置不是我们可以控制的，打个比方，在流体域的交界面上，nloc=5时好比是第五个节点，要用到这个节点在固体域一侧的温度梯度，但是在固体域内，nloc=5就不是对应的这个节点了。这种一一对应的关系我们是不能得到的，所以我就把变量拟合了一下，这样的话，只要有坐标就可以得到另一侧计算域对应同样节点上的变量的。

40534427

可能Fluent udf 在处理此类问题上更方便点，希望有会员能够效仿你的方法做出来，看来楼主在CFX 二次开发应用很有造诣，不知道作出这样一个案例用了多长时间。

coolwolf0100

40534427 发表于 2012-11-3 18:56
可能Fluent udf 在处理此类问题上更方便点，希望有会员能够效仿你的方法做出来，看来楼主在CFX 二次开发应 ...

实不相瞒，我做了大概一个月了。每天都是早8点到晚10点！

40534427

楼主精神可嘉！值得学习啊。另外不知道楼主用的是否 Visual 2008+interFortran 11组合，IMSL 7 我也安装了，但是你的子程序一直没有调试成功过，其他的JB我都可以调试成功的，请楼主分享一下如何正确安装编程环境，这个还是有难度的。

coolwolf0100

我用的vs2010，fortran版本是Intel Parallel Studio XE 2011。另外，编译imsl的dll我是用的论坛上‘抛弃cfx5mkext'那个帖子上的方法。你需要把imsl的lib文件以及include路径添加到fortran的路径里面。您如果不清楚的话加我qq吧：547132660

greencker

看看 学学 顶顶

wjc0309

入门者，感觉做的很好！

fanhongshuai1

真是很好啊，谢谢楼主

LOGO213

楼主做的东西非常有借鉴意义，多谢了！

return776

好久没见这么好的帖子了

有几个想问的

1.收敛后，子弹壁面的Y+有多少
2.初始化时的初场 是from入口还是自己设定参数 还是用了udf什么的？
3.fluent算超音速时，经常遇到湍流受限问题，CFX会遇到吗？


还有想问的 先写这三个吧  希望楼主有空能回复下，，