李雪松，程秀生，苗麗穎

(吉林大學(xué) 汽車工程學(xué)院，吉林 長(zhǎng)春，130022)

液力緩速器作為一種輔助制動(dòng)裝置，依靠工作輪內(nèi)液流的作用將車輛的動(dòng)能轉(zhuǎn)化為液體的熱能，再通過冷卻器散熱的方式實(shí)現(xiàn)車輛制動(dòng)[1]。當(dāng)液力緩速器工作在部分充液工況下，其內(nèi)部是復(fù)雜的黏性、不可壓縮、非穩(wěn)態(tài)氣-液兩相流動(dòng)，流場(chǎng)測(cè)試試驗(yàn)較復(fù)雜，且設(shè)備昂貴。隨著計(jì)算機(jī)技術(shù)和計(jì)算流體動(dòng)力學(xué)的迅速發(fā)展，對(duì)液力緩速器內(nèi)流場(chǎng)進(jìn)行數(shù)值模擬的方法比較經(jīng)濟(jì)可行。為全面揭示液力緩速器流場(chǎng)的內(nèi)部特性以及內(nèi)外特性之間的變化規(guī)律，目前國(guó)內(nèi)學(xué)者常用雷諾平均法(RNS)進(jìn)行流場(chǎng)的數(shù)值計(jì)算[2-4]，其弊端是忽略了湍流的若干細(xì)節(jié)，所以，對(duì)封閉的雷諾時(shí)均方程進(jìn)行求解并不能有效地模擬流動(dòng)中分離、渦旋以及擴(kuò)散等現(xiàn)象[5]。本文作者借助CFD軟件平臺(tái)，采用更為準(zhǔn)確的大渦模擬法對(duì)液力緩速器部分充液流場(chǎng)進(jìn)行數(shù)值模擬，獲取液力緩速器內(nèi)部流場(chǎng)速度和壓力分布特性，以及不同充液率下內(nèi)流場(chǎng)結(jié)構(gòu)的變化和兩相體積分布情況，并基于流場(chǎng)數(shù)值解計(jì)算液力緩速器的外特性。

1 控制方程及計(jì)算模型

大渦模擬(LES)是介于直接數(shù)值模擬(DNS)與雷諾平均法(RANS)之間的一種湍流數(shù)值模擬方法,是目前CFD研究的熱點(diǎn)之一[6-8]，其基本假設(shè)是：動(dòng)量、能量、質(zhì)量方程及其他標(biāo)量主要由大渦輸運(yùn)；流動(dòng)的幾何和邊界條件決定了大渦的特性，而流動(dòng)特性主要在大渦中體現(xiàn)；小尺度渦旋受幾何合邊界條件影響較小，并且各向同性，直接求解大渦，而小尺度渦對(duì)大尺度渦運(yùn)動(dòng)的影響則通過一定的模型在針對(duì)大尺度渦的瞬時(shí)N-S方程體現(xiàn)出來(lái)[9]。

LES的控制方程是對(duì)N-S方程在波數(shù)空間或者物理空間進(jìn)行過濾得到的，通過去掉比過濾寬度或者給定物理寬度小的渦旋，得到其控制方程：

為使控制方程組封閉，必須建立亞格子尺度模型(SGS模型)，目前應(yīng)用最廣泛的漩渦黏性模型方程為：

式中：δij是“Kronecker delta”符號(hào)(當(dāng)i=j時(shí)，δij=1；當(dāng) i≠j時(shí)為亞格子尺度的湍動(dòng)黏度，根據(jù)Smagorisnsky-Lily模型定義，其中LS為網(wǎng)格濾波寬度。

式中：k為 VonKarman常數(shù)，k＝0.418 7；CS為Samagorin常數(shù)，CS=0.1；d為到最近的壁面的距離；V為計(jì)算單元體積。

根據(jù)液力緩速器樣機(jī)建立的三維模型以及提取的計(jì)算區(qū)域網(wǎng)格模型分別如圖1和圖2所示。定子和轉(zhuǎn)子循環(huán)圓有效直徑分別為293 mm和296 mm，其葉片數(shù)量分別為34和36，均為40°前傾葉片?？紤]到液力緩速器定子內(nèi)有進(jìn)出油口，定子葉片厚度不同，工作流道幾何結(jié)構(gòu)不是周期對(duì)稱，因此，選取整個(gè)流道空間作為計(jì)算模型，網(wǎng)格單元總數(shù)為794 640個(gè)。

圖1 液力緩速器CAD模型Fig.1 Hydraulic retarder CAD model

圖2 液力緩速器計(jì)算網(wǎng)格模型Fig.2 Computational mesh models of hydrodynamic retarder

2 數(shù)值計(jì)算方法

2.1 控制方程的分離式求解法

流場(chǎng)數(shù)值計(jì)算的本質(zhì)就是求解離散后的控制方程組，其基本過程是在空間上將計(jì)算域離散成許多小的體積單元，在每個(gè)體積單元上對(duì)離散后的控制方程組進(jìn)行求解。求解方法可分為分離式求解法和耦合式求解法。分離式求解法通常用于不可壓縮流體數(shù)值計(jì)算，其求解時(shí)同時(shí)考慮所有控制單元，順序、逐個(gè)地求解各變量代數(shù)方程組。目前工程上使用最廣泛的壓力修正法，其單個(gè)時(shí)間步的求解過程如圖3所示。本文運(yùn)用大渦模擬法對(duì)液力緩速器內(nèi)不可壓縮的三維瞬態(tài)流場(chǎng)進(jìn)行模擬，也必須采用分離式求解法求解其控制方程組[10-11]。

圖3 分離求解法計(jì)算流程圖Fig.3 Flow chart for segregated method

2.2 可動(dòng)區(qū)域耦合計(jì)算的滑動(dòng)網(wǎng)格法

液力緩速器的瞬態(tài)特性主要體現(xiàn)在其工作腔內(nèi)液體的流動(dòng)狀態(tài)為湍流非穩(wěn)態(tài)流動(dòng)，當(dāng)轉(zhuǎn)子轉(zhuǎn)速一定、工作腔內(nèi)充液率保持不變時(shí)，轉(zhuǎn)子和定子之間相對(duì)位置的不斷變化引起工作流道內(nèi)的流場(chǎng)特性呈現(xiàn)周期性的變化，當(dāng)緩速器工作腔內(nèi)的充液率及轉(zhuǎn)子轉(zhuǎn)速不斷變化時(shí)，工作流道內(nèi)呈現(xiàn)復(fù)雜的非穩(wěn)態(tài)流動(dòng)形態(tài)。

為了對(duì)轉(zhuǎn)子和定子工作流道進(jìn)行統(tǒng)一計(jì)算，對(duì)轉(zhuǎn)子和定子的接合面采用了滑動(dòng)網(wǎng)格理論來(lái)傳遞兩個(gè)區(qū)域之間的參數(shù)。滑動(dòng)網(wǎng)格法是模擬多移動(dòng)參考系流場(chǎng)的最精確的瞬態(tài)計(jì)算方法，其基本原理為：將計(jì)算區(qū)域中相鄰子單元的分界面互相耦合形成“網(wǎng)格分界面”，數(shù)值計(jì)算中單元區(qū)域在離散步驟中沿著網(wǎng)格分界面相互滑動(dòng)(旋轉(zhuǎn)或平移)，同時(shí)分界面也隨時(shí)間發(fā)生變化(如圖 4所示)[12]。滑動(dòng)網(wǎng)格模擬的瞬態(tài)問題大部分是時(shí)間周期性的，即計(jì)算區(qū)域的速度是周期復(fù)現(xiàn)的。

設(shè)T為瞬態(tài)計(jì)算的周期，在計(jì)算區(qū)域的一些流動(dòng)特性函數(shù)Φ為

轉(zhuǎn)子在一個(gè)周期T內(nèi)運(yùn)動(dòng)2π，又其轉(zhuǎn)速為1 200 r/min，則計(jì)算周期T可用下式計(jì)算：

圖4 網(wǎng)格滑移前后葉輪之間相對(duì)位置Fig.4 Relative position of impellers for sliding mesh theory

單位時(shí)間步長(zhǎng)以不大于單位網(wǎng)格單元大小為選取原則，文中選取的單位時(shí)間步t=0.12 ms。

液力緩速器工作在部分充液狀態(tài)時(shí)，氣相與液相之間相互混摻，將混合模型與歐拉模型交替運(yùn)用在其兩相流數(shù)值計(jì)算中，在多相流模型中首先采用混合模型獲得初始解，并以此作為起點(diǎn)，啟動(dòng)歐拉模型繼續(xù)計(jì)算。在制動(dòng)特性研究中往往關(guān)心液力緩速器內(nèi)外特性的關(guān)系，即轉(zhuǎn)子轉(zhuǎn)速、充液率以及制動(dòng)轉(zhuǎn)矩的變化規(guī)律，而這主要是由液相起支配作用，因此設(shè)定液相為主相，氣相為附加相。

假設(shè)液力緩速器內(nèi)的循環(huán)液體作等溫流動(dòng)且沒有泄露，認(rèn)為其在流動(dòng)過程中為不可壓縮黏性流體，且流體與葉輪間的流固耦合作用不致引起流道的變形。流道內(nèi)壁與葉片表面近壁處速度場(chǎng)計(jì)算采用速度無(wú)滑移邊界條件，進(jìn)出油口采用壓力邊界條件。

3 數(shù)值計(jì)算結(jié)果分析

為了準(zhǔn)確掌握液力緩速器內(nèi)部分充液時(shí)的流場(chǎng)分布特性，對(duì)轉(zhuǎn)子轉(zhuǎn)速為1 200 r/min，充液率f分別為25%和75%的數(shù)值計(jì)算結(jié)果進(jìn)行對(duì)比分析。

3.1 節(jié)面流場(chǎng)特性分析

液力緩速器內(nèi)不同充液率下葉片壓力面的靜壓力分布如圖5所示。較高的充液率能導(dǎo)致較大的靜壓力，因?yàn)楫?dāng)轉(zhuǎn)子轉(zhuǎn)速相同時(shí)，充液率越大，流體密度也越大，則引起的靜壓力越大。將不同流動(dòng)單元中壓力面的靜壓分布(圖 5中(Ⅰ)～(Ⅳ))對(duì)比分析發(fā)現(xiàn)：其分布趨勢(shì)相同，即從循環(huán)圓中心到外環(huán)，靜壓力遞增，入口處出現(xiàn)局部高壓區(qū)，主要是由于定-轉(zhuǎn)子之間過大的轉(zhuǎn)速差導(dǎo)致此處所受液流沖擊最大所致；在循環(huán)圓中心，渦旋和二次流現(xiàn)象頻繁，因此存在環(huán)形低壓區(qū)。充液率較高時(shí)，葉片內(nèi)的壓差較大，液力損失也大。

不同充液率下吸力面的靜壓力分布如圖6所示。由于入油口分布在不同流動(dòng)單元內(nèi)葉片吸力面上，因此不同流動(dòng)單元內(nèi)吸力面靜壓分布差別較大(圖 6中(Ⅰ)～(Ⅳ))。轉(zhuǎn)子葉片中，液流離心力作用使得葉片根部出現(xiàn)高壓區(qū)；定子吸力面內(nèi)沒有大面積高壓區(qū)，低壓區(qū)則是由于沖向壓力面的液流在吸力面產(chǎn)生了尾流造成的。由圖6中(Ⅲ)和(Ⅳ)可見：定子葉片的進(jìn)油口附近壓力較低，正是為了保證液流以高速進(jìn)入工作腔。

3.2 弦面流場(chǎng)特性分析

圖7所示為不同弦面相對(duì)速度矢量分布。中間弦面的渦旋及二次流動(dòng)現(xiàn)象幾乎均布在各流動(dòng)單元內(nèi)兩葉輪的交接處，這主要是由于壓力面和吸力面之間存在較大壓差。此外，該弦面內(nèi)流體質(zhì)點(diǎn)間的相對(duì)劇烈運(yùn)動(dòng)也是導(dǎo)致渦旋形成的重要因素，充液率越高，渦旋在各流動(dòng)單元內(nèi)分布越均勻[13-14]。如圖中A處所示，充液率較低時(shí)，在進(jìn)油口較集中的部分流動(dòng)單元內(nèi)出現(xiàn)了逆流以及小渦旋。

圖8所示為相對(duì)應(yīng)的弦面靜壓力分布圖。從整體來(lái)看，不同充液率時(shí)壓力分布趨勢(shì)相同。充液率為75%時(shí)，整個(gè)弦面內(nèi)的壓力差以及吸力面與壓力面之間的壓差較25%充液時(shí)大。壓力差是使液流偏離循環(huán)流動(dòng)方向形成二次流的重要因素[15]，與圖7對(duì)比可發(fā)現(xiàn)逆流及二次流較明顯的部分流動(dòng)單元內(nèi)，其靜壓值也偏低。

圖5 壓力面靜壓分布Fig.5 Static pressure distribution of pressure surface

圖6 吸力面靜壓分布Fig.6 Static pressure distribution of suction surface

圖7 弦面相對(duì)速度矢量分布Fig.7 Velocity vector distribution of chord surface

圖8 弦面靜壓分布Fig.8 Static pressure distribution of chord surface

3.3 液相體積分布

圖9所示為不同充液率時(shí)的液相體積分布。由于液體在工作腔內(nèi)的循環(huán)流動(dòng)，使液相主要集中在流道外環(huán)，而氣相分布在內(nèi)環(huán)。充液率為75%時(shí)，工作流道的大部分區(qū)域被液相占據(jù)，而在循環(huán)圓中心區(qū)域形成氣環(huán)；充液率為25%時(shí)，氣相占據(jù)了工作區(qū)域的大部分空間。

4 性能預(yù)測(cè)與試驗(yàn)對(duì)比

通過對(duì)液力緩速器內(nèi)流場(chǎng)進(jìn)行瞬態(tài)數(shù)值計(jì)算，可以得到轉(zhuǎn)子葉片表面的速度和壓力梯度，據(jù)此獲得轉(zhuǎn)子葉片表面各單元的壓力和剪切力，再通過對(duì)旋轉(zhuǎn)軸求矩可得到該單元對(duì)旋轉(zhuǎn)軸的作用力矩。對(duì)所有轉(zhuǎn)子葉片表面單元進(jìn)行積分求和即可得到流體作用在轉(zhuǎn)子上的總力矩?；瑒?dòng)網(wǎng)格模擬的瞬態(tài)問題與時(shí)間t有關(guān)，經(jīng)過若干周期，轉(zhuǎn)矩T已經(jīng)隨時(shí)間進(jìn)行周期性變化，因此可以進(jìn)行數(shù)據(jù)采集。

表1和表2所示分別為數(shù)值計(jì)算結(jié)果以及在汽車動(dòng)態(tài)模擬國(guó)家實(shí)驗(yàn)室液力傳動(dòng)試驗(yàn)臺(tái)獲得的試驗(yàn)數(shù)據(jù)。由表1，2和圖10可知：基于三維瞬態(tài)流場(chǎng)數(shù)值解的轉(zhuǎn)矩計(jì)算結(jié)果與試驗(yàn)測(cè)試結(jié)果吻合較好，最大相對(duì)誤差在8%以內(nèi)，說明采用的數(shù)值計(jì)算方法可以較準(zhǔn)確地預(yù)測(cè)液力緩速器的性能。

圖9 液相體積分布Fig.9 Volume distribution of liquid-phase

表1 數(shù)值計(jì)算結(jié)果Table 1 Numerical simulation results

表2 試驗(yàn)數(shù)據(jù)結(jié)果Table 2 Experimental data results

圖10 數(shù)值計(jì)算結(jié)果與試驗(yàn)結(jié)果對(duì)比Fig.10 Comparison of numerical results with experiments

5 結(jié)論

(1) 通過定性對(duì)比發(fā)現(xiàn)，采用大渦模擬法，并利用滑動(dòng)網(wǎng)格理論處理定-轉(zhuǎn)子之間的流動(dòng)參數(shù)傳遞，可以較準(zhǔn)確地模擬液力緩速器內(nèi)氣-液兩相流場(chǎng)的真實(shí)流動(dòng)狀態(tài)，并預(yù)測(cè)其性能。

(2) 外特性數(shù)值計(jì)算結(jié)果與試驗(yàn)結(jié)果吻合良好，因此可以用軟件計(jì)算來(lái)代替部分試驗(yàn)工作，其結(jié)果可用于液力緩速器的設(shè)計(jì)及其結(jié)構(gòu)優(yōu)化。

[1]何仁. 汽車輔助制動(dòng)裝置[M]. 北京: 化學(xué)工業(yè)出版社, 2005:19-21.HE Ren. Automotive auxiliary braking system[M]. Beijing:Chemical Industry Press, 2005: 19-21.

[2]何延?xùn)|, 馬文星, 劉春寶. 液力偶合器部分充液流場(chǎng)數(shù)值模擬與特性計(jì)算[J]. 農(nóng)業(yè)機(jī)械學(xué)報(bào), 2009, 40(5): 24-28.HE Yan-dong, MA Wen-xing, LIU Chun-bao. Numerical simulation and characteristic calculation of hydrodynamic coupling[J]. Transaction of CSAM, 2009, 40(5): 24-28.

[3]嚴(yán)軍, 何仁, 魯明. 液力緩速器變?nèi)~片數(shù)的三維數(shù)值模擬[J].江蘇大學(xué)學(xué)報(bào): 自然科學(xué)版, 2009, 30(1): 27-31.YAN Jun, HE Ren, LU Ming. Numerical simulation of hydraulic retarder with different blade number[J]. Journal of Jiangsu University: Natural Science Edition, 2009, 30(1): 27-31.

[4]張德生, 趙繼云, 劉立寶, 等. 基于 CFD 的桃形腔偶合器流場(chǎng)分析及結(jié)構(gòu)優(yōu)化[J]. 中國(guó)礦業(yè)大學(xué)學(xué)報(bào), 2010, 39(5):687-692.ZHANG De-sheng, ZHAO Ji-yun, LIU Li-bao, et al. Flow field analysis and structure optimization of peach shaped chamber hydrodynamic coupling based on CFD[J]. Journal o f China University of Mining & Technology, 2010, 39(5): 687-692.

[5]才委, 馬文星, 褚亞旭, 等. 液力變矩器導(dǎo)輪流場(chǎng)數(shù)值模擬與試驗(yàn)[J]. 農(nóng)業(yè)機(jī)械學(xué)報(bào), 2007, 38(8): 11-14.CAI Wei, MA Wen-xing, CHU Ya-xu, et al. Numerical simulation and experimental research on flow field in the stator of a torque converter[J]. Transactions of the Chinese Society for Agricultural Machinery, 2007, 38(8): 11-14.

[6]Brandt T T. Usability of explicit filtering in large eddy simulation with a low-order numerical scheme and different subgrid-scale models[J]. International Journal for Numerical Methods in Fluids, 2008, 57(7): 905-928.

[7]ZHANG Yan-hong, YANG Chao, MAO Zai-sha. Large eddy simulation of the gas-liquid flow in a stirred tank[J]. AICHE Journal, 2008, 54(8): 1963-1974.

[8]Grigoriadis D G E, Bartzis J B, Goulas A. Efficient treatment of complex geometries for large eddy simulations of turbulent flows[J]. Computers and Fluid, 2004, 33(2): 201-222.

[9]王福軍. 計(jì)算流體動(dòng)力學(xué)分析—CFD軟件原理與應(yīng)用[M]. 北京: 清華大學(xué)出版社, 2004: 139-142.WANG Fu-jun. Computational fluid dynamics analysis:Principle and application of CFD software[M]. Beijing: Tsinghua University Press, 2004: 139-142.

[10]Felten F, Fautrelle Y, Du Terrail Y, et al. Numerical modeling of electrognetically-riven turbulent flows using LES methods[J].Applied Mathematical Modelling, 2004, 28(1): 15-27.

[11]Feiz A A, Ould-Rouis M, Lauriat G. Large eddy simulation of turbulent flow in a rotating pipe[J]. International Journal of Heat and Fluid Flow, 2003, 243(3): 412-420.

[12]FLUENT Inc. FLUENT user’s guide[M]. Fluent Inc, 2003:239-244.

[13]Mary I, Sagaut P. Large eddy simulation of flow around an airfoil near stall[J]. AIAA Journal, 2002, 40(6): 1139-1145.

[14]Julian R E, Smolarkiewicz K. Eddy resolving simulations of turbulent solar convection[J]. International Journal for Numerical Method in Fluids, 2002, 39(9): 855-864.

[15]Bai L, Fieblg M, Mitra N K. Numerical analysis of turbulent flow in fluid couplings[J]. ASME Journal of Fluids Engineering,1997, 119(3): 569-576.