刁 璇,楊海青
(南京航空航天大學(xué)能源與動力學(xué)院,南京 210016)
對二沖程發(fā)動機而言,進氣方式可分為對稱進氣和非對稱進氣兩種。對稱進氣方式雖然結(jié)構(gòu)簡單、可靠,但存在高速和低速性能矛盾的缺點。由簧片閥控制的非對稱進氣方式在中低速不會產(chǎn)生反噴,穩(wěn)定性和經(jīng)濟性較好[1]。
本文的研究對象是二沖程發(fā)動機進氣系統(tǒng)中關(guān)鍵節(jié)流部件簧片閥。國內(nèi)外很多關(guān)于簧片閥的研究主要側(cè)重于以實驗方式測量在發(fā)動機工作過程中簧片閥升程隨曲軸轉(zhuǎn)角的變化關(guān)系,以及采用數(shù)學(xué)理論推導(dǎo)的方式分析閥體內(nèi)的氣體運動以及簧片的振動特性[2-5]。前者必須在簧片閥結(jié)構(gòu)生產(chǎn)加工出來后才能進行,若以此種方式進行簧片閥的預(yù)研工作,則成本高且效率低;后者可作為研究簧片閥內(nèi)部流體與簧片運動機理的手段,但難以直接運用于工程實踐中。
采用數(shù)值模擬的方法可以有效地預(yù)測簧片的運動,極大地縮短設(shè)計開發(fā)周期,并能獲得一些采用實驗研究法所無法得到的參數(shù),如簧片閥體內(nèi)的流場與壓力分布等[6-7]。在很多場合下,數(shù)值模擬的研究方法具有無法取代的優(yōu)勢。因此,研究適合簧片閥的仿真方法很有必要。
對于二沖程發(fā)動機中簧片閥的仿真工作,國內(nèi)外相關(guān)的研究并不是很多。國內(nèi)有采用一維仿真方法計算簧片閥的升程和流量的研究;國外有學(xué)者進行過一維數(shù)值模擬、二維計算流體力學(xué)(CFD)的研究,并在此二維模型的基礎(chǔ)上進行了流固耦合(FSI)研究[8],以及采用Fluent軟件中用戶自定義功能(UDF)控制簧片運動的軌跡進行三維 CFD 計算研究[9]。
本文采用基于Ansys Workbench的雙向流固耦合仿真研究方法。它不同于借助第三方軟件(如MPCCi)進行流體域與固體域數(shù)據(jù)交換,所有流體域與固體域的計算及數(shù)據(jù)交換都在Workbench平臺中完成。該方法也不同于前文中提到的采用Fluent軟件中的UDF控制簧片運動軌跡實現(xiàn)流體與固體的“同時計算”。每個瞬間作用在簧片上每個位置的流體力都由軟件計算流體運動得到,而簧片的變形對流體域的影響將通過Fluent軟件中動網(wǎng)格功能實現(xiàn),在每個時間步長點上進行流體域與固體域的耦合計算。只要時間步長取得足夠小就能無限接近每個瞬時的耦合計算,即能無限接近實際的簧片閥內(nèi)運動。
簧片閥安裝在二沖程發(fā)動機進氣系統(tǒng)中,位于曲軸箱前面,控制發(fā)動機的整個進氣過程。在閥類型上屬于單向閥,可防止在發(fā)動機進氣過程中流過閥體的氣體出現(xiàn)反噴現(xiàn)象,有效確保了進入曲軸箱中的進氣量。閥體主要由閥殼、閥芯、簧片、限位板組成。簧片有6塊,每3片一組,分別位于上下兩側(cè),通過螺釘將簧片與限位板緊固在閥芯上?;善y模型見圖1,物理參數(shù)見表1。
圖1 簧片閥模型
表1 簧片閥物理參數(shù)
簧片閥位于進氣管與曲軸箱之間,在發(fā)動機進氣過程中,流過簧片閥氣體的運動過程十分復(fù)雜。這是由于在進氣過程中進氣管中的氣體壓力是不斷波動的,且曲軸箱受活塞上下運動的影響不斷向缸內(nèi)掃氣的過程導(dǎo)致了曲軸箱中氣體壓力也在不斷地變化。在進氣過程中,氣體存在熱交換的過程,使得氣體的溫度也會發(fā)生變化。因此,需要對實際流動過程做必要的簡化及假定[10-11],具體如下:
1)只考慮空氣在簧片閥內(nèi)的流動,忽略進氣道噴射型(PFI)發(fā)動機進氣系統(tǒng)中的燃油成分。
2)氣體在簧片閥內(nèi)的流動沒有熱交換,認為是絕熱流動過程,比熱容保持不變。
3)認為進氣壓力是不變的,不考慮進氣管中的壓力波動;考慮到閥內(nèi)的氣體在流通方向的流量及速度遠大于另兩個方向,將流動過程當(dāng)做一元定常等熵過程處理。
4)將在簧片閥內(nèi)流動的氣體視為理想氣體,即氣體嚴(yán)格遵從理想氣體狀態(tài)方程。
圖2為簧片閥內(nèi)氣體流動的情況。進口處氣體的壓力、溫度、密度分別為 P0,T0,ρ0;氣體的流速為V0;簧片閥出口處的氣體壓力、溫度、密度分比為 P1,T1和 ρ1;氣體流速為 V1。
圖2 簧片閥內(nèi)流動情況
根據(jù)理想氣體的一元定常流動能量方程[12-13]:
氣體流過簧片閥的一元等熵流動伯努利方程為
對于理想氣體則有
式(1)~(4)中各物理量分別為:P為壓強(Pa);ρ為密度(kg/m3);H為焓(J/kg);K為等熵指數(shù);Cp為定壓比熱容(J/kg·k)。
由式(3)與(4)得到
將式(2)代入式(5)得到
由氣體連續(xù)性方程得到通過簧片閥的氣體質(zhì)量流量微分方程為
式(7)中:m1為通過簧片閥的氣體質(zhì)量(kg);?1為簧片閥流量系數(shù);A1為簧片閥出口截面積(m3)。
將式(6)代入(7)中得到通過簧片閥氣體質(zhì)量流量微分方程:
在簧片運動過程中,當(dāng)外力的頻率達到簧片的固有頻率時,簧片的振幅會達到最大,此時簧片的應(yīng)力與應(yīng)變會非常大,對限位板和閥芯的撞擊也最為劇烈,十分容易損傷簧片。因此,應(yīng)保證在發(fā)動機工作過程中,氣體推力的頻率應(yīng)遠離簧片的固有頻率。對本文研究的簧片閥所安裝的發(fā)動機來說,設(shè)計工況點為6300 r/min,二沖程發(fā)動機曲軸每轉(zhuǎn)一圈,發(fā)動機完成一次工作循環(huán),則簧片閥需要開啟一次。在發(fā)動機處于設(shè)計工況點時,氣體推力的頻率為105 Hz,應(yīng)保證在設(shè)計工況點氣體的推力頻率小于簧片的固有頻率。
根據(jù)簧片的結(jié)構(gòu)和安裝情況,可將簧片視為等截面的懸臂梁來研究其振動特性。需要說明的是,本研究不考慮簧片的振動阻尼,只研究簧片在主彎曲平面里的彎曲振動,不考慮剪切力的影響[11]。圖3為懸臂梁形狀的簧片模型。
圖3 懸臂梁結(jié)構(gòu)簧片模型
本研究中:Q為截面切力;Z為彎曲位移;M為截面彎矩;ρ為簧片密度;EI為彎曲剛度;A為截面積。由于這里研究簧片的固有特性,故令作用在簧片上的氣體推力為0。
對照模型圖3,根據(jù)達朗貝爾原理建立Z方向上的力平衡方程:
再建立任意微元段力矩平衡方程:
忽略dx的二階小項得到
整理式(9)與(11)得到簧片自由振動方程:
代入懸臂梁結(jié)構(gòu)所對應(yīng)的邊界條件:x=0,Z(0)=Z'(0)=0;x=l,Z″(l)=Z?(l)=0。
求解式(12)可以得到簧片振動的頻率方程:
利用Matlab軟件求解方程(13),發(fā)現(xiàn)該方程有無窮多的解,對應(yīng)了簧片振動的多個固有頻率。取第1階固有頻率與發(fā)動機氣體推力頻率對比,更高階的固有頻率遠大于發(fā)動機中氣體的推力頻率,不需要考慮。計算得出k1l=1.8751。再由k與ω的關(guān)系以及表(1)中簧片的尺寸、物理參數(shù)求出第1階固有頻率ω1/2π約為120 Hz。發(fā)動機在設(shè)計工況點的氣體推力頻率為105 Hz,小于簧片的1階固有頻率,在該轉(zhuǎn)速附近不會發(fā)生共振,保證了簧片閥在最常用工況下的使用安全。但當(dāng)發(fā)動機轉(zhuǎn)速接近7200 r/min時,氣體推力的頻率就會接近簧片的1階固有頻率,應(yīng)避免發(fā)動機長時間工作在該轉(zhuǎn)速附近。
圖4為簧片閥雙向流固耦合計算原理。
圖4 簧片閥雙向流固耦合計算原理
如圖4所示,進行簧片閥雙向流固耦合計算時,首先要分別對流體域以及固體域進行網(wǎng)格劃分、設(shè)置邊界條件、選定材料屬性;然后初始化模型;先在流體域進行CFD計算,當(dāng)計算收斂后,流體軟件將計算得到的流體信息傳遞給固體計算軟件,固體計算軟件根據(jù)流體信息計算出固體域的變形信息;變形后的固體會產(chǎn)生新的流場幾何邊界,流體軟件通過動網(wǎng)格的光順和重構(gòu)功能得到新的流體域網(wǎng)格,并進行下一時間步的流體CFD計算;當(dāng)設(shè)定的所有時間步計算完成后,就可以導(dǎo)出流體域和固體域的計算結(jié)果進行后處理。
在Ansys Workbench14.0版本中建立簧片閥流固耦合仿真模型,如圖5所示。幾何建模采用UG軟件,然后倒入Workbench的幾何模塊中。流體部分采用Fluent軟件計算,固體部分采用Transient Structure模塊計算,流體域與固體域的數(shù)據(jù)交換采用System Coupling模塊完成。
2.2.1 流體域模型建立
建立模型時使用Workbench中自帶的Meshing畫網(wǎng)格工具。由于存在動網(wǎng)格區(qū)域,所以網(wǎng)格類型選擇四面體,網(wǎng)格量為154584。湍流模型為k-ε兩方程模型,近壁區(qū)采用標(biāo)準(zhǔn)壁面模型。采用基于壓力、絕對速度的SIMPLE求解方式。對流差分格式為1階迎風(fēng)格式。動網(wǎng)格方法采用Smoothing和Remeshing,將流體域和固體域需要進行數(shù)據(jù)交換的面設(shè)置為System Coupling類型的動網(wǎng)格區(qū)域。流體域網(wǎng)格模型如圖6所示。
圖5 Workbench中仿真模型
圖6 流體域網(wǎng)格模型
2.2.2 固體域模型建立
固體域為簧片閥的6塊簧片。該型簧片閥的6塊簧片是相互獨立的,在根部沒有互相連接,不需要對簧片做分離的簡化??紤]到流體域是劃分的四面體網(wǎng)格,為保持流體域與固體域在接觸面上網(wǎng)格類型的一致性,在Transient Structure模塊中采用Meshing工具劃分四面體網(wǎng)格,網(wǎng)格量為38946。對固體域需要設(shè)置約束條件,將簧片被螺釘固定在閥芯上的一端設(shè)置為Fixed Support約束類型。其他5個面均與流體接觸,受到流體力的影響,設(shè)置為Fluid Solid interface。約束類型如圖7所示。
圖7 簧片約束類型
2.2.3 流體域與固體域數(shù)據(jù)交換模塊設(shè)置
System Coupling模塊負責(zé)處理流體與固體的數(shù)據(jù)交換,需要在該模塊中設(shè)置耦合計算的時間信息,且要保證與流體域及固體域中設(shè)置的計算時間信息一致。將流體與固體需要進行數(shù)據(jù)交換的面設(shè)置成Data Transfer類型。還需要說明流體與固體哪一個首先進行計算。對于流體驅(qū)動固體運動的問題,應(yīng)將流體設(shè)置為首先進行計算。
保持簧片閥仿真模型的進出口壓力不變,研究在固定壓差下,簧片在氣體壓差力作用下的升程。本文研究了在5組不同的氣體壓差力(以大氣壓為單位)下簧片的尖端升程曲線和簧片閥中氣體流場變化。表2列出了在5組氣體壓差下的簧片閥穩(wěn)定升程值。圖8為各組壓差下簧片升程曲線圖。
表2 不同氣體壓差對應(yīng)的簧片升程值
圖8 不同壓差下簧片升程曲線
從圖8可以看出:當(dāng)簧片閥前后氣體壓差保持不變時,簧片的升程是從0開始,最終穩(wěn)定在一個較小區(qū)間內(nèi)的過程。當(dāng)計算時間足夠長時,簧片的升程會最終穩(wěn)定在區(qū)間內(nèi)的某個值,此時簧片的升程就是理論上在某一壓差下的穩(wěn)定升程值。但考慮到計算機性能有限,計算雙向流固耦合問題需要的時間很長,所以當(dāng)簧片升程波動足夠小時,取波動區(qū)間的中間值作為簧片升程最終的穩(wěn)定值。
分析圖8中升程曲線的變化規(guī)律可以發(fā)現(xiàn):在簧片閥進出口氣體壓差力作用下,在簧片變形初期,簧片的升程會到達一個較大的波峰。此時,簧片的升程值為整個變化過程中的最大值,但簧片無法一直保持在該位置。在該波峰過后,簧片的升程值會迅速下降,下降到一定值之后,簧片的升程值會再次上升,但幅度遠低于首次上升時的值。之后,簧片升程值會不斷的上下波動,試圖尋求一個平衡位置。
出現(xiàn)上述變化規(guī)律的原因可以從使簧片發(fā)生變形的作用力來考慮。作用在簧片上的力有兩方面組成:一是簧片前后空氣的壓差力;二是簧片自身的回彈力?;善纳套兓?guī)律其實正反映了這兩種力的變化過程。當(dāng)氣體壓差力大于簧片的回彈力時,簧片的升程會增大;當(dāng)回彈力大于氣體壓差力時,升程會減小。直到氣體壓差力與簧片的變形程度處于相對平衡時,簧片的升程才會穩(wěn)定在某一位置。圖9可見簧片中最小應(yīng)力變化,以此來說明簧片的回彈力變化過程。
圖9 不同壓差下簧片尖端應(yīng)力值
對比圖8與圖9可以發(fā)現(xiàn):簧片尖端應(yīng)力變化過程與簧片升程曲線的總體趨勢保持一致,符合實際的變化規(guī)律。由于簧片的變形過程是非線性的,簧片的升程越大,則簧片內(nèi)產(chǎn)生的應(yīng)力也越大,產(chǎn)生相同的升程變化量所需的氣體壓差力也越大。從表2就可以看出這一點。圖10為0.1 atm氣體壓差下簧片閥截面總壓云圖,設(shè)置的操作工況為1 atm,反映了氣體通過壓差逐步打開簧片的過程。
圖10 0.1atm氣體壓差下簧片閥內(nèi)總壓云圖
從圖10也可以看出:在簧片的運動過程中,雖然簧片閥進出口面上壓力保持不變,但簧片周圍的氣體壓力是在不斷變化的,它與簧片的回彈力一起控制著簧片的變形程度。圖11為在不同壓差下通過簧片閥的氣體體積流量。
圖11 不同壓差下通過簧片閥氣體體積流量
非定常邊界條件需要控制進出口面上的壓力隨時間發(fā)生變化,利用Fluent軟件中用戶自定義功能(UDF)來實現(xiàn),入口面上的壓力變化規(guī)律如圖12所示,出口面壓力為99 Pa。
圖12 入口面壓力變化圖
進出口面上的氣體壓差最小為1 Pa,最大為6 Pa,最大壓差出現(xiàn)在0.01 s這一時刻。通過流固耦合計算來驗證簧片的升程曲線能否響應(yīng)邊界壓力的變化規(guī)律。圖13為計算得到的簧片升程曲線。
圖13 非定常邊界條件下簧片升程曲線
圖13表明了簧片升程曲線符合邊界壓力的變化規(guī)律:升程曲線隨入口面壓力增大而增大、減小而減小,并在0.01 s時達到了最大值。因為設(shè)置了最小壓差1 Pa,所以在簧片運動的初期,升程上升得較快,一段時間后升程上升速度放緩,出現(xiàn)了一個波峰,這類似于定常邊界條件下簧片升程曲線上第一個波峰的形成原因;在計算結(jié)束時,簧片的升程比1 Pa固定壓差下的穩(wěn)定升程要大,原因是受到了上游較大氣體壓差的影響。說明簧片升程能準(zhǔn)確地響應(yīng)邊界面上的壓力變化,則可以完善對二沖程發(fā)動機的三維CFD數(shù)值模擬工作。此前對二沖程發(fā)動機的CFD模擬工作大都是以掃氣道入口面作為整個模型的入口邊界,取曲軸箱中氣體的平均壓力作為入口面上的壓力邊界條件。這種方法忽略了整個進氣系統(tǒng)對發(fā)動機進氣過程的影響,與發(fā)動機實際的工作情況存在較大差距。
1)分別對簧片閥內(nèi)氣體流動情況和簧片振動特性做了說明,得到了通過簧片閥的氣體質(zhì)量流量微分方程,并計算得出簧片的固有頻率,分析了簧片的共振情況。
2)對簧片閥分別進行了定常邊界和非定常邊界條件下的流固耦合仿真計算,得到了定常邊界條件下簧片的升程變化曲線、閥體內(nèi)部氣體壓力云圖和通過簧片閥的氣體流量;證明了在非定常邊界條件下,簧片升程能夠響應(yīng)邊界面上壓力的變化,從而可以進一步完善二沖程發(fā)動機整機三維CFD數(shù)值模擬的工作。
本文對簧片運動過程的研究沒有涉及簧片撞擊限位板和閥芯的過程。對于不同邊界條件下,簧片撞擊限位板和閥芯會產(chǎn)生多大的沖擊應(yīng)力及簧片的彈跳高度等問題,需要在今后的工作中做進一步研究。對于二沖程發(fā)動機整機的CFD模擬工作,由于模型的網(wǎng)格量較大,且存在模擬活塞運動的動網(wǎng)格區(qū)域和簧片閥流固耦合區(qū)域,對計算機配置要求較高,故最好采用高性能計算機進行計算。
[1]楊光興,葉茂炎,程善斌.摩托車發(fā)動機原理與設(shè)計[M].武漢:武漢測繪科技大學(xué)出版社,1993.
[2]黃偉.NF50摩托車發(fā)動機簧片閥座組合的試驗研究[J].摩托車技術(shù),1998,3:7 -8.
[3]由麗娜,毛華永,戴永圣.小型二沖程汽油機進氣系統(tǒng)的試驗研究[J].摩托車技術(shù),1997(6):3-8.
[4]Mitianiec W,Bogusz A.Theoretical and Experimental Study of Gas Flow through Reed Valve in a Two-Stroke Engine[J].SAE Paper,1996(1):1802.
[5]錢興華.無升程閥限制器的條狀簧片閥開發(fā)[J].流體工程,1992,20(10):22 -26.
[6]高孝洪.內(nèi)燃機工作過程數(shù)值計算[M].北京:國防工業(yè)出版社,1986.
[7]朱訪君,吳堅.內(nèi)燃機工作過程數(shù)值計算及優(yōu)化[M].北京:國防工業(yè)出版社,1997.
[8]Dave A,Siddiqui A,Probst D,et al.Development of a Reed Valve Model for Engine Simulations for Two-Stroke Engines[J].SAE Paper,2004(1):1455.
[9]Battistoni M,Grimaldi C N,Baudille R,et al.Development of a Model for the Simulation of a Reed Valve Based Secondary Air Injection System for SI Engines[J].SAE Paper,2005(1):0224.
[10]吳丹青,叢敬同.壓縮機簧片閥的數(shù)學(xué)模擬與設(shè)計[M].北京:機械工業(yè)出版社,1993.
[11]鄧志勇.K157發(fā)動機進氣系統(tǒng)單向閥的研究[D].濟南:山東大學(xué),2006.
[12]陳卓如,金朝銘,王洪杰,等.工程流體力學(xué)[M].北京:高等教育出版社,2004.
[13]劉錚,張楊軍.內(nèi)燃機一維非定常流動[M].北京:清華大學(xué)出版社,2007.