謝東升,劉 禎,李 正
(1.湖北文理學院 純電動汽車動力系統設計與測試湖北省重點實驗室,湖北 襄陽 441053;2.湖北文理學院 汽車與交通工程學院, 湖北 襄陽 441053)
隨著膨脹機應用范圍越來越廣和應用要求不斷提高,就透平式膨脹機和活塞式膨脹機而言,渦旋膨脹機作為一種新型容積式膨脹機,具有結構緊湊、體積小、重量輕、能量轉換效率高以及可處理小氣量膨脹工質等優點[1-2],更能滿足人們對其更高的要求,因此近年來在空調制冷、壓縮空氣儲能技術、有機朗肯循環、低溫余熱回收等方面得到廣泛應用[3-6]。在渦旋膨脹機工作過程中,高溫高壓氣體經吸氣口進入推動渦旋盤運動,不斷增大進氣腔與膨脹腔的體積,從而實現氣體體積的膨脹,并將氣體的壓力能轉化為機械能[7],在實現整個吸氣、膨脹、排氣的周期性運作中,膨脹機動靜渦旋盤工作環境復雜,不僅受到2個渦旋盤之間的摩擦,還承受著較高的熱負荷和氣體力,這使得渦旋盤尤其是渦旋齒極易發生變形。變形會使軸向間隙變大,也會對2個渦旋齒的嚙合效果以及徑向間隙的密封產生影響,從而影響膨脹機的工作性能和可靠性,因此,渦旋膨脹機渦旋盤的變形研究十分重要。
目前,關于渦旋盤的變形研究大都集中在對渦旋盤進行氣體力載荷和熱載荷分析,得到相應的應力和變形,且研究對象多為渦旋壓縮機。郭守寧等[8]通過ANSYS分析軟件對氣體力作用下的渦旋盤的應力與變形規律作了探討;楊廣衍等[9-10]對溫度場下動靜渦旋盤產生的熱應力和熱變形進行了仿真研究;劉振全等[11]對動渦旋盤進行了有限元劃分,對其受力下的剛度和強度進行了分析;劉濤等[12]對變截面動渦旋盤的熱應力和熱變形進行了研究。后來,不少學者又通過耦合的方法對渦旋盤多場載荷下的變形進行了研究,如劉國平等[13-15]對變化的多場載荷進行耦合分析,探究渦旋動盤的應力變形情況。以上研究雖對渦旋盤設計優化有一定的指導作用,但渦旋膨脹機與渦旋壓縮機的工作過程以及工作條件不同,兩者之間承受的溫度載荷和氣體力載荷差異較大[16],所以不能用渦旋壓縮機的結果來指導渦旋膨脹機的設計優化。
此外,目前對渦旋膨脹機的研究大多為膨脹機內部流場的數值模擬研究[17-19],或者非穩態流動特性分析[20-22],對渦旋盤的變形研究較少。此前,劉禎等[16]通過對比2種不同型線下的流場特性,探討了不同型線參數對所受溫度和壓力載荷的影響規律,得到了不同型線膨脹機渦旋齒的變形特性,但并未通過有限元方法對膨脹機渦旋盤的變形進行系統分析。本文通過有限元分析方法,對渦旋膨脹機動渦旋盤進行2種熱邊界載荷的施加,得到其溫度場分布及熱變形和熱應力結果,并進行對比分析;對所受氣體力作用下的變形分析采用基于熱力學模型計算和基于流場結果2種方式進行對比研究,最后通過耦合的方法,對膨脹機渦旋盤溫度場與氣體力共同作用下的變形進行了仿真分析。通過溫度場下的對比分析、氣體力作用下的對比分析以及耦合作用下的對比分析,更全面的得到渦旋膨脹機渦旋盤形變分布以及應力分布規律,對渦旋膨脹機的強度設計、可靠性研究以及優化改進提供了一定的理論參考。
渦旋盤幾何模型建立的幾何參數,主要包括:渦旋齒齒高h、齒厚t、基圓半徑R、節距Pt以及渦旋齒圈數N、渦旋盤直徑D等。計算公式為:
t=2αR
(1)
Pt=2πR
(2)
(3)
(4)
其中:α為圓漸開線發生角;Ror為主軸偏心距;φe為最終展角;δ為邊緣間隙。
本研究中渦旋盤幾何參數值如表1所示。

表1 主要幾何參數
渦旋型線采用圓漸開線型線,齒頭部分采用雙圓弧修正法進行修正,依據渦旋型線參數方程建立動渦旋盤的三維模型如圖1,其漸開線方程表達式如式(5)所示:

(5)

圖1 動渦旋盤三維模型
本研究中渦旋盤材料為鋁合金4032,其基本參數如表2所示。
在渦旋膨脹機實際工作過程中,渦旋盤受氣體力作用,主軸頂部緊貼軸承孔頂面,限制了軸向位移,且由于十字滑環,限制了動渦盤發生自轉。所以,其位移約束可簡化為:
1) 端板周圍側壁Z方向的位移為零;
2) 動渦旋盤軸承孔內壁面X、Y方向的位移為零;
3) 動渦旋盤軸承孔頂部Z方向位移為零。

表2 材料屬性
在對渦旋盤網格劃分時,由于膨脹機實際工況下,端板變形很小,渦旋齒變形較大,所以對渦旋齒采用較為精細的正六面體網格劃分,端板采用自動網格劃分,設置全局網格尺寸2.5 mm。劃分網格后的有限元模型共有23 040個單元,67 823個節點,如圖2所示。

圖2 有限元模型
2.1.1基于線性變化的溫度場設置
根據某車用發動機有機朗肯余熱回收系統實際工況,設定進氣溫度為132 ℃,出口溫度為47 ℃,環境溫度為25 ℃。在ANSYS Workbench中建立穩態熱分析模塊,由渦旋盤幾何特性和溫度分布特性分析,可在柱坐標系下通過端板半徑的變化來表達溫度的分布情況。于是,在圓柱坐標系下,可將渦旋盤的溫度簡化為沿端板半徑方向呈線性變化的溫度場分布,其變化規律為:
T=47+85(1-x/79)
(6)
式中:x為柱坐標下渦旋盤端板變化半徑(mm);79為渦旋盤端板最大半徑(mm);85為進口溫度與出口溫度的差值。
2.1.2基于流場溫度的溫度場設置
由該渦旋膨脹機的流場數值模擬結果,以渦旋齒始端為起點,沿渦旋齒壁面連續選取相鄰間隔相同展開角閾值的位置的點。本文以間隔π/6選取,探查其溫度值,可獲得膨脹腔內渦旋盤壁面的部分流場溫度值,通過對溫度值進行平均化處理,得到動渦旋盤壁面的邊界熱載荷。在進行溫度場設置時,利用Workbench環境中的表格數組進行加載,可得到渦旋齒的溫度場分布;此外,因渦旋齒受溫度影響嚴重,渦旋盤端板影響很小,為方便載荷施加,對渦旋盤端面仍施加分布規律為:T=47+85(1-x/79)的線性溫度,亦能體現端面的分布情況,對整個渦旋盤的溫度分布影響很小,能有效的模擬渦旋盤的溫度場分布
本研究中選取吸氣結束時,對動渦旋盤在2種溫度邊界下進行熱應力和熱變形分析,分別得到線性溫度場下的結果云圖3和基于流場溫度下的結果云圖4。由2種溫度邊界下的分析云圖可知,高段溫度均出現在渦旋齒齒頭部分,且由中心向外圍逐漸降低,但基于流場下的溫度場分布更符合實際情況。在熱載荷作用下,最大變形發生在渦旋齒齒頭部分,且頂部位置更加突出,并呈現由齒頭向齒尾非線性逐漸減小的趨勢。

圖3 線性溫度場下的熱分析

圖4 流場溫度場下的熱分析
為了更直觀地觀察渦旋齒的變形情況,以渦旋齒齒頭頂端為原點,型線弧長為x坐標,變形量為y坐標,繪制渦旋齒頂部變形曲線和以渦旋齒齒頭與端板接觸始端為起點,型線弧長為x坐標,應力值為y坐標,繪制渦旋齒與端板接觸處的應力分布曲線,分別如圖5和圖6所示,圖中M代表線性溫度場下的應力和變形結果,N代表基于流場溫度下的應力和變形結果。

圖5 渦旋齒頂部變形曲線

圖6 渦旋齒與端板接觸處的應力分布曲線
通過對比分析,在熱變形上,由于分析求解時,在渦旋盤端板上施加了Z方向的約束,兩者最大變形均位于齒頭部分,但后者最大熱變形116.3 μm大于前者111.9 μm。由變形曲線可知,在弧長為0~100 mm范圍內,變形量下降速度較快,而后比較平緩,在齒尾處又有所上升,這是因為渦旋齒前端膨脹腔體積隨主軸轉動變化較大,且由于齒頭部分經修正后齒壁較厚,相對于弧長100 mm以后的渦旋齒,溫度梯度變化較大,變形量下降速度較快,后因熱脹冷縮,齒尾呈現外凸的趨勢,局部變形量有所上升。另外,不難發現后者渦旋齒整體變形分布大于前者,這是由于實際情況下膨脹腔腔內溫度分布是不均勻造成的,文獻[17-21]中對膨脹機內部流場特性展開了數值模擬研究,得到了流場下的溫度分布的復雜性和不均勻性,并且文獻[16]中已經驗證了數值模擬結果相對于性能試驗結果的有效性,所以基于流場結果得到的渦旋膨脹機運行參數更符合其實際工況,溫度場分布更合理,其造成的熱變形更能有效的模擬真實變形情況;對于熱應力,由于在動渦盤端板下部的軸承孔內壁面以及頂部設置了位移約束,所以應力分布表現在軸承孔與端板交界處和渦旋齒與端板接觸位置,通過兩者熱應力的比較可以看出,應力分布趨勢大致相同,最大熱應力集中在渦旋齒齒頭與端板接觸部分,且前者61 MPa小于后者69 MPa,在弧長為200~400 mm,兩者應力分布差異較大。
由渦旋膨脹機工作機理可知,在吸氣結束時,膨脹機內部氣量最大,渦旋盤所受氣體力最大[7],因此本研究中選擇主軸轉角為吸氣結束時的狀態進行分析。此狀態下,由于渦旋圈數小,此時僅存在一對中心吸氣腔、一對壓縮腔和一對排氣腔。
3.1.1基于熱力學模型的氣體力載荷
渦旋膨脹機渦旋盤在運動過程中,較高的氣體溫度和氣體力作用是導致渦旋齒變形的主要因素。而渦旋盤所受的氣體作用力是渦旋齒內外壁面的壓力所致,且由于不同位置壓力分布不同,內外側存在一定壓差,氣體徑向力作用在內外側壁面,如圖7所示為壁面壓力分布示意圖。所以氣體徑向壓力載荷可簡化為相鄰膨脹腔的壓力差,作用于內壁表面。

圖7 壓力分布示意圖
而渦旋盤所受的氣體作用力的大小受膨脹機膨脹腔容積大小的影響,膨脹腔容積大小又隨主軸轉角的變化而變化,呈現周期性、時變性的變化特點。通過文獻可知渦旋膨脹機膨脹腔容積計算公式為:

(7)
式中:Vi為各膨脹腔容積編號;P為渦旋體節距;t為渦旋齒厚;h為渦旋齒高。
假設膨脹過程是絕熱進行的,那么通過氣體絕熱過程的P-V關系式可以求出第i個膨脹腔內的氣體壓力,P-V關系表達式如下:
(8)
式中:Ps為吸氣壓力;k為氣體等熵指數,取k為4/3;VS為初始容積;Vi為第i個膨脹腔的容積。
通過上式氣體壓力計算公式可得中心進氣腔理論計算值為1.1 MPa,膨脹腔理論計算值為0.584 MPa,排氣腔腔理論計算值為0.3 MPa。
3.1.2基于流場的氣體力載荷
測定流場數值模擬下各膨脹腔的壓力平均值,形成如表3所示的壓力載荷對照表,再在ANSYS中對渦旋盤施加Fluent仿真結果下的各腔壓力值,將其自動轉換到對應的各節點上。

表3 壓力載荷結果對照表 MPa
對動渦旋盤分別施加基于熱力學模型計算得到氣體力載荷和基于流場Fluent仿真結果得到氣體力載荷,應力和變形結果云圖分別如圖8和圖9所示。拾取與溫度場相同的頂部形變路徑和底部應力分布路徑,得到氣體力作用下的渦旋齒頂部變形曲線如圖10所示和底部應力分布曲線如圖11所示,圖中P代表基于熱力學模型下的應力和變形結果,Q代表基于流場氣體力下的應力和變形結果。

圖8 熱力學模型下的氣體力應力變形分析

圖9 Fluent仿真結果下的氣體力應力變形分析

圖10 氣體力作用下的渦旋齒頂部變形曲線

圖11 氣體力作用下的渦旋齒與端板接觸處的應力分布曲線
通過對比分析發現,2種氣體力作用下的最大變形均在渦旋齒齒頭處,而且變形量降速很快,此后有所減緩,且都有沿外圈減小的趨勢;與溫度場作用下的變形曲線對比可知,氣體力作用下的變形遠小于熱載荷作用下的變形;此外,云圖上的B點和變形曲線末端的變化都展現了變形量有所增大的趨勢,這是由于渦旋齒Z軸方向被約束,相對于齒高來說,其齒厚較薄,在載荷作用下產生了類似懸臂梁的效果,應力產生集中,因此形變也相應較大;兩者整體的應力和形變分布基本一致,驗證了計算的正確性;兩者膨脹腔A處的變形均高于兩側,對應了曲線中弧長為280~320 mm中的一點,且在這個范圍內,基于流場氣體力下的渦旋盤形變值高于基于熱力學模型下氣體力造成的形變值,整個形變曲線呈現不同幅度的波動,分析可知是由于齒壁內外側存在壓力差造成的,且相鄰膨脹腔壓差值越大,其波動幅度越大,通過表3可知Fluent仿真結果下的相鄰膨脹腔的壓力差高于熱力學模型計算下的相鄰膨脹腔壓力差,故基于流場結果下的點A處的形變值高于熱力學模型計算下的點A處的形變值。另外,基于流場下的點A處的形變值和點B處的應力值均略高于基于熱力學模型下的A處形變值和B處應力值,一是因為氣體力分布不均,二是因為膨脹機實際運行下,情況更復雜,會受到環境、設備、摩擦等因素的影響從而有所變化。
由于基于熱力學模型下和基于流場下的氣體力作用下的應力變形情況基本相同,故選取其中之一分別與2種形式下的溫度場進行耦合分析。圖12為線性溫度場與流場下氣體力的耦合分析結果;圖13為流場下的溫度場和流場下的氣體力的耦合分析結果;圖14為耦合作用下的渦旋齒頂部變形曲線;圖15為耦合作用下的渦旋齒與端板接觸處的應力分布曲線,圖中S代表線性溫度場和流場下氣體力耦合的結果,R代表流場下的溫度場和流場下氣體力的耦合結果。

圖12 線性溫度場與流場下氣體力的耦合結果

圖13 流場下的溫度場和氣體力的耦合結果

圖14 耦合作用下的渦旋齒頂部變形曲線

圖15 耦合作用下的渦旋齒與端板接觸處的應力分布
通過耦合分析可知,耦合作用下的渦旋盤渦旋齒變形更嚴重,最大變形量分別為125 μm和127 μm,大于任何一個載荷單獨施加下的形變量,且不是單獨施加結果下的簡單線性疊加,其中熱載荷影響較大;耦合作用下的變形趨勢與氣體力、溫度載荷單獨作用施加下的變形趨勢相似,都呈現向外圈逐漸減小的趨勢,最大變形發生在渦旋齒中心齒頭部分,齒頂位置變形更突出;渦旋齒最大應力集中在渦旋齒與端板接觸的位置;耦合過后的變形曲線和應力分布曲線較單獨施加作用下的結果更穩定,波動幅度較小,說明溫度場與氣體力之間在共同作用時相互之間存在一定影響;線性溫度場和基于流場下的氣體力耦合后,渦旋盤變形量低于流場下的溫度場和流場下的氣體力耦合后的變形量,而應力值分布卻略大于后者,這是因為經流場分析得到的溫度場和氣體壓力更符合渦旋膨脹機實際工況,分布更合理。
1) 在單獨施加溫度和單獨施加氣體力下,溫度載荷引起的應力和變形遠大于氣體力引起的應力和變形,所以可以得出溫度載荷是導致渦旋齒變形的主要原因,降低渦旋膨脹機的工作溫度能有效減少變形量和熱應力。
2) 基于線性分布的溫度熱載荷和基于流場得到的熱載荷在某些特殊轉動位置下造成的熱變形和熱應力相差不大,可作為一種備選施加方案,但流場下的熱載荷分布更合理,更符合實際。
3) 理論計算得到的氣體力載荷和基于流場得到的氣體力載荷所得的變形效果很接近,可互相作為一種驗證方式。
4) 由于各工作腔內氣體壓力不同,相鄰膨脹腔存在壓差,壓差越大,波動越大,形變量越大,在與吸氣腔相連的第一膨脹腔形變最明顯。
5) 各載荷單獨施加作用下和耦合作用下的最大變形均位于渦旋齒齒頭頂部,所受應力集中于渦旋齒和渦旋端板連接處以及軸承孔與端板的交界處,應力和變形分布由渦旋齒中心向外圍逐漸減小。且耦合作用下的變形值不是單獨作用下變形值的線性疊加,耦合下的值大于單獨施加得到的值。