高艷婧,倪藝銘,劉宏偉,林勇剛
(浙江大學(xué) 流體動力與機電系統(tǒng)國家重點實驗室,浙江 杭州,310027)
我國近海岸海流因水土流失導(dǎo)致含沙量較高[1-3]。泥沙中含有大量的固體顆粒和其他雜質(zhì),會對水力機械水動力性能造成嚴(yán)重影響,加速水力機械的磨損,嚴(yán)重影響其工作效率和使用壽命[4]。浙江省杭州灣位于河流入海口處,屬于近海沿岸海域,水域較淺,底部泥沙在水流的曳力作用下懸浮移動,導(dǎo)致海流能機組在含沙的工況下運行[5]。顆粒與流體之間的相互作用是否會改變翼型截面的流場以及影響葉輪的捕獲功率尚且未知,因此,研究兩相流環(huán)境下的葉輪捕獲功率可以為后續(xù)海流能機組選址以及葉片設(shè)計提供參考。
國內(nèi)外大部分的研究均針對海流能機組葉片表面因沙粒沖擊、空化以及海水腐蝕引起的粗糙度增加而造成的功率損失[6-8]展開,對于水沙環(huán)境直接對海流能機組捕獲功率的影響的研究甚少。ORME等[9]利用風(fēng)洞研究了藤壺對機翼升力系數(shù)和阻力系數(shù)的影響,發(fā)現(xiàn)機翼上覆蓋不同尺寸和密度的理想藤壺時,翼型升阻比降低。BATTEN等[10]基于葉素動量理論(BEM)研究了污垢對葉片粗糙度增加的潛在影響,發(fā)現(xiàn)在葉片被污染的情況下阻力系數(shù)增加,而升力系數(shù)沒有改變,基于葉素動量理論預(yù)測在高的葉尖速下功率系數(shù)下降6%~8%。WALKER 等[11]通過風(fēng)洞實驗測試了NACA63618 翼型在清潔和粗糙條件下的升阻力系數(shù),而且通過水池實驗測試了模型葉輪在不同表面條件下的捕獲性能,發(fā)現(xiàn)升力最大下降幅度的平均值為11%,阻力在重度粗糙下能上升150%,捕獲功率最大降低19%。而對于沙粒與流體的相互作用對海流能葉片的捕獲功率的影響尚不明確。不同海域的顆粒直徑、粒徑分布以及顆粒濃度均有差異[12],研究這些因素對于葉輪捕獲功率的影響,可為海流能機組選址和葉片設(shè)計提供參考。
在本研究中,以120 kW 海流能葉片翼型NACA63-為研究對象,采用CFD-DPM 模型(考慮顆粒與流體相互作用)模擬不同顆粒屬性下翼型的升力系數(shù)和阻力系數(shù)。根據(jù)沿同一顆粒編號軌跡運動的顆粒速度和有/無顆粒的流體速度,分析顆粒對流體速度造成的影響,再依據(jù)顆粒偏離流線的程度揭示在流體速度變化的情況下的升阻力系數(shù)的變化。
翼型是海流能葉片的重要組成部分,翼型參數(shù)直接影響海流能葉片的性能。NACA63-翼型是海流能葉片最常用的翼型。水流經(jīng)過翼型前緣附近的壓力最高點即駐點后,分成上、下兩股,分別沿翼型上、下表面流過,在后緣重新匯合,翼型上下表面分別形成變化的壓力。根據(jù)伯努利方程,當(dāng)水流流經(jīng)翼型凸面時,平均流速增大、壓力減小,此面被稱為吸力面。當(dāng)水流流經(jīng)翼型凹面時,平均流速減小、壓力增大,此面被稱為壓力面。葉片功率與翼型升力系數(shù)呈正相關(guān),與翼型阻力系數(shù)呈負(fù)相關(guān)。
翼型升力源于吸/壓力面壓差,吸力面?zhèn)攘黧w速度先增大后減小,壓力面?zhèn)攘黧w速度先減小再增大。同一流線位置的流體速度變化大小可以反映升力系數(shù)大小,即同一流線位置流體速度變化越大,升力系數(shù)越大。在相同的雷諾數(shù)和攻角下,流線距離翼型表面越遠,沿流線流體速度變化越小,如圖1所示。當(dāng)升力系數(shù)變化值相同時,流線距離翼型吸/壓力面越遠,沿同一流線流體速度變化值越小,如圖2所示。由圖2可以推斷,沿同一流線流體速度變化值相同時,流線距離翼型吸/壓力面越遠,升力系數(shù)變化值越大。翼型阻力系數(shù)也與雷諾數(shù)和攻角相關(guān),隨著流體速度增大,翼型阻力系數(shù)減小。
圖1 不同流線位置的流體速度Fig. 1 Fluid velocities at different streamline positions
圖2 升力系數(shù)變化值相同(攻角從6°到8°)時不同流線位置的流體速度變化值Fig. 2 Change values of fluid velocity at different streamline positions when the lift coefficient changes by the same value(angle of attack from 6° to 8°)
考慮到隨著顆粒直徑或顆粒形狀因子的變化,顆粒偏離流線的程度不同,故在研究不同直徑、不同濃度和不同顆粒形狀因子下的顆粒速度和無顆粒流體速度時,以同一顆粒編號的顆粒軌跡為研究對象。由遠離吸力側(cè)指向遠離壓力側(cè)方向,顆粒編號逐漸增大,顆粒編號為30000 和編號為36000的顆粒軌跡位置如圖3所示。
圖3 顆粒編號為30000和編號為36000的顆粒軌跡位置Fig. 3 Particle track positions when particle IDs are 30000 and 36000
以葉尖NACA63-翼型為研究對象,基于葉素理論[13-14],將120 kW 海流能葉片分成為若干翼型葉素。在ICEM軟件中對翼型和流域進行幾何建模和網(wǎng)格劃分。圖4所示為計算域和邊界條件,以翼型尾部為圓心,上游與翼型的距離為弦長c的12.5倍,以減小邊界對翼型附近流體流動的影響;下游與翼型的距離為16 倍弦長,以保證流體能夠充分流動;展向取1倍弦長。顆粒噴射面與翼型的距離為弦長的6倍,顆粒噴射面是在FLUENT中創(chuàng)建的虛擬面,通過設(shè)置節(jié)點數(shù)來保證顆粒噴射面網(wǎng)格均勻[15],顆粒噴射面長度為10 倍弦長,保證升阻力系數(shù)不再隨顆粒噴射面長度的增加而變化。
進口處采用速度入口邊界,出口處采用壓力出口邊界,將翼型定義為無滑移壁面,顆粒在顆粒噴射面以與流體相同的速度進入流場。顆粒在翼型壁面時設(shè)置為“reflect”,在出口處設(shè)置為“escaped”。湍流模型采用SSTk-ω模型,邊界層湍流和自由剪切流均具有較理想的精度與穩(wěn)定性。速度與壓力的耦合選擇SIMPLE算法,動量、湍流動能和耗散率均采用二階迎風(fēng)格式,并采用基于壓力耦合的隱式求解器[16]。
圖4 計算域和邊界條件Fig. 4 Calculation domain and boundary conditions
流場網(wǎng)格劃分采用結(jié)構(gòu)化網(wǎng)格方法,并用C-Block 方法對翼型附近的網(wǎng)格進行加密[17]?;赟STk-ω模型,首層網(wǎng)格量綱一高度y+要小于1[18],首層網(wǎng)格高度為0.02 mm。
為保證求解的準(zhǔn)確性,進行網(wǎng)格無關(guān)性驗證[19]。隨著網(wǎng)格數(shù)量的增加(從網(wǎng)格1 到網(wǎng)格6),翼型下游8倍弦長位置處的速度以及湍流動能基本上不發(fā)生改變(圖5)??紤]到計算精度和時間成本,選擇網(wǎng)格3,翼型表面分布340 個網(wǎng)格節(jié)點,翼型壁面法向分布180個網(wǎng)格節(jié)點,網(wǎng)格增長率為1.2,網(wǎng)格單元總數(shù)為32 萬個左右。圖6 所示為結(jié)構(gòu)化網(wǎng)格劃分后的流體域網(wǎng)格和翼型表面網(wǎng)格。
圖5 不同網(wǎng)格分辨率翼型下游8倍弦長位置處的速度和湍流動能變化曲線Fig. 5 Velocity and turbulent kinetic energy curves at 8 chord lengths downstream of airfoils with different grid resolutions
圖6 計算域網(wǎng)格Fig. 6 Computational meshes used in the simulation
2.3.1 SSTk-ω模型驗證
圖7所示為無顆粒海水條件下的翼型升阻力系數(shù)仿真結(jié)果與通過風(fēng)洞實驗得到的升阻力系數(shù)[20](相同雷諾數(shù))的對比,可見升阻力系數(shù)仿真結(jié)果與實驗結(jié)果[20]具有較高的吻合性。
圖7 升阻力系數(shù)仿真和實驗結(jié)果對比Fig. 7 Comparison of the simulated and experimental results of lift and drag coefficients
2.3.2 DPM模型驗證
將受沙粒沖擊的120 kW 海流能機組葉片的磨損數(shù)據(jù)作為實驗數(shù)據(jù),120 kW 海流能機組葉片參數(shù)如表1 所示,包含30 個翼型葉素。選擇考慮顆粒形狀因子和沖擊角的Ahlert侵蝕模型[21],仿真模型中顆粒屬性和葉片壁面屬性如表2所示。詳細(xì)求解過程如下:首先,通過求解雷諾平均Navier-Stokes方程得到流體流動的解;然后,考慮流體對顆粒的阻力獲得顆粒碰撞信息;最后,依據(jù)顆粒碰撞信息通過Ahlert侵蝕模型預(yù)測磨損特性[22]。
表1 120 kW海流能機組葉片參數(shù)Table 1 Parameters of the 120 kW tidal current turbine blade
表2 材料屬性Table 2 Material properties
圖8(a)所示為實際海況下120 kW 海流能葉片的沙粒沖蝕圖,圖8(b)所示為考慮粒徑分布(表3)的120 kW 海流能葉片仿真沖蝕云圖。圖9 所示為最大侵蝕位置(以翼型前緣點為原點,沿弦長方向沙粒對壁面造成的最大侵蝕距離)的仿真和實驗結(jié)果對比,可見最大侵蝕位置從葉根到葉尖逐漸向前緣移動,與實驗結(jié)果較吻合。
圖8 120 kW海流能葉片的顆粒沖擊沖蝕圖Fig. 8 Erosion diagrams of the 120 kW tidal current turbine blade
表3 舟山海況顆粒粒徑分布Table 3 Particle size distribution in Zhoushan sea condition
圖9 最大侵蝕位置仿真和實驗結(jié)果對比Fig. 9 Comparison of the maximum erosion location of the simulation and the experiment results
本研究利用CFD-DPM 模型計算120 kW 海流能機組葉片各翼型的升力系數(shù)和阻力系數(shù),再根據(jù)葉素動量理論計算葉輪的輸出轉(zhuǎn)矩和功率。不同流速下葉輪功率的仿真結(jié)果與實驗結(jié)果的對比見圖10。由圖10 可以看出,模擬的葉輪功率與實驗結(jié)果較吻合,驗證了DPM模型的準(zhǔn)確性。
仿真算例的運行條件如表4所示,研究不同顆粒屬性下翼型升阻力系數(shù)。在本算例中,顆粒直徑dP為500 μm,質(zhì)量濃度CP為1.5 g/L,形狀因子f為1,弦長L為1 m,攻角α為6°,相對流體速度U為16 m/s。
表4 仿真算例運行條件Table 4 Operating conditions of the simulation cases
3.1.1 顆粒直徑對翼型升力系數(shù)的影響
圖11所示為翼型升力系數(shù)隨顆粒直徑的變化。由圖11可知:當(dāng)顆粒直徑小于100 μm時,翼型升力系數(shù)大于無顆粒時的翼型升力系數(shù);當(dāng)顆粒直徑大于100 μm 時,翼型升力系數(shù)小于無顆粒時的翼型升力系數(shù);隨著顆粒直徑增加,翼型升力系數(shù)在先減小后增大的趨勢中波動。
圖11 翼型升力系數(shù)隨顆粒直徑的變化Fig. 11 Change of airfoil lift coefficient with particle diameter
圖12所示為無顆粒時的流線,圖13所示為不同顆粒直徑下的顆粒軌跡。由圖12 和圖13 可知:當(dāng)顆粒直徑為20 μm 時,顆粒軌跡與流線基本重合;隨著顆粒直徑的增加,顆粒軌跡偏離流線的程度增大。
圖12 無顆粒流線Fig. 12 The streamlines without particles
圖13 不同顆粒直徑下的顆粒軌跡Fig. 13 Particle trajectories with different particle diameters
圖14 所示為當(dāng)顆粒直徑為20 μm 時,沿同一顆粒編號軌跡運動的顆粒和流體速度曲線。由于顆粒直徑較小、慣性較小,顆粒軌跡主要受流體阻力控制[23-25]。由圖14 可知:在顆粒噴射面顆粒速度急劇下降,隨流體速度變化而變化,等效于流體密度增加,從而導(dǎo)致翼型升力增加。
圖14 沿5號顆粒軌跡運動的顆粒速度和流體速度(顆粒直徑為20 μm)Fig. 14 Velocities of particle and fluid along the particle trajectory of the No.5 particle(particle diameter is 20 μm)
圖15所示為當(dāng)顆粒直徑大于100 μm時,沿相同顆粒編號軌跡運動的顆粒速度和無顆粒流體速度。由圖15 可知:受翼型影響,壓力面流體速度下降,流體對顆粒的阻力同樣導(dǎo)致顆粒速度下降。但流體速度的減小程度大于顆粒速度的減小程度,即顆粒速度大于流體速度??紤]顆粒對流體的影響,壓力面流體速度減小的程度降低,導(dǎo)致升力系數(shù)減小。
圖15 不同顆粒直徑下沿5號顆粒軌跡運動的顆粒速度和無顆粒流體的速度Fig. 15 Particle velocity and the velocity of particle-free fluid along the trajectory of No.5 particle with different particle diameters
從圖15 還可以看出:隨著顆粒直徑增加,Up-fn增加(Up-fn為沿相同顆粒編號軌跡運動的顆粒速度與無顆粒流體的速度之差),單個顆粒對流體的影響增大,Ufp-fn增加(Ufp-fn為沿相同顆粒編號軌跡運動的下流體速度與無顆粒流體的速度之差);隨著顆粒直徑繼續(xù)增大,顆粒質(zhì)量濃度不變,與流體作用的顆粒數(shù)減少,導(dǎo)致Ufp-fn減小;隨著顆粒直徑增大,Ufp-fn先增大后減小,沿同一顆粒軌跡的流體速度變化值先增大后減小。
隨著顆粒直徑的增大,顆粒的慣性增大,具有相同顆粒編號的顆粒遠離壓力面,如圖13所示。從圖13 可見:當(dāng)沿同一顆粒編號軌跡運動的流體速度變化值相同時,大直徑顆粒會使升力系數(shù)(相對于無顆粒時的升力系數(shù))減小更多,即升力系數(shù)更?。挥捎谕痪幪栴w粒的位置不同導(dǎo)致流體速度變化值與升力系數(shù)的變化值不具有統(tǒng)一性,因此,隨著直徑增大,升力系數(shù)在先減小后增大的趨勢中波動;當(dāng)顆粒直徑為2 500 μm 時,升力系數(shù)變化最大,和無顆粒相比升力系數(shù)減小了0.141%。
3.1.2 顆粒直徑對翼型阻力系數(shù)的影響
圖16所示為翼型阻力系數(shù)隨顆粒直徑的變化。由圖16 可知:隨著顆粒直徑增大,翼型阻力系數(shù)先增大后減小。
圖16 翼型阻力系數(shù)隨顆粒直徑的變化Fig. 16 Change of airfoil drag coefficient with particle diameter
圖17 所示為當(dāng)顆粒直徑為1 500 μm 時,沿23000 號和45000 號顆粒軌跡運動的顆粒流體速度和無顆粒流體速度,其中,路徑長度以顆粒噴射面作為起點。由圖17 可知:翼型阻力系數(shù)隨速度的減小而增大,與無顆粒流體速度相比,等效流體速度的減小程度越大,翼型阻力系數(shù)越大;隨著粒徑的增大,Ufp-fn先增大后減小,等效流體速度減小程度先增大后減小,阻力系數(shù)先增大后減小;當(dāng)顆粒直徑為1 500 μm時阻力系數(shù)變化最大,與無顆粒流體相比阻力系數(shù)增加了1.872%。
圖17 沿45000號和23000號顆粒軌跡運動的無顆粒流體速度與顆粒作用下的流體速度Fig. 17 Velocities of fluid with or without particles influence along the trajectories of No.23000 and No.45000 particles
3.2.1 顆粒濃度對翼型升力系數(shù)的影響
翼型升力系數(shù)隨顆粒質(zhì)量濃度的變化如圖18所示。由圖18 可知:翼型升力系數(shù)與顆粒質(zhì)量濃度呈線性關(guān)系。圖19 所示為不同質(zhì)量濃度下沿5號顆粒軌跡運動的顆粒速度和無顆粒流體的速度。由圖19 可知:隨著顆粒質(zhì)量濃度的增加,顆粒與無顆粒流體速度差保持不變,即單個顆粒對流體速度變化的影響不變。圖20 所示為不同顆粒質(zhì)量濃度下顆粒軌跡。從圖20 可以看出:隨著顆粒質(zhì)量濃度增加,顆粒偏離流線的程度不變,單個顆粒對翼型升力系數(shù)變化的影響保持不變;與流場相互作用的顆粒數(shù)與顆粒質(zhì)量濃度成正比,即對于相同的粒徑,當(dāng)顆粒質(zhì)量濃度增加值相同時,與流體相互作用的顆粒數(shù)增加值相同。因此,翼型升力系數(shù)與顆粒質(zhì)量濃度呈線性關(guān)系。
圖18 翼型升力系數(shù)隨顆粒質(zhì)量濃度的變化Fig. 18 Change of airfoil lift coefficient versus particle mass concentration
當(dāng)顆粒質(zhì)量濃度為2 g/L,粒徑為2 500 μm時,升力系數(shù)降低了0.188%??紤]到舟山海況中的最大顆粒質(zhì)量濃度為5.522 g/L,當(dāng)粒徑為2 500 μm時,升力系數(shù)降低了0.478%。
3.2.2 顆粒質(zhì)量濃度對翼型阻力系數(shù)的影響
圖21 所示為翼型阻力系數(shù)隨顆粒質(zhì)量濃度的變化,由圖21 可知:隨著顆粒質(zhì)量濃度增加,翼型阻力系數(shù)呈線性增大;當(dāng)粒徑為1 500 μm,顆粒質(zhì)量濃度分別為2.000 g/L 和5.522 g/L 時,與無顆粒流體相比阻力系數(shù)分別增加了2.418% 和6.676%。
圖21 翼型阻力系數(shù)隨顆粒質(zhì)量濃度的變化Fig. 21 Change of airfoil drag coefficient with particle mass concentration
3.3.1 顆粒形狀因子對翼型升力系數(shù)的影響
顆粒形狀是造成海流能機組磨損嚴(yán)重的重要原因之一,下面對顆粒形狀因子對翼型升阻系數(shù)的影響進行研究。顆粒的形狀因子f指與顆粒體積相同的球體表面積與顆粒實際表面積的比值,其定義如下:
式中:a為與非球形顆粒體積相同的球形顆粒表面積;As為非球形顆粒的表面積。
圖22 所示為不同顆粒形狀因子下翼型升力系數(shù)隨粒徑的變化。由圖22 可知:隨著顆粒形狀因子的增大,使顆粒翼型升力系數(shù)大于無顆粒時流體翼型升力系數(shù)的臨界顆粒直徑減小。小形狀因子的顆粒具有較大的表面積,受到流體更大的阻力,促使顆粒跟隨流體移動。不同形狀因子下的顆粒軌跡如圖23所示。由圖23可知:隨著顆粒形狀因子增大,顆粒偏離流線的程度增大。顆粒能夠與流體保持相對靜止運動的臨界顆粒直徑減小,因此,能夠使顆粒翼型升力系數(shù)大于無顆粒時流體翼型升力系數(shù)的臨界顆粒直徑減小。
圖22 不同顆粒形狀因子下翼型升力系數(shù)隨粒徑的變化Fig. 22 Change of airfoil lift coefficient with particle diameter under different particle shape factors
圖23 不同顆粒形狀因子下顆粒軌跡Fig. 23 Particle trajectories under different particle shape factors
圖24 所示為翼型升力系數(shù)隨顆粒形狀因子的變化,由圖24 可知:隨著顆粒形狀因子增加,翼型升力系數(shù)減小。圖25 所示為不同顆粒形狀因子下沿相同編號顆粒軌跡運動的顆粒速度與無顆粒流體速度。由圖25 可知:隨顆粒形狀因子增加,顆粒偏離流體流線的程度增大,Up-fn增大。圖26所示為不同顆粒形狀因子下沿相同顆粒編號軌跡運動的顆粒速度和無顆粒流體速度之差。由圖26可知:顆粒對流體的影響增強,沿同一顆粒軌跡的流體速度變化程度減小,Ufp-fn增大。
圖24 翼型升力系數(shù)隨顆粒形狀因子的變化Fig. 24 Change of airfoil lift coefficient with particle shape factor
圖25 不同顆粒形狀因子下沿相同編號顆粒軌跡運動的顆粒速度與無顆粒流體速度Fig. 25 Particle velocity and the particle-free fluid velocity along the same particle ID trajectory under different particle shape factors
經(jīng)分析可得,隨著顆粒形狀因子增加,相同顆粒編號的顆粒遠離壓力表面。在沿同一編號顆粒軌跡運動的流體速度變化值相同的情況下,形狀因子大的顆粒升力系數(shù)(與無顆粒升力系數(shù)相比)降低程度更大,升力系數(shù)較小。而且隨著顆粒形狀因子增加,Ufp-fn增大,沿同一編號顆粒軌跡運動的流體速度變化值增大。故隨顆粒形狀因子增大,升力系數(shù)減小。當(dāng)形狀因子為0.9、顆粒直徑為2 500 μm 時,升力系數(shù)減小幅度最大,與無顆粒時相比減少了0.18%。
3.3.2 顆粒形狀因子對翼型阻力系數(shù)的影響
圖27所示為翼型阻力系數(shù)隨形狀因子的變化。由圖27 可知:翼型阻力系數(shù)隨顆粒形狀因子的增大而增大;隨著顆粒形狀因子的增加,Ufp-fn增大,等效于相對流體速度的減少程度增大,導(dǎo)致阻力系數(shù)增大;當(dāng)形狀因子為0.9,顆粒直徑為1 500 μm時,阻力系數(shù)變化最大,與無顆粒時相比增加了1.315%。
1) 當(dāng)顆粒直徑小于100 μm 時,顆粒慣性小,顆粒和流體可以作為單相,相當(dāng)于增加了流體密度,升力系數(shù)大于無顆粒時的升力系數(shù);當(dāng)顆粒直徑大于100 μm 時,顆粒偏離流線,升力系數(shù)小于無顆粒時的升力系數(shù)。
2) 隨著顆粒直徑增大,顆粒慣性增強,顆粒偏離流線,升力系數(shù)在先減小后增大的趨勢中波動。
3) 隨著顆粒形狀因子增大,顆粒表面積減小,顆粒受流體阻力減小,偏離流線程度增大,升力系數(shù)減小,阻力系數(shù)增大。
4) 顆粒質(zhì)量濃度對翼型升阻力系數(shù)影響最大,隨顆粒質(zhì)量濃度增大,升力系數(shù)呈線性增大,阻力系數(shù)呈線性減小,當(dāng)顆粒質(zhì)量濃度為5.522 g/L時,相比于無顆粒時升力系數(shù)減小了0.478%,阻力系數(shù)增大了6.676%。