曹月昊 ，孫峪鵬 ，尉洋，黃羽鵬，周億莉，姚松

(1.中南大學(xué) 交通運(yùn)輸工程學(xué)院，軌道交通安全教育部重點(diǎn)實(shí)驗(yàn)室，湖南 長(zhǎng)沙，410075；2.中南大學(xué) 軌道交通安全關(guān)鍵技術(shù)國(guó)際合作聯(lián)合實(shí)驗(yàn)室，湖南 長(zhǎng)沙，410075；3.中南大學(xué) 軌道交通列車安全保障技術(shù)國(guó)家地方聯(lián)合工程研究中心，湖南 長(zhǎng)沙，410075)

在城市公共交通體系中，地鐵因具有安全環(huán)保、運(yùn)量大等特點(diǎn)成為了主流的交通出行方式。新形勢(shì)要求地鐵各系統(tǒng)應(yīng)滿足人性化設(shè)計(jì)，乘客舒適性已經(jīng)成為地鐵車輛性能評(píng)估的主要指標(biāo)之一。在地鐵各系統(tǒng)中，地鐵空調(diào)通風(fēng)系統(tǒng)調(diào)節(jié)列車內(nèi)外空氣，且送風(fēng)風(fēng)道作為其重要部件，負(fù)責(zé)引導(dǎo)氣流流向，控制送風(fēng)出口送風(fēng)狀態(tài)，影響著乘客的整體舒適度。地鐵在運(yùn)營(yíng)過(guò)程中常出現(xiàn)送風(fēng)量不均勻、風(fēng)向直吹、冷熱變化過(guò)大等問(wèn)題，為合理組織車廂內(nèi)氣流分布，控制車內(nèi)的微風(fēng)速，許多學(xué)者對(duì)地鐵送風(fēng)風(fēng)道結(jié)構(gòu)進(jìn)行優(yōu)化，以改善風(fēng)道送風(fēng)的均勻性。

王斌[1]通過(guò)流場(chǎng)仿真觀察到靠近空調(diào)機(jī)組的送風(fēng)風(fēng)道出口存在回流現(xiàn)象，基于經(jīng)驗(yàn)在送風(fēng)風(fēng)道進(jìn)口兩側(cè)增加穿孔擋板使得此現(xiàn)象得到緩解。劉家林等[2]采取在送風(fēng)風(fēng)道內(nèi)部增加相應(yīng)隔板、封堵一些送風(fēng)口、加裝導(dǎo)流管等內(nèi)部擾流措施，使得風(fēng)道送風(fēng)更為均勻。尤立偉等[3]利用不均勻系數(shù)對(duì)風(fēng)道內(nèi)的氣流組織進(jìn)行評(píng)價(jià)，并通過(guò)調(diào)整孔板參數(shù)對(duì)非變截面風(fēng)道進(jìn)行結(jié)構(gòu)優(yōu)化。胡滋艷[4]研究了地鐵送風(fēng)風(fēng)道出口導(dǎo)流板的送風(fēng)角度對(duì)溫度場(chǎng)、速度場(chǎng)、流場(chǎng)渦結(jié)構(gòu)及PMV-PPD 評(píng)價(jià)指標(biāo)的影響，通過(guò)參數(shù)優(yōu)化得到了最優(yōu)的送風(fēng)角度。KIM等[5]利用2k因子設(shè)計(jì)方法對(duì)頂風(fēng)風(fēng)道系統(tǒng)的設(shè)計(jì)參數(shù)進(jìn)行有效性評(píng)價(jià)(其中，k為影響因素?cái)?shù)量，2表示影響因素的兩個(gè)方面)，并采用響應(yīng)面法設(shè)計(jì)了具有改善流動(dòng)均勻性的優(yōu)化模型，有效地抑制初始模型中存在的非均勻流分量。夏春晶等[6]在分析風(fēng)道通風(fēng)性能的基礎(chǔ)上對(duì)風(fēng)道主體、送風(fēng)孔縫、送風(fēng)格柵等結(jié)構(gòu)進(jìn)行優(yōu)化，有效地改善了送風(fēng)均勻性。

以上地鐵送風(fēng)風(fēng)道結(jié)構(gòu)優(yōu)化方法的核心是借助CFD 技術(shù)，依據(jù)經(jīng)驗(yàn)在送風(fēng)風(fēng)道中增添構(gòu)型簡(jiǎn)單的導(dǎo)流結(jié)構(gòu)，并采用尺寸優(yōu)化方法[7]尋求此初始構(gòu)型限制下的局部最優(yōu)解。此類方法的優(yōu)化潛力嚴(yán)重依賴于導(dǎo)流結(jié)構(gòu)的初始構(gòu)型，優(yōu)化設(shè)計(jì)的理念局限于原有結(jié)構(gòu)的改進(jìn)，無(wú)法實(shí)現(xiàn)最優(yōu)的效果。本文則將拓?fù)鋬?yōu)化方法應(yīng)用于地鐵車廂送風(fēng)風(fēng)道結(jié)構(gòu)優(yōu)化中，以風(fēng)道出口處氣流均勻性偏差最小化與進(jìn)出口平均壓差最小化為目標(biāo)，從拓?fù)錁?gòu)型尋求當(dāng)前邊界條件下固體與流體區(qū)域的最優(yōu)分布，優(yōu)化設(shè)計(jì)具有良好流場(chǎng)性能的新型導(dǎo)流結(jié)構(gòu)，送風(fēng)風(fēng)道出口處速度均勻性與進(jìn)出口平均壓差都明顯降低。

1 初始送風(fēng)風(fēng)道流場(chǎng)性能分析

1.1 幾何模型

現(xiàn)有地鐵送風(fēng)風(fēng)道分布在車廂頂部，且每節(jié)車廂布置2 個(gè)空調(diào)機(jī)與2 組送風(fēng)風(fēng)道。2 組送風(fēng)風(fēng)道相對(duì)車廂中心呈對(duì)稱分布，分別控制車廂橫向兩側(cè)送風(fēng)狀態(tài)。提取單側(cè)送風(fēng)風(fēng)道作為流場(chǎng)分析的研究對(duì)象，并在簡(jiǎn)化模型的過(guò)程中保留出口導(dǎo)流板、風(fēng)道固定內(nèi)筋等關(guān)鍵結(jié)構(gòu)。圖1所示該地鐵送風(fēng)風(fēng)道的初始結(jié)構(gòu)，風(fēng)道內(nèi)腔高為0.29 m，風(fēng)道進(jìn)口寬為0.85 m，風(fēng)道出口較為狹長(zhǎng)，長(zhǎng)為11.44 m，且部分出口區(qū)域靠近進(jìn)風(fēng)口，部分出口區(qū)域靠近風(fēng)道尾部尖端。

圖1 地鐵車廂送風(fēng)風(fēng)道的布局與結(jié)構(gòu)Fig.1 Layout and structure of air supply duct for subway cars

1.2 物性參數(shù)與邊界條件

此送風(fēng)風(fēng)道標(biāo)準(zhǔn)送風(fēng)量為4 300 m3/h，內(nèi)部空氣密度取1.225 kg/m3，動(dòng)力黏度取1.86×10-5kg/(m·s)。風(fēng)道進(jìn)口設(shè)為速度進(jìn)口邊界條件，進(jìn)口速度根據(jù)標(biāo)準(zhǔn)參數(shù)換算為6.862 8 m/s；風(fēng)道出口設(shè)為壓力出口邊界條件；風(fēng)道壁面設(shè)為無(wú)滑移壁面邊界條件。

1.3 計(jì)算網(wǎng)格與湍流模型

以六面體網(wǎng)格為主對(duì)網(wǎng)格進(jìn)行劃分，網(wǎng)格的基礎(chǔ)邊長(zhǎng)設(shè)為0.02 m，并在風(fēng)道壁面添加3層邊界層，總體網(wǎng)格數(shù)約為80 萬(wàn)個(gè)。具體網(wǎng)格如圖2所示。

圖2 風(fēng)道網(wǎng)格劃分Fig.2 Grid division of air duct

地鐵送風(fēng)風(fēng)道內(nèi)流動(dòng)為湍流流動(dòng)，且馬赫數(shù)低于0.3，流體的密度變化很小，故采用定常、黏性、不可壓縮N-S方程與標(biāo)準(zhǔn)k-ε湍流模型[8]進(jìn)行仿真模擬。

1.4 結(jié)果分析

圖3(a)所示為送風(fēng)風(fēng)道出口表面的速度分布情況。由于風(fēng)道出口構(gòu)型較狹長(zhǎng)，且進(jìn)口處額定流速較快，出口表面處的速度分布產(chǎn)生了明顯的分層，靠近進(jìn)口側(cè)的出口區(qū)域速度較小，而遠(yuǎn)離進(jìn)口側(cè)的出口尖端區(qū)域則產(chǎn)生了較大的速度，部分區(qū)域的速度達(dá)6 m/s。圖3(b)所示為風(fēng)道出口表面法向速度矢量的分布情況?？拷M(jìn)口側(cè)的出口區(qū)域產(chǎn)生了回流現(xiàn)象，而遠(yuǎn)離進(jìn)口側(cè)的出口尖端區(qū)域則出現(xiàn)了法向風(fēng)速過(guò)高情況。

圖3 風(fēng)道出口的速度分布情況Fig.3 Velocity distribution at the outlet of air duct

影響送風(fēng)風(fēng)道流場(chǎng)性能的主要因素為送風(fēng)均勻性與進(jìn)出口壓差，而工程上需在一定進(jìn)出口壓差限制下盡可能提高出口均勻性，以保證客室內(nèi)氣流分布的合理性。送風(fēng)均勻性評(píng)價(jià)指標(biāo)為速度均勻性系數(shù)γm(速度完全均勻時(shí)，γm=1)，進(jìn)出口壓差定義為進(jìn)出口壓強(qiáng)之差Δp，其計(jì)算公式分別為：

其中：n為出口平面網(wǎng)格個(gè)數(shù)；vi為每個(gè)網(wǎng)格處的平均風(fēng)速；vmean為整個(gè)出口處的平均風(fēng)速；pin為進(jìn)口處平均壓強(qiáng)；pout為出口處平均壓強(qiáng)。為評(píng)估送風(fēng)風(fēng)道出口的送風(fēng)均勻性，分別將出口表面速度以及法向速度矢量代入式(1)，以求得送風(fēng)風(fēng)道出口速度均勻性系數(shù)γt與法向速度均勻性系數(shù)γv。

根據(jù)式(1)和(2)可求得初始送風(fēng)風(fēng)道出口速度均勻性系數(shù)γt為0.7912，法向速度均勻性系數(shù)γv為0.585 9，進(jìn)出口壓差Δp為23.879 22 Pa。由于送風(fēng)風(fēng)道出口具有“橫向窄、縱向長(zhǎng)”的特性，出口處送風(fēng)均勻性主要取決于出口處法向速度分布，因此，減少出口近端的回流現(xiàn)象與遠(yuǎn)處尖端的風(fēng)速過(guò)高情況是亟待解決的問(wèn)題。

2 拓?fù)鋬?yōu)化理論

2.1 拓?fù)鋬?yōu)化問(wèn)題設(shè)計(jì)

拓?fù)鋬?yōu)化求解的一般形式[9]為：

其中：γ為設(shè)計(jì)變量γ組成的向量，流體拓?fù)鋬?yōu)化中γ定義為材料的滲透能力(反滲透的固體材料單元γ為0，流體單元γ為1，中間單元?jiǎng)t考慮為不同孔隙率的多孔介質(zhì)材料)；C為優(yōu)化問(wèn)題的目標(biāo)函數(shù)；R為優(yōu)化問(wèn)題的控制方程；U為隨γ變化的狀態(tài)變量；g為優(yōu)化問(wèn)題的約束方程。

2.2 適應(yīng)拓?fù)鋬?yōu)化的修正k-ε模型

固體結(jié)構(gòu)拓?fù)鋬?yōu)化常采用變密度法對(duì)彈性模量進(jìn)行插值[10]，流體拓?fù)鋬?yōu)化則需要在流動(dòng)控制方程中添加達(dá)西插值項(xiàng)[11]，使得γ實(shí)現(xiàn)反滲透固體與流體單元之間的相互轉(zhuǎn)化。拓?fù)鋬?yōu)化應(yīng)用于層流流動(dòng)，只需要對(duì)N-S 方程中的動(dòng)量模型進(jìn)行修正[12-15]；而拓?fù)鋬?yōu)化應(yīng)用于湍流流動(dòng)，在對(duì)標(biāo)準(zhǔn)k-ε模型中動(dòng)量方程進(jìn)行修正的同時(shí)，還需考慮湍流流動(dòng)中拓?fù)渥冃我鸬耐牧髂芰孔兓?，將類似的插值?xiàng)加入k與ε的傳輸方程[16]中，使得固體區(qū)域邊界速度、湍流動(dòng)能以及湍流耗散能均設(shè)定為0。

修正后的不可壓縮流標(biāo)準(zhǔn)k-ε模型為：

在懲罰項(xiàng)中，λ為布林克曼懲罰系數(shù)(將優(yōu)化過(guò)程中產(chǎn)生的多孔介質(zhì)單元推向反滲透固體或流體單元，取1 000)；α(γ)函數(shù)為實(shí)現(xiàn)固體與流體空間變化的插值函數(shù)；q為達(dá)西罰函數(shù)系數(shù)(控制插值函數(shù)的曲率，取0.1)。

2.3 拓?fù)鋬?yōu)化算法流程

為實(shí)現(xiàn)拓?fù)浞植疾粩嘞蚋鼉?yōu)構(gòu)型迭代，且考慮到求解復(fù)雜度問(wèn)題，采用伴隨法[17]對(duì)優(yōu)化問(wèn)題的靈敏度進(jìn)行計(jì)算，選擇移動(dòng)漸近線方法[18]求解優(yōu)化模型，同時(shí)采用亥姆霍茲方程[19]對(duì)材料體積因子進(jìn)行過(guò)濾，并將過(guò)濾的結(jié)果通過(guò)雙曲正切投影[20]進(jìn)行處理以獲得清晰的拓?fù)溥吔?。整個(gè)拓?fù)鋬?yōu)化算法流程如圖4所示。

圖4 拓?fù)鋬?yōu)化流程Fig.4 Flow of turbulent topology optimization

3 送風(fēng)風(fēng)道結(jié)構(gòu)拓?fù)鋬?yōu)化設(shè)計(jì)

在上述流場(chǎng)分析中，風(fēng)道出口處氣流速度差異主要在于風(fēng)道結(jié)構(gòu)的變截面設(shè)計(jì)?？紤]到三維湍流有限元計(jì)算與伴隨靈敏度求解所需時(shí)間較長(zhǎng)，為提高優(yōu)化效率，將優(yōu)化區(qū)域簡(jiǎn)化為風(fēng)道的水平橫截面，保留初始風(fēng)道的變截面特性，采用二維簡(jiǎn)化模型對(duì)風(fēng)道進(jìn)行拓?fù)鋬?yōu)化設(shè)計(jì)。與三維初始風(fēng)道的邊界條件相同，二維風(fēng)道進(jìn)口風(fēng)速為6.862 8 m/s。

在本次優(yōu)化設(shè)計(jì)中，二維流道出口區(qū)域被劃分為50 個(gè)等長(zhǎng)區(qū)域(即n=50)，設(shè)每個(gè)小出口的平均速度為vj，整個(gè)出口平均速度為vmean，進(jìn)口處平均壓強(qiáng)為pin，出口處平均壓強(qiáng)為pout。優(yōu)化的目標(biāo)是盡可能減小風(fēng)道出口均勻性偏差，同時(shí)避免風(fēng)道進(jìn)出口壓差過(guò)大。優(yōu)化中無(wú)約束，目標(biāo)函數(shù)C包括出口均勻性偏差與進(jìn)出口平均壓差兩部分[21]，其中出口均勻性偏差為1-γm，進(jìn)出口平均壓差為Δp。計(jì)算當(dāng)前區(qū)域下1-γm與Δp的初始值，結(jié)合均勻性目標(biāo)在風(fēng)道優(yōu)化中的重要性，且考慮到壓差目標(biāo)的作用是在實(shí)現(xiàn)均勻性目標(biāo)過(guò)程中避免風(fēng)道阻力急劇增加，經(jīng)過(guò)反復(fù)試算后，最終調(diào)節(jié)目標(biāo)函數(shù)中權(quán)重因子ω1與ω2分別為400 和1。整體目標(biāo)函數(shù)公式為

由于流體拓?fù)鋬?yōu)化中存在流體與固體的轉(zhuǎn)化關(guān)系，為防止優(yōu)化進(jìn)程中進(jìn)口或出口處出現(xiàn)固體阻塞從而影響優(yōu)化迭代，將二維風(fēng)道主體內(nèi)腔作為設(shè)計(jì)區(qū)域，進(jìn)口區(qū)域與出口區(qū)域作為非設(shè)計(jì)區(qū)域，優(yōu)化的初始設(shè)計(jì)布局如圖5所示。

圖5 風(fēng)道優(yōu)化問(wèn)題的設(shè)計(jì)模型Fig.5 Design model of air duct optimization problem

4 送風(fēng)風(fēng)道優(yōu)化結(jié)果分析

4.1 優(yōu)化結(jié)果

計(jì)算使用高性能工作站(2 個(gè)E5-2650 v4 處理器，24 核，128 G 內(nèi)存)，歷時(shí)約4 h，迭代每步耗時(shí)約140 s。伴隨著優(yōu)化進(jìn)程，在二維設(shè)計(jì)域腔內(nèi)逐漸生成邊界較為清晰的固體區(qū)域，約經(jīng)過(guò)100次迭代后目標(biāo)函數(shù)收斂。優(yōu)化進(jìn)程如圖6所示。

圖6 優(yōu)化進(jìn)程Fig.6 Optimization processes

優(yōu)化得到的構(gòu)型在邊緣處仍然存在少許中間滲透率單元，為得到清晰的幾何邊界，選取γ≤0.5的區(qū)域作為固體域，γ＞0.5作為流體域，最終得到的拓?fù)錁?gòu)型如圖7所示。

圖7 優(yōu)化后拓?fù)錁?gòu)型Fig.7 Layout of the optimized topology

結(jié)合上述二維設(shè)計(jì)空間的優(yōu)化結(jié)果，提取優(yōu)化所得的二維固體域，拉伸0.29 m 得到三維拓?fù)鋵?dǎo)流結(jié)構(gòu)，使其高度與初始風(fēng)道內(nèi)腔高度相同，并將此固體域賦予到初始風(fēng)道內(nèi)腔的相應(yīng)位置上，重建后可得新型三維送風(fēng)風(fēng)道模型。新型送風(fēng)風(fēng)道結(jié)構(gòu)如圖8所示。

圖8 新型送風(fēng)風(fēng)道結(jié)構(gòu)Fig.8 Three-dimensional mapping model of the optimized configuration

4.2 優(yōu)化后風(fēng)道出口速度分布

為驗(yàn)證新型三維送風(fēng)風(fēng)道內(nèi)流場(chǎng)性能，將重建后的三維送風(fēng)風(fēng)道模型進(jìn)行仿真分析，其物性參數(shù)、進(jìn)出口邊界條件與初始風(fēng)道的物性參數(shù)保持一致，內(nèi)部拓?fù)鋵?dǎo)流結(jié)構(gòu)與風(fēng)道外壁設(shè)置為非滑移壁面條件。數(shù)值模擬收斂后得到風(fēng)道出口處速度分布情況以及法向速度矢量分布情況，如圖9所示。

圖9 優(yōu)化后風(fēng)道出口的速度分布情況Fig.9 Velocity distribution at the outlet of the optimized air duct

風(fēng)道出口處的速度分布較優(yōu)化前有明顯改善，進(jìn)口側(cè)出口處的低風(fēng)速區(qū)域變小，尖端側(cè)出口處的高風(fēng)速區(qū)域消失。從出口處法向速度矢量分布可知，進(jìn)口側(cè)的出口區(qū)域風(fēng)速有所增加，且回流現(xiàn)象基本消除；尖端側(cè)的出口區(qū)域風(fēng)速過(guò)高情況減弱，風(fēng)速分布更均勻。

根據(jù)式(1)與(2)，優(yōu)化后的送風(fēng)風(fēng)道出口速度均勻性系數(shù)γt為0.889 5，法向速度均勻性系數(shù)為0.785 4，進(jìn)出口平均壓差為20.909 1 Pa。優(yōu)化后的風(fēng)道送風(fēng)性能提升顯著，兩項(xiàng)均勻性系數(shù)分別增加12%和34%，同時(shí)，進(jìn)出口壓差減小12%，優(yōu)化前后對(duì)比如表1所示。

表1 初始風(fēng)道與優(yōu)化后風(fēng)道送風(fēng)性能對(duì)比Table 1 Comparison of the air supply performance between the initial air duct and the optimized air duct

4.3 拓?fù)浣Y(jié)果評(píng)價(jià)

經(jīng)過(guò)優(yōu)化后所得的新型送風(fēng)風(fēng)道并沒(méi)有改變幾何外形，繼續(xù)沿用初始送風(fēng)風(fēng)道的外部構(gòu)型，只是將拓?fù)鋬?yōu)化所得的導(dǎo)流結(jié)構(gòu)添加到風(fēng)道內(nèi)腔中。此導(dǎo)流結(jié)構(gòu)外形并不規(guī)則，相對(duì)傳統(tǒng)的導(dǎo)流板制造更復(fù)雜，但其使得送風(fēng)風(fēng)道的流場(chǎng)性能大幅度提升，在提高風(fēng)道送風(fēng)均勻性的同時(shí)，降低了進(jìn)出口平均壓差，對(duì)于地鐵送風(fēng)風(fēng)道優(yōu)化設(shè)計(jì)具有較好的參考價(jià)值。

5 結(jié)論

1) 地鐵送風(fēng)風(fēng)道出口較為狹長(zhǎng)，在額定送風(fēng)條件下，初始送風(fēng)風(fēng)道在進(jìn)口側(cè)處的出口區(qū)域存在回流現(xiàn)象，在尾部尖端的出口區(qū)域存在風(fēng)速過(guò)高情況。

2) 本文將變密度拓?fù)鋬?yōu)化理論應(yīng)用于湍流流場(chǎng)中，在k-ε湍流模型的動(dòng)量方程、湍流動(dòng)能方程以及湍流耗散率方程中添加類似的達(dá)西插值項(xiàng)來(lái)實(shí)現(xiàn)反滲流固體與流體單元的相互轉(zhuǎn)化，并采用伴隨靈敏度法求解目標(biāo)梯度，提出了適用于湍流流動(dòng)的流體拓?fù)鋬?yōu)化方法。

3) 將此拓?fù)鋬?yōu)化方法應(yīng)用于地鐵送風(fēng)風(fēng)道的優(yōu)化設(shè)計(jì)，目標(biāo)函數(shù)采用多目標(biāo)函數(shù)加權(quán)處理的方式，考慮風(fēng)道出口處氣流均勻性偏差與進(jìn)出口平均壓差兩部分。優(yōu)化后的地鐵送風(fēng)風(fēng)道出口均勻性顯著提高，初始風(fēng)道出口處存在的回流現(xiàn)象、尖端風(fēng)速過(guò)高情況基本消除。風(fēng)道出口處截面速度均勻性系數(shù)增加12%，法向速度均勻性系數(shù)增加34%，進(jìn)出口平均壓差減小12%。