基于COMSOL的線性壓縮機動力學(xué)仿真與實驗驗證

孔家煊，趙欽宇，晁春陽，周文杰，王 博，程 路，甘智華*，徐 軍

（1.浙江大學(xué)制冷與低溫研究所 浙江省制冷與低溫技術(shù)重點實驗室，杭州 310027；2.杭州電子科技大學(xué)能量利用與自動化研究所，杭州 310018；3.日照華斯特林科技有限公司，山東 日照 276800）

0 引言

近年來，由于紅外技術(shù)和超導(dǎo)領(lǐng)域的持續(xù)發(fā)展，低溫制冷機的需求日益增長[1?2]。斯特林制冷機作為低溫制冷機的重要分支，具有本征效率高、穩(wěn)定性好、結(jié)構(gòu)緊湊和使用壽命長等優(yōu)點。斯特林制冷機的驅(qū)動源是高頻的交變壓力，壓縮活塞和氣缸間的側(cè)向力不利于延長壽命和減小噪聲，因此使用直線電機驅(qū)動壓縮機，不但可以省去運動轉(zhuǎn)換裝置；同時能通過調(diào)節(jié)電壓來調(diào)節(jié)活塞行程[3]。

建立合理的數(shù)值模型對斯特林制冷機以及線性壓縮機的設(shè)計非常重要，許多研究者已經(jīng)開展了相關(guān)工作。Clark將線性壓縮機的機械常量、電參數(shù)、動子質(zhì)量等效到一個電路等效模型中，以匹配一臺斯特林制冷機為例，通過非線性靜磁場的有限元方法，以最小線圈電阻損耗和最小設(shè)備體積為目標，對線性壓縮機進行了設(shè)計[4]。Almajri等[5]用COMSOL建立了斯特林制冷機的計算流體力學(xué)（CFD）模型，通過實驗驗證后使用該模型研究了相角、孔隙率、回?zé)崞鏖L度對制冷量的影響。上海技術(shù)物理研究所的鄧偉峰等[6]采用Maxwell建立了自由活塞式線性壓縮機電磁?機械?聲耦合模型，模擬了不同負載條件下壓縮機電磁參數(shù)和動力特性的瞬態(tài)變化規(guī)律。

目前，在回?zé)崾降蜏刂评錂C設(shè)計領(lǐng)域，一維數(shù)值計算軟件SAGE在制冷機端較為完善，而在壓縮機端僅提供了幾個簡單的接口，簡化了電磁場的計算[7]，無法對壓縮機的設(shè)計起到良好的指導(dǎo)作用。為了彌補SAGE在壓縮機建模上的不足，本文采用COMSOL對線性壓縮機進行建模仿真，并通過實驗對仿真模型進行驗證，以期為線性壓縮機的設(shè)計和優(yōu)化改進提供一種簡易有效的方案。

1 建模

高頻回?zé)崾降蜏叵到y(tǒng)中各狀態(tài)量在時域中都可看作交變的正弦波[8]，基本的表示形式為：

利用歐拉公式改寫式中動態(tài)項Xd：

其中：表示交變的物理量，可根據(jù)正弦波的定義將其分解為Xˉ與動態(tài)Xd。對于制冷系統(tǒng)，運行過程中存在多個交變量，在時域中無法直觀的表現(xiàn)各物理信號之間的關(guān)系。一種更常用且直觀的做法是將時域上的交變信號映射到頻域中。在復(fù)平面中，動態(tài)項映射的結(jié)果是以幅值Xa為長度，ωt+φ為方向的向量。

基于以上原理，在復(fù)平面內(nèi)建立壓縮活塞的質(zhì)量?彈簧?阻尼模型，該模型在復(fù)平面中某一時刻的受力分析如圖 1 所示[9]，其中點折線箭頭表示此時刻位移、速度和加速度的相位，實線箭頭表示壓縮活塞在工作狀態(tài)下所受的力，表示電機推力，表示彈性力，表示慣性力，表示阻尼力，表示氣體力。

圖1 復(fù)平面活塞受力分析圖Fig.1 Force analysis of piston in the complex plane

式中：pa和xa分別為壓力和位移在時域中的幅值，A代表活塞截面積；ω為線性壓縮機的角頻率；θ為回?zé)崞鳠岫速|(zhì)量流領(lǐng)先壓力波的相位角，對于斯特林制冷機，實際設(shè)計中應(yīng)將其控制在40°至45°，從而減小回?zé)崞鲹p失；kp和Rp為合并常數(shù)項后的結(jié)果，分別代表氣體彈性剛度和氣體阻尼系數(shù)。二者表征了制冷機側(cè)對線性壓縮機的影響，通過改變壓力幅值和相位角可以模擬壓縮機在連接不同工況制冷機時的運行情況。

COMSOL是一款基于有限元方法的多物理場仿真軟件，通過其內(nèi)置的模塊可以便捷地實現(xiàn)對電磁學(xué)、流體流動、傳熱學(xué)等領(lǐng)域的仿真和對不同物理場的耦合。本文采用COMSOL對一款商用Cryo S 100斯特林制冷機動磁式線性壓縮機工作過程中的電磁場及壓縮機在電磁場中的運動狀況進行了動態(tài)分析。

線性壓縮機的裝配模型與仿真模型的對應(yīng)關(guān)系如圖2所示，運動活塞由板彈簧支撐，與氣缸內(nèi)壁形成間隙密封；永磁體2是活塞的一部分，位于內(nèi)定子1和齒形磁軛3之間的氣隙內(nèi)。當(dāng)齒形磁軛包裹的線圈4中通入正弦電流時，永磁體會受到交變的軸向電磁力從而帶動活塞做往復(fù)直線運動。

圖2 線性壓縮機裝配模型與仿真模型對照圖Fig.2 The assembly model and simulation model of linear compressor

建立仿真模型的參數(shù)如表1所示。為了提高計算效率，對裝配模型進行以下的假設(shè)和簡化：

（1）忽略模型在徑向上的運動；

（2）簡化活塞的完整形狀，將所有的質(zhì)量和受力集中到永磁體上；

（3）假定板彈簧對活塞提供的軸向力與位移是線性相關(guān)的，且摩擦因素被簡化為固定的機械阻尼系數(shù)；

（4）使用相對磁導(dǎo)率為1的氣體介質(zhì)填充模型組件以外的區(qū)域。

表1 線性電機的基本參數(shù)Tab.1 Primary structure parameter of linear motor

使用不同的磁性本構(gòu)關(guān)系描述磁場中軟磁材料和永磁體的磁性特征。對于內(nèi)定子和磁軛使用的Somaloy軟磁材料，通過導(dǎo)入圖3所示的外部的材料磁化屬性表，對其磁性能進行了描述。對于永磁體，指定了剩余磁感應(yīng)強度關(guān)系式，如式（4）所示，其中表示剩余磁感應(yīng)強度，μ0表示真空磁導(dǎo)率，μr表示永磁體的相對磁導(dǎo)率，表示徑向矯頑力。對應(yīng)該系統(tǒng)中使用的永磁體，μr=1.02，

圖3 軟磁材料B-H曲線Fig.3 B-H curve of Somaloy

根據(jù)圖1，使用“全局常微分和微分代數(shù)方程”接口導(dǎo)入帶氣體力的運動方程描述活塞的機械運動。電磁場和運動方程通過位移量x和電磁力FI相互耦合，在線圈內(nèi)指定輸入電流，可以求解電磁力、位移量和速度量隨時間的變化關(guān)系。

為仿真模型添加移動網(wǎng)格時，考慮到運動部件只有活塞，可以將模型如圖4所示分為內(nèi)定子、永磁體和包裹線圈的磁軛三個部分。僅在中間的永磁體區(qū)域使用移動網(wǎng)格接口。在該區(qū)域內(nèi)，將圖4中紅色線框所包裹的永磁體和其周圍的一部分氣體作為整體，指定其軸向位移為微分方程接口求解得到的位移量x，上下方的氣體域設(shè)置為被動變形的自由變形域，并使用了映射網(wǎng)格來減少由于網(wǎng)格劃分帶來的計算誤差，防止發(fā)散。

圖4 模型網(wǎng)格剖分Fig.4 model meshed geometry

電磁場的模擬必須使用麥克斯韋應(yīng)力張量對活塞受到的電磁力進行精確的求解，因此在永磁體和氣隙部分繪制了更為精密的網(wǎng)格。同時為了更好地保證移動域和非移動域之間物理量的連續(xù)性，在二者相交的邊界上，移動域一側(cè)必須進行更精細的邊界單元劃分。最終網(wǎng)格的平均質(zhì)量為0.82。經(jīng)無關(guān)性驗證，進一步加密網(wǎng)格對求解結(jié)果的影響不超過5%，求解器步長設(shè)置為0.001 s。

2 實驗流程

2.1 空載實驗裝置

圖5為測試空載工況活塞位移的裝置示意圖，線性壓縮機由變頻電源驅(qū)動，可直接讀出輸入電壓、輸出電流和有效功率；壓縮機活塞板簧軸向上方裝有激光位移傳感器，能對活塞板彈簧的運動進行實時測量，活塞板彈簧的中心孔通過緊固螺母和活塞相連接，二者的運動始終具有一致性；在空載條件下設(shè)置輸出電壓和頻率為80 Hz、7 V，記錄此時壓縮機活塞的位移結(jié)果，與仿真結(jié)果進行對比。

圖5 空載工況下活塞位移測試裝置實物圖和示意圖Fig.5 Physical diagram and schematic diagram of piston displacement test device under on-load condition

2.2 比推力實驗裝置

將壓縮機直線電機平穩(wěn)放置在拉壓力測試儀的托盤中心。直線電機由直流電源驅(qū)動，可以直接讀出輸入電壓和輸入電流；拉壓力測試儀在壓縮活塞軸向上方連接有壓力觸頭，可以通過手輪控制的限位機構(gòu)精確限制壓縮活塞的位置，并對活塞的推力進行測量。

測量時，設(shè)置直流電源輸出電流依次為1 A、2 A、3 A、4 A、5 A；記錄活塞在不同限位條件下推力與電流的變化關(guān)系，從而對仿真結(jié)果進行對比驗證。

3 討論

3.1 空載實驗

在模型中輸入與實驗相同的空載工況，對活塞的位移情況進行了計算。采用測試裝置對多臺壓縮機進行了空載測試，測試結(jié)果如圖6所示。在隨機抽取的50臺線性壓縮機中，有四臺壓縮機的活塞位移與計算值出現(xiàn)較大偏差，其余樣機的測試結(jié)果平均分布在計算值附近。進一步比較后發(fā)現(xiàn)，出現(xiàn)較大偏差的四臺線性壓縮機的活塞與排出器之間的位移相位差遠大于其余線性壓縮機的位移相位差，推測是在裝配過程中零件安裝不當(dāng)產(chǎn)生的影響。剔除錯誤樣本后，測試值與計算值的平均偏差為3.49%。

圖6 空載工況下活塞位移的結(jié)果分析圖Fig.6 The result of displacement in no-load condition

3.2 比推力測試

在SAGE模型中，壓縮機端的影響是通過指定比推力引入的，為驗證壓縮機模型的合理性，使用前文所述測量方式，依次測量了活塞在距離平衡位置0 mm、1 mm、2 mm、3 mm、4 mm處的比推力。同時在模型中將壓縮活塞的永磁體固定在對應(yīng)的位置，設(shè)定線圈的電流參數(shù)，從穩(wěn)態(tài)的角度對壓縮活塞在不同位置上的比推力進行求解，結(jié)果如圖7所示，實測值略大于計算值，最大偏差出現(xiàn)在離初始位置1 mm處，為7.2%；最小偏差出現(xiàn)在離初始位置4 mm處，為4.37%，無論是實測值還是計算值，比推力的最大變化不超過0.4 N∕A。在模型中進一步增大與初始點的距離后，發(fā)現(xiàn)當(dāng)壓縮活塞位移超過6 mm時，活塞會由于逐漸脫離線圈生成的磁場，導(dǎo)致比推力迅速衰減。

3.3 負載條件下與Sage模型比較

該模型中，制冷機側(cè)的影響被簡化為氣體力，為了驗證這種方法的可行性，使用SAGE模型計算了該款制冷機在100 K、193 K、253 K和293 K四個工況點下的動態(tài)壓力幅值、壓縮機出口的壓力?質(zhì)量流相位、阻抗等參數(shù)，將它們換算為氣體力的形式輸入COMSOL模型進行計算，并對比了用兩種模型計算的活塞位移幅值。該制冷機使用表1所列的線性電機。計算結(jié)果如圖8和表2所示，在不同的冷頭極限溫度下，對壓縮活塞位移的計算結(jié)果均吻合的較好，最大偏差不超過3%。

圖7 比推力隨位移的變化關(guān)系比較曲線Fig.7 The measurement and simulated result of specific thrust vs.displacement

圖8 兩種模型計算結(jié)果對比曲線Fig.8 Comparison of displacement in two models

表2 基于壓力幅值的位移幅值比較結(jié)果Tab.2 Comparison of displacement base on pressure amplitude

4 結(jié)論

本文基于斯特林制冷機電?機械?聲的耦合原理，將聲阻抗拆分為位移和速度方向分量，建立了帶負載條件的直線電機COMSOL模型，并進行了空載測試、比推力測試和與SAGE模型的對比驗證。其中空載測試的50臺電機中，測試值與計算值的平均偏差為3.49%；對位于磁場不同位置處活塞的比推力測試表明，測試值和計算值的最大偏差出現(xiàn)在距離平衡位置1 mm處，為7.2%，最小偏差出現(xiàn)在距離平衡位置4 mm處，為4.37%。計算結(jié)果表明，當(dāng)壓縮活塞距離平衡位置6 mm以上時，比推力將迅速衰減；在不同負載條件下，該模型的活塞位移計算結(jié)果和SAGE的計算結(jié)果均吻合良好，最大偏差不超過3%。