平流层浮空器具有广泛的应用前景，可以用于对地观测、战区侦查以及远程通信。日本ISAS组织在2002年放飞了总体积达到600 000 m³的高空气球^[1]。2005年，美国洛克希德·马丁公司设计了长152.4 m、直径45.7 m、总体积105 000 m³的高空飞艇^[2]，并且于2011年7月进行了缩比验证艇的首次试飞，缩比验证艇体积为14 150 m^3[3]。韩国设计了长50 m、最大直径为12.5 m、总体积为409 m³、载重量为100 kg的飞艇^[4]。
浮空器高空飞行时，囊体内部氦气温度不断变化，白天太阳辐射囊体导致其内部温度增高，夜间囊体内部温度降低，对于长航时浮空器来说，昼夜温差导致囊体产生内外压差，且昼夜温差越大囊体内外压差越大，昼夜温差越小内外压差越小。因此，减小囊体昼夜温差至关重要。
文献[5]对比了囊体吸收率、发射率和吸收发射比等辐射参数对浮空器热性能的影响；文献[6]列出了飞艇热特性模型建立的难点问题，并提出了建模的一些建议；文献[7]分析了囊体蒙皮热特性和云层对浮空器热特性的影响；文献[8]建立了飞艇的传热数学模型，仿真分析了蒙皮的温度以及艇内气体的温度；文献[9]建立了浮空器的热动力模型，分析了浮空器在上升和下降过程中的热特性；文献[10]建立了浮空器热力学模型，并且通过地面试验验证了热力学模型的正确性；文献[11]通过地面试验分析了浮空器在不同太阳光辐射和空气对流情况下的热特性，证明了太阳辐射对浮空器热特性起到了关键性的作用；文献[12]建立了浮空器热力学模型，得出囊体内部不同部位的温度相差20 K；文献[13]运用传热学基本原理分析了浮空器热状况的形成机制与影响因素；文献[14]对浮空器上升过程中超冷和下降过程中的超热问题进行了研究；文献[15]采用了25 m和35 m长度的飞艇进行了热力学测试，并成功预测了飞艇的上升过程和飞行轨迹。
Kriging优化方法已经应用于各领域^[16-17]。文献[18-19]将Kriging优化方法应用于气动优化设计；文献[20]将Kriging优化方法应用于结构优化设计，得到了金属加筋板的优化布局；文献[21]将Kriging优化方法应用于机械零件的形状优化，结合有限元分析，得到最优的零件形状。
本文将囊体分为48个部分，每个部分采用不同吸收率和发射率的涂层材料。建立浮空器囊体热力学模型，采用拉丁超立方体方法选取样本点，对样本点进行热力学分析得到响应值，采用Kriging模型进行模拟，并求出囊体最优的材料分布。
1 Kriging模型 Kriging模型利用样本点的响应可以得出最优解。若已知样本点X=[x¹ x² … x^m]^T，xⁱ为n维行向量，其维数为设计变量的个数，对应的响应为Y=[y¹ y² … y^m]^T为q维行向量。第i维响应yⁱ可以表示为

式中：f(β_i, X)为回归模型, β_i为回归系数；z_i(X)为随机过程。
Kriging模型要求预测模型的方差最小，可以得出待测点的响应估计值为

式中：β^*、γ^*为权重系数，样本点一定时，β^*、γ^*为固定值。计算出f(x)和r(x)，即可求出待测点的响应值。
2 浮空器囊体热力学模型 2.1 大气模型 根据国际标准大气模型，不同高度H下的大气参数计算公式为

式中：T为大气温度，K；P为大气压强，Pa；ρ为大气密度，kg/m³。
浮空器飞行在20 km高空，大气环境温度为216.5 K，大气压力为5 475 Pa，大气密度为0.089 kg/m³。
2.2 太阳辐射模型 太阳方位通过太阳高度角α、方位角θ表示，如图 1所示。太阳光线的方向计算公式为



太阳高度角α计算公式为

式中：φ为当地纬度；δ为太阳赤纬角；ω为时角，上午为负，中午为0，下午为正。
太阳赤纬角和时角计算公式为


式中：n为当天在一年中的顺序数，如元旦为n=1；t_s为真太阳时。
太阳方位角θ的计算公式为

太阳对地面的辐射强度Q_D计算公式为

式中：S为大气上方的太阳辐射，根据世界气象组织测量数据，S=1 367 W/m²；B为大气消光系数。
大气散射辐射强度Q_Atm计算公式为

式中：C为常数；ψ为散射面与水平面之间的夹角。
地面反射辐射强度Q_Ref计算公式为

式中：ρ_g为地面的平均反射率，地面的反射率与地区纬度、地面情况等因素有关。北京地区草地地面反射率测量数据为0.174~0.219，混凝土地面反射率可达0.33~0.37，因此本文中计算时选取地面反射率ρ_g=0.2。
2.3 云遮系数 浮空器飞行过程中，空中云的数量直接影响飞艇表面受到的太阳辐射、大气散射和地面反射辐射。云遮系数C_F定义为天空被云层全部遮挡的比例。云遮系数可以通过天空云量计算，本文中分析时假设云遮系数C_F=0.9。
云遮系数对太阳辐射与地面反射辐射的影响如表 1所示。
表 1 云遮系数对太阳辐射的影响 Table 1 Influence of cloud cover coefficient on solar radiation


表选项






浮空器飞行高度为20 km，云层在对流层高空，高度小于20 km，因此本文分析时取H＞H_c。云层中水滴和冰晶的散射，使云层表面形成了比较强的反射面，不同类型的云层反射率均不同，云层的平均反射率为60%^[22]，本文选取ρ_c=0.6。
2.4 浮空器蒙皮辐射
2.4.1 浮空器蒙皮外表面辐射 浮空器蒙皮外表面辐射包括蒙皮外表面与地面和大气之间的长波辐射。
浮空器蒙皮外表面与地面的长波辐射强度计算公式为

式中：ε为蒙皮发射率；T_a为蒙皮温度；T_g为地球表面温度，地球表面平均温度为288 K^[22]；σ为斯蒂芬-玻尔兹曼常数，σ=5.67×10^-8 W/(m²·K⁴)；?_a为接受辐射的表面对地面的角系数; τ_a为红外透射率。

其中：p_a为浮空器周围大气压力；p₀为地球表面的大气压力。
浮空器蒙皮外表面与大气的长波辐射强度计算公式为

式中：φ_a为接受辐射的表面对天空的角系数；T_s为天空的有效温度，计算公式为

其中：T_Atm为大气的温度。

2.4.2 浮空器蒙皮内表面辐射 假设浮空器由N个灰表面组成的封闭腔体，内表面是漫射的。若每个表面的温度T_i已知，任意表面的辐射换热计算公式为

式中：


其中：J_i为离开表面i的辐射，包括自身热辐射和反射；G_i为投入到表面i的辐射；X_{i, j}为表面j对表面i的角系数。
2.5 浮空器蒙皮对流换热
2.5.1 浮空器蒙皮外表面强迫对流换热 浮空器蒙皮外表面对流换热计算公式为

式中：T_{f, i}为蒙皮第i个单元的温度；A_{f, i}为蒙皮第i个单元的面积；h_C为对流换热系数，其计算公式为

其中：λ_air为空气的热传导系数；Pr_air为空气的普朗特数；d为浮空器的直径；Re为空气的雷诺数。

2.5.2 浮空器蒙皮内表面自然对流换热 浮空器内表面与囊体内氦气产生自然对流换热，对流换热计算公式为

式中：T_He为氦气的温度; h_CI为内表面对流换热系数, 其计算公式为

其中：λ_He为氦气的热传导系数；Pr_He为氦气的普朗特数。

2.5.3 氦气温度计算 内部氦气的平均温度计算公式为^[23]

式中：c_v为氦气的定容比热容；m_He为氦气的质量；Q_{CI, He}为氦气与蒙皮交换热量，其计算公式为

将式(23)代入式(22)，可得氦气温度与蒙皮温度的关系表达式为

式中：T_f为蒙皮的平均温度；A为蒙皮的面积。
2.6 浮空器蒙皮热平衡模型 将浮空器囊体划分为N个微元，第i个囊体微元的热平衡方程为

式中：M_{f, i}为第i个蒙皮微元的质量；c_f为材料的比热容；Q_{D, i}为太阳直接辐射；Q_{Atm, i}为大气散射辐射；Q_{Ref, i}为地面反射辐射；Q_{B→g, i}为蒙皮外表面与地面之间的长波辐射；Q_{B→a, i}为蒙皮外表面与大气之间的长波辐射。
3 优化方法 3.1 试验设计方法 囊体热力学优化设计变量为囊体蒙皮涂层材料的吸收率和发射率。优化目标为囊体内氦气昼夜温差最小。
优化(优化流程如图 2所示)分为2个部分：第1部分是囊体热力学计算，第2部分是Kriging模型优化。囊体热力学计算在CFD/FLUENT软件中完成，计算出蒙皮的温度分布和囊体内部气体的平均温度。Kriging模型优化在ISIGHT软件中完成，ISIGHT软件通过样本点的响应值建立Kriging模型，得出最优解。


3.2 优化流程 囊体内部氦气温度与蒙皮温度计算如式(24)所示。可以看出，氦气温度最终与蒙皮温度相等。为了提高计算效率，本文以蒙皮昼夜温差最小为目标，优化完成后，计算氦气昼夜温差。
假设T_D表示N个蒙皮微元白天的平均温度，T_N表示N个蒙皮微元晚上的平均温度。白天T_D的计算公式为

假设?T=T_D-T_N，优化目标?T为最小，即

4 算例分析 4.1 算例描述 算例为正球形囊体，具体尺寸和飞行参数如表 2所示，浮升气体为氦气。
表 2 浮空器参数 Table 2 Parameters of aerostat


表选项






常用蒙皮材料吸收率η、发射率ε的参数值如表 3所示。
表 3 常用材料的热辐射特性参数 Table 3 Thermal radiation characteristic parameters of common materials


表选项






根据浮空器热平衡模型，分析不同吸收率、发射率特性蒙皮材料对氦气平衡温度的影响。分析对象特征如下：取蒙皮材料发射率ε=η+0.5(0.1≤η≤0.4)。
将囊体划分为48个部分，如图 3所示。


4.2 算例仿真 首先检验CFD/FLUENT囊体热特性分析的有效性，根据文献[5]，蒙皮材料吸收率为0.1，发射率为0.6时囊体的昼夜温差最小，囊体的温度分布如图 4所示。


分析可得囊体内部氦气昼夜温差为31 K，与文献[5]相近。证明热力学模型的准确性。
采用拉丁超立方体方法选取样本点500个，取样本点1，其48个面的参数如表 4所示，其囊体温度分布如图 5所示，昼夜温差为54 K。
表 4 样本点1的吸收率 Table 4 Absorption rate of sample point 1


表选项








经过优化后，囊体白天蒙皮温度分布如图 6所示，夜间囊体温度与样本点1相差不大，昼夜温差为28.6 K。优化前，囊体内部氦气的昼夜温差为31 K, 因此经过优化，囊体内部氦气昼夜温差减少7.7%。


5 结论 1) 优化后氦气昼夜温差减小至28.6 K，比采用同一种蒙皮材料加工的囊体昼夜温差减少7.7%。
2) 优化结果为囊体背光部位采用发射率为0.9的材料, 其他部位采用吸收率为0.1的材料。
囊体迎光面和背光面采用不同热特性材料时可降低氦气昼夜温差。囊体背光面与浮空器飞行地点的纬度、太阳赤纬角、时角和囊体姿态有关，在进行囊体实际设计中，可以根据浮空器飞行地点和飞行时间判断背光面的位置。浮空器飞行时可以调节囊体姿态，使得囊体背光面始终背对太阳光。因此，本文提出的优化与设计方法对工程应用有一定的指导意义。

参考文献