穆浩然馬小舟毛艷軍高 翔

(1.大連理工大學(xué) 建設(shè)工程學(xué)部,遼寧 大連116024;2.海岸和近海工程國家重點(diǎn)實(shí)驗(yàn)室,遼寧 大連116024)

輻射應(yīng)力概念最初是由Longuest-Higgins和Stewart[1]提出并給出了適用于計(jì)算微幅波輻射應(yīng)力值的簡化公式。在近岸波浪的傳播運(yùn)動(dòng)中,由于海水深度變淺導(dǎo)致了波浪產(chǎn)生嚴(yán)重變形甚至破碎,從而引起輻射應(yīng)力的強(qiáng)烈改變,輻射應(yīng)力變化對(duì)沿岸波生流、波浪增水(set-up)、波浪減水(set-down)有著重要影響,也是低頻波浪生成的重要驅(qū)動(dòng)力。Longuest-Higgins和Stewart[2]認(rèn)為波浪輻射應(yīng)力的增加將促使波群中自由長波的釋放;Symonds等[3]建立了控制方程中含有輻射應(yīng)力受迫項(xiàng)的碎波拍模型,認(rèn)為波群中波浪破碎點(diǎn)前后移動(dòng)引起輻射應(yīng)力在時(shí)間及空間上的改變產(chǎn)生了低頻波浪;Kostense[4]定性地驗(yàn)證了短波群波浪破碎點(diǎn)移動(dòng)導(dǎo)致低頻波浪產(chǎn)生的理論。因此,為了更加準(zhǔn)確地了解輻射應(yīng)力這一波浪特性,眾多學(xué)者通過理論分析,數(shù)值和實(shí)驗(yàn)?zāi)M的手段進(jìn)行輻射應(yīng)力的計(jì)算。對(duì)于常水深Graham[5]提出用有限差分法計(jì)算輻射應(yīng)力張量;鄭金海和嚴(yán)以新[6]利用線性波理論研究波浪輻射應(yīng)力張量隨深度變化的分布規(guī)律,并給出適用于任意波向角的輻射應(yīng)力的表達(dá)式;Xia等[7]提出了輻射應(yīng)力垂向剖面的定義和計(jì)算公式,并利用斜坡上波生流的實(shí)驗(yàn)結(jié)果證明了模型的準(zhǔn)確性;Stive和Wind[8]利用實(shí)驗(yàn)數(shù)據(jù),通過線性外推法近似求得湍流波浪速度剖面來計(jì)算輻射應(yīng)力沿程變化;溫秀媛等[9]利用改進(jìn)色散關(guān)系的Boussinesq方程推導(dǎo)出了全新的輻射應(yīng)力計(jì)算公式,探究輻射應(yīng)力對(duì)波浪增減水的影響。然而,上述輻射應(yīng)力計(jì)算方法都難以獲得碎浪帶內(nèi)詳細(xì)的流場信息,并存在著一定程度的條件假設(shè)和公式簡化,因此計(jì)算所得輻射應(yīng)力結(jié)果誤差較大。

基于OpenFOAM[10]軟件開源程序建立數(shù)值水槽則是模擬波浪破碎區(qū)流場信息的一個(gè)較為可行的手段,可以詳細(xì)研究完整計(jì)算域內(nèi)波浪輻射應(yīng)力隨波浪傳播而產(chǎn)生的變化過程。查晶晶和萬德成[11]基于OpenFOAM軟件進(jìn)行的數(shù)值造波實(shí)驗(yàn)驗(yàn)證了此數(shù)值方法造波和消波方式的可靠性;王東旭等[12]利用此軟件較好地模擬了孤立波在潛礁上的波浪破碎及水躍現(xiàn)象;毛艷軍等[13]應(yīng)用Open FOAM軟件求解了箱式浮式防波堤的水動(dòng)力性能,并分析了其作為垂蕩浮子式波能轉(zhuǎn)換裝置的能量轉(zhuǎn)換性能;姚宇等[14]運(yùn)用Open-FOAM軟件模擬孤立波在島礁地形的傳播,數(shù)值輸出的破碎區(qū)內(nèi)波形及流速與實(shí)驗(yàn)結(jié)果對(duì)比良好;張陳浩和鄭茜[15]利用OpenFOAM軟件模擬規(guī)則波在淺灘上的破碎變形,模擬結(jié)果與實(shí)驗(yàn)結(jié)果相吻合,且輸出了能較好體現(xiàn)波浪破碎整體過程的流場信息。

本文擬應(yīng)用OpenFOAM軟件來模擬波浪在潛堤地形下的傳播破碎,在進(jìn)行平底水槽線性波的數(shù)值模擬和數(shù)值模型模擬波浪的準(zhǔn)確性及計(jì)算輻射應(yīng)力方法的有效性驗(yàn)證后,進(jìn)行規(guī)則波在潛堤地形下傳播破碎的模擬研究,在此基礎(chǔ)上對(duì)輻射應(yīng)力沿程變化的特點(diǎn)及其對(duì)波浪增減水的影響開展了研究,并探討了平底地形下輻射應(yīng)力與波幅偏移量之間的關(guān)系。

1 理論基礎(chǔ)

1.1 控制方程

本文采用基于Open FOAM軟件的造波模塊waves2Foam[16]開展波浪傳播變形的數(shù)值模擬。坐標(biāo)系統(tǒng)為笛卡爾坐標(biāo)系,x方向?yàn)椴ɡ说娜肷浞较?y方向?yàn)榇瓜?z方向?yàn)檠厮鄣膶挾确较?基于不可壓縮流體的運(yùn)動(dòng)特點(diǎn),控制方程為:

式(1)～式(5)中,ρ為混合流體的密度;u為流體的速度場;t為時(shí)間;μ為動(dòng)力黏性系數(shù);g為重力加速度;x為位置矢量;fσ為自由表面的張力;σ為張力系數(shù);κ為自由表面的平均曲率;n為界面單位法向量;α為流體體積分?jǐn)?shù)。

對(duì)于自由表面的處理,根據(jù)VOF法將2種流體看作一種混合流體,實(shí)現(xiàn)每個(gè)單元相界面的追蹤,引入相函數(shù):

在求解時(shí),式(2)中的混合流體的密度ρ可用液體密度ρ1和氣體密度ρ2表示,動(dòng)力黏性系數(shù)μ可用液體黏性系數(shù)μ1和氣體黏性系數(shù)μ2來表示:

同時(shí),α滿足輸運(yùn)方程:

傳統(tǒng)的VOF法需要反復(fù)地進(jìn)行自由表面的重構(gòu),降低了計(jì)算效率,因此Weller[17]引入了額外的人工壓縮項(xiàng)來實(shí)現(xiàn)對(duì)自由水面的精確捕捉,式(9)可變?yōu)?/p>

式中:uc為適用于壓縮自由面及調(diào)整自由面尖銳程度的速度場,新增加的人工壓縮項(xiàng)只對(duì)界面過渡區(qū)起作用,不會(huì)影響到周圍流場。

1.2 輻射應(yīng)力公式

將作用于單位面積水柱體的總動(dòng)量流的時(shí)均值減去沒有波浪作用時(shí)的靜水壓力定義為波浪輻射應(yīng)力。Longest-Higgins和Stewart[18]給出的水平x方向輻射應(yīng)力(S xx)的定義公式為

式中:u為質(zhì)點(diǎn)速度在x方向的分量;h為水深;η為水面高程;i為平均水面高程;p為總壓。

對(duì)于微幅波,輻射應(yīng)力計(jì)算公式可簡化為

2 模型驗(yàn)證

2.1 波形驗(yàn)證

表1 不同工況下的網(wǎng)格尺寸Table 1 Sensitivity analysis of grid size

圖1 不同網(wǎng)格尺寸下波面變化與解析波面對(duì)比Fig.1 Comparison diagram of wave surface between experimental results and numerical simulation with different grid size

將模型作用于模擬平底水槽中線性波的傳播,其模擬工況為:波高0.02 m,波周期2 s,水深1 m,計(jì)算域長30 m,高1.1 m,數(shù)值水槽前端和后端各設(shè)置6 m松弛區(qū)(消波),模擬計(jì)算時(shí)長30 s,采用自動(dòng)調(diào)整計(jì)算時(shí)間步長的方法,設(shè)置輸出時(shí)間步長Δt=0.02 s。模型應(yīng)用Case1,Case2和Case3三種不同大小的網(wǎng)格設(shè)置來進(jìn)行敏感性分析,具體網(wǎng)格尺寸見表1。不同網(wǎng)格尺寸下波面對(duì)比如圖1所示,由圖可見,Case1網(wǎng)格尺度下模擬的波浪波峰值明顯低于解析結(jié)果,由此產(chǎn)生的誤差將影響到之后的數(shù)值計(jì)算,而Case2和Case3網(wǎng)格尺度下模擬的波浪則與理論解對(duì)應(yīng)很好。因此,為保證波浪輻射應(yīng)力計(jì)算結(jié)果的準(zhǔn)確性并且能夠節(jié)約計(jì)算資源,本文采用Case2網(wǎng)格尺寸進(jìn)行模擬計(jì)算。得到了較好的數(shù)值模擬結(jié)果,其波面與理論波面的運(yùn)動(dòng)趨勢及波峰波谷值均吻合較好,誤差較小,由此說明本文所用數(shù)值模型可以很好地模擬波浪運(yùn)動(dòng)。

2.2 輻射應(yīng)力計(jì)算結(jié)果驗(yàn)證

將上述工況中數(shù)值模擬計(jì)算的輻射應(yīng)力沿程變化結(jié)果與微幅波理論下輻射應(yīng)力結(jié)果進(jìn)行比對(duì),對(duì)比結(jié)果如圖2所示。由式12和式13計(jì)算得到微輻波理論下的輻射應(yīng)力值為0.4709 N/m,數(shù)值模擬結(jié)果中輻射應(yīng)力值有較小幅度的振蕩,可能是由波浪沿水槽傳播過程中波幅產(chǎn)生微小差異所致,但整體而言,本模型模擬結(jié)果與理論結(jié)果對(duì)應(yīng)得較好,說明本數(shù)值模型模擬波浪運(yùn)動(dòng)輸出的流場、壓力場較準(zhǔn)確,輻射應(yīng)力計(jì)算方法也是可行的。

圖2 平底水槽中數(shù)值模擬輻射應(yīng)力值及微輻波理論輻射應(yīng)力結(jié)果對(duì)比Fig.2 Comparison of numerical simulation and airy wavetheory results of radiation stress at flat tank

3 潛堤地形下輻射應(yīng)力計(jì)算及波浪增減水模擬

潛堤地形的實(shí)驗(yàn)案例采用Beji和Battjes[19]的實(shí)驗(yàn)地形及工況,數(shù)值模擬Case A和CaseB兩種工況下波浪的傳播,并計(jì)算輻射應(yīng)力及波浪增減水的沿程變化情況。本文采用層流模型來進(jìn)行數(shù)值計(jì)算,利用相對(duì)較小的網(wǎng)格尺度進(jìn)行波浪破碎的求解,能夠得到較好的波浪模擬結(jié)果,且與湍流模型RANS相比,層流模型計(jì)算所得的破波區(qū)內(nèi)波浪運(yùn)動(dòng)與實(shí)際波浪運(yùn)動(dòng)更為相似,能更好地模擬出波浪的破碎形態(tài)和波浪破碎后波面的起伏振蕩現(xiàn)象,由于篇幅有限,在此不作更多的說明。數(shù)值模型布置如圖3所示,堤頂高度為0.3 m,波浪沿水槽正向傳播。數(shù)值模擬計(jì)算域長40.0 m,前后各設(shè)置8.0 m松弛區(qū),計(jì)算時(shí)長為50 s,CaseA,CaseB均采用自動(dòng)調(diào)整計(jì)算時(shí)間步長的方法,輸出時(shí)間步長分別為Δt=0.02 s,Δt=0.025 s。總網(wǎng)格量均為80萬,x方向網(wǎng)格尺寸為0.01 m,y方向網(wǎng)格尺寸為0.0025 m。數(shù)值模擬實(shí)驗(yàn)工況如表2所示。

圖3 潛堤地形布置(m)Fig.3 Layout of the submerged bars(m)

表2 潛堤地形實(shí)驗(yàn)工況列表Table 2 List of experimental cases of submerged bars topography

3.1 波形對(duì)比

當(dāng)波浪傳播到潛堤前坡坡頂位置時(shí),非線性效應(yīng)達(dá)到最大,波浪發(fā)生明顯的破碎。將堤頂G2,G3和G4三處波面的時(shí)間變化曲線與實(shí)驗(yàn)結(jié)果對(duì)比,結(jié)果見圖4,可以看出,在波浪初破階段,Case A和CaseB工況的G2,G3處數(shù)模波峰波谷值與實(shí)驗(yàn)結(jié)果對(duì)應(yīng)較好;在波浪完全破碎后(G4),2個(gè)工況的數(shù)值模擬計(jì)算結(jié)果均出現(xiàn)了低估峰值和谷值的現(xiàn)象。但是數(shù)值模擬的兩工況的整體運(yùn)動(dòng)趨勢,特別是破碎后波浪振蕩趨勢與實(shí)驗(yàn)現(xiàn)象基本一致,說明該數(shù)值模型可以較好地反映實(shí)驗(yàn)中波浪破碎后的實(shí)際運(yùn)動(dòng),模擬的波浪破碎具有較大的可信度。圖5給出了潛堤堤頂位置(24～26 m處)波浪破碎時(shí)的流場信息,可以清楚地看到,當(dāng)破碎發(fā)生時(shí),破碎點(diǎn)周圍水質(zhì)點(diǎn)的速度方向均指向破碎點(diǎn)處,越靠近破碎點(diǎn)的水質(zhì)點(diǎn)速度值越大,在破碎波面處水質(zhì)點(diǎn)運(yùn)動(dòng)最為劇烈,速度值達(dá)到最大。而距離破碎點(diǎn)較遠(yuǎn)的流場中,水質(zhì)點(diǎn)速度較小,水體流動(dòng)較為溫和,但是受前一破碎波浪的影響,水體中水質(zhì)點(diǎn)移動(dòng)較為混亂,方向具有不確定性,這與波浪運(yùn)動(dòng)的實(shí)際情況較為符合,因此本數(shù)值方法計(jì)算所得流場信息較為可信。

圖4 CaseA和CaseB兩種工況下G2、G3、G4處波面對(duì)比Fig.4 Comparison of wave surface at G2,G3 and G4 under CaseA and CaseB

圖5 Case A工況堤頂位置流場分布Fig.5 The flow field distribution at the top of the dam under Case A

3.2 輻射應(yīng)力與波浪增減水

當(dāng)波浪傳播到潛堤區(qū)域時(shí),會(huì)發(fā)生淺水變形并最終產(chǎn)生破碎,繼而波高驟降,波浪動(dòng)量減小,輻射應(yīng)力將沿程發(fā)生變化,進(jìn)而影響平均水平面的改變,這一波浪平均水平面的變化便稱為波浪的增減水。當(dāng)規(guī)則波沿x方向正向入射時(shí),依據(jù)波浪增減水方程可知:

當(dāng)輻射應(yīng)力(S xx)沿程增長,即時(shí),則波浪將產(chǎn)生減水現(xiàn)象;當(dāng)輻射應(yīng)力沿程降低,即時(shí),則波浪將產(chǎn)生增水現(xiàn)象。波浪增減水的值即為波面升高的時(shí)均值與靜水面位置的差值。本文根據(jù)數(shù)值模型輸出的波面數(shù)據(jù),進(jìn)行周期內(nèi)時(shí)間平均,得到Case A和CaseB兩個(gè)工況下的波浪增減水變化。

圖6和圖7給出了潛堤地形下2種工況的波浪輻射應(yīng)力和波浪增減水的沿程變化曲線。2個(gè)工況均能明顯看出,當(dāng)波浪傳播在堤前平底時(shí),輻射應(yīng)力基本呈穩(wěn)定值,到潛堤前坡坡底附近其值開始增大,直至堤頂處達(dá)到最大值,隨后產(chǎn)生下降的趨勢,在潛堤后坡面處由于水深的增加,輻射應(yīng)力值再次增大,直至堤后平底水槽處趨于平穩(wěn)。而波浪增減水的沿程變化趨勢恰好與之相反。由于波浪完全破碎,波浪成分不穩(wěn)定,在堤后坡面至平底處的輻射應(yīng)力和平均水平面均有明顯震蕩,但是整體變化趨勢依舊能較好地滿足波浪增減水方程。

圖6 CaseA波浪輻射應(yīng)力與增減水沿程變化情況Fig.6 Variation of wave radiation stress and wave set-up,wave set-down under Case A

圖7 CaseB波浪輻射應(yīng)力與增減水沿程變化情況Fig.7 Variation of wave radiation stress and wave set-up,wave set-down under CaseB

另外,在堤前平底處,模型穩(wěn)定后的輻射應(yīng)力結(jié)果較小,這是由于隨著波浪的傳播,堤前反射波使波幅隨著時(shí)間推移產(chǎn)生下沉現(xiàn)象,從而導(dǎo)致輻射應(yīng)力定義公式中的速度積分項(xiàng)和動(dòng)壓積分項(xiàng)變小,最終導(dǎo)致了輻射應(yīng)力結(jié)果偏低。由圖8可見,Case A工況下15 m水槽處波面過程曲線及波幅偏移值隨時(shí)間變化的曲線,將波幅偏移值量化為波浪上下幅值之和,并將其進(jìn)行無量綱化,其中au表示波浪的上幅值,ad表示波浪的下幅值。波幅的偏移值在30 s后基本達(dá)到穩(wěn)定狀態(tài),同時(shí)由圖9可見數(shù)值水槽前端平底15 m處不同周期內(nèi)的波幅偏移值與其對(duì)應(yīng)的輻射應(yīng)力值關(guān)系,較明顯地反映出潛堤工況下的平底水槽部分輻射應(yīng)力結(jié)果與對(duì)應(yīng)周期內(nèi)波浪幅值的偏移值成正相關(guān),且受波幅變動(dòng)影響較為敏感,波幅偏移值越小,輻射應(yīng)力值越低。結(jié)合圖8和圖9來看,圖8中當(dāng)波浪波幅偏移值達(dá)到穩(wěn)定狀態(tài)時(shí),其值對(duì)應(yīng)于圖9中輻射應(yīng)力值大概為0.5 N/m處的點(diǎn),較理論預(yù)估值低,而波浪傳播初期的波浪偏移值對(duì)應(yīng)的同周期內(nèi)的輻射應(yīng)力結(jié)果約為3 N/m,這與理論預(yù)估值基本一致,由此可見平底水槽處波浪穩(wěn)定時(shí)輻射應(yīng)力結(jié)果偏低是波幅偏移值變小造成的。而造成波浪波幅產(chǎn)生振蕩的原因還需進(jìn)一步地分析探討。

圖8 CaseA工況中水槽15 m處波面變化及波幅偏移量變化Fig.8 Wave surface change and amplitude offset change at 15 m on the tank under Case A

圖9 CaseA工況水槽15 m處不同周期內(nèi)波浪波幅偏移量與其對(duì)應(yīng)輻射應(yīng)力值關(guān)系Fig.9 The relation figure between wave amplitude offset and corresponding value of radiation stress at 15 m of tank under CaseA

4 結(jié) 語

本文應(yīng)用Open FOAM數(shù)值模型模擬有破碎情況下潛堤地形的波浪傳播,詳細(xì)地輸出了碎浪區(qū)波浪運(yùn)動(dòng)的流場信息,完整地計(jì)算了波浪在堤前平底、堤身及堤后平底處整個(gè)沿程的輻射應(yīng)力和波浪增減水變化,并探究了堤前平底處輻射應(yīng)力值受波幅偏移值的影響。

結(jié)果表明,數(shù)值模擬波浪在平底水槽中傳播的波形及輻射應(yīng)力結(jié)果與解析結(jié)果吻合較好;在含潛堤地形的水槽中,得到了與實(shí)驗(yàn)結(jié)果相吻合的破波區(qū)波面對(duì)比圖,并計(jì)算輸出了整個(gè)水槽內(nèi)包括碎浪帶及堤后平底水槽處的波浪輻射應(yīng)力變化,發(fā)現(xiàn)波浪完全破碎后堤后坡面及平底處輻射應(yīng)力呈先增大后趨于平穩(wěn)的趨勢,且堤后平底處較堤前平底處輻射應(yīng)力值更大。

輻射應(yīng)力結(jié)果與波浪增減水變化趨勢在碎浪帶內(nèi)也能較好地對(duì)應(yīng)波浪增減水方程;并且由于波浪反射造成的波幅的上浮與下沉,會(huì)間接導(dǎo)致輻射應(yīng)力結(jié)果的增大與減小,且輻射應(yīng)力值受波幅偏移值影響較為敏感。