吉 然，陈超飞，何佳豪，崔梦祺，陈代芬，刘少俊,3

(1.江苏科技大学能源与动力学院，江苏 镇江 212003)(2.中船重工第七〇三研究所无锡分部，江苏 无锡 214026)(3.浙江大学能源清洁利用国家重点实验室国家环境保护燃煤大气污染控制工程技术中心，浙江 杭州 310027)

制冷系统在舰船航行过程中为人员正常工作和物资储存提供基本保障，制冷剂是制冷系统中必不可少的一部分。舰船在复杂的冲击环境中，不可避免地会发生冷媒的泄漏，以一体式空调系统为例，每发生15%的制冷剂设计量泄漏[1]，能效比(coefficient of performance，COP)就降低5%，从而影响其正常工作，因此需要及时补充制冷剂。由此可知，舰船制冷剂储存容器在设计时，必须考虑航行中的船体颠簸使一些部件或设备产生过大应力造成永久断裂或者形变[2]，从而出现制冷剂泄漏的情况。本文利用动态设计分析方法(dynamic design analysis method，DDAM)和动态时域仿真，对参照压力容器[3]标准设计出的制冷剂储罐进行研究，考察最大应力响应，确定其强度是否满足要求。

1 制冷剂储罐参数计算

1.1 储罐结构基本参数

本文设计可重复充装制冷剂R404A的14 L舰船用制冷剂储罐，其工作压力为0.5 MPa，设计压力为2.0 MPa，工作温度和设计温度都为-10～40 ℃，其主要设计参数：储罐的圆筒壁外径D0为200 mm，工作容积V为14 L，高度L为450 mm。储罐为单极孔的合金制储罐。

1.2 储罐的材料和强度

对于船用设备，应首先采用强度高、延伸率大于10%的材料，避免使用脆性以及对缺口敏感的材料。本文选取的材料为TC4钛合金，TC4强度高、密度小，强度超过了大部分合金结构钢的强度。热处理后按照表1所示物理性能参数取值来校核储罐罐身屈服强度和抗拉强度值。

表1 TC4钛合金物性参数

1.3 结构简图

储罐结构简图如图1所示，δe=δh=0.8 mm，L=450 mm，D0=200 mm+2δh，H=75 mm+δh，h=25 mm。其中，δe为筒体罐身厚度，δh为储罐封头厚度，H为封头高度与直边段之和，h为直边段高度。

图1 储罐结构简图

2 有限元模型的建立及分析

2.1 储罐模型建立

储罐模型具有结构对称性，其完整的三维模型可以看作由1/2的模型绕储罐中轴旋转得到。因为其余结构在物理性质、约束施加以及载荷分布方面基本一致，在添加边界条件的时候，只需施加具有连续性的位移约束[4]就可以使用一半的模型进行模拟仿真，极大减少了计算量。利用ANSYS Workbench建立的模型如图2所示。

图2 制冷剂储罐1/2结构示意图

2.2 网格划分

本文采用Adaptive的方法进行网格划分：先从边开始划分网格，再在曲率较大的地方细化边网格。划分完成后得到如图3所示的四面体网格，节点数为87 428个，单元数为43 100个。

图3 模型整体网格示意图

2.3 施加边界条件

有限元模型的边界条件是由储罐的实际约束条件和加载条件决定的，本文根据储罐模型底部的结构特征和受力特点，对模型施加边界约束条件：上截面施加y轴固定约束，对称截面上采用完全约束；瓶底内壁面按分析需要施加均匀载荷(压力) 以及用R404A液体密度定义的静水压力载荷。

2.4 网格无关性验证

网格划分的方法不同以及数量差异会对数值模拟结果的计算产生一定影响，从而对分析的准确性产生干扰[5]。下面对采用相同方法划分的5类单元数不一致的网格模型进行无关性检测，网格的单元数见表2。

表2 网格无关性检验结果

由表可知，网格单元数对综合应力最大值的影响不可忽视，当网格单元数从低于一万的A类单元数增大至C类单元数时，综合应力最大值上升50 MPa，占A类网格最大值的15.7%；而当划分类型超过C类网格到达E类时，网格单元数增加近200%，综合应力最大值仅上升10 MPa，占C类网格最大值的2%。由此可以认为，单元数增大到一定值后，计算结果趋于稳定，综合考虑后，确定E类网格单元数可满足无关性要求。

3 求解与分析

3.1 3D模型静态强度分析

依据胡金伯格方程按照3D封头椭圆实体形状建模：

(1)

(2)

式中：σφ为径向应力；σθ为环向应力；R2为薄壁封头上各点的曲率半径；P为容器内压；t为容器壁厚；a,b为容器剖面椭圆长、短轴；x为封头剖面各点至中轴线的距离。σφ在X-Y坐标系中的方向如图4所示，而σθ始终在K点处与封头剖面保持垂直状态。

图4 薄壁结构应力分析示意图

由公式(1)、(2)可知椭球壳应力与P和t有关，与长轴与短轴之比a/b有关。当a=b时，球壳趋于正圆球壳，最大应力为圆筒壳中周向应力的一半；而a/b增大时，椭球壳中应力持续增大，直到a/b比值达到1.414时，椭圆封头赤道处的σθ将从拉应力变为压应力。

3.2 ANSYS计算结果与胡金伯格方程解的对比

通过有限元分析软件中的Normal Stress后处理模块对椭圆球壳与直边段的交界处，即封头赤道处任意一点沿X和Y轴方向的受力进行了数值模拟分析，结果如图5，6所示。

图6 无冲击径向应力云图

在封头赤道截面处任一点，取x=a，将内压P与壁厚t代入方程，依据建模数据取a/b为2，计算结果对比见表3。表中RATIO为仿真计算结果与方程计算结果的比值。

表3 理论计算与模拟仿真计算结果的对比

由表3可见，无载荷加速度时，赤道处一点上的径向应力、环向应力的理论计算结果在数值和方向上与软件模拟结果相互吻合，且最大数值差距不超过12%，可见理论方法与有限元方法均具有一定的参考性。

3.3 基于DDAM的储罐模态分析

ANSYS中设定最大分析阶数为8，根据相关规范[6]，保留3个冲击方向上模态质量大于10%的所有模态，结果见表4。在通过模态扩展之后，表中仅保留了主要振型的有效质量，对应的主要模态振型如图7所示。

根据标准GJB 1060.1—91[7]可得关于水面舰船的冲击设计值，见表5，表中的v0为基准速度，m/s。制冷剂储罐为甲板部位安装设备，其设计冲击载荷谱由式(3)确定。

表4 不同冲击方向上的主要模态信息

图7 各阶振型图

(3)

式中：aa为设计加速度，m/s2；a0为基准加速度，m/s2；ma为设备有效模态质量。

表5 船体安装部位的设计冲击载荷谱

将式(3)得到的设计加速度依次代入表5， 可以得到不同冲击方向上的主要载荷谱，见表6。表中aD为最终设计加速度。

表6 不同冲击方向上的设计载荷谱

由图8可以看出，在横向和纵向冲击响应的应力最大值出现在椭圆封头与直边段的交界处，分别为380.46 MPa和381.19 MPa，而垂向冲击响应的应力最大值出现在了储罐接口附近，达到了791.76 MPa。因此对R404A制冷剂储罐来说，在8.42 Hz的1阶垂向冲击环境下的响应虽然符合强度规范，但却是最为危险的冲击环境，需要引起重视；而2阶横向和6阶纵向其次，表现较为安全。

图8 3个冲击方向储罐综合应力响应云图

3.4 储罐时域动态冲击分析

时域动态冲击分析，即对3D模型进行短时间内的瞬态动力学分析，通过输入数据作为时间函数的载荷，研究结构在随时间任意变化的载荷作用下，模型非线性响应的过程[8]。

根据前联邦德国国防军舰艇建造规范BV/0430冲击安全性[9]以及相关文献资料[10]，将冲击设计谱值转化成加速度时程曲线函数，对3D模型施加0.03 s内Y方向上的垂向冲击加速度，如图9所示。

图9 动态冲击加速度曲线图

为了便于得到各变量的时间历程曲线，可以将Solution结果链接至Mechanical APDL中利用ANSYS经典中的命令流对结果进行时间历程后处理。本文选取1489号节点作为研究对象，位置如图10所示。

本文利用ANSYS经典中的POST26命令对结果进行时间历程后处理，通过对该节点微分得到该节点的速度信息，再次微分得到加速度信息，最后结果如图11所示。

图10 1489号节点位置示意图

图11 X向、Y向、Z向加速度时域曲线

从图中可以看到3个方向的加速度中，Y向的加速度变化最大，加速度响应幅值最大达到1 233 m/s2，X向在0.005 s附近达到响应极限，而Z向的加速度响应幅值波动较小。

1489号节点位于上封头内部与筒体交接附近的应力带区域，从图12可看到该节点的应力在300 MPa到350 MPa之间上下浮动，与图13数值基本接近。

图12 节点最大合成应力时域曲线

由图9可见，在BV/0430中5.3 ms正向冲击载荷达到58g的最高值，而综合应力最大值仅为371 MPa，低于节点单向冲击考察方法峰值约62.9%。

4 结论

本文主要对理论设计的舰船制冷剂储罐进行动态设计分析方法的有限元仿真和动态时域冲击分析，通过得到的最大应力响应值与分布来分析、考察在设计工况下的冲击响应是否符合要求。

图13 上封头内部应力云图

1)经过分析发现，在8.42 Hz的垂向冲击下，局部最大应力达到了791.76 MPa，制冷剂储罐结构处于较为危险的状态，但符合强度规范。

2)考察点单方向冲击输入会影响其余方向的应力响应，产生小幅波动。从考察点综合应力来看，距离屈服强度依旧剩余61.4%的安全裕量。

3)动态设计分析方法显示纵向和横向的综合应力最大值均出现在容器内壁对接的焊缝融合区域，在产品成型加工过程中，往往会在这些区域产生残余拉伸应力。而该区域往往也是介质压力引起最大拉伸应力的地方，会严重影响结构的疲劳寿命。分析可知，若要提高该类储罐的使用期限，在加工过程中需要通过热处理工艺减少相关区域的残余应力。