杜平安，盧 燕，王 玨，凌明祥

（1. 電子科技大學 機械電子工程學院，成都 611731；2. 中國工程物理研究院 總體工程研究所，綿陽 621900）

精密離心機流-熱-固多場耦合計算方法

杜平安1，盧 燕1，王 玨2，凌明祥2

（1. 電子科技大學 機械電子工程學院，成都 611731；2. 中國工程物理研究院 總體工程研究所，綿陽 621900）

離心機是慣性導航系統(tǒng)加速度計的標定設備，轉(zhuǎn)盤變形將嚴重影響標定精度?；隈詈侠碚撗芯苛穗x心機流-熱-固多場耦合問題的數(shù)值計算方法，分析了此類問題的全耦合數(shù)學模型，研究了求解數(shù)學方程的強、弱耦合算法及其應用。根據(jù)耦合方式不同定義了三種分析方法：直接法、順序迭代法和順序耦合法，建立了它們與強、弱耦合算法的聯(lián)系。驗證方法正確性后，對比離心機溫度場強、弱耦合算法的計算結果發(fā)現(xiàn)，強耦合算法計算時間長，但模擬細節(jié)好，更接近真實情況，弱耦合算法計算時間短，但精度較低，故采用強耦合算法計算離心機流場和溫度場。三種分析方法計算熱-結構耦合的結果相近，故選擇順序耦合法計算離心機結構變形。

多場耦合；計算方法；離心機；流場；溫度場；結構場

離心機是慣性導航系統(tǒng)加速度計的標定和校準設備[1]，離心機轉(zhuǎn)盤加速度計安裝位置的半徑誤差是影響標定精度的主要因素。對于轉(zhuǎn)盤半徑為1.1 m的離心機，1 μm 的半徑變化引起的加速度相對誤差為0.91×l0-6，因此對于標定精度要求為 10-6級的精密離心機，其半徑誤差需控制在微米級。45號鋼線膨脹系數(shù)為11.59×l0-6/℃，線性尺寸為1 m的結構在單位溫度變化下線性尺寸變化達11.59 μm。因此，準確計算轉(zhuǎn)盤溫度及熱變形對控制精密離心機的標定精度具有重要意義。

本文研究離心機流-熱-固多場耦合計算方法。首先介紹涉及共軛傳熱的流-熱-固多場耦合數(shù)學模型，包括流體流動、熱能傳輸及結構變形相互耦合過程的數(shù)學描述；然后研究求解耦合方程的強、弱耦合算法，并提出在軟件上應用這兩種算法的多場耦合分析方法，定義了直接法、順序耦合法和順序迭代耦合法；最后進行了離心機的多場耦合計算，對比了共軛傳熱強、弱耦合算法的計算結果，并分析了原因，最終計算了離心機轉(zhuǎn)盤的熱變形。

1 離心機流-熱-固多場耦合現(xiàn)象

離心機由主軸、轉(zhuǎn)盤、止推/徑向氣體軸承、機箱等組成，圖1為結構簡圖。電機驅(qū)動轉(zhuǎn)盤旋轉(zhuǎn)，帶動周圍空氣氣團運動；氣團由于湍流耗散、摩擦將動能轉(zhuǎn)化為內(nèi)能，導致流場溫升；流場熱量一部分被機箱吸收并傳遞到外界，一部分通過對流換熱傳給轉(zhuǎn)盤，導致轉(zhuǎn)盤發(fā)生熱變形。上述物理現(xiàn)象涉及流-熱-固三個物理場的耦合。

圖1 離心機結構示意圖Fig.1 Structure diagram of centrifuge

2 流-熱-固多場耦合的數(shù)學模型

本文建立的數(shù)學模型為流場、溫度場和結構場的全耦合模型。流體流動質(zhì)量守恒方程和動量守恒方程的張量形式為

式中：ρ為流體密度；t為時間；xi為坐標，i=（1,2,3）；u為速度，下標對應x、y、z三個方向的速度分量；Sm為質(zhì)量源項；p為流體壓力；τij為作用在與i方向垂直的平面上 j方向的應力；gi為 i方向重力加速度；Fi為作用外力 i方向分量；μ為流體動力粘度；Sij為流體變形率張量；δij為克羅內(nèi)爾符號； u′為湍流脈動速度；“ ˉ ”表示平均。

流體熱能傳輸方程為

式中：h為流體比焓，包括內(nèi)能和壓力能；Sh為能量源項；cp為流體定壓比熱容；T為流體溫度，T0為流體參考溫度，比焓可由式(5)計算。流體速度和壓力分布將影響溫度場。

共軛傳熱指流體和固體界面的熱邊界條件由熱能傳輸過程動態(tài)決定，而不能預先給出。共軛傳熱時，固體域溫度場控制方程為

式中：u為固體自身運動速度；h為固體的焓；T為固體溫度；′為固體內(nèi)部體積熱源。其中右邊兩項分別表示熱傳導引起的熱流以及固體本身的體積熱源，左邊第二項表示固體自身的運動而引起的對流傳熱。該方程與流體熱能傳輸方程形式一致，只是默認固體動力粘度為無窮大后，式(4)變成了式(6)。

熱-結構耦合涉及熱彈性理論問題，求解此類問題關鍵在于找出滿足給定邊界條件且適合下述兩個方程的解[2-3]。以應變、溫差表示應力的廣義Hooke定律張量形式為

考慮表面力和體積力熱彈性力學平衡微分方程為

式中：σ為應力；G為材料剪切彈性模量；ε為應變；δ為泊松比；e為體積應變，e=εx+εy+εz；β為材料線膨脹系數(shù)，ΔT為溫度變化量；λ為拉梅常數(shù)，λ=2Gδ/(1-2δ)；Fi為體積力的i方向分量。

流體運動會影響流體和固體的溫度分布，若流體動力粘度μ是溫度的函數(shù)，溫度場也將影響流體運動，溫升會導致結構熱變形，變形的結構將改變流體域大小，因此三個物理場相互影響，變量間相互作用，故在相應的方程中都有體現(xiàn)對方的變量。求解這樣的流-熱-固耦合問題需要聯(lián)立求解上述耦合控制方程。

3 耦合方程的求解算法

3.1 強耦合算法

求解耦合控制方程有強耦合算法和弱耦合算法。強耦合法將各場控制方程耦合到同一個方程矩陣中，直接聯(lián)立求解，在同一個時間步內(nèi)同時求解所有場的全部未知量，能一次性求解出所有場的響應[4]。

理論上強耦合現(xiàn)象需采用強耦合算法才能得到準確的解，但由于流體和結構的剛度不同，強耦合算法在求解耦合方程時容易出現(xiàn)病態(tài)矩陣，故不能充分利用現(xiàn)有計算軟件，難以解決一般的耦合問題。

3.2 弱耦合算法

弱耦合算法是在同一時間步內(nèi)依次對各場控制方程求解，通過交換耦合面的數(shù)據(jù)得到收斂結果，從而實現(xiàn)耦合求解。該類算法可以充分利用現(xiàn)有軟件，可避免強耦合算法無法或難以求解的問題。

弱耦合算法由于兩個場分開計算，其離散方法不同導致網(wǎng)格不一致，故需要在兩類網(wǎng)格間建立對應關系再傳遞界面數(shù)據(jù)，因此弱耦合算法包含搜索算法和映射算法。

3.2.1 搜索算法

用來搜索單元和點的算法稱為搜索算法[5]，目前常用的有強力搜索和桶式搜索。強力搜索即當前節(jié)點在另一個網(wǎng)格所有的單元上執(zhí)行一個循環(huán)，根據(jù)最小距離準則在搜索的結果中選擇最佳匹配，見圖 2(a)。節(jié)點N1位于單元e1和e2中，根據(jù)準則會映射到最小距離的單元e1上。當節(jié)點不能映射到任何單元時，就根據(jù)最近位置準則選擇最佳匹配，見圖2(b)，節(jié)點N1沒有映射到任一單元，就映射到最近節(jié)點N1?上。

以流體單元搜索對應的結構節(jié)點為例，桶式搜索過程見圖3。首先把流體單元分為三個桶，桶內(nèi)不止一個單元，且包含周圍結構節(jié)點；然后找到結構節(jié)點 N1，并確認其在桶Box3內(nèi)，再使用強力搜索算法依次查找桶內(nèi)三個單元：e4、e5、e6，計算得到節(jié)點與單元e5的距離最近，根據(jù)最小距離準則把N1映射到單元e5上。

3.2.2 映射算法

數(shù)據(jù)映射從數(shù)學上可看作數(shù)據(jù)在不同網(wǎng)格間的相互插值問題，因此研究映射算法即研究插值算法。在建立節(jié)點-單元間對應關系后，需要在映射點完成插值計算。在耦合界面上需滿足變形協(xié)調(diào)、能量守恒及作用力平衡條件。目前有局部插值法和全局插值法。圖4為不同網(wǎng)格間數(shù)據(jù)傳遞示意圖。

全局插值法將發(fā)送端的每一個節(jié)點映射到接收端的單元上，然后在所有節(jié)點數(shù)據(jù)（已知）的基礎上構造一個整體上的數(shù)據(jù)與坐標的插值函數(shù)，再根據(jù)此數(shù)據(jù)函數(shù)可直接求得插值點的數(shù)據(jù)。全局插值法可避免節(jié)點單元間的搜索過程，相對局部插值法對耦合界面的網(wǎng)格形狀及質(zhì)量要求均較低，如徑向基函數(shù)是常見的全局插值法。局部插值法是接收端的每一個節(jié)點映射到發(fā)送端的單元上，然后只選取在插值點附近的一些節(jié)點或者單元上的數(shù)據(jù)（已知）進行插值。局部插值法易理解，應用方便，且只需少量單元數(shù)據(jù)，但容易導致局部總力不平衡，整體位移協(xié)調(diào)性較差。

圖2 強力搜索算法Fig.2 Powerful search algorithm

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

圖4 不同網(wǎng)格間數(shù)據(jù)傳遞Fig.4 Data transmission between different grids

3.3 建立分析方法與耦合算法聯(lián)系

求解耦合方程有強、弱耦合算法，而耦合算法在軟件上的應用可通過多種分析方法實現(xiàn)，因此不同分析方法應用的耦合算法也不盡相同。本文將某種分析方法應用的耦合算法描述為某分析方法是強/弱耦合。

直接法只包含一次分析，分析計算結束后能得到多個場的響應，它通過計算包含所有未知量的單元矩陣或單元載荷向量來實現(xiàn)耦合。直接法不一定是強耦合，當所有場的全部變量是在同一時間步內(nèi)求解，它是強耦合，否則為弱耦合。如用Fluent求解共軛換熱問題時，采用耦合求解器則是強耦合，采用分離求解器則是弱耦合。

順序迭代法通常不止包含一次分析，在每一時間步交換耦合界面信息，把上場計算結果作為載荷傳遞到下一個場，在界面迭代直至得到收斂結果。順序耦合法是先計算一個場，計算結束后直接將結果作為載荷傳遞到下一個場，其根本是只對單個場的一次求解，數(shù)據(jù)傳遞無完整迭代收斂過程。

順序迭代法和順序耦合法是弱耦合，前者較后者更接近真實情況，能應用于耦合較為強烈的場合，后者只能應用于耦合效應不明顯的場合。

4 精密離心機流-熱-固多場耦合計算

目前沒有軟件能采用強耦合算法直接求解流-熱-固耦合問題，對于涉及共軛傳熱的流-熱-固耦合通常分為共軛傳熱求解和熱-結構耦合求解。離心機的多場耦合計算只考慮流場、溫度場和結構場的單向耦合。

4.1 本文方法驗證

首先將流場和溫度場的測量值與仿真值對比，以驗證本文采用的仿真方法。離心機轉(zhuǎn)盤直徑2.2 m，流場直徑 2.5658 m，機箱高度 0.514 m。其流體域采用MRF方法模擬，流場在機箱底部為壓力出口，轉(zhuǎn)速為270 r/min，環(huán)境溫度為20℃，機箱與外界自然對流。

圖5為流場監(jiān)測點位置示意圖，其中：點A距機箱底部62 mm，距主軸中心990 mm；點B位于負載面（加速度計安裝位置）正上方，距主軸中心726.8 mm，距機箱頂部62 mm；點C與點B處于同一平面，在接口的正上方，距主軸中心550 mm。

表1為流場測量值與仿真計算值的對比，可以看到三個點的計算誤差都在10%以內(nèi)。

另外，還監(jiān)測了負載面加速度計安裝位置的溫度，監(jiān)測點溫度值為20.50℃，溫升0.5℃，仿真計算溫升值為0.53℃，相對誤差6%。因此認為強耦合算法計算離心機流場、溫度場方法可行。

圖5 流場監(jiān)測點位置示意圖Fig.5 The sketch map of flow field monitoring

表1 流場監(jiān)測點實驗值與仿真值對比Tab.1 Monitoring value contrast of flow field

4.2 共軛傳熱強、弱耦合算法計算

分別采用強、弱耦合算法計算離心機300 r/min工況下的共軛換熱。其中，強耦合算法利用有限體積法計算共軛傳熱問題，弱耦合算法采用順序耦合分析法實現(xiàn)，采用有限體積法和有限元法結合。

兩種算法的計算網(wǎng)格都采用四面體非結構化網(wǎng)格，并加密耦合界面。強耦合法網(wǎng)格見圖6(a)，網(wǎng)格長寬比9.2、縱橫比0.77、數(shù)量483萬；弱耦合法流體域網(wǎng)格見圖6(b)，網(wǎng)格長寬比8.95、縱橫比0.78、數(shù)量396萬。

圖6 離心機流場、溫度場網(wǎng)格劃分結果Fig. 6 Grid of centrifuge flow field and temperature field

圖 7(a)為離心機流場溫度分布。強耦合法計算的溫度較高，原因是考慮了固體，熱阻增大，功耗不變的條件下溫升會增大；圖7(b)為轉(zhuǎn)盤溫度分布及點（0, -1.206, -0.01）到點（0, 1.206, -0.01）溫度分布曲線，弱耦合法轉(zhuǎn)盤最大溫升為 0.655℃，強耦合法為0.792℃，二者相差0.137℃，結合45號鋼的線膨脹系數(shù)，其產(chǎn)生的熱變形為1.6 μm，這對于滿足10-6級的標定精度要求很勉強，因此，弱耦合算法不能滿足離心機仿真計算要求。兩者溫度分布規(guī)律一致，最高溫出現(xiàn)在整流罩背風面，強耦合法模擬的溫度稍高，這是由于其流場溫度也較高。圖7(c)為整流罩溫度分布，強耦合法在整流罩整個背風面溫度都較高，而弱耦合法溫度梯度較大，可見后者模擬細節(jié)也不如強耦合法。

綜上，離心機溫流場、溫度場計算可采用強耦合算法，弱耦合算法（應用順序耦合分析法）則不適用。

圖7 離心機溫度場計算結果對比Fig.7 Contrast on calculation results of centrifuge temperature field

圖8 環(huán)肋散熱器結構示意圖Fig.8 Structure of ring rib radiator

4.3 熱-結構耦合計算

本節(jié)首先以散熱器為例，分析了直接法、順序迭代法和順序耦合法計算熱應力的精度。

圖8為散熱器結構示意，取圖中陰影部分為研究對象。環(huán)肋及管道均為不銹鋼，熱導率25.96 W/(m?℃)，彈性模量1.93×105 MPa，管內(nèi)為熱流體，溫度250℃，對流換熱系數(shù)249.23 W/(m2?℃)，管內(nèi)壓力6.89 MPa；管外為空氣，溫度 39℃，對流換熱系數(shù) 62.3 W/(m2?℃)。如表2所示，可見三種方法的計算結果相差不大?？紤]到計算時間、資源等，本文采用順序耦合法計算離心機的熱變形。

以強耦合算法計算出的溫度作為載荷傳遞給轉(zhuǎn)盤，計算得到離心機柱坐標下的熱變形如圖9所示。轉(zhuǎn)盤熱變形量越遠離旋轉(zhuǎn)中心越大，最大變形發(fā)生在整流罩上，為6.44 μm，這是一個變形累積結果，加速度計安裝位置熱變形量為4.9 μm。

離心機轉(zhuǎn)盤的變形實際上是溫升、離心載荷、靜不平衡量、回轉(zhuǎn)誤差共同作用結果。圖9為離心機動態(tài)半徑測試示意圖[9]，采用定位局部平均測試法，在轉(zhuǎn)盤邊緣設計定位槽，利用電容傳感器連續(xù)測試轉(zhuǎn)盤邊緣半徑變化量。在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 動態(tài)半徑挖槽測試示意圖Fig.10 Excavated test of dynamic radius

5 結 論

本文研究了流-熱-固多物理場耦合的數(shù)學模型及求解算法。強耦合算法直接求解多場耦合方程，弱耦合算法則將上場計算結果作為載荷傳遞給下一物理場，精度較前者低，但適用范圍更廣。分析了直接法、順序迭代法和順序耦合法的特點，直接法能實現(xiàn)強耦合計算，其余只能實現(xiàn)弱耦合計算。

對比離心機溫度場強、弱耦合算法的計算結果，得出強耦合算法模擬流-熱耦合細節(jié)更好，結果更真實。通過三種分析方法計算環(huán)肋散熱器，表明三種方法計算熱-固耦合相差不大，考慮到計算資源等問題，選擇順序耦合法計算離心機的熱變形。

本文利用熱-流-固耦合方法計算了離心機由于氣流運動的湍流耗散、摩擦引起的轉(zhuǎn)盤熱變形。轉(zhuǎn)盤變形是多因素作用的結果。計算轉(zhuǎn)盤綜合變形后便可獲得動態(tài)半徑相對誤差，進而可根據(jù)該誤差在測量系統(tǒng)中進行補償，以滿足精密離心機的測量精度要求。

（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] 孔繁余, 夏斌, 張慧, 等. 基于熱-結構耦合的熱水循環(huán)泵結構強度研究[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] 李勇, 林緬, 張召彬. 流-熱-固耦合滲流的數(shù)學模型及其應用[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] 孫學功. 多物理場耦合界面數(shù)據(jù)傳遞的仿真研究[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] 成永博. 基于多學科統(tǒng)一建模的精密離心機動態(tài)半徑誤差分離仿真研究[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] 李德才, 魏艷勇, 韋俊新, 等. 基于流固耦合理論的長壽命陀螺儀浮子內(nèi)部流場及溫度場分析[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

國家重大科學儀器設備開發(fā)專項（2011YQ130047）

杜平安（1962—），男，教授，博士生導師，從事機電系統(tǒng)多場耦合仿真研究。E-mail: dupingan@uestc.edu.cn