一种基于ANSYS-APDL语言开发的复杂热环境下球形光学头罩瞬态热-结构耦合分析方法

摘要

本发明公开了一种基于ANSYS-APDL语言开发的复杂热环境下球形光学头罩瞬态热-结构耦合分析方法，包括以下步骤：一、基于APDL语言的任意球形光学头罩有限元模型的建立；二、基于APDL语言的球形光学头罩计算时间参数以及材料物理参数的确定；三、基于APDL语言的随时间变化的热流密度以及气动压力的换热以及力边界条件的加载；四、基于APDL语言的求解控制以及结果后处理。相比较直接采用用户界面进行分析的过程，本发明避免了同一类问题多次进行加载费事、费力、易错等缺点，而且通过宏的使用成功实现了头罩外表面热流密度和气动压力随时间变化情况下头罩气动热响应的动态计算。

说明书

技术领域

本发明属于气动光学效应研究领域，涉及一种基于ANSYS-APDL语言开发的复杂热环境下球形光学头罩瞬态热-结构耦合分析方法。

背景技术

气动热环境下球形光学头罩由于热传导和热应力等效应导致头罩材料的外形和光学性质发生变化，进而影响来自目标的红外辐射光线传输，使得高速飞行器成像探测系统对目标图像产生像模糊、抖动、偏移和能量衰减，最终将严重影响高速飞行器的探测性能。飞行器在大气层中高速飞行时其头罩与来流之间发生剧烈的相互作用，头罩温度升高，使光学成像探测系统处于复杂的气动热环境中，这种效应称为飞行器光学头罩气动热效应。球形头罩气动热效应中离不开热量的传递，因此在换热、力学边界条件已知的情况下，利用有限元分析软件ANSYS对头罩进行热-结构耦合计算，可得到气动热环境下头罩的温度场、热应力场、热应变场和热形变场等结果数据文件。后续我们可以根据已得到的结果数据利用热光效应和弹光效应原理，计算头罩在气动热环境下的折射率场，为气动热环境下高速飞行器外流场和头罩光传输仿真分析和头罩热辐射计算提供数据支持。

发明内容

本发明的目的是提供一种基于ANSYS-APDL语言开发的复杂热环境下球形光学头罩瞬态热-结构耦合分析方法。在本发明中，用ANSYS-APDL语言开发了一种外表面热流密度以及气动压力随时间变化的瞬态仿真计算法，通过开发的APDL命令流可以一次性计算得到球形光学头罩在气动热环境下的热响应结果(包括光学头罩在气动热环境下的温度场、热应力场、热应变场和热形变场)，避免了之前同一类问题多次使用用户界面进行加载费事、费力、易错等缺点，而且通过宏的使用成功实现了头罩外表面热流密度和气动压力随时间变化情况下头罩气动热响应的动态计算。

本发明的目的是通过以下技术方案实现的：

一种基于ANSYS-APDL语言开发的复杂热环境下球形光学头罩瞬态热-结构耦合分析方法，包括以下步骤：

一、基于APDL语言的任意球形光学头罩有限元模型的建立

在对复杂热环境下球形光学头罩进行瞬态热-结构耦合分析时，首先应依据分析的类型选定对应的有限元模型，在APDL中依靠命令ET进行有限元模型的定义。由于本发明要对复杂热环境下球形头罩进行热-结构分析，所以在所编写的APDL语言中所采用的类型为plane13与solid5。在选定有限元类型并进行定义后，进行仿真模型的建立过程。依据球形头罩的内外半径r’与r以及头罩顶端到底端的距离d三个几何参数，进行球形头罩计算模型的建立。在使用APDL语言建立部分圆环模型时，需要将头罩顶端到底端的距离d转化为部 分圆环所对的圆心角θ，利用APDL中命令中PCIRC命令生成对应角度的计算模型。在本发明中，所建立的部分圆环结构如图1所示，图示中的圆心角为83.475度。

接下来通过编写APDL命令对上述建立的部分圆环的四条边界线进行网格大小以及网格划分类型的控制，可以按照要求精度的高低自由限制划分网格的大小，网格划分类型采用二维平铺式。然后采用APDL中的AMESH命令划分网格，在本发明中采用上述方法划分网格后结果如图2所示。

在得到如图2所示的有限元模型后，我们通过APDL中VROTAT命令控制该模型绕坐标轴或是过两个特定关键点的直线旋转一定的度数，进而生成最终的有限元分析模型。具体旋转的度数应按照实际模型来确定，在本发明中图3-4中所示的球形头罩内外侧有限元模型是由图2所示的部分平面圆环旋转360度形成的。如果仅分析其中的一部分，可以选择旋转0°到360°的任意角度来建立模型。

二、基于APDL语言的球形光学头罩计算时间参数以及材料物理参数的确定

在使用APDL进行计算时，需要对计算时间总长、计算时间步长、估计熔化时间、飞行时间间隔以及子载荷的步数进行合理的设定。由于在后续研究中要利用本发明所得的数据结果结合热光效应与弹光效应计算气动热环境下高速飞行器光学头罩的折射率，所以在对材料的物理特性输入时，需要考虑材料的非线性特性。

在本发明中，使用APDL载入的头罩光学晶体材料的物理特性参 数包括：该材料的弹光系数、导热系数、弹性模量、热膨胀系数、屈服强度、熔点、切变弹性模量、o光折射率、e光折射率、o光热光系数、e光热光系数、比热容、密度和泊松比。在球形光学头罩材料确定的情况下，以上数据均可查阅相关资料得到。在得到以上所提的材料物理参数后，使用APDL中的MP，TB，TBTEMP，TBDATA等命令对相应的物理量进行定义与输入。

三、基于APDL语言的随时间变化的热流密度以及气动压力的换热以及力边界条件的加载

在对球形光学头罩进行热-结构耦合计算时需要输入换热与力边界条件。本发明在利用ANSYS对球形光学头罩进行热分析时，所加载的换热边界条件包括：头罩外侧气流温度、头罩内侧气流温度、头罩内侧对流换热系数、头罩外表面气动热流密度分布数据文件(flux_out.dat)、头罩初始温度和头罩热应力计算参考温度。其中头罩内侧温度分布假设为均匀分布，取为300K，头罩内侧对流换热系数可依据头罩内外温度差以及头罩材料查阅资料得到，本发明中球形光学头罩内侧对流换热系数为10。

flux_out.dat文件的数据结构如图5所示：文件中第1行，第14行，…第1+13*(i-1)行分别表示第i个时间间隔表示的计算时间，文件第2-13行，第2+13(i-1)-2+13*i行分别表示第i个时间间隔内头罩外表面第1-12个区域内气动热流密度，因此本发明可以计算当头罩外表面热流密度随时间变化情况下头罩气动热响应。

对应的力学边界条件包括：头罩内侧均一分布气压以及头罩外表 面气动压力分布数据文件pres_out.dat，pres_out.dat文件的数据如附图6所示：文件中第1行，第14行，…第1+13*(i-1)行分别表示第i个时间间隔表示的计算时间，文件第2-13行，第2+13(i-1)-2+13*i行分别表示第i个时间间隔内头罩外表面第1-12个区域内气动压力，因此本发明可以计算当头罩外表面气动压力随时间变化情况下头罩气动热响应。

四、基于APDL语言的求解控制以及结果后处理

根据表1，这部分内容对应于表1中的运行2与运行3。本发明在进行求解计算时，充分利用了APDL语言的二次开发功能，采用宏进行多层循环的控制。在本发明中采用*if-then-else，*endif命令在每个载荷步进行之前比较计算所得的最大温度与材料熔点的关系，来判断此载荷步是否进行。采用*do，*enddo命令来控制循环的次数以及计算的时间总长。当计算达到时间总长时，循环结束，然后可以根据实际要求计算得到时间总长范围内不同时刻的热响应结果。

另外，在本发明中通过编写APDL命令控制ANSYS通用后处理器，经过ANSYS求解后，对球形头罩所有节点的温度、热应变、热应力、热形变数据利用*get，*do，*enddo，*cfopen，*cfclos等命令按照要求的格式输出数据文件。在本发明中，采用以下APDL语言获得了86406个节点沿X，Y，Z轴方向的变形数据，其他结果数据均可采用类似方法进行输出。

*do，j，1，npoints

*get，ux，node，j，u，x

*get，uy，node，j，u，y

*get，uz，node，j，u，z

*vwrite，j，ux，uy，uz

(f15.1，f25.10，f25.10，f25.10)

*enddo

另外为了数据更加直观，在仿真过程中可以根据需求输出温度、不同轴方向的热应变、热应力、热形变等值线分布图，在本发明中给出了计算时刻为1s时的沿X轴方向的热形变与沿Y轴方向的热应力等值分布图结果如图8-9所示。

在本发明中，采用了ANSYS批处理(Batch)方式，首先通过使用APDL命令流，可灵活建立任意尺寸的球形头罩，且头罩顶端到底端距离也可依据用户要求做出改动。接下来在本发明中依靠APDL语言对所建立的计算模型进行网格划分，并将划分网格后节点的坐标数据写入文件。然后对前边已经建立的有限元模型进行换热与力边界条件载荷的施加，并进行求解。在换热边界条件的施加方式上，本发明可以实现随时间动态变化的热流的施加，进而可以得到总时间范围内任意时刻球形光学头罩的温度场、热应力场、热应变场和热形变场。最后通过APDL控制ANSYS后处理器按要求输出头罩温度、热形变、热应变、热应力数据以及对应的等值线分布图。其次，可以根据实际情况对所编写的APDL命令进行简单的修改，使其能拥有更更加广泛的应用范围，不仅仅局限于球形光学头罩的分析。

本发明通过APDL命令可以实现任意尺寸球形光学头罩模型的 建立，通过控制头罩底端到顶端的距离与部分圆环所对的圆心角的转化关系，达到了灵活建立模型的目的。另外相比较直接采用用户界面进行分析的过程，本发明避免了同一类问题多次进行加载费事、费力、易错等缺点，而且通过宏的使用成功实现了头罩外表面热流密度和气动压力随时间变化情况下头罩气动热响应的动态计算。

附图说明

图1为利用APDL语言建立的对应角度的部分圆环结构图；

图2为部分圆环结构二维网格划分结果示意图；

图3为由APDL从面旋转到体建立的球形头罩内侧有限元模型；

图4为由APDL从面旋转到体建立的球形头罩外侧有限元模型；

图5为球形头罩外表面气动热流密度数据存放格式示意图；

图6为球形头罩外表面气动压力数据存放格式示意图；

图7为球形头罩所划分12个区域的示意图；

图8为计算时刻为1s时头罩沿X轴方向的热形变等值分布图；

图9为计算时刻为1s时头罩沿Y轴方向的热应力等值分布图；

图10为采用APDL实现球形光学头罩瞬态热-结构耦合分析的整体流程图。

具体实施方式

下面结合附图对本发明的技术方案作进一步的说明，但并不局限于此，凡是对本发明技术方案进行修改或者等同替换，而不脱离本发明技术方案的精神和范围，均应涵盖在本发明的保护范围中。

本发明提供了一种基于ANSYS-APDL语言开发的复杂热环境下 球形光学头罩瞬态热-结构耦合分析方法，如图10及表1所示，具体内容介绍如下：

一、基于APDL语言的任意球形光学头罩有限元模型的建立，包括有限元类型的选择、头罩尺寸在APDL中的的正确转化、网格大小的控制以及旋转关键点的选择。

本发明在建立球形头罩模型时采用由面旋转到体的过程。如图1所示，在确定的直角坐标系下，通过控制所建立的部分圆环的圆心角θ来控制所建立头罩顶端到底端的距离d，二者满足的数学关系如下：

$\theta =arccos\left(\frac{r-d}{d}\right);$

式中：r——圆环的外径。

这样依据实际仿真的头罩的尺寸以及底端到顶端的距离结合上边数学关系，即可建立正确的计算模型。

在确定好圆环尺寸后，对圆环进行网格划分，因为采用的有限元类型为plane 13，所以本发明在平面所建网格为四边形网格，划分完网格后的示意图如图2所示。由于旋转出现新的面，需要提前对其单元类型进行定义，定义其有限元类型为solid 5，在设定好网格划分的个数后，对前边具有特定圆心角的部分圆环沿着确定的两个关键点进行旋转，得到的旋转后的有限元模型如图3-4所示。依据本发明中的APDL语言，可以实现将球形光学头罩进行合理分区分块。

二、基于APDL语言的球形光学头罩计算时间参数以及材料物理参数的确定

在使用APDL进行计算时，需要对计算时间总长、计算时间步 长、估计熔化时间、飞行时间间隔以及子载荷的步数进行合理的设定。由于在后续研究中要利用本发明所得的数据结果结合热光效应与弹光效应计算气动热环境下高速飞行器光学头罩的折射率，所以在对材料的物理特性输入时，需要考虑材料的非线性特性。在本发明中，使用APDL载入的头罩光学晶体材料的物理特性参数包括：该材料的弹光系数、导热系数、弹性模量、热膨胀系数、屈服强度、熔点、切变弹性模量、o光折射率、e光折射率、o光热光系数、e光热光系数、比热容、密度和泊松比。

三、基于APDL语言的随时间变化的热流密度以及气动压力的换热以及力边界条件的加载

本发明在利用ANSYS对球形光学头罩进行热分析时，所加载的换热边界条件包括：头罩外侧气流温度、头罩内侧气流温度、头罩内侧对流换热系数、头罩外表面气动热流密度分布数据文件(flux_out.dat)、头罩初始温度和头罩热应力计算参考温度。

flux_out.dat文件的数据结构如图5所示：文件中第1行，第14行，…第1+13*(i-1)行分别表示第i个时间间隔表示的计算时间，文件第2-13行，第2+13(i-1)-2+13*i行分别表示第i个时间间隔内头罩外表面第1-12个区域内气动热流密度，因此本发明可以计算当头罩外表面热流密度随时间变化情况下头罩气动热响应。

对应的力学边界条件包括：头罩内侧均一分布气压以及头罩外表面气动压力分布数据文件pres_out.dat，pres_out.dat文件的数据如附图6所示：文件中第1行，第14行，…第1+13*(i-1)行分别表示第i 个时间间隔表示的计算时间，文件第2-13行，第2+13(i-1)-2+13*i行分别表示第i个时间间隔内头罩外表面第1-12个区域内气动压力，因此本发明可以计算当头罩外表面气动压力随时间变化情况下头罩气动热响应。

球形头罩12个分区的结构如图7所示。

四、基于APDL语言的求解控制以及结果后处理

本发明在进行求解计算时，充分利用了APDL语言的二次开发功能，采用宏进行多层循环的控制，在每个载荷步进行之前比较计算所得的最大温度与材料熔点的关系，来判断此载荷步是否进行。当计算达到时间总长时，循环结束，然后可以根据实际要求得到计算时间总长范围内不同时刻的热响应结果。

在本发明中通过编写APDL命令控制ANSYS通用后处理器，经过ANSYS求解后，对球形头罩所有节点的温度、热应变、热应力、热形变坐标数据按照APDL控制的格式输出数据文件。为了数据更加直观，在仿真过程中可以根据需求输出温度、不同轴方向的热应变、热应力、热形变等值线分布图，在本发明中给出了计算时刻为1_s时的沿X轴方向的热形变与沿Y轴方向的热应力等值分布图结果如图8-9所示。

表1 APDL运行控制