水沙環境下海流能機組翼型水動性能

高艷婧，倪藝銘，劉宏偉，林勇剛

(浙江大學 流體動力與機電系統國家重點實驗室，浙江 杭州，310027)

我國近海岸海流因水土流失導致含沙量較高[1-3]。泥沙中含有大量的固體顆粒和其他雜質，會對水力機械水動力性能造成嚴重影響，加速水力機械的磨損，嚴重影響其工作效率和使用壽命[4]。浙江省杭州灣位于河流入海口處，屬于近海沿岸海域，水域較淺，底部泥沙在水流的曳力作用下懸浮移動，導致海流能機組在含沙的工況下運行[5]。顆粒與流體之間的相互作用是否會改變翼型截面的流場以及影響葉輪的捕獲功率尚且未知，因此，研究兩相流環境下的葉輪捕獲功率可以為后續海流能機組選址以及葉片設計提供參考。

國內外大部分的研究均針對海流能機組葉片表面因沙粒沖擊、空化以及海水腐蝕引起的粗糙度增加而造成的功率損失[6-8]展開，對于水沙環境直接對海流能機組捕獲功率的影響的研究甚少。ORME等[9]利用風洞研究了藤壺對機翼升力系數和阻力系數的影響，發現機翼上覆蓋不同尺寸和密度的理想藤壺時，翼型升阻比降低。BATTEN等[10]基于葉素動量理論(BEM)研究了污垢對葉片粗糙度增加的潛在影響，發現在葉片被污染的情況下阻力系數增加，而升力系數沒有改變，基于葉素動量理論預測在高的葉尖速下功率系數下降6%～8%。WALKER 等[11]通過風洞實驗測試了NACA63618 翼型在清潔和粗糙條件下的升阻力系數，而且通過水池實驗測試了模型葉輪在不同表面條件下的捕獲性能，發現升力最大下降幅度的平均值為11%，阻力在重度粗糙下能上升150%，捕獲功率最大降低19%。而對于沙粒與流體的相互作用對海流能葉片的捕獲功率的影響尚不明確。不同海域的顆粒直徑、粒徑分布以及顆粒濃度均有差異[12]，研究這些因素對于葉輪捕獲功率的影響，可為海流能機組選址和葉片設計提供參考。

在本研究中，以120 kW 海流能葉片翼型NACA63-為研究對象，采用CFD-DPM 模型(考慮顆粒與流體相互作用)模擬不同顆粒屬性下翼型的升力系數和阻力系數。根據沿同一顆粒編號軌跡運動的顆粒速度和有/無顆粒的流體速度，分析顆粒對流體速度造成的影響，再依據顆粒偏離流線的程度揭示在流體速度變化的情況下的升阻力系數的變化。

1 理論基礎

翼型是海流能葉片的重要組成部分，翼型參數直接影響海流能葉片的性能。NACA63-翼型是海流能葉片最常用的翼型。水流經過翼型前緣附近的壓力最高點即駐點后，分成上、下兩股，分別沿翼型上、下表面流過，在后緣重新匯合，翼型上下表面分別形成變化的壓力。根據伯努利方程，當水流流經翼型凸面時，平均流速增大、壓力減小，此面被稱為吸力面。當水流流經翼型凹面時，平均流速減小、壓力增大，此面被稱為壓力面。葉片功率與翼型升力系數呈正相關，與翼型阻力系數呈負相關。

翼型升力源于吸/壓力面壓差，吸力面側流體速度先增大后減小，壓力面側流體速度先減小再增大。同一流線位置的流體速度變化大小可以反映升力系數大小，即同一流線位置流體速度變化越大，升力系數越大。在相同的雷諾數和攻角下，流線距離翼型表面越遠，沿流線流體速度變化越小，如圖1所示。當升力系數變化值相同時，流線距離翼型吸/壓力面越遠，沿同一流線流體速度變化值越小，如圖2所示。由圖2可以推斷，沿同一流線流體速度變化值相同時，流線距離翼型吸/壓力面越遠，升力系數變化值越大。翼型阻力系數也與雷諾數和攻角相關，隨著流體速度增大，翼型阻力系數減小。

圖1 不同流線位置的流體速度Fig. 1 Fluid velocities at different streamline positions

圖2 升力系數變化值相同(攻角從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°)

考慮到隨著顆粒直徑或顆粒形狀因子的變化，顆粒偏離流線的程度不同，故在研究不同直徑、不同濃度和不同顆粒形狀因子下的顆粒速度和無顆粒流體速度時，以同一顆粒編號的顆粒軌跡為研究對象。由遠離吸力側指向遠離壓力側方向，顆粒編號逐漸增大，顆粒編號為30000 和編號為36000的顆粒軌跡位置如圖3所示。

圖3 顆粒編號為30000和編號為36000的顆粒軌跡位置Fig. 3 Particle track positions when particle IDs are 30000 and 36000

2 計算模型

2.1 計算域和邊界條件

以葉尖NACA63-翼型為研究對象，基于葉素理論[13-14]，將120 kW 海流能葉片分成為若干翼型葉素。在ICEM軟件中對翼型和流域進行幾何建模和網格劃分。圖4所示為計算域和邊界條件，以翼型尾部為圓心，上游與翼型的距離為弦長c的12.5倍，以減小邊界對翼型附近流體流動的影響；下游與翼型的距離為16 倍弦長，以保證流體能夠充分流動；展向取1倍弦長。顆粒噴射面與翼型的距離為弦長的6倍，顆粒噴射面是在FLUENT中創建的虛擬面，通過設置節點數來保證顆粒噴射面網格均勻[15]，顆粒噴射面長度為10 倍弦長，保證升阻力系數不再隨顆粒噴射面長度的增加而變化。

進口處采用速度入口邊界，出口處采用壓力出口邊界，將翼型定義為無滑移壁面，顆粒在顆粒噴射面以與流體相同的速度進入流場。顆粒在翼型壁面時設置為“reflect”，在出口處設置為“escaped”。湍流模型采用SSTk-ω模型，邊界層湍流和自由剪切流均具有較理想的精度與穩定性。速度與壓力的耦合選擇SIMPLE算法，動量、湍流動能和耗散率均采用二階迎風格式，并采用基于壓力耦合的隱式求解器[16]。

圖4 計算域和邊界條件Fig. 4 Calculation domain and boundary conditions

2.2 網格劃分

流場網格劃分采用結構化網格方法，并用C-Block 方法對翼型附近的網格進行加密[17]。基于SSTk-ω模型，首層網格量綱一高度y+要小于1[18]，首層網格高度為0.02 mm。

為保證求解的準確性，進行網格無關性驗證[19]。隨著網格數量的增加(從網格1 到網格6)，翼型下游8倍弦長位置處的速度以及湍流動能基本上不發生改變(圖5)。考慮到計算精度和時間成本，選擇網格3，翼型表面分布340 個網格節點，翼型壁面法向分布180個網格節點，網格增長率為1.2，網格單元總數為32 萬個左右。圖6 所示為結構化網格劃分后的流體域網格和翼型表面網格。

圖5 不同網格分辨率翼型下游8倍弦長位置處的速度和湍流動能變化曲線Fig. 5 Velocity and turbulent kinetic energy curves at 8 chord lengths downstream of airfoils with different grid resolutions

圖6 計算域網格Fig. 6 Computational meshes used in the simulation

2.3 模擬方法驗證

2.3.1 SSTk-ω模型驗證

圖7所示為無顆粒海水條件下的翼型升阻力系數仿真結果與通過風洞實驗得到的升阻力系數[20](相同雷諾數)的對比，可見升阻力系數仿真結果與實驗結果[20]具有較高的吻合性。

圖7 升阻力系數仿真和實驗結果對比Fig. 7 Comparison of the simulated and experimental results of lift and drag coefficients

2.3.2 DPM模型驗證

將受沙粒沖擊的120 kW 海流能機組葉片的磨損數據作為實驗數據，120 kW 海流能機組葉片參數如表1 所示，包含30 個翼型葉素。選擇考慮顆粒形狀因子和沖擊角的Ahlert侵蝕模型[21]，仿真模型中顆粒屬性和葉片壁面屬性如表2所示。詳細求解過程如下：首先，通過求解雷諾平均Navier-Stokes方程得到流體流動的解；然后，考慮流體對顆粒的阻力獲得顆粒碰撞信息；最后，依據顆粒碰撞信息通過Ahlert侵蝕模型預測磨損特性[22]。

表1 120 kW海流能機組葉片參數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 所示為最大侵蝕位置(以翼型前緣點為原點，沿弦長方向沙粒對壁面造成的最大侵蝕距離)的仿真和實驗結果對比，可見最大侵蝕位置從葉根到葉尖逐漸向前緣移動，與實驗結果較吻合。

圖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 最大侵蝕位置仿真和實驗結果對比Fig. 9 Comparison of the maximum erosion location of the simulation and the experiment results

本研究利用CFD-DPM 模型計算120 kW 海流能機組葉片各翼型的升力系數和阻力系數，再根據葉素動量理論計算葉輪的輸出轉矩和功率。不同流速下葉輪功率的仿真結果與實驗結果的對比見圖10。由圖10 可以看出，模擬的葉輪功率與實驗結果較吻合，驗證了DPM模型的準確性。

3 結果及分析

仿真算例的運行條件如表4所示，研究不同顆粒屬性下翼型升阻力系數。在本算例中，顆粒直徑dP為500 μm，質量濃度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 顆粒直徑

3.1.1 顆粒直徑對翼型升力系數的影響

圖11所示為翼型升力系數隨顆粒直徑的變化。由圖11可知：當顆粒直徑小于100 μm時，翼型升力系數大于無顆粒時的翼型升力系數；當顆粒直徑大于100 μm 時，翼型升力系數小于無顆粒時的翼型升力系數；隨著顆粒直徑增加，翼型升力系數在先減小后增大的趨勢中波動。

圖11 翼型升力系數隨顆粒直徑的變化Fig. 11 Change of airfoil lift coefficient with particle diameter

圖12所示為無顆粒時的流線，圖13所示為不同顆粒直徑下的顆粒軌跡。由圖12 和圖13 可知：當顆粒直徑為20 μm 時，顆粒軌跡與流線基本重合；隨著顆粒直徑的增加，顆粒軌跡偏離流線的程度增大。

圖12 無顆粒流線Fig. 12 The streamlines without particles

圖13 不同顆粒直徑下的顆粒軌跡Fig. 13 Particle trajectories with different particle diameters

圖14 所示為當顆粒直徑為20 μm 時，沿同一顆粒編號軌跡運動的顆粒和流體速度曲線。由于顆粒直徑較小、慣性較小，顆粒軌跡主要受流體阻力控制[23-25]。由圖14 可知：在顆粒噴射面顆粒速度急劇下降，隨流體速度變化而變化，等效于流體密度增加，從而導致翼型升力增加。

圖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所示為當顆粒直徑大于100 μm時，沿相同顆粒編號軌跡運動的顆粒速度和無顆粒流體速度。由圖15 可知：受翼型影響，壓力面流體速度下降，流體對顆粒的阻力同樣導致顆粒速度下降。但流體速度的減小程度大于顆粒速度的減小程度，即顆粒速度大于流體速度。考慮顆粒對流體的影響，壓力面流體速度減小的程度降低，導致升力系數減小。

圖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為沿相同顆粒編號軌跡運動的下流體速度與無顆粒流體的速度之差)；隨著顆粒直徑繼續增大，顆粒質量濃度不變，與流體作用的顆粒數減少，導致Ufp-fn減小；隨著顆粒直徑增大，Ufp-fn先增大后減小，沿同一顆粒軌跡的流體速度變化值先增大后減小。

隨著顆粒直徑的增大，顆粒的慣性增大，具有相同顆粒編號的顆粒遠離壓力面，如圖13所示。從圖13 可見：當沿同一顆粒編號軌跡運動的流體速度變化值相同時，大直徑顆粒會使升力系數(相對于無顆粒時的升力系數)減小更多，即升力系數更小；由于同一編號顆粒的位置不同導致流體速度變化值與升力系數的變化值不具有統一性，因此，隨著直徑增大，升力系數在先減小后增大的趨勢中波動；當顆粒直徑為2 500 μm 時，升力系數變化最大，和無顆粒相比升力系數減小了0.141%。

3.1.2 顆粒直徑對翼型阻力系數的影響

圖16所示為翼型阻力系數隨顆粒直徑的變化。由圖16 可知：隨著顆粒直徑增大，翼型阻力系數先增大后減小。

圖16 翼型阻力系數隨顆粒直徑的變化Fig. 16 Change of airfoil drag coefficient with particle diameter

圖17 所示為當顆粒直徑為1 500 μm 時，沿23000 號和45000 號顆粒軌跡運動的顆粒流體速度和無顆粒流體速度，其中，路徑長度以顆粒噴射面作為起點。由圖17 可知：翼型阻力系數隨速度的減小而增大，與無顆粒流體速度相比，等效流體速度的減小程度越大，翼型阻力系數越大；隨著粒徑的增大，Ufp-fn先增大后減小，等效流體速度減小程度先增大后減小，阻力系數先增大后減小；當顆粒直徑為1 500 μm時阻力系數變化最大，與無顆粒流體相比阻力系數增加了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 顆粒質量濃度的影響

3.2.1 顆粒濃度對翼型升力系數的影響

翼型升力系數隨顆粒質量濃度的變化如圖18所示。由圖18 可知：翼型升力系數與顆粒質量濃度呈線性關系。圖19 所示為不同質量濃度下沿5號顆粒軌跡運動的顆粒速度和無顆粒流體的速度。由圖19 可知：隨著顆粒質量濃度的增加，顆粒與無顆粒流體速度差保持不變，即單個顆粒對流體速度變化的影響不變。圖20 所示為不同顆粒質量濃度下顆粒軌跡。從圖20 可以看出：隨著顆粒質量濃度增加，顆粒偏離流線的程度不變，單個顆粒對翼型升力系數變化的影響保持不變；與流場相互作用的顆粒數與顆粒質量濃度成正比，即對于相同的粒徑，當顆粒質量濃度增加值相同時，與流體相互作用的顆粒數增加值相同。因此，翼型升力系數與顆粒質量濃度呈線性關系。

圖18 翼型升力系數隨顆粒質量濃度的變化Fig. 18 Change of airfoil lift coefficient versus particle mass concentration

當顆粒質量濃度為2 g/L，粒徑為2 500 μm時，升力系數降低了0.188%。考慮到舟山海況中的最大顆粒質量濃度為5.522 g/L，當粒徑為2 500 μm時，升力系數降低了0.478%。

3.2.2 顆粒質量濃度對翼型阻力系數的影響

圖21 所示為翼型阻力系數隨顆粒質量濃度的變化，由圖21 可知：隨著顆粒質量濃度增加，翼型阻力系數呈線性增大；當粒徑為1 500 μm，顆粒質量濃度分別為2.000 g/L 和5.522 g/L 時，與無顆粒流體相比阻力系數分別增加了2.418% 和6.676%。

圖21 翼型阻力系數隨顆粒質量濃度的變化Fig. 21 Change of airfoil drag coefficient with particle mass concentration

3.3 顆粒形狀因子的影響

3.3.1 顆粒形狀因子對翼型升力系數的影響

顆粒形狀是造成海流能機組磨損嚴重的重要原因之一，下面對顆粒形狀因子對翼型升阻系數的影響進行研究。顆粒的形狀因子f指與顆粒體積相同的球體表面積與顆粒實際表面積的比值，其定義如下：

式中：a為與非球形顆粒體積相同的球形顆粒表面積；As為非球形顆粒的表面積。

圖22 所示為不同顆粒形狀因子下翼型升力系數隨粒徑的變化。由圖22 可知：隨著顆粒形狀因子的增大，使顆粒翼型升力系數大于無顆粒時流體翼型升力系數的臨界顆粒直徑減小。小形狀因子的顆粒具有較大的表面積，受到流體更大的阻力，促使顆粒跟隨流體移動。不同形狀因子下的顆粒軌跡如圖23所示。由圖23可知：隨著顆粒形狀因子增大，顆粒偏離流線的程度增大。顆粒能夠與流體保持相對靜止運動的臨界顆粒直徑減小，因此，能夠使顆粒翼型升力系數大于無顆粒時流體翼型升力系數的臨界顆粒直徑減小。

圖22 不同顆粒形狀因子下翼型升力系數隨粒徑的變化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 所示為翼型升力系數隨顆粒形狀因子的變化，由圖24 可知：隨著顆粒形狀因子增加，翼型升力系數減小。圖25 所示為不同顆粒形狀因子下沿相同編號顆粒軌跡運動的顆粒速度與無顆粒流體速度。由圖25 可知：隨顆粒形狀因子增加，顆粒偏離流體流線的程度增大，Up-fn增大。圖26所示為不同顆粒形狀因子下沿相同顆粒編號軌跡運動的顆粒速度和無顆粒流體速度之差。由圖26可知：顆粒對流體的影響增強，沿同一顆粒軌跡的流體速度變化程度減小，Ufp-fn增大。

圖24 翼型升力系數隨顆粒形狀因子的變化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

經分析可得，隨著顆粒形狀因子增加，相同顆粒編號的顆粒遠離壓力表面。在沿同一編號顆粒軌跡運動的流體速度變化值相同的情況下，形狀因子大的顆粒升力系數(與無顆粒升力系數相比)降低程度更大，升力系數較小。而且隨著顆粒形狀因子增加，Ufp-fn增大，沿同一編號顆粒軌跡運動的流體速度變化值增大。故隨顆粒形狀因子增大，升力系數減小。當形狀因子為0.9、顆粒直徑為2 500 μm 時，升力系數減小幅度最大，與無顆粒時相比減少了0.18%。

3.3.2 顆粒形狀因子對翼型阻力系數的影響

圖27所示為翼型阻力系數隨形狀因子的變化。由圖27 可知：翼型阻力系數隨顆粒形狀因子的增大而增大；隨著顆粒形狀因子的增加，Ufp-fn增大，等效于相對流體速度的減少程度增大，導致阻力系數增大；當形狀因子為0.9，顆粒直徑為1 500 μm時，阻力系數變化最大，與無顆粒時相比增加了1.315%。

4 結論

1) 當顆粒直徑小于100 μm 時，顆粒慣性小，顆粒和流體可以作為單相，相當于增加了流體密度，升力系數大于無顆粒時的升力系數；當顆粒直徑大于100 μm 時，顆粒偏離流線，升力系數小于無顆粒時的升力系數。

2) 隨著顆粒直徑增大，顆粒慣性增強，顆粒偏離流線，升力系數在先減小后增大的趨勢中波動。

3) 隨著顆粒形狀因子增大，顆粒表面積減小，顆粒受流體阻力減小，偏離流線程度增大，升力系數減小，阻力系數增大。

4) 顆粒質量濃度對翼型升阻力系數影響最大，隨顆粒質量濃度增大，升力系數呈線性增大，阻力系數呈線性減小，當顆粒質量濃度為5.522 g/L時，相比于無顆粒時升力系數減小了0.478%，阻力系數增大了6.676%。