一种柔性浮空器流热结构耦合计算方法

本发明涉及浮空器设计，具体为一种柔性浮空器流热结构耦合计算方法。

背景技术：

1、软式平流层飞行器具有体积重量巨大，质量、惯量变化显著，热力、动力耦合严重等特点，其内部气体的热力学状态对飞行状态和运行安全具有重要影响。飞行器运行过程中涉及太阳辐射、地球反射、红外辐射等热辐射，外部大气风场的强迫对流换热，内部气体的自然对流换热，热力膨胀和压缩产生的致冷和致热，以及蒙皮复合材料的导热等综合物理场作用，其热特性的仿真和计算较为复杂。

2、软式平流层飞行器的结构一般为轻质柔性蒙皮材料，依靠内部超压形成一定的外形包络，在内部超压或外部压力变化时，飞行器结构容易发生非线性变形，对飞行器的安全运行带来一定影响。飞行器的正常飞行期间，内部主要受到浮升气体超压以及流动作用，其中超压主要来源于浮升气体的昼夜温度变化，流动则主要来自内部气体自然对流；外部则受到不同风速带来的气动压力作用。因此飞行器的形变是一种充气膜结构的流动-传热-膜结构流热固耦合变形。

3、平流层飞行器流-热-结构耦合特性对飞行器循环能源系统、浮升能力、飞行控制特性以及蒙皮材料都有不同程度的影响，开展平流层飞行器流-热-结构耦合特性计算方法研究，从而有效预测飞行器在工作过程中的状态，能够促进平流层飞行器设计和控制的发展，保证其稳定性和可控性，为今后更多的平流层飞行器技术发展奠定基础。研究这些问题不仅是掌握平流层飞行器技术的重要基础，而且对丰富和发展传热学的研究内涵具有重要意义。

4、目前的文献中一般只考虑平流层飞行器流-热-结构单向耦合，实际上，流-热-结构之间相互影响，单向耦合的计算方法由于未考虑高度变化引起的内外压差变化，因此误差较大，而传统的双向耦合的计算方法过于复杂，难以快速求解。

技术实现思路

1、本发明的目的在于提供一种柔性浮空器流热结构耦合计算方法，以解决上述背景技术中提出的目前柔性浮空器流热结构耦合计算过程中由于存在流-热-结构之间相互影响，而导致单向耦合的计算方法误差较大，双向耦合的计算方法又过于复杂，难以快速求解的问题。

2、为实现上述目的，本发明提供如下技术方案：一种柔性浮空器流热结构耦合计算方法，包括以下步骤：

3、步骤s1：建立平流层环境热模型；

4、来自太阳的直接辐射id可表示为：

5、id＝τatm·isun

6、其中，τatm为太阳直接辐射在地球大气环境中的透过率，isun表示地球大气外太阳直接辐射强度；

7、地球外太阳直接辐射强度isun为：

8、

9、式中，i0为太阳辐射常数，i0＝1367w/m2，eearth为地球偏心率，eearth＝0.016708，为地球真近点角，可表示为：

10、

11、其中日角θday可表示为：

12、θday＝2π·(n-n0)/365.2422

13、其中，n表示当前日期在一年中的天数，n0表示天数的闰年修正项；

14、

15、式中，year表示计算日期对应的年份；

16、太阳直接辐射在地球大气环境中的透过率τatm为：

17、

18、式中，ph和p0分别表示高度为h处的大气压强和海平面大气压强，clow和chigh分别表示低空校准因子和高空校准因子，λam表示太阳辐射穿过大气层时的空气质量比，空气质量比λam与飞艇高度角θele有关，则不同高度角θele对应的空气质量比如下式所示：

19、

20、在大地坐标系中，太阳光直接辐射的方向是太阳高度角θele和太阳方位角θazi的函数。

21、

22、其中太阳高度角和方位角可表示为：

23、θele＝arcsin(sin(θdec)·sin(φ)+cos(θdec)·cos(φ)·cos(θhour))；

24、

25、其中，θdec表示太阳赤纬角，φ表示当地纬度，θhour表示太阳时角；

26、θdec＝[0.3723+23.2567·sin(θday)+0.1149·sin(2θday)-0.1712·sin(3θday)-0.758·cos(θday)+0.3656·cos(2θday)+0.0201·cos(3θday)]·π/180

27、θhour＝15(time+et/60-12.0)·π/180

28、et＝0.0028-1.9857·sin(θday)+9.9059·sin(2θday)-7.0924·cos(θday)-0.6882·cos(2θday)

29、天空散射辐射强度可表示为：

30、

31、对于任意倾角为的平面，可近似认为其所受到的天空散射辐射强度为：

32、

33、此外，当太阳光经过大气时，还有一部分能量会被大气中的水蒸气、二氧化碳和臭氧等气体吸收，引起大气温度上升，进而产生天空红外辐射，一般情况下，假设大气为黑体，根据stefan-boltzman定律，天空红外辐射功率可表示为：

34、iir_sky＝σ·tsky

35、其中，tsky为天空等效温度，可表示为天空发射率εsky和大气温度tatm的函数：

36、

37、地面红外辐射主要是来自地面的长波辐射，辐射强度与地球表面发射率εg和地球表面温度tg有关，地面处的辐射强度为：

38、

39、其中，εg为地球表面的发射率，与地表属性有关，沙漠地区的发射率为0.85左右，地球表面平均发射率约为0.95，有积雪的情况下取0.98；te0为地球表面温度，地表辐射的方向垂直于地面，在向高空辐射的过程中同样存在损耗，在高度为h的空间，地表辐射强度可表示为

40、iir_gh＝τir_gh·iir_g

41、其中τir_gh表示高度h处地表辐射在大气中的透过率，其表达式为：

42、τir_gh＝1.716-0.5(exp(-0.65ph/p0)+exp(-0.95ph/p0))；

43、步骤s2：建立平流层飞艇热模型；

44、临近空间飞艇热环境可以分为外部热环境和内部热环境两部分，其中外部热环境由太阳直接辐射、大气散射、大气红外辐射、地球红外辐射以及气球与外部环境之间的对流换热组成，内部热环境则由内部气体自然对流和内部气体与蒙皮之间的红外辐射组成；

45、蒙皮吸收的太阳直接辐射、大气散射和地球反射能量可表示为：

46、qd＝α·id·sprojected·(1+τ/(1-r))

47、qs＝α·is·ssurf·(1+τ/(1-r))

48、qr＝α·ir·sprojected·(1+τ/(1-r))

49、来自大气红外辐射、地球红外辐射以及内部气球的红外辐射可表示为：

50、qir_atm＝αir·iir_atm·ssurf·(1+τir/(1-rir))

51、qir_ear＝αir·iir_ear·sprojected·(1+τir/(1-rir))

52、qir_he＝αir·iir_he·ssurf_he·(1/(1-rir))

53、其中，α、τ和r分别表示蒙皮的可见光平均吸收率、透过率和反射率，αir、τir和rir分别表示蒙皮的红外平均吸收率、透过率和反射率，ssurf和sprojected分别表示囊体的表面积和投影面积；

54、蒙皮与外界环境之间的对流换热可表示为：

55、qconv_ex＝hex·ssurf·(tatm-tf)

56、hex为外部对流换热系数，可表示为：

57、

58、其中，hex_force为外部强制对流换热系数，hex_natural为外部自然对流换热系数，外部强制对流换热系数hex_force可表示为：

59、

60、其中，reatm表示大气雷诺数，katmh表示高度为h处的大气热导率，l0表示飞艇特征长度；

61、

62、

63、

64、其中，d表示飞艇特征直径，υatm表示空气流速，ρatmh、μatmh和tatmh分别表示高度h处大气密度、大气粘性系数和大气温度；

65、外部自然对流换热系数hex_natural可表示为

66、

67、其中，nuatm_natural为外部大气努塞尔数，l0为飞艇特征长度；

68、

69、式中，gratmh和pratmh分别表示高度为h处的大气格拉晓夫数和普朗特数；

70、

71、pratmh＝0.804-3.25·10-4·tatmh；

72、步骤s3：结合飞艇外部风场的风速和风向，基于上述平流层环境热模型和平流层飞艇热模型，采用有限元方法求解飞艇内部氦气温度场和氦气平均温度；

73、步骤s4：建立飞艇外部流场模型，求解飞艇外部流场在飞艇表面的气动压力；基于步骤s3中所建立的飞艇仿真模型，采用fluent软件计算外部流场在飞艇表面的气动压力，计算域前端选择速度入口，后端选择压力出口，采用稳态压力基求解连续方程和rans方程，湍流模型选择realizable模型，压力和速度耦合采用simple算法，并采用二阶迎风格式进行方程离散，计算残差设置为10-4；

74、步骤s5：建立临近空间飞艇结构变形仿真模型；将步骤s3中建立的飞艇几何模型导入abaqus软件中，根据蒙皮材料实际力学性能设置蒙皮的弹性模量和泊松比，同时设置约束条件为：飞艇头部约束x、y、z方向的位移，飞艇尾部仅约束y和z方向的位移，x方向的位移不作约束；

75、步骤s6：按照如下步骤求解飞艇流-热-结构耦合效应下的结构变形：

76、1)基于步骤s3中计算得到的飞艇内部氦气平均温度t，根据理想气体状态方程计算飞艇内压：

77、

78、其中，mhe为飞艇内部氦气质量，r为氦气的摩尔气体常数，v为飞艇体积；

79、2)结合飞艇当前的飞行高度h，基于美国标准大气模型计算得到飞艇外

80、部的大气压力为：

81、

82、根据飞艇内部压力pin和外部大气压力patmh，得到内外压差δp1：

83、

84、3)基于将飞艇的膨胀/收缩过程视为恒温膨胀/收缩的假设，将飞艇内部压力pin、外界大气压力patmh、以及步骤s4中计算得到的飞艇表面气动压力分布pw导入abaqus软件并加载到飞艇蒙皮表面，利用abaqus软件计算飞艇结构变形以及形变后的体积v。

85、4)根据浮重平衡计算平衡后高度h，计算方法如下：

86、飞艇所受到的浮力可表示为b＝ρatm(h)gv，其中ρatm(h)为高度h处的大气密度，g为重力加速度，根据美国标准大气模型，大气温度tatm和气压patm可表示为

87、

88、

89、因此大气密度ρatm(h)可表示为：

90、

91、其中rm_atm为大气的摩尔气体常数；飞艇所受重力为g＝mg，其中m为飞艇的质量；

92、当浮重平衡时，有ρatm(h)gv＝mg，可得高度即密度关于高度的函数的反函数，其中

93、

94、联立上述公式，即可计算得到高度h；

95、5)根据步骤1)更新飞艇形变后的内压pin；

96、6)根据步骤2)更新大气压力patmh，并计算内外压差δp2＝|pin-patmh|；

97、7)假设气动压力分布pw不随飞艇结构变形而变化，设置收敛条件ε，一般可取ε＝5pa，在此前提下，进行如下判定：

98、若|δp1-δp2|≤ε，则认为飞艇形变达到稳定，对应的结构变形云图、蒙皮应力分布云图及变形后的体积v即为最终稳定状态；

99、若|δp1-δp2|＞ε，则飞艇形变未达到稳定状态，利用abaqus软件导出该非稳定状态的飞艇体积v，并重新设定内外压差为重复步骤3)～7)，直至|δp1-δp2|≤ε，此时认为飞艇形变达到稳定，对应的结构变形云图、蒙皮应力分布云图及变形后的体积v即为最终稳定状态。

100、作为一种优选的技术方案，在步骤s3中，采用有限元方法分析飞艇内部氦气温度场和氦气平均温度：在space claim软件建立飞艇几何模型，其中飞艇坐标系采用右手坐标系，飞艇轴向为x方向，高度方向为z方向，利用fluent meshing进行网格划分，网格划分时采用poly-hexcore体网格；基于步骤s1中建立的平流层环境热模型和步骤s2中建立的飞艇热模型，在fluent软件中编写飞艇热模型udf程序，利用udf程序准确加载步骤s2中建立的飞艇热环境模型以及当前飞艇高度处对应的风场风速，基于有限体积法求解质量、动量和能量控制方程。

101、作为一种优选的技术方案，有限体积法的流场控制方程如下：

102、流体域质量控制方程为：

103、

104、流体域动量控制方程为：

105、

106、流体域能量控制方程为：

107、

108、边界条件为：

109、质量和动量守恒方程：假设飞艇蒙皮不存在氦气泄漏，因此飞艇与外部环境之间不存在质量转移，另外，假设蒙皮的运动速度为零。

110、能量守恒方程：在飞艇与外界环境之间存在两种形式的热交换，即热对流和热辐射，因此蒙皮的能量方程可表示为

111、stf＝qd+qs+qr+qir_atm+qir_ear+qir_grid+qir_in+qconv_ex+qconv_in-qpower。

112、作为一种优选的技术方案，利用fluent软件及udf程序进行热分析后，得到氦气温度场，采用体积平均法获得氦气的平均温度，其中热分析时，部分参数选择情况如下：

113、(1)设置工作条件为操作压力5829pa，工作温度230k，假设飞艇内部氦气为不可压缩气体，选择基于压力法的求解器pressure-based，时间类型选择稳态，速度格式为绝对速度；

114、(2)打开能量方程；

115、(3)湍流模型采用标准k-ε模型；

116、(4)辐射模型选择s2s模型；

117、(5)边界条件选择：蒙皮和太阳能电池的壁面条件选择混合边界条件；

118、(6)求解方法选择：求解方法选择coupled方法，压力选择presto，动量方程、湍流动能方程、湍流耗散率方程和能量方程均选择二阶迎风格式；

119、(7)松弛因子选择：选择压力松弛因子0.5，动量松弛因子0.5，密度松弛因子1，体积力松弛因子1，湍流动能松弛因子0.75，湍流耗散率松弛因子0.75，湍流粘度松弛因子1，能量松弛因子0.75；

120、(8)选择标准初始化，参考温度220k。

121、与现有技术相比，本发明的有益效果是：

122、(1)飞艇持续受到外界热环境作用，其温度变化不大，近似为准稳态过程，因此可将飞艇的膨胀/收缩过程视为恒温膨胀/收缩。此外，飞艇外部风场环境主要影响飞艇的热特性。外部风场环境引起的飞艇表面气动压力对结构变形的影响较小，同时飞艇结构变形之后引起的气动压力分布的变化也较小，因此可以将气动压力与结构变形之间的双向耦合简化为单向耦合关系。基于上述假设，可以将平流层飞艇流-热-结构双向耦合简化为流-热-结构单向耦合，从而大大简化飞艇的流-热-结构耦合效应的计算过程，在保证一定计算精度的前提下提高计算效率。

123、(2)目前文献中对于飞艇流-热-结构耦合的计算方法，一般只计算到步骤s6(3)，实际情况是，该步骤中得到的结构变形以及形变后的体积v并非稳定后的状态，飞艇的结构变形与飞行高度强烈耦合，具体体现在：飞艇体积发生变化，体积变化导致飞行高度变化，由于不同高度下大气压力不同，因此飞艇的内外压差随高度发生变化，进而又导致飞艇的体积发生变化，浮重失衡导致体积和高度进一步变化。此过程为结构变形与飞行高度的双向耦合过程。

124、本发明将上述双向耦合过程解耦，通过步骤s6(4)～步骤s6(7)的迭代计算，可以获得较为准确的平流层飞艇流-热-结构耦合特性，同时大大提升了计算效率。