基于熵產(chǎn)理論的三維渦激振動(dòng)低速水流能俘能分析

收稿日期：2022-06-25

基金項(xiàng)目：國(guó)家自然科學(xué)基金（520609010；51969009）

通信作者：羅竹梅（1979—），女，博士、教授、博士生導(dǎo)師，主要從事流固耦合及工程應(yīng)用理論方面的研究。112374102@qq.com

DOI：10.19912/j.0254-0096.tynxb.2022-0932 文章編號(hào)：0254-0096（2023）06-0045-08

摘 要：以剛性圓柱作為俘能結(jié)構(gòu)，通過(guò)強(qiáng)流固耦合數(shù)值方法模擬俘能結(jié)構(gòu)在低速水流環(huán)境下的三維渦激振動(dòng)，對(duì)共振約化速度（[Ur=2～12]）范圍內(nèi)進(jìn)行俘能結(jié)構(gòu)的俘能大小及俘能效率分析。為獲得俘能結(jié)構(gòu)三維尾流對(duì)俘能大小的影響，采用熵產(chǎn)理論并結(jié)合尾流特性以捕捉俘能結(jié)構(gòu)尾渦的能量損失來(lái)源及分布，獲得不同振動(dòng)分支下俘能大小與尾流耗散間的關(guān)系。結(jié)果表明，基于熵產(chǎn)理論的能量損失分析可準(zhǔn)確捕捉到俘能結(jié)構(gòu)的能量耗散趨勢(shì)。黏性熵產(chǎn)主要發(fā)生在圓柱表面，湍流熵發(fā)生在耗散率較大的尾流區(qū)。隨著約化速度的增加，上端分支內(nèi)，當(dāng)[Ur=6]時(shí)俘能結(jié)構(gòu)的俘能效率最大達(dá)到30.5%，俘能結(jié)構(gòu)俘獲的水動(dòng)能增加，圓柱表面的黏性熵?fù)p耗和圓柱周?chē)拔擦鲄^(qū)湍流熵產(chǎn)損失也增加，俘能結(jié)構(gòu)表面的黏性熵產(chǎn)損失和尾流區(qū)湍流熵產(chǎn)損失增大成為俘能效率下降的主要原因。

關(guān)鍵詞：海洋能；俘能結(jié)構(gòu)；渦激振動(dòng)；流固耦合；能量耗散；熵產(chǎn)理論

中圖分類(lèi)號(hào)：TK79；TK730.2 文獻(xiàn)標(biāo)志碼：A

0 引 言

實(shí)現(xiàn)碳達(dá)峰、碳中和目標(biāo)愿景，發(fā)展清潔可再生能源是重中之重。作為可再生能源之一的海洋能，自20世紀(jì)70年代開(kāi)始，其開(kāi)發(fā)就備受各沿海國(guó)家的重視。在海流流速穩(wěn)定的區(qū)域，可獲得穩(wěn)定可靠的海流能。目前人們對(duì)海流能的利用主要集中在潮流能的開(kāi)發(fā)利用方面，即通過(guò)利用水下的渦輪發(fā)電設(shè)備可從有較大流速的潮流中獲取能量，但基于渦激振動(dòng)從低速海流或河流中獲取能量方面的研究較少［1］。由于水下渦輪機(jī)受水流速度的限制，在低速水流中無(wú)法實(shí)現(xiàn)有效能量轉(zhuǎn)換，基于渦激振動(dòng)驅(qū)動(dòng)俘獲低速海流能或河流能成為清潔能源利用的一種新選擇［2］。

利用渦激振動(dòng)從低速海流中俘獲能量的研究目前還停留在理論和實(shí)驗(yàn)階段。相關(guān)文獻(xiàn)主要來(lái)源于Michael Bernitsas教授研究團(tuán)隊(duì)，該團(tuán)隊(duì)對(duì)從低速海流中俘獲海流能的能量轉(zhuǎn)換裝置VIVACE進(jìn)行了一系列實(shí)驗(yàn)。Garcia［3］根據(jù)能量現(xiàn)象學(xué)思想建立出一個(gè)數(shù)學(xué)方程模型，將由Bernitsas團(tuán)隊(duì)對(duì)單個(gè)光滑的圓柱進(jìn)行的渦激振動(dòng)的實(shí)驗(yàn)數(shù)據(jù)代入方程，得到方程中的參數(shù)值。該數(shù)學(xué)模型還捕捉到渦激振動(dòng)響應(yīng)的初始分支、上端分支和下端分支，但未捕捉到符合實(shí)驗(yàn)結(jié)果的振幅，且模型中較小的能量變化都會(huì)對(duì)振幅產(chǎn)生較大影響。Bernitsas等［4］和Lee等［5］在高雷諾數(shù)和高阻尼條件下進(jìn)行渦激振動(dòng)實(shí)驗(yàn)，通過(guò)優(yōu)化俘能裝置的阻尼比，得到最大振幅比為1.78和最大能量轉(zhuǎn)換效率為30.8%以及平均能量轉(zhuǎn)換效率為22%。

利用數(shù)值方法進(jìn)行渦激振動(dòng)俘能研究是Bernitsas團(tuán)隊(duì)的另一種研究手段，能輔助獲取俘能過(guò)程中流場(chǎng)變化信息與俘能參數(shù)和俘能大小間的關(guān)系，便于進(jìn)行試驗(yàn)條件下難以實(shí)現(xiàn)的高雷諾數(shù)下的渦激振動(dòng)響應(yīng)結(jié)果的研究。孫海等［6］采用二維非定常雷諾平均法（RANS），通過(guò)被動(dòng)湍流控制對(duì)粗糙表面的單個(gè)彈性支撐圓柱的馳振現(xiàn)象進(jìn)行模擬，得到與實(shí)驗(yàn)基本一致的響應(yīng)結(jié)果及渦脫落模式。丁琳等［7］采用OpenFOAM數(shù)值求解了兩個(gè)彈性支撐圓柱的橫向流激振動(dòng)，采用被動(dòng)湍流控制觀(guān)察到圓柱的馳振現(xiàn)象。振幅響應(yīng)除了初始分支、上端分支外，還有馳振分支。丁文俊等［8］采用相同的數(shù)值方法通過(guò)被動(dòng)湍流控制方法研究了高雷諾數(shù)下的流激振動(dòng)響應(yīng)，獲得較高的馳振振幅響應(yīng)。

僅有極少數(shù)的國(guó)外學(xué)者已對(duì)三維渦激振動(dòng)進(jìn)行系統(tǒng)分析研究，Kondo［9］采用Petrov-Galerkin有限元法來(lái)求解三維的Navier-Stokes方程，模擬三維圓柱體的雙自由度渦激振動(dòng)。Mittal［10］通過(guò)對(duì)質(zhì)量比[m*=10]和雷諾數(shù)[Re=1000]的圓柱體模型進(jìn)行三維數(shù)值模擬，重點(diǎn)研究約化速度[Ur]為4.2～6.2時(shí)上下分支間的滯后過(guò)渡。目前，對(duì)圓柱體中的渦激振動(dòng)效應(yīng)進(jìn)行了大量系統(tǒng)的動(dòng)力學(xué)實(shí)驗(yàn)分析和計(jì)算數(shù)值研究，但大多數(shù)研究集中在抑制二維圓柱渦激振動(dòng)響應(yīng)、減少工程破壞方面，對(duì)三維渦激振動(dòng)俘能研究較少。

能量轉(zhuǎn)換過(guò)程中，不可避免地存在耗散效應(yīng)，局部熵產(chǎn)率（local entropy production rate，LEPR）正受到越來(lái)越多的關(guān)注。在整個(gè)流動(dòng)區(qū)域上通過(guò)LEPR積分得到的產(chǎn)生熵與流動(dòng)損失密切相關(guān)。Kock等［11］對(duì)管道流動(dòng)進(jìn)行分析，指出熵產(chǎn)可作為一個(gè)有效的指標(biāo)，以確定傳熱的湍流過(guò)程的效率。此外，Herwig等［12］利用管道層流和紊流的實(shí)驗(yàn)數(shù)據(jù)，驗(yàn)證了熵產(chǎn)理論的正確性。通過(guò)對(duì)3種不同的二維粗糙度類(lèi)型進(jìn)行模擬，計(jì)算任意形狀的壁面粗糙度所要求的等效沙粒粗糙度。在確定產(chǎn)生熵的基礎(chǔ)上，驗(yàn)證了耗散模型適用于層流和湍流流動(dòng)。

熵產(chǎn)理論普遍用于可壓縮流體中，主要用于評(píng)估熱力傳動(dòng)的能量損失，而在不可壓縮過(guò)流部件中的應(yīng)用非常少。針對(duì)渦激振動(dòng)現(xiàn)象，大多數(shù)研究者僅局限于二維渦激振動(dòng)的水動(dòng)力特性研究，真實(shí)低速海流環(huán)境下的三維渦激振動(dòng)俘能研究較少，且基于熵產(chǎn)理論對(duì)渦激振動(dòng)進(jìn)行能量捕捉及能量損失分析目前較少研究者涉及。為了研究渦激振動(dòng)三維尾流效應(yīng)對(duì)俘能大小的影響、準(zhǔn)確捕捉圓柱四周及尾渦的能量變化規(guī)律，本研究采用流固耦合數(shù)值方法對(duì)均勻低速來(lái)流下的俘能結(jié)構(gòu)進(jìn)行三維渦激振動(dòng)模擬，結(jié)合熵產(chǎn)理論，從能量損失的角度分析不同振動(dòng)分支下俘能大小與尾流特性間的關(guān)系，尋找各分支下影響俘能能量耗散的主要原因。

1 渦激振動(dòng)俘能模型

1.1 幾何模型

低速水流中渦激振動(dòng)驅(qū)動(dòng)發(fā)電的俘能結(jié)構(gòu)如圖1所示。整個(gè)俘能結(jié)構(gòu)為單圓柱構(gòu)成的質(zhì)量-彈簧-阻尼振動(dòng)系統(tǒng)，結(jié)構(gòu)模型中振動(dòng)圓柱直徑[D=0.1] m，長(zhǎng)[L=0.5] m，質(zhì)量比[m*=2.4]，質(zhì)量阻尼比[m*ζ=0.013]。將圓柱置于均勻來(lái)流中，其在水中的固有頻率[fn=1.46 ]Hz，為模擬真實(shí)低速海流或河流環(huán)境，實(shí)驗(yàn)流速[U=0.28～1.68] m/s，約化速度[Ur=2～12]。

1.2 數(shù)值計(jì)算模型

1.2.1 流體控制方程

均勻來(lái)流作用條件下俘能結(jié)構(gòu)的渦激振動(dòng)屬流固耦合問(wèn)題，選取RANS方程作為流體動(dòng)力學(xué)控制方程，結(jié)合SST [k]-ω湍流模型對(duì)流場(chǎng)進(jìn)行基本理論方法求解，采用有限體積法對(duì)控制方程進(jìn)行離散，其流體控制方程為：

[?ui?xi=0] （1）

[?ρui?t+?ρuiuj?xj=-?p?xi+μ▽2ui-?ρui′uj′?xj] （2）

[-ρui′uj′=μt?ui?xj+?uj?xi-23ρktδij] （3）

式中：[ui]和[ui′]——對(duì)應(yīng)的[i]方向速度分量和速度變化；[xi]——[x]方向的笛卡爾坐標(biāo)系；[ρ]和[μ]——流體的密度和動(dòng)力黏度；[t]——時(shí)間；[p]——壓力；[-ρui′uj′]——雷諾應(yīng)力；[μt]——湍流渦流黏度；[kt]——湍動(dòng)能，從湍流振蕩過(guò)程產(chǎn)生的速度場(chǎng)中獲得；[δij]——Kronecker函數(shù)。

SST [k]-ω湍流模型當(dāng)中包括湍動(dòng)能[k]和湍流耗散率[ω]的輸運(yùn)方程。其中湍動(dòng)能[k]的輸運(yùn)方程為：

[?ρk?t+?ρujk?xj=P-β*ρωk+??xjμ+σkμi?k?xj] （4）

湍流耗散率ω的輸運(yùn)方程為：

[?ρω?t+?ρujω?xj=γvtP-βρω2+" " " " " " ??xjμ+σωμt?ω?xj+21-F1ρσωω·?k?xj·?ω?xj] （5）

式中：[P]——雷諾應(yīng)力和平均運(yùn)動(dòng)變形率張量的二重標(biāo)量積；[β*]——熱膨脹系數(shù)；[σk]——湍動(dòng)能Prandtl數(shù)；[γ、][β、][σω]——經(jīng)驗(yàn)系數(shù)；[vt]——渦黏性系數(shù)；[F1]——混合函數(shù)，目的是使輸運(yùn)方程從近壁面的[k]-ω模型過(guò)渡到外部的[k-ε]模型，其定義為：

[?=F1?1+1-F1?2] （6）

式中：[?1]——近壁面的[k]-ω模型，[γ1=β1/β*-][σω1κ2β*]；[?2]——外部的[k-ε]模型，[γ2=β2/β*][-σω2κ2β*]。其中，[σω1、][σω2、][β1、][β2]——在近壁面[k]-ω模型和外部[k-ε]模型中求解[k]和ω方程時(shí)的經(jīng)驗(yàn)系數(shù)；[κ]——Von Karman常數(shù)。

1.2.2 固體控制方程

基于數(shù)值仿真方法，對(duì)自由表面下的圓柱振動(dòng)方程進(jìn)行離散求解。固體運(yùn)動(dòng)方程采用四階龍格庫(kù)塔法，基于用戶(hù)自定義程序UDF求解每個(gè)時(shí)間步長(zhǎng)的圓柱運(yùn)動(dòng)方程，其計(jì)算結(jié)果傳遞給流體方程，實(shí)現(xiàn)雙向流固耦合計(jì)算，該數(shù)值方法不僅計(jì)算效率較高，計(jì)算精度也較為滿(mǎn)意。四階龍格庫(kù)塔法求解圓柱運(yùn)動(dòng)方程如式（7）所示。

[fv，y=FfluidtM-2ξω0yt-ω02yt] （7）

采用四階龍格庫(kù)塔法將方程離散：

[vtn+1=vtn+Δt6×K1+2K2+2K3+K4] （8）

[ytn+1=ytn+Δt×vtn+Δt26K1+K2+K3] （9）

式中：[Δt]——時(shí)間步長(zhǎng)；[K1]、[K2]、[K3]、[K4]——四階龍格庫(kù)塔的轉(zhuǎn)換因子。

如果確定了[tn]時(shí)刻的柱體結(jié)構(gòu)的位移[y（tn）]和速度[v（tn）]，可通過(guò)式（7）～式（9）來(lái)求解[tn+1]時(shí)刻的位移[y（tn+1）]。

1.3 計(jì)算域及網(wǎng)格劃分

為進(jìn)一步實(shí)現(xiàn)在真實(shí)來(lái)流條件作用下俘能結(jié)構(gòu)的三維渦激振動(dòng)，建立出如圖2所示的流固耦合的計(jì)算模型。模型中俘能結(jié)構(gòu)為剛性單圓柱質(zhì)量-彈簧-阻尼系統(tǒng)，圓柱直徑為[D]，流體域長(zhǎng)為[20D]，高為[5D]，寬為[10D]。圓柱中心距離上游進(jìn)口面為[5D]，距離出口面為[15D]。置于流體域中的俘能結(jié)構(gòu)在均勻來(lái)流激勵(lì)下，可捕捉低速水流環(huán)境下自由渦激振動(dòng)的精細(xì)流場(chǎng)信息。

采用四邊形網(wǎng)格對(duì)流體的計(jì)算域進(jìn)行網(wǎng)格分布的精確劃分，沿水流流動(dòng)方向，網(wǎng)格呈非均勻分布，為獲得圓柱表面及尾流精細(xì)的流場(chǎng)信息，對(duì)圓柱周?chē)饔虻木W(wǎng)格進(jìn)行加密，如圖3所示。流體域中的上下兩個(gè)邊界均為固壁界面，圓柱體與流體的交接界面處設(shè)置為流固耦合（fluid-structure-interaction，F(xiàn)SI）邊界。結(jié)構(gòu)表面采用無(wú)滑移邊界條件，即圓柱體表面的流體速度與圓柱振動(dòng)速度相同，入口速度邊界設(shè)置為自由來(lái)流速度。下游邊界條件設(shè)定為壓力出口，參考?jí)毫榱恪?/p>

俘能結(jié)構(gòu)渦激振動(dòng)產(chǎn)生較大位移，結(jié)構(gòu)的運(yùn)動(dòng)必然引起網(wǎng)格變形，本論文選取彈簧光順模型中的擴(kuò)散光順?lè)ǎ╠iffusion），通過(guò)求解式（10）的擴(kuò)散方程來(lái)設(shè)置網(wǎng)格節(jié)點(diǎn)：

[▽?duì)扫寀=0] （10）

式中：[ι]——擴(kuò)散系數(shù)；[u]——網(wǎng)格的運(yùn)動(dòng)速度。

通過(guò)UDF自定義程序?qū)⑦\(yùn)動(dòng)圓柱設(shè)為剛性體，在流場(chǎng)中可以進(jìn)行整體運(yùn)動(dòng)，而不會(huì)發(fā)生彎曲扭轉(zhuǎn)和局部形變。

1.4 網(wǎng)格無(wú)關(guān)性驗(yàn)證

為量化數(shù)值計(jì)算結(jié)果對(duì)網(wǎng)格密度的依賴(lài)性，對(duì)約化速度[Ur=5]時(shí)4種不同網(wǎng)格密度下俘能結(jié)構(gòu)的渦激振動(dòng)進(jìn)行數(shù)值計(jì)算，結(jié)果見(jiàn)表1。由表1可看出，隨著圓柱表面徑向網(wǎng)格尺寸[Δr]的減少和總節(jié)點(diǎn)數(shù)[Nnode]的增加，振幅比[A*]和頻率比[f/fn]（[fn]為結(jié)構(gòu)在水中的固有頻率）結(jié)果逐漸趨于穩(wěn)定，采用第3種網(wǎng)格劃分方式既能滿(mǎn)足計(jì)算精度的要求，又能提高計(jì)算效率。

2 熵產(chǎn)理論

使用傳統(tǒng)的壓差評(píng)估方法無(wú)法準(zhǔn)確確定能量耗散的主要位置。圓柱的渦激振動(dòng)過(guò)程中，由于工作介質(zhì)中的黏度和產(chǎn)生的雷諾應(yīng)力，機(jī)械能被不可逆地轉(zhuǎn)化為了內(nèi)能?；跓崃W(xué)第二定律，將熵產(chǎn)分析結(jié)合流體力學(xué)和傳熱理論來(lái)計(jì)算實(shí)際的能量損失，可準(zhǔn)確測(cè)量和捕捉能量損失?？紤]到水的比熱容較高，三維圓柱渦激振動(dòng)的數(shù)值計(jì)算中，認(rèn)為溫度為恒定，不考慮傳熱引起的熵產(chǎn)。在忽略熱傳導(dǎo)的情況下，邊界層中的黏性力將不可逆地轉(zhuǎn)化熱傳導(dǎo)流體的動(dòng)能和勢(shì)能為內(nèi)能，從而引起熵的增加。

熵產(chǎn)主要包括黏度熵產(chǎn)、湍流熵產(chǎn)和溫度熵產(chǎn)。由于熵產(chǎn)是不可壓縮流體的狀態(tài)參數(shù)，引入熵產(chǎn)率來(lái)表示不可壓縮流體的輸運(yùn)方程，如式（11）所示。

[ρ?s?t+u?s?x+v?s?y+w?s?z=-divqT+ΦT+ΦΘT2] （11）

式中：[s]——熵產(chǎn)；[u、v]和[w]——3個(gè)方向（[x、y]和[z]）的速度；[q]——熱流密度；[Φ]——黏性耗散產(chǎn)生的熵；[ΦΘ]——溫度變化產(chǎn)生的熵產(chǎn)。

不可壓縮流體的黏性耗散可通過(guò)式（12）計(jì)算。

[Φ=μ2?u?x2+?v?y2+?w?y2+" " " "?u?y+?v?x2+?u?z+?w?x2+?v?z+?w?y2] （12）

湍流中，局部熵產(chǎn)率（LEPR）由以下兩部分組成，其中一部分是時(shí)間平均運(yùn)動(dòng)所導(dǎo)致的熵產(chǎn)率即[S?D]，而另一部分是由速度波動(dòng)導(dǎo)致的熵產(chǎn)率即[S?D′]?？傡禺a(chǎn)率計(jì)算公式為：

[S?D=S?D+S?D′] （13）

[S?D=2μeffT?u1?x12+?u2?x22+?u3?x32+μeffT" " " " " ?u2?x1+?u1?x22+?u3?x1+?u1?x32+?u2?x3+?u3?x22] （14）

[S?D′=2μeffT?u1′?x12+?u2′?x22+?u3′?x32+μeffT" " " " " ?u2′?x1+?u1′?x22+?u3′?x1+?u1′?x32+?u2′?x3+?u3′?x22] （15）

式中：[T]——溫度；[ui]——時(shí)間平均速度；[i]——笛卡爾坐標(biāo)系的3個(gè)方向；有效動(dòng)態(tài)黏度用[μeff]表示，計(jì)算公式為：

[μeff=μ+μt] （16）

然而，在雷諾平均Navier-Stokes方法中，速度方向的[u1′，][u2′，][u3′]的分量不是直接計(jì)算得到的，而是通過(guò)求解[k]和ω方程間接得到的。根據(jù)Kock和Mathieu在[k-ω]湍流模式下的模擬計(jì)算，由速度波動(dòng)引起的局部產(chǎn)生熵可近似計(jì)算如式（17）所示：

[S?D′=βρωkT] （17）

式中：[β]=0.09，湍流耗散速率[ε=βωk]，可對(duì)尾渦處的湍流耗散速率進(jìn)行分析。

因此，通過(guò)對(duì)整個(gè)流動(dòng)區(qū)域的積分，可得到總產(chǎn)生熵。表達(dá)方式如式（18）～式（20）所示：

[SD=VS?DdV] （18）

[SD′=VS?D′dV] （19）

[SD=VS?DdV=SD+SD′] （20）

式中：[SD]——由時(shí)間平均運(yùn)動(dòng)引起的產(chǎn)生熵；[V]——流體區(qū)域的體積；[SD′]——由速度波動(dòng)引起的產(chǎn)生熵；[SD]——全部產(chǎn)生熵。

3 結(jié)果分析

3.1 振動(dòng)響應(yīng)結(jié)果

圖4為約化速度[Ur=2～12]范圍內(nèi)俘能結(jié)構(gòu)渦激振動(dòng)響應(yīng)的振幅比[A*]和頻率比[f/fn]，采用與文獻(xiàn)［13］、文獻(xiàn)［14］數(shù)值實(shí)驗(yàn)、文獻(xiàn)［15］物理實(shí)驗(yàn)完全相同的低質(zhì)量比[m*=2.4]和低質(zhì)量阻尼比[m*ζ=0.013]，并將實(shí)驗(yàn)?zāi)M結(jié)果與文獻(xiàn)結(jié)果加以比較。

從圖4可看出，本文采用的數(shù)值方法在各約化速度下的振幅比[A*]與文獻(xiàn)［15］的物理實(shí)驗(yàn)結(jié)果基本吻合，但上端分支更寬，其中最大振幅比[A*]的值和物理實(shí)驗(yàn)結(jié)果非常接近。當(dāng)頻率比[f/fn]在1.1～1.4之間時(shí)，鎖定現(xiàn)象發(fā)生，俘能結(jié)構(gòu)振幅較大，對(duì)應(yīng)于響應(yīng)的上端分支。此外，本數(shù)值實(shí)驗(yàn)得到的鎖定范圍與文獻(xiàn)［14］、文獻(xiàn)［13］數(shù)值實(shí)驗(yàn)結(jié)果基本一致，但振幅更高。鎖定范圍[Ur=4～10]內(nèi)，頻率比[f/fn]幾乎重疊。

3.2 三維尾流特性分析

渦激振動(dòng)驅(qū)動(dòng)的俘能結(jié)構(gòu)其尾流特性與其俘能結(jié)果有著必然聯(lián)系。捕捉俘能結(jié)構(gòu)的三維尾流特性可為結(jié)構(gòu)俘能中的能量分布及耗散規(guī)律提供流場(chǎng)分析依據(jù)。

圖5為不同約化速度下俘能結(jié)構(gòu)振動(dòng)在振幅比最大位置時(shí)的三維尾渦脫落結(jié)構(gòu)，俘能裝置的尾渦結(jié)構(gòu)與來(lái)流速度和響應(yīng)幅值的相關(guān)性較大。[Ur=3]時(shí)，來(lái)流速度[U]僅為0.28 m/s，此時(shí)振幅響應(yīng)幅值較小，尾流區(qū)未出現(xiàn)明顯的脫落渦對(duì)，見(jiàn)圖5a。在渦激振動(dòng)響應(yīng)鎖定范圍[Ur=4～10]內(nèi)，隨著約化速度的增加，流體與俘能結(jié)構(gòu)的耦合作用增強(qiáng)，圓柱結(jié)構(gòu)后出現(xiàn)明顯的脫落渦對(duì)。此時(shí)渦動(dòng)力更大，響應(yīng)幅值較大，較大的響應(yīng)振幅導(dǎo)致出現(xiàn)更寬的尾流。由于響應(yīng)頻率[f]隨著約化速度增加而增加，尾流渦脫落頻率加快，尾流中脫落渦強(qiáng)度更強(qiáng)，相應(yīng)能量耗散也增加，如圖5b～圖5e所示。當(dāng)響應(yīng)進(jìn)入下端分支時(shí)，振動(dòng)頻率偏離俘能結(jié)構(gòu)的固有頻率，響應(yīng)偏離鎖定范圍，導(dǎo)致渦激振動(dòng)振幅減弱。在較大來(lái)流速度[U]和較低響應(yīng)振幅下，尾流區(qū)內(nèi)出現(xiàn)狹長(zhǎng)的正負(fù)渦對(duì)，尾流寬度變窄，見(jiàn)圖5f。

截取不同約化速度下x-z平面內(nèi)的流場(chǎng)渦量圖，以清楚獲取尾渦模式隨約化速度的變化關(guān)系，如圖6所示。[Ur=3]時(shí)，響應(yīng)處于初始分支，此時(shí)出現(xiàn)一個(gè)經(jīng)典的2S尾流脫落模式，即在單個(gè)振動(dòng)周期內(nèi)的尾流出現(xiàn)一對(duì)渦脫落，這與文獻(xiàn)［15］物理實(shí)驗(yàn)中初始分支的渦脫落模式一致。上端分支內(nèi)，尾渦脫落模式從初始分支內(nèi)的2S模式變?yōu)?P模式，單個(gè)周期內(nèi)圓柱后出現(xiàn)一個(gè)較大渦對(duì)和一個(gè)較小渦對(duì)，隨著尾流向下游的發(fā)展，較小渦對(duì)逐漸耗散，見(jiàn)圖6b。上端分支內(nèi)隨著來(lái)流速度的增加，尾渦脫落頻率增強(qiáng)，流體結(jié)構(gòu)相互作用持續(xù)增強(qiáng)，在[Ur=9]時(shí)尾流中出現(xiàn)[P+T]尾渦模式，單個(gè)周期內(nèi)更多的渦對(duì)脫落。在下端分支中，渦位發(fā)生錯(cuò)亂，展向渦排列非常不規(guī)則。尾渦模式變成P+S模式，圓柱上側(cè)脫離出兩個(gè)正渦，圓柱下側(cè)剛好脫落出一個(gè)完整的負(fù)渦。脫落的渦細(xì)長(zhǎng)，占據(jù)下游較長(zhǎng)的尾流區(qū)。

3.3 基于熵產(chǎn)理論的俘能分析

3.3.1 渦強(qiáng)度和耗散率分析

在圓柱渦激振動(dòng)的過(guò)程中，產(chǎn)生的尾渦一部分轉(zhuǎn)換成圓柱渦激振動(dòng)的升力，另一部分能量卻因?yàn)闇u與渦之間的影響產(chǎn)生湍流能量耗散［16］。為了探究鎖定區(qū)內(nèi)俘能圓柱尾渦耗散規(guī)律，將[Ur=5]時(shí)俘能結(jié)構(gòu)尾流區(qū)捕捉到的4個(gè)脫落正渦中心[b1、b2、b3、b4]設(shè)為監(jiān)測(cè)點(diǎn)，提取4個(gè)監(jiān)測(cè)點(diǎn)的渦強(qiáng)度，計(jì)算由速度波動(dòng)引起的局部產(chǎn)生熵[S?D′]中的[ε]，得到渦強(qiáng)度和尾渦耗散率[ε]隨時(shí)間變化的關(guān)系圖，如圖7所示。尾流渦主要由流向渦對(duì)控制，兩個(gè)相鄰渦對(duì)之間的距離約為1個(gè)圓柱直徑。從曲線(xiàn)圖中可清晰地看出，當(dāng)結(jié)構(gòu)產(chǎn)生較大振動(dòng)幅值時(shí)，不斷有較強(qiáng)的渦對(duì)從圓柱上下表面脫落。由于振動(dòng)位移隨時(shí)間發(fā)生周期性變化，第1對(duì)脫落的渦強(qiáng)度和耗散率隨時(shí)間也發(fā)生變化。當(dāng)脫落渦到監(jiān)測(cè)點(diǎn)[b1]處，此時(shí)渦強(qiáng)度和湍流耗散率分別為0.023和0.144，隨著脫落渦向下游不斷演變，渦強(qiáng)度不斷降低，湍流耗散率也越來(lái)越小，在[b4]點(diǎn)達(dá)到最低，此時(shí)渦強(qiáng)度和湍流耗散率分別降低了97%和93%。隨著時(shí)間的增加，越遠(yuǎn)離結(jié)構(gòu)的尾流區(qū)其尾渦強(qiáng)度和耗散率ε越小，渦耗散越快，流場(chǎng)受結(jié)構(gòu)的影響越小。

圖8顯示在4個(gè)約化速度（Ur=3、5、9、11）下第1個(gè)脫落正渦中心處的渦強(qiáng)度和耗散率三維曲線(xiàn)圖。結(jié)合圖4a和圖9可看出，渦強(qiáng)度及渦耗散率與不同約化速度下的振幅及俘能大小相關(guān)。當(dāng)俘能結(jié)構(gòu)處在較低振幅的初始分支（Ur=3）時(shí)，來(lái)流速度低，渦強(qiáng)度較小，為0.071，相應(yīng)的渦耗散率為0.063。隨著約化速度的增加，結(jié)構(gòu)受到的來(lái)流速度增加，在上端分支鎖定區(qū)間（Ur=5和Ur=9），流體結(jié)構(gòu)耦合作用持續(xù)增強(qiáng)，渦強(qiáng)度更大，俘獲的水動(dòng)能[PVIV]增加。由于結(jié)構(gòu)振幅較大，尾流區(qū)流場(chǎng)受結(jié)構(gòu)干擾影響較大，相對(duì)于初始分支，尾渦耗散率增加了10倍。響應(yīng)進(jìn)入下端分支（Ur=11）時(shí)，結(jié)構(gòu)振動(dòng)頻率偏離固有頻率，響應(yīng)進(jìn)入非鎖定區(qū)間，振幅降低，渦強(qiáng)度相應(yīng)也降低。此時(shí)大部分水動(dòng)能都耗散在尾渦中，結(jié)構(gòu)從流體中俘獲的能量[PVIV]較少。

3.3.2 俘能結(jié)果分析

圓柱在均勻來(lái)流中俘獲的水動(dòng)能采用式（21）［17］進(jìn)行計(jì)算。

[PVIV=1Tosc0ToscCtotaly2dt] （21）

式中：[Tosc]——單個(gè)周期振動(dòng)時(shí)間；[Ctotal]——俘能結(jié)構(gòu)阻尼。

計(jì)算俘能結(jié)構(gòu)俘獲的的水動(dòng)能[PVIV]和流體能量[Pfluid]的比值可獲得各約化速度下結(jié)構(gòu)的俘能效率[η]。對(duì)約化速度[Ur=2～12]范圍內(nèi)的俘能結(jié)構(gòu)渦激振動(dòng)俘獲的水動(dòng)能[PVIV]和俘能效率[η]進(jìn)行計(jì)算，得到如圖9所示的俘能結(jié)果和能量轉(zhuǎn)換效率隨約化速度變化關(guān)系圖，針對(duì)[Ur=3、][5、][9、][11]這4個(gè)典型約化速度（對(duì)應(yīng)圖中A、B、C、D），給出了圓周表面的渦黏性熵產(chǎn)損失圖。由圖9可知，初始分支中（[Ur=2～4]），來(lái)流速度較低，流體結(jié)構(gòu)干擾作用較弱，圓柱受到較小的升力作用，導(dǎo)致結(jié)構(gòu)從水流中獲取的水動(dòng)能較少，為0.120 W。隨著約化速度增加，響應(yīng)進(jìn)入上端分支（[Ur=5～10]），此時(shí)鎖定區(qū)間振幅明顯增加，隨著渦動(dòng)力的逐漸增強(qiáng)，流體結(jié)構(gòu)耦合作用增強(qiáng)，[Ur=10]時(shí)結(jié)構(gòu)從水流中獲取的能量達(dá)到最大為[6.486 W]。進(jìn)入下端分支后（[Ur=11～12]），此時(shí)來(lái)流速度更大，響應(yīng)處于非鎖定狀態(tài)，獲取的能量出現(xiàn)跳躍性減小，與上端分支相比降低了約80%，導(dǎo)致尾流場(chǎng)紊亂，圓柱表面因黏性耗散的熵產(chǎn)和尾流中因湍動(dòng)能耗散的熵產(chǎn)都增加，水動(dòng)能損失較大，俘能結(jié)構(gòu)僅獲取較少的水動(dòng)能[PVIV]。對(duì)比能量轉(zhuǎn)換效率的曲線(xiàn)圖可知，隨著俘能結(jié)構(gòu)俘獲的水動(dòng)能[PVIV]的增加，俘能效率[η]不斷增加。上端分支中，俘能最大效率達(dá)到30.5%。進(jìn)入下端分支，水流速度的增加使流體能量[Pfluid]隨之增加，由于[Pfluid]的增加速率高于俘能結(jié)構(gòu)在水流中俘獲水動(dòng)能的增加速率，渦表面的黏性熵產(chǎn)損失也隨著流速的增加逐漸增大，導(dǎo)致能量轉(zhuǎn)換效率開(kāi)始下降，圓柱渦激振動(dòng)的黏性熵?fù)p耗增加，水動(dòng)能損失增加，俘能結(jié)構(gòu)獲取的水動(dòng)能減少，能量轉(zhuǎn)換效率也隨之減少。

3.3.3 熵產(chǎn)損失分布及能量耗散分析

能量耗散主要是因?yàn)闇u動(dòng)能和湍動(dòng)能的耗散，這個(gè)過(guò)程常伴隨著速度梯度的巨大變化，從而導(dǎo)致尾流區(qū)不穩(wěn)定的流動(dòng)。從圖8基于熵理論得出的尾渦強(qiáng)度和耗散率圖可知，水流撞擊圓柱表面時(shí)產(chǎn)生的能量損耗，黏性熵?fù)p耗和溫度熵?fù)p耗主要集中在圓柱與水的來(lái)流方向接觸面，尾流區(qū)有湍流熵?fù)p耗的產(chǎn)生。黏性熵產(chǎn)主要是由不可逆黏性耗散引起的，主要發(fā)生在圓柱靠近入口方向。由于均勻來(lái)流受圓柱附近渦的影響，改變了流動(dòng)方向，形成渦流和沖擊損失，圓柱表面仍有小面積黏性耗散。與黏滯產(chǎn)生熵相比，湍流產(chǎn)生熵更大。

圖10將熵產(chǎn)理論用于評(píng)估三維渦激振動(dòng)能量損失，給出了約化速度[Ur=3、][5、][9、][11]時(shí)流場(chǎng)的熵產(chǎn)分布圖，以分析各約化速度下渦激振動(dòng)響應(yīng)與俘能結(jié)構(gòu)表面及尾流能量耗散間的關(guān)系。從圖10可看出，隨著約化速度的不斷增加，圓柱表面的黏性熵?fù)p耗和周?chē)耐牧黛負(fù)p耗逐漸增加。當(dāng)俘能結(jié)構(gòu)響應(yīng)處于上端分支時(shí)（[Ur=5，Ur=9]），從圓柱后脫落的尾流旋轉(zhuǎn)角速度增大，旋渦強(qiáng)度隨切向速度的增大而增大。此時(shí)流動(dòng)受到強(qiáng)烈干擾，黏度熵產(chǎn)和湍流熵產(chǎn)都比較大。能量轉(zhuǎn)換效率與圓柱渦激振動(dòng)響應(yīng)振幅、頻率等息息相關(guān)，結(jié)合圖9能量轉(zhuǎn)換效率的曲線(xiàn)可知，[Ur=5]的時(shí)候能量轉(zhuǎn)換效率高，熵?fù)p耗小。由于尾渦處流體之間的動(dòng)能交換，在俘能結(jié)構(gòu)的尾流區(qū)可觀(guān)察到數(shù)個(gè)湍流產(chǎn)生熵。

4 結(jié) 論

本文通過(guò)雙向流固耦合數(shù)值方法對(duì)均勻來(lái)流下三維圓柱渦激振動(dòng)俘獲低速水流能進(jìn)行研究，利用熵產(chǎn)理論獲得了俘能結(jié)構(gòu)渦激振動(dòng)響應(yīng)，俘能大小和俘能效率與能量損失分布及耗散間的關(guān)系，在本研究實(shí)驗(yàn)參數(shù)范圍內(nèi)，得到如下主要結(jié)論：

1）三維效應(yīng)下，渦激振動(dòng)的振幅頻率響應(yīng)結(jié)果與實(shí)驗(yàn)結(jié)果接近，隨著約化速度的增加，俘能結(jié)構(gòu)尾流的流向渦對(duì)展向渦的控制作用更強(qiáng)。振幅越大，尾流區(qū)脫落渦橫向?qū)挾葘?huì)變更大，尾渦模式也從初始分支的2S模式變?yōu)樯隙朔种У?P和P+T模式，進(jìn)入下端分支，發(fā)展為P+S模式。

2）流體中的水動(dòng)能其增加速率高于俘能結(jié)構(gòu)俘獲水動(dòng)能的增加速率，隨著約化速度的增加，上端分支內(nèi)，當(dāng)[Ur=6]時(shí)俘能結(jié)構(gòu)的俘能效率最大達(dá)到30.5%，俘能結(jié)構(gòu)俘獲的水動(dòng)能增加，而渦表面的黏性熵產(chǎn)損失增大，導(dǎo)致能量轉(zhuǎn)換效率逐漸下降。然而進(jìn)入下端分支后，圓柱渦激振動(dòng)的湍流耗散熵增加，水動(dòng)能損失增加較快，俘能結(jié)構(gòu)獲取的水動(dòng)能驟減80%。

3）利用熵產(chǎn)理論準(zhǔn)確捕捉到能量損失的位置。隨著約化速度的增加，圓柱表面的黏性熵?fù)p耗和圓柱周?chē)拔擦鲄^(qū)湍流熵?fù)p耗逐漸增加，且黏性熵產(chǎn)主要發(fā)生在圓柱表面，湍流熵發(fā)生在耗散率較大的尾流區(qū)。

引用熵產(chǎn)理論可有效預(yù)測(cè)和分析俘能結(jié)構(gòu)的能量損失，為俘能結(jié)構(gòu)設(shè)計(jì)和性能優(yōu)化提供參考，促進(jìn)俘能結(jié)構(gòu)在低速水流能中的應(yīng)用。

［參考文獻(xiàn)］

［1］ BARRERO-GIL" A，" ALONSO" G，" SANZ-ANDERS" A.Energy harvesting from transverse galloping［J］. Journal of sound and vibration， 2010， 329（14）： 2873-2883.

［2］ 羅竹梅， 聶聰， 郭濤. 渦激振動(dòng)驅(qū)動(dòng)的柱群結(jié)構(gòu)集中俘獲海流能研究［J］. 太陽(yáng)能學(xué)報(bào)， 2021， 42（4）： 89-94.

LUO Z M， NIE C， GUO T. Study on centralized harvesting ocean current energy with column-group structure by VIV ［J］. Acta energiae solaris sinica， 2021， 42（4）： 89-94.

［3］ GARCIA E M H. Prediction by energy phenomenology for harnessing" "hydrokinetic" "energy" "using" "vortex-induced vibrations［D］. Michigan： University of Michigan， 2008.

［4］ BERNITSAS M M， BEN-SIMON Y， RAGHAVAN K， et al. The VIVACE converter：model tests at high damping and Reynolds number around 105［J］. Journal of offshore mechanics" "and" "arctic" "engineering，" "2009，" "131（1）： 011102.

［5］ LEE" J" H，" BERNITSAS" M" M." High-damping，" high-Reynolds VIV tests for energy harnessing using the VIVACE" converter［J］." Ocean" engineering，" 2011，" 38（16）： 1697-1712.

［6］ SUN" "H，" "BERNITSAS" "M" nbsp;M." "Bio-inspired" "adaptive damping in hydrokinetic energy harnessing using flow-induced oscillations［J］. Energy， 2019， 176： 940-960.

［7］ DING L， BERNITSAS M M， KIM E S. 2-D URANS vs. experiments of flow induced motions of two circular cylinders in tandem with passive turbulence control for 30， 000＜Re＜105，000［J］. Ocean engineering， 2013， 72： 429-440.

［8］ DING W J， SUN H， XU W H， et al. Numerical investigation on interactive FIO of two-tandem cylinders for hydrokinetic" "energy" "harnessing［J］. Ocean engineering， 2019， 187： 106215.

［9］ KONDO" "N." Three-dimensional" "computation" "for" "flow-induced vibrations in in-line and cross-flow directions of a circular" cylinder［J］. International" journal" for" numerical methods in fluids， 2012， 70（2）： 158-185.

［10］ MITTAL S. The critical mass phenomenon" in" vortex-induced vibration at low Re［J］. Journal of fluid mechanics， 2017， 820： 159-186.

［11］ KOCK F， HERWIG H. Local entropy production in turbulent shear flows： a high-Reynolds number model with wall functions［J］. International journal of heat and mass transfer， 2004， 47（10-11）： 2205-2215.

［12］ HERWIG H， GLOSS D， WENTERODT T. A new approach to understanding and modelling the influence of wall roughness on friction factors for pipe and channel flows［J］. Journal of fluid mechanics， 2008， 613： 35-53.

［13］ GUILMINEAU E， QUEUTEY P. Numerical simulation of vortex-induced vibration of a circular cylinder with low mass-damping in a turbulent flow［J］. Journal of fluids and structures， 2004， 19（4）： 449-466.

［14］ PAN Z Y， CUI W C， MIAO Q M. Numerical simulation of vortex-induced vibration of a circular cylinder at low mass-damping" "using" "RANS" "code［J］." Journal" of" "fluids" and structures， 2007， 23（1）： 23-37.

［15］ KHALAK A， WILLIAMSON C H K. Motions， forces and mode transitions in vortex-induced vibration at low mass-damping［J］. Journal of fluids and structures， 1999， 13（7-8）： 813-851.

［16］ 譚俊哲， 王保振， 袁鵬， 等. 湍流強(qiáng)度對(duì)渦激振動(dòng)潮流能轉(zhuǎn)換裝置振子響應(yīng)影響研究［J］. 太陽(yáng)能學(xué)報(bào)， 2020， 41（10）： 20-26.

TAN J Z， WANG B Z， YUAN P， et al. Study on turbulence intensity influence on cylindrical oscillator response of VIV tidal energy conversion device［J］. Acta energiae solaris sinica， 2020， 41（10）： 20-26.

［17］ 羅竹梅， 張立翔. 影響從渦激振動(dòng)中獲取能量的參數(shù)研究［J］. 振動(dòng)與沖擊， 2014， 33（9）： 12-15.

LUO Z M， ZHANG L X. Influence of parameters on extracting energy from vortex-induced vibration［J］. Journal of vibration and shock， 2014， 33（9）： 12-15.

ENERGY CAPTURE ANALYSIS OF THREE-DIMENSIONAL VORTEX-INDUCED VIBRATION OF LOW-SPEED WATER FLOW BASED ON ENTROPY PRODUCTION THEORY

Li Jun1，Luo Zhumei1，Guo Tao2，Yang Tao1，Gao Suoming3

（1. School of Metallurgy and Energy Engineering， Kunming University of Science and Technology， Kunming 650093， China；

2. School of Engineering Mechanics， Kunming University of Science and Technology， Kunming 650500， China；

3. Kunming Branch of Shenzhen Water Planning and Design Co.， Ltd.， Kunming 650032， China）

Abstract：In this paper， the rigid cylinder is used as the energy-capture structure， and the three-dimensional vortex-induced vibration of the energy-captured system in the low-speed water flow environment is simulated by the strong fluid-solid coupling numerical method. The captured energy quantity and the energy capture efficiency of the capturing structure is analyzed" in the range of resonance reduction velocity（Ur=2-12）. To obtain the effect of the three-dimensional wake of the energy captured structure on the captured energy， the entropy production theory and the wake characteristics are used to capture the energy loss source and distribution of the wake vortex of the energy capture structure， and the relationship between capture energy and wake dissipation under different vibration branches are obtained. The results show that the energy loss analysis based on the entropy production theory can accurately capture the energy dissipation trend of the energy-harvesting structure. The viscous entropy production mainly occurs on the cylindrical surface， and the turbulent entropy occurs in the wake region， where the dissipation rate is more significant. With the increase of the reduced velocity， the kinetic energy of water captured by the energy capture structure in the upper branch increases， the viscous entropy loss on the cylinder surface and the turbulent entropy production loss around the cylinder and in the wake region also increases， and the viscous entropy production loss on the surface of the energy capture structure increases. The increase of turbulent entropy production loss in the wake region is the main reason for the decrease in energy capture efficiency.

Keywords：ocean energy； energy capture structure； vortex-induced vibration； fluid-structure interaction； energy dissipation； entropy generation theory