毛艷軍馬小舟

(大連理工大學 海岸和近海工程國家重點實驗室,遼寧 大連166024)

在可再生能源的開發(fā)探索中,波浪能發(fā)電裝置(Wave Energy Converter,WEC)得到了廣泛的研究和應(yīng)用,但也面臨著建設(shè)成本較高的問題。針對這一問題,近些年相關(guān)學者提出了成本共享的思路,將波浪能發(fā)電裝置與現(xiàn)有海上結(jié)構(gòu)物集成,一方面可以降低建設(shè)成本,另一方面也可以增加裝置的穩(wěn)定性。其中將波浪能發(fā)電裝置與近岸防波堤進行集成具有較高的可行性和適用性[1-3],此類集成裝置的設(shè)計方案可應(yīng)用于現(xiàn)有防波堤的改造以及新型防波堤的結(jié)構(gòu)概念設(shè)計。He和Huang[4]將振蕩水柱式波浪能轉(zhuǎn)換(Oscillating Water Column,OWC)裝置與浮式防波堤進行集成并進行了實驗研究,結(jié)果顯示裝置在一定條件下可以獲得較好的能量捕獲效率和消波效果。寧德志等[1]對一個由垂蕩浮子式波浪能發(fā)電裝置(Heave Oscillating Buny,HOB)與有樁柱約束的浮式防波堤組成的集成系統(tǒng)進行了實驗研究,在合適的波頻范圍和PTO(Power Take Off)系統(tǒng)阻尼條件下可以獲得較高的能量捕獲率和較好的消波效果。針對此類形式,趙玄烈等分別針對波浪能發(fā)電裝置集成于有樁柱約束的浮式防波堤[2]和波浪能發(fā)電裝置集成于固定的箱式防波堤[3]兩種不同的集成方式,進行了基于線性勢流理論的頻域分析。前一種集成方式采用對稱結(jié)構(gòu),存在極限波能捕獲寬度比0.5。后一種集成方式是將浮子作為附屬結(jié)構(gòu)裝配在固定的透空箱式防波堤前,有很好的適用性和可擴展性。同時,理論分析顯示其可以很好地利用反射波能量,波能捕獲寬度比可達0.8甚至更高,因此后者具有更為優(yōu)異的表現(xiàn)。然而上述理論研究是基于線性勢流理論進行,波浪非線性、流體粘性和流動分離的影響均未考慮,且僅能夠考慮簡化為線性阻尼模型的PTO 系統(tǒng)。在實際海況下,波浪的非線性、流體的黏性和流體分離對于共振浮子的運動響應(yīng)往往具有很強的影響,且波浪能發(fā)電裝置的PTO 系統(tǒng)的阻尼模型也不僅局限于線性的形式[5],其中液壓型PTO 也是常用形式之一,液壓型PTO 采用高壓和低壓蓄壓器,可以光滑能量輸出,其PTO 阻尼力形式上在半周期近似滑動摩擦阻尼力[6],因此可以用Coulomb阻尼模型近似液壓型PTO 阻尼模型。針對以上不足,本文應(yīng)用基于黏性流的計算流體力學軟件OpenFOAM對該模型進行模擬研究,研究線性PTO 阻尼模型Fpto=Cpto·V約束下裝置的能量轉(zhuǎn)換性能和消波性能。并與理論分析結(jié)果進行對比,分析在共振區(qū)間上黏性耗散帶來的影響。同時也應(yīng)用2 種二次平方型,和Coulomb阻尼型非線性的PTO 阻尼模型,研究對比非線性PTO 裝置與線性PTO 的區(qū)別和聯(lián)系。同時為改進此類集成裝置的能量捕獲效率,針對浮子形狀進行優(yōu)化。從而得以獲得更高的能量捕獲效率。

1 理論基礎(chǔ)

1.1 控制方程

本文模型采用OpenFOAM 中的overInter Dy MFoam 求解器求解,流體計算采用有限體積法求解不可壓縮的兩相流N-S方程,不可壓縮連續(xù)性方程及動量方程分別為

式中,u為流體的速度;ug為動網(wǎng)格更新時的網(wǎng)格修正速度;ρ為兩相流的混合密度;p為壓力;S為黏性應(yīng)力張量;fb為有重力作用的體積力和表面張力項;sρu為源項,可用于阻尼區(qū)消波,其中s為阻尼系數(shù)。

采用VOF 方法處理兩相流界面,VOF 方法的原理是定義體積分數(shù)α[0,1],其中0 代表空氣相,1 代表水相,兩者之間的界面采用0 到1 的過渡值來表示,從而實現(xiàn)界面的捕捉。造波方法采用速度入口造波,通過給定速度和體積分數(shù)的入口邊界條件而實現(xiàn)數(shù)值造波。消波方法采用阻尼區(qū)消波的方式,通過在動量方程中添加阻尼項,在消波區(qū)內(nèi)設(shè)置阻尼系數(shù)大于零,從而實現(xiàn)數(shù)值消波。

1.2 重疊網(wǎng)格

重疊網(wǎng)格的基本原理是將不同的部件分別畫在不同的網(wǎng)格上,然后融合到一個背景網(wǎng)格上,其中各部件的網(wǎng)格與背景之間網(wǎng)格有重疊區(qū)域,利用這些重疊區(qū)域網(wǎng)格,通過插值進行信息的雙向傳遞,從而在重疊網(wǎng)格的基礎(chǔ)上實現(xiàn)計算域的更新。重疊網(wǎng)格主要包括計算單元(calculated),插值單元(interpolated),洞單元(hole)三類網(wǎng)格(圖1)：洞單元網(wǎng)格作為圖上黑色網(wǎng)格,為結(jié)構(gòu)內(nèi)部網(wǎng)格,在計算時間步內(nèi)的物理量值凍結(jié),不隨計算更新;插值單元網(wǎng)格作為圖上中灰色網(wǎng)格,為信息傳遞的主要部分,在子部件和背景網(wǎng)格上分別有一圈插值單元用于雙向信息的傳遞和插值;計算單元為正常參與流場控制方程計算部分的網(wǎng)格,圖上灰色網(wǎng)格。可參考沈志榮[7]具體網(wǎng)格插值方法和實現(xiàn)。

1.3 PTO 系統(tǒng)

PTO 阻尼系統(tǒng)采用線性阻尼模型和非線性Coulomb阻尼模型。線性阻尼系統(tǒng)可表示為

式中,λpto為線性阻尼系數(shù),其中最優(yōu)阻尼系數(shù)(λoptimal)可根據(jù)勢流理論頻域下的解析解獲得。λoptimal的表達式為

式中,M為質(zhì)量,μ為附加質(zhì)量力系數(shù),λ為輻射阻尼,K為水靜力回復(fù)剛度。

非線性Coulomb阻尼模型可表示為

式中,Coulomb阻尼模型因為其階躍特性容易引起數(shù)值振蕩,因此可以做近似的線性化處理如下：

式中,Cpto為預(yù)設(shè)Coulomb阻尼型PTO 阻尼力,Gpto為線性化位置PTO 阻尼系數(shù)。

PTO 阻尼系統(tǒng)采用線性阻尼模型,非線性Coulomb阻尼模型和線性化處理的Coulomb阻尼模型三種阻尼模型如圖2 所示。

圖1 重疊網(wǎng)格Fig.1 A diagram of the overlapping grid

圖2 線性PTO 阻尼模型,Coulomb阻尼模型和線性化的Coulomb阻尼模型Fig.2 Damping models of linear PTO,Coulomb and linearized Coulomb

2 數(shù)值水槽構(gòu)建和模型驗證

2.1 模型構(gòu)建

參考趙玄烈等[3]給出的結(jié)構(gòu)參數(shù)構(gòu)建模型,模型基本布置如圖3,水深h=10 m,浮子寬度a1=2 m,后方固定式防波堤寬度a2=6 m,浮子與防波堤間距離D=1 m,浮子吃水深度d1=1 m,防波堤吃水深度d2=2.5 m。此模型的模擬工作難點在于動網(wǎng)格的處理方法,因附屬浮子與固定箱式防波堤的距離較近,且共振時浮子運動響應(yīng)較大。Open FOAM 中傳統(tǒng)的隨體網(wǎng)格變形不適用于此類在靠近固定邊界處的局部網(wǎng)格變形過大的模擬工作。Open FOAM-v1906版本中的重疊網(wǎng)格的功能,可適用于海洋工程浮式結(jié)構(gòu)物的模擬工作,因此本文應(yīng)用此種網(wǎng)格更新方式。同時可以利用Open FOAM-v1906中的速度入口造波邊界進行數(shù)值波浪水槽的構(gòu)建。入口U 邊界條件為wave Velocity,相體積分數(shù)邊界條件為wave Alpha,壓力邊界條件為fixed Flux Pressure。重疊網(wǎng)格區(qū)域邊界類型為overset。浮子邊界為固壁邊界條件。自由液面和重疊網(wǎng)格區(qū)域進行如圖4的加密處理以滿足計算精度要求,網(wǎng)格總數(shù)為136 050,近壁面網(wǎng)格大小為0.1 m。時間步長控制采用自動調(diào)整時間步長,最大庫朗數(shù)為0.5。

圖3 附屬結(jié)構(gòu)式集成裝置Fig.3 Layout of the integrated device with a subordinate structure

圖4 附屬結(jié)構(gòu)形式裝置附近重疊網(wǎng)格劃分情況Fig.4 Division of the overlapping grid around the integrated device with a subordinate structure

2.2 模型驗證

本節(jié)驗證模型選取一在波浪作用下垂蕩方向無約束運動的浮式方箱進行模擬,背景網(wǎng)格和重疊網(wǎng)格如圖5所示,近壁面網(wǎng)格大小0.02 m,模型幾何參數(shù)詳以及數(shù)值造波的精度驗證詳見已發(fā)表工作[8]。本節(jié)應(yīng)用重疊網(wǎng)格方法進行了浮子在周期T=1.79 s,入射波高H=0.2 m 下的運動響應(yīng)模擬,并與Isaacson的實驗的RAO 值[9]和應(yīng)用隨體變形網(wǎng)格方法模擬的結(jié)果[8]進行對比。垂蕩方向運動響應(yīng)歷時線如圖6,表明應(yīng)用重疊網(wǎng)格方法進行浮式結(jié)構(gòu)物計算可以獲得穩(wěn)定且可靠的計算結(jié)果,進而驗證重疊網(wǎng)格方法的有效性。

圖5 方箱周圍重疊網(wǎng)格網(wǎng)格劃分Fig.5 Division of the overlapping grid around the permeable box

圖6 重疊網(wǎng)格法模擬的樁式約束方箱垂蕩方向運動響應(yīng)時程與隨體變形網(wǎng)格以及實驗值的對比Fig.6 The response time history of heave motion of the floating pile-restraint box simulated by overlapping grid method comparing with those obtained by body deformation grid method and by experiments

3 PTO 阻尼模型與浮子形狀優(yōu)化

3.1 線性PTO 阻尼模型

趙玄烈等[3]文中線性PTO 阻尼力選取為單個浮子在線性頻域解下的理論最優(yōu)PTO 阻尼系數(shù),而在考慮黏性的CFD 模擬中,該最優(yōu)阻尼系數(shù)可能不為最優(yōu)值,且集成系統(tǒng)結(jié)構(gòu)的非對稱性也使得附加質(zhì)量和輻射阻尼發(fā)生變化,也導(dǎo)致了集成系統(tǒng)的最優(yōu)阻尼系數(shù)不再與單一浮子相同。因此本節(jié)采用線性PTO 模型,從頻域下的理論最優(yōu)阻尼系數(shù)出發(fā),研究黏性流下此結(jié)構(gòu)的最優(yōu)PTO 阻尼的特性。結(jié)合理論分析結(jié)果,選取模擬工況：T=4.74,4.45,3.88,3.66,3.45,3.30,3.18,3.04,2.84,2.57,2.49和2.39s,對應(yīng)的無量綱波數(shù)(kh)的范圍為7.073～1.882,水深和波高設(shè)為恒定值(h=10 m,H i=0.5 m)進行模擬研究。

圖7為各波況下附加浮子在理論最優(yōu)PTO 阻尼λoptimal作用下的能量捕獲寬度比(Capture Width Ratio,CWR)與理論分析結(jié)果的對比。研究發(fā)現(xiàn)相對比于單浮子結(jié)構(gòu)形式[8]約0.30的最大CWR值,能量捕獲寬度比明顯提升,其最大CWR值可達0.7。但是相比于理論分析結(jié)果,在共振區(qū)間kh=2.0～5.0上,CWR值明顯下降,理論最大值接近1.0,模擬值相較于理論值下降約30%,可以明顯看到黏性效應(yīng)對于此類裝置能量捕獲效率的影響。對于kh=6.089 97和kh=7.07兩個短波波況,CWR的極小值和第2個峰值也能夠被有效地模擬獲得。選取kh=3.72波況,對PTO 阻尼系數(shù)進行了放縮,選取放縮系數(shù)0.75,1.25,1.5。模擬發(fā)現(xiàn),1.25λoptimal和1.5λoptimal的CWR值都比1.0λoptimal要略大,關(guān)于在考慮黏性影響下理論最有阻尼值的討論詳見文獻[8]。圖8,圖9是反射系數(shù)(K r)和透射系數(shù)(K t)與解析結(jié)果的對比。反射系數(shù)在共振位置附近出現(xiàn)了下降,在共振位置kh=3.72反射系數(shù)達到極小值,意味著反射波波能被吸收。kh=6.09波況反射系數(shù)偏小,因為在模擬中此波況下波浪仍然可以引起浮子的振蕩。由圖7可見,kh=6時CWR不完全是0,可能是由于非線性原因,導(dǎo)致波浪并不是完全的滿足了此處的kh值,從而能量也沒有被完全反射。

圖7 1.0倍最優(yōu)PTO 阻尼約束下CWR 值與解析解的對比值的對比Fig.7 Comparison between the CWR values and the analytical solutions under the constraint of 1.0 times optimal PTO damping

圖8 1.0倍最優(yōu)PTO 阻尼約束下反射系數(shù)與解析解的對比Fig.8 Comparison between the reflection coefficient(K r),and the analytical solution under the constraint of 1.0 times optimal PTO damping

圖9 1.0倍最優(yōu)PTO 阻尼約束下透射系數(shù)與解析解的對比Fig.9 The comparison between the transmission coefficient Kt and the analytical solution under the constraint of 1.0 times optimal PTO damping

3.2 非線性PTO 阻尼模型

上節(jié)所述線性PTO 模型在微幅低速運動中可以獲得較好的近似效果,但當非線性較強,運動響應(yīng)較大時更適合采用非線性PTO 模型,可選用二次非線性阻尼模型。液壓型PTO 裝置因為有光滑能量輸出的效果,常用于WEC 裝置的PTO 系統(tǒng)中,其中Coulomb阻尼模型可作為液壓型PTO 的有效近似模型。因此本節(jié)選用二次非線性阻尼模型和Coulomb阻尼模型作為集成系統(tǒng)的PTO 裝置進行模擬,其中在原型比尺的計算中采用理想的Coulomb阻尼模型未發(fā)現(xiàn)明顯的數(shù)值振蕩,可以獲得穩(wěn)定的計算結(jié)果,原因可能為原型比尺中模型的流體阻尼和靜水回復(fù)力可以抵消數(shù)值振蕩。阻尼系數(shù)取λquadratic=λlinear,λcoulomb=λlinearVmax,其中λlinear對應(yīng)線性PTO 系統(tǒng)在kh=3.01時的理論最優(yōu)阻尼λlinear=5 720.91 kg/s。Vmax=0.48 m/s為線性PTO 作用下浮子的最大運動速度。如圖10為附加浮子在3種PTO 模型作用下,在共振波況T=3.663 s,H=0.5 m 下的運動響應(yīng)的對比。圖11為3種不同PTO 模型的PTO 阻尼力的對比,可以發(fā)現(xiàn),因阻尼系數(shù)選取相同,模擬速度低于1.0 m/s,所以二次非線性阻尼型PTO 系統(tǒng)的運動響應(yīng)比線性PTO 裝置運動響應(yīng)偏大。PTO 阻尼力偏小,并且在PTO 力為0 kg/s的附近表現(xiàn)出二次拋物線形的PTO 力的時程,證明二次非線性阻尼模型可以很好地應(yīng)用于PTO 系統(tǒng)的模擬。Coulomb阻尼模型因其在半周期上保持恒定最大PTO 力,因此其約束的浮子運動響應(yīng)最小。線性阻尼模型,Coulomb阻尼模型以及二次非線性PTO 阻尼模型下集成裝置的能量捕獲寬度比CWR結(jié)果如表1所示?？梢钥吹皆诖朔N條件下,線性PTO 阻尼模型獲得了最大的CWR值。這里其他兩種非線性PTO 阻尼模型表現(xiàn)較差。但是也有待于進一步的優(yōu)化PTO 阻尼系數(shù)的選取,以獲得更大的CWR值。

圖10 3種不同形式PTO 阻尼約束下的浮子垂蕩運動響應(yīng)Fig.10 The heave motion response of the floater under the restraints of three kinds of PTO damping

圖11 3種不同形式PTO 阻尼力的時程Fig.11 Time history of the three kinds of PTO damping forces

表1 3種阻尼模型下裝置的CWR 值Table 1 The CWR values of the device under three kinds of PTO pumping model

3.3 浮子形狀優(yōu)化

為進一步優(yōu)化裝置的能量捕獲效率和消波性能,浮子形狀是影響裝置能量捕獲寬度比的重要參數(shù)之一[10],因此本節(jié)主要就附屬結(jié)構(gòu)中,前方浮子的形狀進行優(yōu)化。由毛艷軍等[8]研究可知,方箱型浮子在底部迎浪測和背浪測兩個直角處產(chǎn)生了明顯的漩渦,從而導(dǎo)致流場粘性效應(yīng)的增加,因此本節(jié)選用2種避免這種結(jié)構(gòu)形式突變的浮子形狀,修圓角形式和圓底形式。模擬后,由圖12可以看到三者的速度矢量圖,方箱形浮子在底部兩個拐角處產(chǎn)生了明顯的漩渦,漩渦從邊壁上分離向外發(fā)展并逐漸耗散。修圓角形浮子底部也有漩渦產(chǎn)生,但是可以看到漩渦的尺度較小,且并未脫離浮子的壁面。圓底形浮子底部流場則光滑很多,流體沿著底部壁面發(fā)展,對整個流場的擾動較小,從而也對應(yīng)流體粘性和湍流耗散影響較小。

圖12 3種形狀浮子周圍的流場變化情況Fig.12 The flow field variations around the floaters with three kinds of shapes

3種形狀浮子的運動時程如圖13 所示,可以看到由于結(jié)構(gòu)水線面寬度相同,排水體積較為接近。從運動響應(yīng)可以看出,修圓角的浮子和方箱形浮子的運動相位基本上一致,運動響應(yīng)幅值算子(Response Amplitude Operador,RAO)略微增大。圓底形的浮子在運動相位上與方箱形浮子產(chǎn)生了一定的相位偏移,同時RAO也有所增加。如表2所示,為3種浮子形狀算例的透射系數(shù)(K t),反射系數(shù)(K r),運動響應(yīng)幅值算子(RAO),捕獲寬度比(CWR)模擬性結(jié)果對比與方箱形浮子對比,修圓角形浮子和圓底形浮子可以保持透射系數(shù)基本上不變,反射系數(shù)下降,圓底形浮子的反射系數(shù)最小,對應(yīng)的運動響應(yīng)和能量捕獲寬度比都是最大的,其中圓底形浮子的CWR達到0.808,相比于方箱形浮子,提升了近0.1,效率提升明顯。

圖13 3種形狀浮子的垂蕩運動響應(yīng)時程Fig.13 The time history of the heave motion response of the floaters with three kinds of shapes

表2 3種形狀浮子的模擬結(jié)果對比Table 2 Comparison of the results(K t,K r,RAO and CWR)simulated for the floaters with three kinds of shapes

4 結(jié) 論

應(yīng)用OpenFOAM 進行了附屬結(jié)構(gòu)形式的波浪能與防波堤的集成系統(tǒng)的數(shù)值模擬工作,分析了考慮黏性和流動分離影響下集成裝置的運動響應(yīng),能量捕獲寬度比和消波性能。通過3種不同形式PTO 系統(tǒng)的對比,研究了PTO 阻尼系統(tǒng)對集成裝置的影響。進一步優(yōu)化浮子形狀,研究了浮子形狀對裝置能量捕獲效率的影響。研究發(fā)現(xiàn)：

1)此類結(jié)構(gòu)在考慮黏性影響的模擬中仍可以獲得較高的能量捕獲寬度比,CWR值可達0.7。反射波能量得到很好地利用,反射系數(shù)(Kr)在共振區(qū)間上降低,同時透射系數(shù)(Kt)在kh＞1.8的較寬波頻范圍內(nèi)都可以維持在0.5以下,保證良好的消波效果。

2)非線性PTO 的開發(fā)和應(yīng)用,進一步擴展了波浪能發(fā)電裝置的PTO 模擬工作。2種PTO 阻尼模型可以維持較好的數(shù)值穩(wěn)定性,在本文選取的阻尼系數(shù)條件下,非線性PTO 阻尼模型作用的浮子能量捕獲寬度比小于線性PTO 阻尼模型。

3)形狀參數(shù)的優(yōu)化顯示圓底形浮子可以減小流體黏性和流動分離的影響,從而可以進一步提高裝置的能量捕獲效率。