丁廣佳，及春寧，徐曉黎，許棟，呂迎雪

（1.中交天津港灣工程研究院有限公司，天津 300222；2.天津大學(xué)水利工程仿真與安全國家重點實驗室，天津 300072；3.中交第一航務(wù)工程局有限公司，天津 300456；4.中國交建海岸工程水動力重點實驗室，天津 300222）

1 研究背景

我國沿岸海流資源豐富，可為社會經(jīng)濟(jì)發(fā)展提供充足且廉價的清潔能源[1]。然而，我國沿岸海流流速普遍偏低（＜1.5 m/s），使得起動流速較高（＞2.0 m/s）的傳統(tǒng)海流能利用裝置（如渦輪機(jī)、水下風(fēng)車等）的應(yīng)用受到限制。新型VIVACE（Vortex Induced Vibration Aquatic Clean Energy）海流能利用裝置[2]利用柱體在海流作用下的渦激振動俘獲海流能，與波浪能利用裝置[3]相比，具有能量可利用時間長的特點，并且其起動流速低（0.25 m/s）、能量密度高、結(jié)構(gòu)簡單可靠、維護(hù)費(fèi)用低[2]，具有十分突出的經(jīng)濟(jì)效益和社會價值，符合我國“雙碳目標(biāo)”的發(fā)展戰(zhàn)略。

目前，國內(nèi)外多位學(xué)者針對VIVACE的能量利用效率優(yōu)化開展了許多富有成效的研究工作[4-6]，并應(yīng)用被動流動控制[7-9]和主動流動控制[10-12]顯著提高了VIVACE的振幅和流速響應(yīng)區(qū)間。與被動流動控制相比，主動流動控制可根據(jù)裝置所處的實時水動力環(huán)境動態(tài)調(diào)整控制策略，使得VIVACE海流能利用裝置時刻處于最優(yōu)俘能狀態(tài)。及春寧等[13]應(yīng)用主動流動控制，提出了一種附加旋轉(zhuǎn)圓柱的渦激振動發(fā)電裝置，如圖1所示。該裝置可主動調(diào)整附加圓柱的轉(zhuǎn)動，控制主圓柱周圍的流態(tài)，以增大圓柱的振幅和流速響應(yīng)區(qū)間，提高一級能量利用效率。然而，實際海洋環(huán)境惡劣多變，在優(yōu)化其工作性能的同時，也必須考慮其生存性問題。為此，本文針對附加旋轉(zhuǎn)圓柱的渦激振動發(fā)電裝置的水動力特性進(jìn)行研究，為其合理設(shè)計提供必要的理論指導(dǎo)。

圖1 附加旋轉(zhuǎn)圓柱的主動流動控制渦激振動發(fā)電裝置Fig.1 Vortex-induced vibration energy harvester with active flow control by using attached spinning cylinders

本文對及春寧等[13]提出的發(fā)電裝置進(jìn)行二維簡化，采用浸入邊界法模擬帶有兩個附屬旋轉(zhuǎn)圓柱的振動系統(tǒng)的渦激振動響應(yīng)，在俘能阻尼比ζharness，附加圓柱無量綱轉(zhuǎn)速α和方位角θ構(gòu)成的多參數(shù)空間內(nèi)，分析系統(tǒng)的受力和振動幅值。

2 數(shù)值模擬方法

采用嵌入式迭代的浸入邊界法[14]對流固耦合問題進(jìn)行數(shù)值模擬，控制方程如下：

式中：u為速度矢量；t為時間；p為壓強(qiáng)；ρ為流體密度；ν為運(yùn)動黏滯系數(shù)；為梯度算子；上標(biāo)T為矩陣轉(zhuǎn)置；f為附加體積力矢量，代表流固耦合邊界條件。本文針對層流條件下附加旋轉(zhuǎn)圓柱的渦激振動發(fā)電裝置水動力特性開展研究，因此控制方程中無紊流模型。

對以上控制方程采用二階的Admas-Bashforth時間格式進(jìn)行離散，可得守恒形式如下：

式中：I和D為插值函數(shù)；Vn+1為物面邊界點的速度矢量。

對做橫流向振動的剛性圓柱，其運(yùn)動方程可表示為：

式中：y為圓柱的橫流向振動位移；m為整個振動系統(tǒng)（包含主圓柱和附加圓柱）的質(zhì)量；c為系統(tǒng)阻尼，包含發(fā)電機(jī)俘能阻尼charness和結(jié)構(gòu)阻尼cstructure；k為彈簧剛度系數(shù)；Fy為系統(tǒng)受到的橫向流體力。方程采用標(biāo)準(zhǔn)的Newmark-β法求解。

3 問題描述和模型驗證

如圖2所示，主圓柱的直徑為D，兩個對稱布置的附加圓柱（C1和C2）的直徑均為d，方位角為θ（主圓柱與附加圓柱中心連線與水平方向的夾角），主圓柱與附加圓柱的間隙為g，附加圓柱的無量綱轉(zhuǎn)動速度為α=Usurf/U∞，其中Usurf為附加圓柱表面線速度，U∞為均勻來流流速。規(guī)定：圖2中附加圓柱的旋轉(zhuǎn)方向（內(nèi)向轉(zhuǎn)動）為正，反之（外向轉(zhuǎn)動）為負(fù)。俘能阻尼比ζharness=charness/（4πm fn），其中fn為系統(tǒng)固有頻率。為了獲得更高的能量利用效率，本文取cstructure=0?；谥鲌A柱直徑的雷諾數(shù)為Re=U∞D(zhuǎn)/ν=100，質(zhì)量比為m*=m/mf=2.0，其中mf為系統(tǒng)等體積流體質(zhì)量。

計算域大小為100D×100D，如圖2所示。為保證數(shù)值精度，在圓柱周圍采用加密網(wǎng)格，加密區(qū)域為8D×8D，無量綱網(wǎng)格尺寸為Δx/D=Δy/D=1/64。邊界條件設(shè)置：入口為Dirichlet型邊界（u=U∞，v=0），出口為Neumann型邊界（?u/?x=0，?v/?x=0），上下為可滑移邊界（?u/?y=0，v=0）。

圖2 計算域、邊界條件和圓柱布置Fig.2 Computational domain,boundary conditions and cylinder configuration

為了驗證數(shù)值方法的正確性，本文開展了帶附加圓柱的單圓柱渦激振動數(shù)值模擬，并與Mittal[15]的數(shù)值結(jié)果對比。數(shù)值模擬參數(shù)：雷諾數(shù)為Re=100，附加圓柱與主圓柱的直徑比為d/D=1/20，方位角θ=60°，間隙比g/D=0.075，無量綱轉(zhuǎn)速α=2.0。如表1所示，本文結(jié)果與Mittal[15]的計算結(jié)果吻合較好，兩者差別小于等于5%。

表1 主圓柱的升阻力系數(shù)對比Table 1 Comparison of lift and drag coefficients of the main cylinder

4 結(jié)果與討論

數(shù)值模擬中，保持直徑比（d/D=0.125）、間隙比（g/D=0.125）和折合流速（Ur=U∞/fnD=6）不變，通過改變俘能阻尼比ζharness、附加圓柱的方位角θ和無量綱轉(zhuǎn)速α，共開展432組工況的數(shù)值模擬（見表2），系統(tǒng)分析各參量對系統(tǒng)水動力性能的影響。

表2 數(shù)值模擬參數(shù)取值Table 2 Values of numerical simulation parameters

4.1 升力系數(shù)

圖3 給出了不同俘能阻尼比ζharness條件下，升力均方根值CL，rms在α-θ參數(shù)空間內(nèi)的分布?？偟膩砜?，隨著無量綱轉(zhuǎn)速α的增大和方位角θ的減小，CL，rms呈現(xiàn)增大的趨勢。并且隨著ζharness的增加，CL，rms較大的區(qū)域不斷縮小，等值線變稀疏，說明ζharness越大，α和θ對于CL，rms的影響則略有減弱。此外，CL，rms較大的區(qū)域均出現(xiàn)在α＞0的情況下，說明在附加圓柱內(nèi)向轉(zhuǎn)動的模式下，振動系統(tǒng)受到的升力較外向轉(zhuǎn)動的更大。這與內(nèi)向轉(zhuǎn)動的模式下流動更容易在控制圓柱表面發(fā)生流動分離有關(guān)，即：離開附加圓柱的自由剪切層卷起形成的漩渦經(jīng)過主圓柱的側(cè)面，導(dǎo)致主圓柱受到的升力增大。而在外向轉(zhuǎn)動的模式下，從附加圓柱上分離的剪切層重新附著在主圓柱的表面，并在主圓柱的背流側(cè)發(fā)生二次流動分離，導(dǎo)致主圓柱受到的升力較小。另外，內(nèi)向轉(zhuǎn)動的附加圓柱表面的線速度與來流速度方向相反，流動剪切率較大，生成的漩渦渦量較強(qiáng)，導(dǎo)致主圓柱的升力增大。

圖3 升力均方根值在α-θ參數(shù)空間內(nèi)的分布Fig.3 Distribution of root-mean-square lift coefficient inparameter space

圖3 中方框內(nèi)數(shù)字為CL，rms的最大值和最小值。隨著ζharness的增加，CL，rms的最大值和最小值均不斷減小。升力均方根最大值CL，rms=1.22出現(xiàn)在阻尼比ζharness=0.1、α=1.0和θ=50°的工況下。升力均方根最小值CL，rms=0.298出現(xiàn)在阻尼比ζharness=0.4、α=-1.0和θ=50°的工況下。

4.2 阻力系數(shù)

圖4給出了不同俘能阻尼比ζharness條件下，阻力均值CD，mean在α-θ參數(shù)空間內(nèi)的分布。整體看來，隨著ζharness的增加，CD，mean呈現(xiàn)下降趨勢，最大值阻力均值從阻尼比ζharness=0.1時的CD，mean=3.803下降到阻尼比ζharness=0.4的CD，mean=3.127。最大阻力均值出現(xiàn)的位置也隨著ζharness的增加逐漸由參數(shù)空間的右側(cè)中部（α=1.0，θ=70°）向參數(shù)空間的右下角（α=1.0，θ=55°）?？梢钥闯?，不同俘能阻尼比條件下，CD，mean最大值均出現(xiàn)在α=1.0條件下，這與附加圓柱內(nèi)向轉(zhuǎn)動引起的較寬的尾流有關(guān)。尾流較寬時，主圓柱的基底壓強(qiáng)較低，導(dǎo)致主圓柱受到的阻力較高。

圖4 阻力均值在α-θ參數(shù)空間內(nèi)的分布Fig.4 Distribution of mean drag coefficient in parameter space

圖5給出了不同俘能阻尼比ζharness條件下，阻力均方根CD，rms在α-θ參數(shù)空間內(nèi)的分布?？偟膩碚f，隨著ζharness的增加，CD，rms顯著減小。ζharness=0.1時，最大均方根值CD，rms=1.044出現(xiàn)在（α=-1.0，θ=75°）工況下。ζharness=0.4時，最大均方根值CD，rms=0.553出現(xiàn)在（α=1.0，θ=85°）工況下。隨著ζharness的增加，CD，rms在α-θ參數(shù)空間內(nèi)的分布趨于更加均勻。這說明，高俘能阻尼比條件下，阻力均方根值受附加圓柱的轉(zhuǎn)速和方位角的影響明顯減弱。

圖5 阻力均方根在α-θ參數(shù)空間內(nèi)的分布Fig.5 Distribution of root-mean-square drag coefficient in parameter space

4.3 無量綱振幅

圖6 給出了不同俘能阻尼比ζharness條件下，振幅最大值Ymax/D在α-θ參數(shù)空間內(nèi)的分布?？傮w上，隨著ζharness的增加，Ymax/D不斷減小。ζharness=0.1時，振幅最大值Ymax/D=0.78出現(xiàn)在（α=0.25，θ=55°）工況下。ζharness=0.4時，振幅最大值Ymax/D=0.472出現(xiàn)在（α=1.0，θ=55°）工況下。在α-θ參數(shù)平面內(nèi)，大振幅主要分布在方位角較小、旋轉(zhuǎn)速度大于0的范圍內(nèi)。當(dāng)方位角較大（θ＞70°）時，Ymax/D的等值線大致呈水平分布。這說明，方位角較大時，最大振幅受附加圓柱的轉(zhuǎn)速影響較小。

圖6 無量綱振幅最大值在α-θ參數(shù)空間內(nèi)的分布Fig.6 Distribution of maximum non-dimensional vibration amplitude in parameter space

4.4 組合參數(shù)（α/θ）對無量綱振幅的影響

圖7 給出了無量綱振幅Ymax/D隨組合參數(shù)α/θ的變化趨勢?？梢钥闯?，不同俘能阻尼比ζharness條件下，Ymax/D隨著α/θ的增大均呈現(xiàn)先增后減的變化趨勢。隨著ζharness的增大，Ymax/D明顯降低。通過多項式擬合得到變化趨勢的擬合公式為：

圖7 無量綱振幅最大值隨組合參數(shù)α/θ的變化Fig.7 Variation of maximum non-dimensional vibration amplitude with combined parameter α/θ

擬合公式表現(xiàn)為二次曲線形式。二次項系數(shù)表示曲線開口的方向和大小，即曲線變化的快慢；一次項系數(shù)表示曲線頂點位置的偏移，正值向右偏，負(fù)值向左偏；常數(shù)項表示曲線的截距。從不同ζharness條件下擬合公式系數(shù)的變化規(guī)律來看，二次項系數(shù)均為負(fù)值，隨著ζharness增大而略有減??；一次項系數(shù)均為正值，說明最大的Ymax/D均在α/θ＞0（附加圓柱向內(nèi)旋轉(zhuǎn)）的條件下取得，這與圖6中Ymax/D的變化趨勢相同；常數(shù)項隨ζharness增大而減小。二次項系數(shù)很小，說明擬合曲線隨組合系數(shù)α/θ變化較為平緩。相反，常數(shù)項隨ζharness的變化幅度較大，說明最大振幅受俘能阻尼比的影響最為顯著，這與圖6的結(jié)果一致。

5 結(jié)語

本文應(yīng)用嵌入式迭代浸入邊界法，數(shù)值模擬了主動流動控制渦激振動發(fā)電裝置的水動力特性隨俘能阻尼比ζharness、附加圓柱的方位角θ和無量綱轉(zhuǎn)速α的變化規(guī)律，分析了升力系數(shù)均方根值CL，rms、阻力系數(shù)均值CD，mean和均方根值CD，rms、最大無量綱振幅Ymax/D在α-θ參數(shù)空間內(nèi)的分布，給出了最大無量綱振幅隨組合參數(shù)α/θ變化的擬合公式。具體結(jié)論如下：

1）整體來看，俘能阻尼比ζharness對于升力系數(shù)均方根值CL，rms的影響較小，而對阻力系數(shù)均值CD，mean和均方根值CD，rms及最大無量綱振幅Ymax/D則有明顯的影響。

2）增大附加圓柱轉(zhuǎn)速α?xí)?dǎo)致升力均方根值CL，rms、阻力均值CD，mean明顯增大，而最大無量綱振幅Ymax/D變化較小；增大附加圓柱方位角，升、阻力系數(shù)和振幅均呈現(xiàn)先增后減的趨勢。

3）最大無量綱振幅Ymax/D與組合參數(shù)α/θ呈現(xiàn)二次函數(shù)關(guān)系，二次項系數(shù)均為負(fù)值，隨著ζharness增大而略有減小，一次項系數(shù)均為正值，隨著ζharness增大而略有增加，常數(shù)項隨ζharness的增大而明顯減小。