液浮陀螺仪的流固热仿真及特性分析

古 旻, 张 晶, 樊景松, 苏 岩

（1. 南京理工大学机械工程学院, 南京 210094;2. 北京航天控制仪器研究所, 北京 100039）

0 引言

液浮陀螺仪具有体积小、 精度高、 动态性能好、 过载能力强等特点, 广泛应用于船舶导航、导弹的精确打击以及航天航空领域［1-2］。 在实际工况下, 由于内部热源发热和外部加热片的保温加热等以热传导和对流换热方式向其他部件散发热量, 使液浮陀螺仪内部存在着温度梯度, 陀螺仪内部温度不同导致浮油密度差异, 进而引起浮液流动, 影响了液浮陀螺仪的精度［3］。 随着液浮陀螺仪精度要求的不断提高, 仪表内部温度分布对精度的影响愈发明显, 研究人员越来越重视液浮陀螺仪的温度特性。 在掌握仪表温度特性的基础上,对结构进行合理的改进, 并采用精准的温控系统,便能提高陀螺仪温度分布的均匀性, 减小由温度因素带来的干扰力矩, 提高仪表的精度［4-5］。 就目前而言, 相关技术已经成为制约我国惯导精度提高的关键问题之一, 而有限元分析是研究此类问题的关键工具［6-7］。

针对以上问题, 国内外学者进行了大量研究。杨盛林等［8］利用有限元软件ANSYS 得到了陀螺仪壳体的温度分布。 白永杰等［9］运用有限元方法对液浮陀螺仪浮子进行了热应力分析, 得到了工作环境下浮子结构的最大热应力值, 为陀螺仪设计时材料的选择提供了重要依据, 但其未考虑内部流体对浮子的应力情况。 卜石等［10］使用ANSYS Fluent 有限元分析软件得到陀螺仪整体的温度梯度,但将浮液视为固体, 并未考虑其流场对温度分布的影响。 牛文韬等［11］基于ANSYS/CFX 对液浮陀螺仪进行流固热耦合仿真, 分别计算了浮筒外侧的仪表部件-浮油耦合模型以及浮子-内部气体耦合模型, 得到了陀螺仪的温度场以及内部流体的流场。

本文使用多物理场有限元仿真软件对陀螺仪整表进行计算, 分别考虑了其内部流体及固体的特点, 首先介绍了液浮陀螺仪的工作原理, 其次分析了流体及传热问题有限元仿真计算的理论基础, 然后建立了液浮陀螺仪流固热多物理场耦合的仿真模型, 根据实际工作情况分析或理论计算得到模型的边界条件, 最终计算并分析了内部流体的流场以及仪表整体的温度场特性, 为实现液浮陀螺仪的分区温控以及仪表改进提供依据。

1 液浮陀螺仪工作原理

液浮陀螺仪利用高速旋转刚体具有的定轴性和进动性测量角速度, 本文所研究的液浮陀螺仪马达为动压马达, 采用气浮轴承支持, 马达转子高速旋转运动产生较大角动量, 形成陀螺效应。该陀螺仪由三大组件组成: 壳体组件、 波纹管组件以及浮子组件。 浮子组件包括马达和浮子框架,其内部充满气体。 在工作状态下, 马达转子与定子通过高速的相对运动在两者之间形成气膜, 实现非接触悬浮。 壳体上的法兰环通过固定螺钉与钛合金平台固紧, 壳体与浮子之间的间隙充满浮油, 壳体外周均匀缠绕电阻加热丝, 浮子右转子上同样缠有发热线圈, 通过调节这两热源的电流大小对陀螺仪进行温度控制, 进而使浮油的浮力等于浮子的重力, 最大程度地减小摩擦, 提高陀螺仪精度［12］。 图1 为液浮陀螺仪内部结构示意图,由于陀螺仪内部结构复杂, 固体、 浮油和马达内部气体相互影响, 考虑陀螺仪内部散热的同时又要保证浮子的重力等于浮油的浮力, 因此对液浮陀螺仪进行整体有限元仿真分析是陀螺仪设计过程中至关重要的一步。

图1 液浮陀螺仪内部结构示意图Fig.1 Schematic diagram of liquid floated gyroscope internal structure

2 仿真计算基础

本文的仿真计算主要涉及陀螺仪内部流体的流场以及仪器整体的温度场, 结合多物理场有限元软件仿真过程中的实际情况与陀螺仪的工作原理, 现做以下假设:

1）对结构进行简化, 忽略零件的倒角、 不重要的孔等特征, 将非关键组件如接线柱、 螺钉等与基础零件合为一体, 看成一个整体（假设1）;

2）由于陀螺仪的整体温度范围较低, 不考虑热辐射的影响（假设2）;

3）只考虑流体的重力, 不考虑结构的重力、 热变形和热应力（假设3）。

根据以上假设, 建立如下的控制方程。

2.1 流体流动控制方程

在流动过程中, 流体需要始终满足连续性方程、 动量守恒定律和能量守恒定律, 它们是陀螺仪内部流体进行有限元仿真的理论基础［13-14］。

（1）连续性方程

式（1）中,ρ为流体密度,t为时间,Δ为拉普拉斯算子,U为流体的流速。

（2）动量守恒定律

式（2） 中,ui、uj为速度张量,p为静压力,xi、xj为坐标张量,τij为应力张量,ρgi为重力体积力,Fi为体积力。

（3）能量守恒定律

式（3）中,ht为流体的总能,ΔT为温度差,λf为等效导热系数,Sh为外部热源,τ为黏性流体的剪切应力。

2.2 陀螺仪内部传热方程

传热是指当有温度梯度出现时所发生的能量转移。 在自然界中, 热传导、 热对流和热辐射被称为最基本的传热途径, 任何传热问题都可由其中某种或多种组合而成。 基于假设2 以及陀螺仪实际工作情况, 其内部的传热途径主要包括: 固体与固体之间的热传导、 固体与氦气和浮油之间的对流换热。以下是两种传热途径的数学描述, 它们是陀螺仪内部传热问题仿真计算的理论基础［15-16］。

（1）陀螺仪内部固体间热传导

式（4）中,Q为热流量,k为导热系数,ΔT为温度差,A为传热面积。

（2）对流换热

式（5）中,q为对流换热量,Ch为对流换热系数,Tss为固体表面温度,Tsf为接触面流体温度,c为流体热容,u为近壁面流体流速,T+为接触面无量纲温度［17］。

3 有限元模型的建立

上述方程是仿真计算的理论基础, 对于简单情况可以通过上述公式直接求得, 但对于本文所研究的陀螺仪而言, 其结构不规律, 流体运动与传热过程十分复杂, 无法求得解析解。 因此, 需要借助有限元法对上一节中的方程进行处理, 方可求得结果。 液浮陀螺仪有限元模型的建立过程如下。

3.1 结构建模与网格划分

基于假设1, 在多物理场有限元仿真软件中建立简化后的三维模型, 如图2 所示。 其中, 重力方向为-z方向。

陀螺仪各零组件的材料类型及物理属性直接决定了其温度场, 由前文的流体流动控制方程以及传热方程可知, 在仪表的仿真计算过程中需要知道材料的密度、 动力黏度等参数。 经统计整理后, 计算中需要使用到的材料及其相关参数如表1所示。

表1 材料参数表Table 1 Parameters of material

由于三维模型结构复杂, 因此主要以自由四面体为基本网格单元对整个模型进行网格划分。同时, 对关键部位进行局部网格加密, 使用多物理场有限元仿真软件划分网格, 得到浮油的网格划分模型, 如图3 所示。 得到网格的单元数为27320394 个, 平均网格质量为0.7。

图3 浮油网格划分Fig.3 Schematic diagram of oil slick meshing

3.2 流场及热边界条件设置

陀螺仪内部流场包括两部分: 气体流场和浮油流场。 由于高速旋转的马达定子会带动周围气体运动, 基于流体密度、 动力黏度以及流速等参数通过下式可计算得到气体的流场雷诺数（Re ＞2000）, 因此使用湍流模型进行计算。 浮油主要受温度梯度影响发生流动, 流动速率较小, 通过下式可计算得到浮油的流场雷诺数（Re ＜2000）, 采用层流模型进行计算。

式（6） 中,ρ、μ为流体的密度和动力黏性系数,U、L为流场的特征速度和特征长度。

液浮陀螺仪的热源主要包括三部分: 马达热源、 主动温控热源和黏性损耗。 马达热功耗约为3.5W; 主动温控热源包括马达线圈发热、 浮子右转子线圈发热以及壳体外加热片的保温加热, 各部分的热功率如表2 所示; 黏性损耗产生的热功率只需在有限元仿真软件中打开, 系统会自动计算。

表2 热源功率表Table 2 Power of heat sources

考虑到陀螺仪工作状态下的外罩与大气接触,故可将陀螺仪外罩处的对流换热视为大空间自然换热。 根据传热学相关公式［17］, 可计算得到格拉晓夫数Gr 和努塞尔数Nu

式（8）、 式（9）中,g为重力加速度,a为空气体积膨胀系数, Δt为温度差,v为动力黏度,l为特征尺度,D为特征直径,L为特征长度。

进一步得到陀螺仪表面对流换热系数

式（9）中,ε为空气导热系数,L为特征长度。

3.3 计算条件

设置模型的初始温度为实际工作环境温度328K, 在启动一段时间后, 陀螺仪内部流场及温度场会趋于稳定状态, 因此设置求解类型为冻结转子（稳态）, 同时为保证计算结果的精确程度, 相对容差设为0.001 作为收敛条件。

4 流场与温度场计算结果

计算结束后累计迭代总次数为160 次, 通过温度实验测量得到某点的温度为329K, 通过对多物理场有限元仿真模型的流场及温度场进行稳态计算, 得到该点的温度为330.176K, 误差为0.36%。两者温度相近, 误差在可接受的范围之内, 可认为陀螺仪温度仿真分析的结果合理可信。 进一步针对模型网格进行优化后, 网格细化20%, 陀螺仪整体温度范围如表3 所示, 表明仿真模型的网格单元加密20%后, 仿真温度梯度变化量为0.05%。同时, 选取关键部位, 如图4 所示, 得到其温度曲线。 由图4 可知, 仿真模型在关键部位的50μm 空间间隔内, 温度梯度分辨率达到0.01K。

表3 网格细化前后温度对比Table 3 Temperature comparison before and after grid refinement

图4 关键部位及其温度曲线Fig.4 Schematic diagram of key part and its temperature curve

4.1 气体与浮油流场

由式（5）可知, 固体与流体之间的热传导需要知道流体的流速, 因此有必要对陀螺仪流体的流场进行仿真计算。 计算得到陀螺仪内部流体的流场如图5 和图6 所示。

图5 气体流场流速图Fig.5 Flow velocity diagram of gas flow field

图6 浮油流场流速图Fig.6 Flow velocity diagram of oil slick flow field

由图5 可知, 高速旋转的马达转子带动周围的气体一起运动, 形成涡流。 气体的最大流速为67.40m/s, 在转子距轴线半径最大处。

由图6 可知, 中段浮油即浮筒与壳体之间微小间隙内的浮油运动相对比较稳定, 而左右段浮油在非均匀热场中的运动十分复杂。 气体的最大流速为3.04 ×10-6m/s, 在右端盖过油孔处。

4.2 温度场

对整表进行计算, 得到液浮陀螺仪垂直输入轴方向的温度分布切面, 如图7 所示。 由图7 可知, 仪表在正常工作条件下的温度范围为328K ～341.38K, 温度最高点位于马达定子与转子高速运动时形成的气膜处, 温度最低点位于法兰环与外界相接触的位置。 马达内部气体的温度云图如图8所示, 浮油的温度云图及部分关键位置的温度如图9 所示。 由图8、 图9 可知, 计算得到的气体温度范围为338.43K ～341.32K, 温差为2.89K; 浮油的温度范围为337.23K ～338.84K, 温差为1.61K, 系统自动计算得到流体流动过程中产生的黏性损耗占总热量的13.90%。 若不考虑马达转子与定子之间气体体积变化产生的热量, 陀螺仪整体温度范围为328K ～335.49K, 可知马达内部气体收缩膨胀产生的热量约占总热量的56%。 图9 中的两个曲线图为靠近壳体部位的浮油温度分布图,上面的曲线图位置为远离法兰环挂耳的位置, 法兰环正对点温度为338.43K; 下面的曲线图位置为靠近法兰环挂耳的位置, 法兰环正对点温度为338.29K。 从中可知在关键位置处, 温度精度可达0.01K 以下, 且法兰环挂耳位置对浮油温度梯度影响较大, 靠近挂耳位置与远离挂耳位置的温度相差0.14K。

图7 整表的温度梯度图Fig.7 Temperature gradient diagram of the whole liquid floated gyroscope

图9 浮油温度梯度云图及关键位置温度Fig.9 Nephogram of temperature gradient for oil slick and temperatures at key positions

5 结论

本文通过使用多物理场有限元仿真软件, 建立了液浮陀螺仪整表的流固热耦合多物理场计算模型, 充分考虑了陀螺仪固体与流体之间的相互作用, 结合实际工作环境, 得到模型准确的计算边界条件, 最终计算得到陀螺仪内部温度分布情况以及内部流体的运动传热状态, 为后续陀螺仪温控系统设计提供依据。 具体可得到以下结论:

1）法兰环挂耳对液浮陀螺仪浮筒外侧浮油温度梯度影响显著, 可有针对性地设计控温系统;

2）流体在固体接触表面由于摩擦产生的热量约占总热量的13.90%, 马达内部气体收缩膨胀产生的热量约占总热量的56%, 其余热量为主动热源产生;

3）浮油温度梯度较明显, 最高温与最低温相差1.61K, 会造成浮油密度分布不均, 引起浮油流动, 进而影响浮子外侧的应力分布和陀螺仪精度。