精密离心机流-热-固多场耦合计算方法

杜平安，卢 燕，王 珏，凌明祥

（1. 电子科技大学 机械电子工程学院，成都 611731；2. 中国工程物理研究院 总体工程研究所，绵阳 621900）

精密离心机流-热-固多场耦合计算方法

杜平安1，卢 燕1，王 珏2，凌明祥2

（1. 电子科技大学 机械电子工程学院，成都 611731；2. 中国工程物理研究院 总体工程研究所，绵阳 621900）

离心机是惯性导航系统加速度计的标定设备，转盘变形将严重影响标定精度。基于耦合理论研究了离心机流-热-固多场耦合问题的数值计算方法，分析了此类问题的全耦合数学模型，研究了求解数学方程的强、弱耦合算法及其应用。根据耦合方式不同定义了三种分析方法：直接法、顺序迭代法和顺序耦合法，建立了它们与强、弱耦合算法的联系。验证方法正确性后，对比离心机温度场强、弱耦合算法的计算结果发现，强耦合算法计算时间长，但模拟细节好，更接近真实情况，弱耦合算法计算时间短，但精度较低，故采用强耦合算法计算离心机流场和温度场。三种分析方法计算热-结构耦合的结果相近，故选择顺序耦合法计算离心机结构变形。

多场耦合；计算方法；离心机；流场；温度场；结构场

离心机是惯性导航系统加速度计的标定和校准设备[1]，离心机转盘加速度计安装位置的半径误差是影响标定精度的主要因素。对于转盘半径为1.1 m的离心机，1 μm 的半径变化引起的加速度相对误差为0.91×l0-6，因此对于标定精度要求为 10-6级的精密离心机，其半径误差需控制在微米级。45号钢线膨胀系数为11.59×l0-6/℃，线性尺寸为1 m的结构在单位温度变化下线性尺寸变化达11.59 μm。因此，准确计算转盘温度及热变形对控制精密离心机的标定精度具有重要意义。

本文研究离心机流-热-固多场耦合计算方法。首先介绍涉及共轭传热的流-热-固多场耦合数学模型，包括流体流动、热能传输及结构变形相互耦合过程的数学描述；然后研究求解耦合方程的强、弱耦合算法，并提出在软件上应用这两种算法的多场耦合分析方法，定义了直接法、顺序耦合法和顺序迭代耦合法；最后进行了离心机的多场耦合计算，对比了共轭传热强、弱耦合算法的计算结果，并分析了原因，最终计算了离心机转盘的热变形。

1 离心机流-热-固多场耦合现象

离心机由主轴、转盘、止推/径向气体轴承、机箱等组成，图1为结构简图。电机驱动转盘旋转，带动周围空气气团运动；气团由于湍流耗散、摩擦将动能转化为内能，导致流场温升；流场热量一部分被机箱吸收并传递到外界，一部分通过对流换热传给转盘，导致转盘发生热变形。上述物理现象涉及流-热-固三个物理场的耦合。

图1 离心机结构示意图Fig.1 Structure diagram of centrifuge

2 流-热-固多场耦合的数学模型

本文建立的数学模型为流场、温度场和结构场的全耦合模型。流体流动质量守恒方程和动量守恒方程的张量形式为

式中：ρ为流体密度；t为时间；xi为坐标，i=（1,2,3）；u为速度，下标对应x、y、z三个方向的速度分量；Sm为质量源项；p为流体压力；τij为作用在与i方向垂直的平面上 j方向的应力；gi为 i方向重力加速度；Fi为作用外力 i方向分量；μ为流体动力粘度；Sij为流体变形率张量；δij为克罗内尔符号； u′为湍流脉动速度；“ ¯ ”表示平均。

流体热能传输方程为

式中：h为流体比焓，包括内能和压力能；Sh为能量源项；cp为流体定压比热容；T为流体温度，T0为流体参考温度，比焓可由式(5)计算。流体速度和压力分布将影响温度场。

共轭传热指流体和固体界面的热边界条件由热能传输过程动态决定，而不能预先给出。共轭传热时，固体域温度场控制方程为

式中：u为固体自身运动速度；h为固体的焓；T为固体温度；′为固体内部体积热源。其中右边两项分别表示热传导引起的热流以及固体本身的体积热源，左边第二项表示固体自身的运动而引起的对流传热。该方程与流体热能传输方程形式一致，只是默认固体动力粘度为无穷大后，式(4)变成了式(6)。

热-结构耦合涉及热弹性理论问题，求解此类问题关键在于找出满足给定边界条件且适合下述两个方程的解[2-3]。以应变、温差表示应力的广义Hooke定律张量形式为

考虑表面力和体积力热弹性力学平衡微分方程为

式中：σ为应力；G为材料剪切弹性模量；ε为应变；δ为泊松比；e为体积应变，e=εx+εy+εz；β为材料线膨胀系数，ΔT为温度变化量；λ为拉梅常数，λ=2Gδ/(1-2δ)；Fi为体积力的i方向分量。

流体运动会影响流体和固体的温度分布，若流体动力粘度μ是温度的函数，温度场也将影响流体运动，温升会导致结构热变形，变形的结构将改变流体域大小，因此三个物理场相互影响，变量间相互作用，故在相应的方程中都有体现对方的变量。求解这样的流-热-固耦合问题需要联立求解上述耦合控制方程。

3 耦合方程的求解算法

3.1 强耦合算法

求解耦合控制方程有强耦合算法和弱耦合算法。强耦合法将各场控制方程耦合到同一个方程矩阵中，直接联立求解，在同一个时间步内同时求解所有场的全部未知量，能一次性求解出所有场的响应[4]。

理论上强耦合现象需采用强耦合算法才能得到准确的解，但由于流体和结构的刚度不同，强耦合算法在求解耦合方程时容易出现病态矩阵，故不能充分利用现有计算软件，难以解决一般的耦合问题。

3.2 弱耦合算法

弱耦合算法是在同一时间步内依次对各场控制方程求解，通过交换耦合面的数据得到收敛结果，从而实现耦合求解。该类算法可以充分利用现有软件，可避免强耦合算法无法或难以求解的问题。

弱耦合算法由于两个场分开计算，其离散方法不同导致网格不一致，故需要在两类网格间建立对应关系再传递界面数据，因此弱耦合算法包含搜索算法和映射算法。

3.2.1 搜索算法

用来搜索单元和点的算法称为搜索算法[5]，目前常用的有强力搜索和桶式搜索。强力搜索即当前节点在另一个网格所有的单元上执行一个循环，根据最小距离准则在搜索的结果中选择最佳匹配，见图 2(a)。节点N1位于单元e1和e2中，根据准则会映射到最小距离的单元e1上。当节点不能映射到任何单元时，就根据最近位置准则选择最佳匹配，见图2(b)，节点N1没有映射到任一单元，就映射到最近节点N1́上。

以流体单元搜索对应的结构节点为例，桶式搜索过程见图3。首先把流体单元分为三个桶，桶内不止一个单元，且包含周围结构节点；然后找到结构节点 N1，并确认其在桶Box3内，再使用强力搜索算法依次查找桶内三个单元：e4、e5、e6，计算得到节点与单元e5的距离最近，根据最小距离准则把N1映射到单元e5上。

3.2.2 映射算法

数据映射从数学上可看作数据在不同网格间的相互插值问题，因此研究映射算法即研究插值算法。在建立节点-单元间对应关系后，需要在映射点完成插值计算。在耦合界面上需满足变形协调、能量守恒及作用力平衡条件。目前有局部插值法和全局插值法。图4为不同网格间数据传递示意图。

全局插值法将发送端的每一个节点映射到接收端的单元上，然后在所有节点数据（已知）的基础上构造一个整体上的数据与坐标的插值函数，再根据此数据函数可直接求得插值点的数据。全局插值法可避免节点单元间的搜索过程，相对局部插值法对耦合界面的网格形状及质量要求均较低，如径向基函数是常见的全局插值法。局部插值法是接收端的每一个节点映射到发送端的单元上，然后只选取在插值点附近的一些节点或者单元上的数据（已知）进行插值。局部插值法易理解，应用方便，且只需少量单元数据，但容易导致局部总力不平衡，整体位移协调性较差。

图2 强力搜索算法Fig.2 Powerful search algorithm

图3 桶式搜索算法Fig.3 Bucket-type search algorithm

图4 不同网格间数据传递Fig.4 Data transmission between different grids

3.3 建立分析方法与耦合算法联系

求解耦合方程有强、弱耦合算法，而耦合算法在软件上的应用可通过多种分析方法实现，因此不同分析方法应用的耦合算法也不尽相同。本文将某种分析方法应用的耦合算法描述为某分析方法是强/弱耦合。

直接法只包含一次分析，分析计算结束后能得到多个场的响应，它通过计算包含所有未知量的单元矩阵或单元载荷向量来实现耦合。直接法不一定是强耦合，当所有场的全部变量是在同一时间步内求解，它是强耦合，否则为弱耦合。如用Fluent求解共轭换热问题时，采用耦合求解器则是强耦合，采用分离求解器则是弱耦合。

顺序迭代法通常不止包含一次分析，在每一时间步交换耦合界面信息，把上场计算结果作为载荷传递到下一个场，在界面迭代直至得到收敛结果。顺序耦合法是先计算一个场，计算结束后直接将结果作为载荷传递到下一个场，其根本是只对单个场的一次求解，数据传递无完整迭代收敛过程。

顺序迭代法和顺序耦合法是弱耦合，前者较后者更接近真实情况，能应用于耦合较为强烈的场合，后者只能应用于耦合效应不明显的场合。

4 精密离心机流-热-固多场耦合计算

目前没有软件能采用强耦合算法直接求解流-热-固耦合问题，对于涉及共轭传热的流-热-固耦合通常分为共轭传热求解和热-结构耦合求解。离心机的多场耦合计算只考虑流场、温度场和结构场的单向耦合。

4.1 本文方法验证

首先将流场和温度场的测量值与仿真值对比，以验证本文采用的仿真方法。离心机转盘直径2.2 m，流场直径 2.5658 m，机箱高度 0.514 m。其流体域采用MRF方法模拟，流场在机箱底部为压力出口，转速为270 r/min，环境温度为20℃，机箱与外界自然对流。

图5为流场监测点位置示意图，其中：点A距机箱底部62 mm，距主轴中心990 mm；点B位于负载面（加速度计安装位置）正上方，距主轴中心726.8 mm，距机箱顶部62 mm；点C与点B处于同一平面，在接口的正上方，距主轴中心550 mm。

表1为流场测量值与仿真计算值的对比，可以看到三个点的计算误差都在10%以内。

另外，还监测了负载面加速度计安装位置的温度，监测点温度值为20.50℃，温升0.5℃，仿真计算温升值为0.53℃，相对误差6%。因此认为强耦合算法计算离心机流场、温度场方法可行。

图5 流场监测点位置示意图Fig.5 The sketch map of flow field monitoring

表1 流场监测点实验值与仿真值对比Tab.1 Monitoring value contrast of flow field

4.2 共轭传热强、弱耦合算法计算

分别采用强、弱耦合算法计算离心机300 r/min工况下的共轭换热。其中，强耦合算法利用有限体积法计算共轭传热问题，弱耦合算法采用顺序耦合分析法实现，采用有限体积法和有限元法结合。

两种算法的计算网格都采用四面体非结构化网格，并加密耦合界面。强耦合法网格见图6(a)，网格长宽比9.2、纵横比0.77、数量483万；弱耦合法流体域网格见图6(b)，网格长宽比8.95、纵横比0.78、数量396万。

图6 离心机流场、温度场网格划分结果Fig. 6 Grid of centrifuge flow field and temperature field

图 7(a)为离心机流场温度分布。强耦合法计算的温度较高，原因是考虑了固体，热阻增大，功耗不变的条件下温升会增大；图7(b)为转盘温度分布及点（0, -1.206, -0.01）到点（0, 1.206, -0.01）温度分布曲线，弱耦合法转盘最大温升为 0.655℃，强耦合法为0.792℃，二者相差0.137℃，结合45号钢的线膨胀系数，其产生的热变形为1.6 μm，这对于满足10-6级的标定精度要求很勉强，因此，弱耦合算法不能满足离心机仿真计算要求。两者温度分布规律一致，最高温出现在整流罩背风面，强耦合法模拟的温度稍高，这是由于其流场温度也较高。图7(c)为整流罩温度分布，强耦合法在整流罩整个背风面温度都较高，而弱耦合法温度梯度较大，可见后者模拟细节也不如强耦合法。

综上，离心机温流场、温度场计算可采用强耦合算法，弱耦合算法（应用顺序耦合分析法）则不适用。

图7 离心机温度场计算结果对比Fig.7 Contrast on calculation results of centrifuge temperature field

图8 环肋散热器结构示意图Fig.8 Structure of ring rib radiator

4.3 热-结构耦合计算

本节首先以散热器为例，分析了直接法、顺序迭代法和顺序耦合法计算热应力的精度。

图8为散热器结构示意，取图中阴影部分为研究对象。环肋及管道均为不锈钢，热导率25.96 W/(m•℃)，弹性模量1.93×105 MPa，管内为热流体，温度250℃，对流换热系数249.23 W/(m2•℃)，管内压力6.89 MPa；管外为空气，温度 39℃，对流换热系数 62.3 W/(m2•℃)。如表2所示，可见三种方法的计算结果相差不大。考虑到计算时间、资源等，本文采用顺序耦合法计算离心机的热变形。

以强耦合算法计算出的温度作为载荷传递给转盘，计算得到离心机柱坐标下的热变形如图9所示。转盘热变形量越远离旋转中心越大，最大变形发生在整流罩上，为6.44 μm，这是一个变形累积结果，加速度计安装位置热变形量为4.9 μm。

离心机转盘的变形实际上是温升、离心载荷、静不平衡量、回转误差共同作用结果。图9为离心机动态半径测试示意图[9]，采用定位局部平均测试法，在转盘边缘设计定位槽，利用电容传感器连续测试转盘边缘半径变化量。在300 r/min工况下圆弧面半径变化量为 15.6 μm，不考虑温升时计算出的机械变形量为11.3μm，本文计算的热变形量为4.8 μm，因此综合变形为16.1 μm，相对实验值误差为3.21%。因此可以说明本文热变形多物理场仿真计算是合理的。

表2 三种耦合方法的计算结果对比Tab.2 Contrast on three analysis methods

图9 离心机热变形云图Fig.9 Nephogram of centrifuge’s thermal deformation

图10 动态半径挖槽测试示意图Fig.10 Excavated test of dynamic radius

5 结 论

本文研究了流-热-固多物理场耦合的数学模型及求解算法。强耦合算法直接求解多场耦合方程，弱耦合算法则将上场计算结果作为载荷传递给下一物理场，精度较前者低，但适用范围更广。分析了直接法、顺序迭代法和顺序耦合法的特点，直接法能实现强耦合计算，其余只能实现弱耦合计算。

对比离心机温度场强、弱耦合算法的计算结果，得出强耦合算法模拟流-热耦合细节更好，结果更真实。通过三种分析方法计算环肋散热器，表明三种方法计算热-固耦合相差不大，考虑到计算资源等问题，选择顺序耦合法计算离心机的热变形。

本文利用热-流-固耦合方法计算了离心机由于气流运动的湍流耗散、摩擦引起的转盘热变形。转盘变形是多因素作用的结果。计算转盘综合变形后便可获得动态半径相对误差，进而可根据该误差在测量系统中进行补偿，以满足精密离心机的测量精度要求。

（References）：

[1] Fang Q, Huang S X. UKF for integrated vision and inertial sensors based on three-view geometry[J]. IEEE Sensors Journal, 2013, 13(7): 2711-2719.

[2] Larin V B, Tunik A .A. On correcting the system of inertial navigation[J]. Journal of Automation and Information Sciences, 2010, 42(8): 13-26.

[3] Zhang Qun, Toshiaki H. Studies of the strong coupling and weak coupling methods in FSI analysis[J]. International Journal for Numerical Method in Engineering, 2004, 60: 2013-2029.

[4] 孔繁余, 夏斌, 张慧, 等. 基于热-结构耦合的热水循环泵结构强度研究[J]. 机械工程学报, 2014, 50(18): 155-162. Kong Fan-yu, Xia Bin, Zhang Hui, et al. Research on structural strength for hot water circulation pump based on thermal-structural coupling[J]. Journal of Mechanical Engineering, 2014, 50(18): 155-165.

[5] Malatip A, Wansophark N, Dechaumphai P. Fractional four-step finite element method for analysis of thermally coupled fluid-solid interaction problems[J]. Applied Mathematics and Mechanics, 2012, 33(1): 99-116.

[6] Deng W, Xi W B. Multi-physics coupling simulation of resistive heat exchanger in EAST.[J]. Atomic Energy Science and Technology, 2015, 49(7): 1330-1337.

[7] 李勇, 林缅, 张召彬. 流-热-固耦合渗流的数学模型及其应用[J]. 水动力学研究与进展, 2015, 30(1): 56-63. Li Yong, Lin Mian, Zhang Zhao-bin. Numerical model of thermal-hydrological-mechanical coupling and its application [J]. Journal of Hydrodynamics, 2015, 30(1): 56-63.

[8] 孙学功. 多物理场耦合界面数据传递的仿真研究[J].计算机仿真, 2015, 32(1): 23-28. Sun Xue-gong. Simulation of data exchange of multiphysics coupling surfaces[J]. Computer Simulation, 2015, 32(1): 23-28.

[9] 成永博. 基于多学科统一建模的精密离心机动态半径误差分离仿真研究[D]. 成都: 中国工程物理研究院, 2013. Cheng Yong-bo. Simulation research on dynamic radius error separation techniques of precision centrifuge based on the general multi-discipline modeling[D]. Chengdu, China Academy of Engineering Physics, 2013.

[10] Vlad G, Umesh J, Nick H. Coupled fluid-structure transient thermal analysis of a gas turbine internal air system with multiple cavities[J]. Journal of Engineering for Gas Turbines and Power, 2012, 134(10): 102508.

[11] 李德才, 魏艳勇, 韦俊新, 等. 基于流固耦合理论的长寿命陀螺仪浮子内部流场及温度场分析[J]. 中国惯性技术学报, 2007, 15(6): 721-724. Li De-cai Wei Yan-yong Wei Jun-xin, et al. Analysis on temperature field and inner flow field of long-life gyro′s floater with Fluid-Structure coupling theory[J]. Journal of Chinese Inertial Technology, 2007, 15(6): 721-724.

Numerical calculation methods of multi-physics coupling of precision centrifuge

DU Ping-an1, LU Yan1, WANG Jue2, LING Ming-xiang2
(1. School of Mechatronics Engineering, University of Electronic Science and Technology of China, Chengdu 611731, China; 2. Institute of Systems Engineering, China Academy of Engineering Physics, Mianyang 621900, China)

A centrifuge is a calibrating device for the accelerometer in inertial navigation system, and the deformation of the rotator could seriously affect the calibrating accuracy. The numerical methods of fluid-heatsolid physics coupling of precision centrifuge are studied based on the coupling theory. The mathematical models of multi-physics coupling are presented, and the strong and weak coupling algorithms to solve the coupling-equations are introduced. Three analysis methods (i.e. the direct method, the sequence iteration method and the sequential coupling method) are defined according to different coupling modes, and the relationship with strong /weak coupling algorithm is given. After validating the proposed methods, the results of centrifuge temperature field by the two algorithms are compared, and it is found that the computation time of the strong coupling algorithm is longer, but the simulation results are better and are closer to the real situation. The calculation time of the weak coupling algorithm is shorter, but its accuracy is lower. Thus the strong coupling algorithm is adopted to simulate the centrifuge’s fluid field and temperature field. Due to the results from three analysis methods are close to each other, the sequence coupling method is selected to calculate the thermal deformation.

multi-physics coupling; calculation method; centrifuge; fluid field; temperature field; structure field

U666.1

A

1005-6734(2016)02-0275-06

10.13695/j.cnki.12-1222/o3.2016.02.025

2015-11-25

2016-03-17

国家重大科学仪器设备开发专项（2011YQ130047）

杜平安（1962—），男，教授，博士生导师，从事机电系统多场耦合仿真研究。E-mail: dupingan@uestc.edu.cn