高 瑞,王一帆,張倍銘,吳賀賀,林 東,葉建設(shè)
(1 中國水產(chǎn)科學(xué)研究院漁業(yè)機械儀器研究所,上海 200092;2 青島海洋科技中心,青島 266237;3 北京跟蹤與通信技術(shù)研究所,北京 100094;4 海鵬海洋工程(上海)有限公司,上海 201803)
高海況打撈技術(shù)是中國海洋捕撈技術(shù)的拓展與延伸[1]。高海況打撈系統(tǒng)主要由打撈船、被打撈浮體及打撈網(wǎng)具組成。被打撈浮體與打撈船在高海況下易受風(fēng)浪影響,且相對于船舶,浮體尺寸很小,其運動還易受到船舶輻射擾流干擾;打撈網(wǎng)具是由目腳、結(jié)節(jié)、沉子、浮子等構(gòu)件所組成的復(fù)雜柔性結(jié)構(gòu)體,各構(gòu)件因其型式、材質(zhì)等不同,所受波浪力、慣性力、拖曳力等各不相同。因此,打撈船、網(wǎng)具與浮體在高海況下受到的流體作用力具有強非線性特點,三者之間的任一運動狀態(tài)的改變,都會影響到整個系統(tǒng)的運動響應(yīng)。
單獨針對船舶或浮體的運動響應(yīng)研究方面,許勇等[2]采用三維勢流理論及多體動力學(xué)理論建立了波浪中近距離并行航行多船波浪作用力及運動響應(yīng)的計算模型,并通過模型試驗進(jìn)行了驗證。Zhou等[3]基于SST k-Ω湍流模型,比較分析了單船與兩船并行航行的水動力差異。鄭平宇等[4]利用AQWA軟件開展了補給船與接收船并行補給過程中的耐波性研究,獲取的結(jié)果與譜分析結(jié)果進(jìn)行了比較。謝楠等[5]利用三維線性勢流理論計算了兩個浮體的水動力相互作用,并與試驗結(jié)果進(jìn)行了對比。杜一豪等[6]采用統(tǒng)計學(xué)方法深入探討了不規(guī)則波浪作用下Wigley型船的運動響應(yīng)問題,結(jié)果表明船舶橫搖方向與升沉和縱搖方向隨機運動的響應(yīng)特征有顯著差異。陳京普等[7-8]研究了10萬噸級油輪在長峰和短峰不規(guī)則波中的運動響應(yīng)問題。張杰等[9]對國際集裝箱標(biāo)模在迎浪規(guī)則波中的運動響應(yīng)和波浪增阻等進(jìn)行數(shù)值計算,計算得到的船體運動幅值能與試驗結(jié)果相一致。呂向琪等[10]采用頻域Rankine面元法計算了Wigley I、Wigley III和S175在有航速時船舶的水動力系數(shù)、波浪激勵力和運動響應(yīng)。李冬琴等[11]采用單參數(shù)Lagrangian支持向量回歸算法用于訓(xùn)練并構(gòu)建了代理模型以預(yù)報船舶耐波性能。
單獨針對柔性網(wǎng)具的研究方面,Lee等[12]建立了水下柔性結(jié)構(gòu)數(shù)值模型,系統(tǒng)地考慮了材料的剛度和數(shù)值穩(wěn)定性,并且提出了網(wǎng)目合并處理辦法。孫霄峰等[13]建立了單船拖網(wǎng)網(wǎng)具的集中質(zhì)量模型,在建立操縱性方程時將船體力、推進(jìn)器、舵力、風(fēng)載荷等因素考慮其中。Tsukrov[14]采用等效網(wǎng)單元對網(wǎng)衣系統(tǒng)進(jìn)行等效合并從而建立了網(wǎng)衣系統(tǒng)的連續(xù)有限單元法。Fredheim[15]采用了有限元辦法對流力作用下的錐形網(wǎng)衣變形特性,結(jié)合模型試驗結(jié)果研究了目腳和結(jié)節(jié)水動力系數(shù)的取值范圍。Wan等[16]考慮了目腳的彈性變形和水動力的耦合作用,獲取了水下網(wǎng)具系統(tǒng)的平衡狀態(tài)和張力分布。林禮群等[17]開展了靜水條件下網(wǎng)具與浮體的打撈試驗,研究了船舶速度、網(wǎng)口高度、打撈物體吃水深度等參數(shù)對網(wǎng)具的影響。陳昌平等[18]基于大渦模擬開展了金屬網(wǎng)衣多孔小直徑網(wǎng)狀結(jié)構(gòu)的水動力響應(yīng)特性仿真分析,得到了不同目腳尺寸和網(wǎng)線直徑組合條件下鋅鋁合金網(wǎng)衣的平面受力計算結(jié)果。陳天華等[19]基于集中質(zhì)量點法和網(wǎng)目群化技術(shù)研究了樁柱式圍網(wǎng)單元網(wǎng)片在水流作用下的水動力特性。倪震宇等[20]使用數(shù)值計算方法對單錨張綱張網(wǎng)漁具的水動力特性和形狀進(jìn)行模擬,獲取了網(wǎng)衣和綱索的張力載荷分布及隨流速的變化情況。
針對高海況下的柔性網(wǎng)具與浮體及打撈船耦合影響的研究較為有限,孫洪波等[21]、王飛等[22]、苑志江等[23]、朱克強等[24]、朱軍等[25]對船-拖纜或船/纜/體系統(tǒng)開展了操縱性運動仿真分析,但是拖纜為單個纜繩而非網(wǎng)具。
本研究通過數(shù)值仿真技術(shù),開展了浮體與打撈船的耦合水動力研究,將浮體運動響應(yīng)及浮體與船舶相對位置作為設(shè)計輸入,并開展了浮體與柔性網(wǎng)具的耦合水動力特性研究。
1.1.1 控制方程
連續(xù)性方程與動量方程如下所示:
(1)
(2)
湍流模式采用Realizablek-ε兩方程模型,具體表達(dá)如下:
(3)
式中:μt為流體渦粘系數(shù),Gk為速度梯度產(chǎn)生的湍動能,Gb為浮力產(chǎn)生的湍動能,Ym為可壓縮湍流中波動膨脹的貢獻(xiàn),αk與αε分別為湍動能k與湍流耗散率ε的普朗特數(shù),C2、C1ε和C3ε為常量,Sk與Sε為用戶定義的源項。C1為時均應(yīng)變率的函數(shù),表示為:
(4)
自由液面捕捉采用VOF法,其輸運方程如下:
(5)
1.1.2 運動方程
將打撈船與浮體的運動視作剛體運動,本研究主要考慮縱傾與升沉兩個自由度,運動方程如下:
(6)
(7)
網(wǎng)具主要由細(xì)長網(wǎng)線、浮子、沉子等構(gòu)成,可將網(wǎng)具離散為通過無質(zhì)量桿件連接的質(zhì)量點的集合。浮子、沉子等視作球體,質(zhì)量集中于各自形心且其受力作用于質(zhì)量點上,網(wǎng)線受力及屬性(質(zhì)量、浮力、曳力等)被平均分配給該線段的兩個端點結(jié)節(jié)處,線自身被視為無質(zhì)量桿件,只傳遞拉力。網(wǎng)具運動方程如下所示[27,29-30]:
(8)
式中:mi為第i個質(zhì)量點的質(zhì)量;Δmi為第i個質(zhì)量點的附加質(zhì)量;(ti,bi,ni)為第i個質(zhì)量點的空間位置坐標(biāo);T為質(zhì)量點受到的拉力,下標(biāo)x表示受到的x方向上的力,其余類推;F為質(zhì)量點受到的水動力;W為質(zhì)量點水中重量,即重力與浮力的合力。
質(zhì)量點j對質(zhì)量點i所受拉力由下式表達(dá):
(9)
式中:E為材料彈性模量;d為網(wǎng)線直徑;lij為質(zhì)量點i、j之間的網(wǎng)線實際長度;l0為質(zhì)量點i、j之間的網(wǎng)線未伸長長度。作用在質(zhì)量點i上的拉力在空間上的分量由下式表達(dá):
(10)
式中:N表示與質(zhì)量點i相連的所有質(zhì)量點的集合。
由于網(wǎng)線、浮子、沉子等網(wǎng)具構(gòu)件直徑相對于波長、波高為小量,且由于其特征尺寸與入射波波長之比小于0.2,可采用莫里森方程計算網(wǎng)具在波浪載荷作用下的水動力,公式如下:
(11)
式中:Cd為阻力系數(shù);A為投影面積;vcx、vcy及vcz為流速在3個坐標(biāo)軸上的分量。
針對船體及浮體,采用重疊網(wǎng)格技術(shù)實現(xiàn)其運動姿態(tài)的模擬。運用重疊網(wǎng)格技術(shù)時,一般將計算域分為背景區(qū)域和重疊區(qū)域,重疊區(qū)域即為運動物體的周邊區(qū)域。各區(qū)域網(wǎng)格單獨生成,流場數(shù)據(jù)通過重疊網(wǎng)格邊界條件實現(xiàn)耦合。重疊網(wǎng)格示意圖如圖1所示,船體周邊立方體即為重疊區(qū)域。
圖1 重疊網(wǎng)格示意
計算模型縮尺比為1∶10,計算用某打撈船如圖2所示,由于計算不涉及船體上層建筑,因而為簡化模型信息、減少計算時間,這里不考慮針對上層建筑進(jìn)行建模工作。
圖2 打撈船三維模型
該打撈船主尺度無量綱化后如表1所示。
表1 打撈船主尺度無量綱值
計算用浮體如圖3所示,主尺度參數(shù)如表2所示。
表2 浮體主尺度無量綱值
圖3 浮體三維模型
2.2.1 計算域與邊界條件
計算域的選取應(yīng)確保流體在邊界處不發(fā)生回流以干擾到浮體及打撈船周圍區(qū)域,同時需綜合考慮邊界選取對流動特性和計算成本的影響,本研究計算域及邊界條件設(shè)置如表3所示。
表3 計算域及邊界條件設(shè)置
2.2.2 網(wǎng)格劃分
網(wǎng)格劃分具體見圖4。
圖4 網(wǎng)格劃分
維系波浪在整個計算域中不隨空間和時間發(fā)生數(shù)值耗散是波浪條件下數(shù)值模擬的重要前提。一般而言,在一個完整波浪的波高方向上,應(yīng)不少于20個網(wǎng)格數(shù)量;在一個完整波浪的波長方向上,應(yīng)不少于40個網(wǎng)格數(shù)量。為將波浪完全包裹,在波浪區(qū)域劃分網(wǎng)格時,還需將加密網(wǎng)格區(qū)域的高度設(shè)置為1.2倍的波高。為盡可能減少計算所需時間成本,加速收斂,同時避免波浪在出口處產(chǎn)生回流以影響到船舶及浮體周圍流場,在距離船舶2倍波長的下游處,通過設(shè)置粗網(wǎng)格對波浪進(jìn)行數(shù)值消波,粗網(wǎng)格在波浪范圍的設(shè)置基本為波浪區(qū)域網(wǎng)格設(shè)置的一半左右。
2.2.3 求解設(shè)置
設(shè)置時間及空間差分均為二階格式,求解方式采用基于壓力的瞬態(tài)求解,時間步長為0.005 s,設(shè)置時間步長時應(yīng)保證庫朗數(shù)小于0.5。庫朗數(shù)計算公式為:
Co=vC·Δt/min(Δx)
(12)
式中:Co為庫朗數(shù);vC為流速;Δt為時間步長;min(Δx)為最小網(wǎng)格尺寸。
2.2.4 計算工況
為獲取浮體與船舶之間的最佳相對位置,以便于打撈船對浮體進(jìn)行打撈作業(yè),根據(jù)實際打撈經(jīng)驗給出9種計算工況,見表4。
表4 計算工況
造波算例驗證如圖5所示。
圖5 造波算例驗證
根據(jù)上述網(wǎng)格及設(shè)置策略,以波高0.1 m、波長1.5 m,以規(guī)則波均勻流場作為初始值開展造波算例驗證,圖5a顯示了計算20 s后的波形圖,圖5b顯示了波高監(jiān)測情況。其中,左側(cè)為進(jìn)口,右側(cè)為出口,出口處往前至2倍波長處設(shè)置為消波區(qū)。由圖5可知,消波區(qū)域有效降低了波幅,減少了波的傳輸能量,防止波在計算域中的回流,造波區(qū)域波幅沒有明顯的衰減。
3.1 網(wǎng)具等效簡化
為反映網(wǎng)具變形及受力情況,建模過程中,僅考慮網(wǎng)線的軸向彈性力,其他諸如質(zhì)量、浮力等參數(shù)都平均分配至網(wǎng)線兩端的節(jié)點處,波浪載荷同樣也加載至節(jié)點處。此外,考慮到網(wǎng)具精確建模的計算成本過高,參考國內(nèi)外相關(guān)研究[26-28]對其進(jìn)行等效簡化。
具體而言,以r為等效比,可得到如下關(guān)系式:
(13)
deq·Leqtotal=d·Ltotal
(14)
式中:m為橫向網(wǎng)目數(shù)量;n為縱向網(wǎng)目數(shù)量;l為目腳長度;d為目腳直徑;Ltotal為網(wǎng)線總長度,下標(biāo)eq為等效后對應(yīng)數(shù)據(jù)。為保持等效后網(wǎng)片的質(zhì)量相同,目腳等效橫截面積由下式計算:
(15)
計算用模型如圖6所示。
圖6 網(wǎng)具示意圖
經(jīng)簡化后的網(wǎng)具各構(gòu)件參數(shù)見表5。其中,浮子總浮力2.5 kN,平均分配至圖6b中的浮子上,沉子總質(zhì)量100 kg,平均分配至圖6b中的沉子上。
表5 簡化后網(wǎng)具各構(gòu)件參數(shù)
為確保計算模型在動態(tài)時域分析中能夠?qū)崿F(xiàn)收斂,對網(wǎng)囊和浮體相對位置進(jìn)行了限定,以確保重物有合理的邊界條件約束。分析時,時間步長取0.005 s。
以工況3為例,打撈船與浮體在高海況下的升沉、縱搖時歷曲線如圖7所示。打撈船在高海況下的運動響應(yīng)呈現(xiàn)較強的周期性變化,且響應(yīng)幅值也較小。相比打撈船,浮體由于主尺度遠(yuǎn)小于船舶,且處于長波浪中,波浪力載荷對浮體的影響遠(yuǎn)大于船舶,體現(xiàn)為波浪中浮體的運動幅值大于船舶運動幅值。
圖7 打撈船與浮體運動時歷曲線
此外,由于打撈船在開展捕撈作業(yè)時與浮體距離較近,浮體除受到波浪載荷外,還受到船體的輻射擾流影響,導(dǎo)致浮體的縱搖響應(yīng)出現(xiàn)較強的不規(guī)則性,圖8中打撈船與浮體在波浪中的運動姿態(tài)變化體現(xiàn)了這點。
圖8 打撈船與浮體運動姿態(tài)變化
為分析浮體與打撈船沿船寬及船長方向不同位置受力情況,圖9給出不同計算工況下浮體橫向受力情況,其中計算工況1~5是浮體在相對打撈船同一縱向位置、不同橫向位置下受到的橫向力,計算工況6~9是浮體在相對打撈船同一橫向位置、不同縱向位置下受到的橫向力。
圖9 浮體在打撈船不同位置下的橫向受力
針對不同橫向位置處,浮體受到的橫向力不存在明顯的相位差。且其大小與浮體及打撈船之間的間距呈現(xiàn)負(fù)相關(guān),即當(dāng)浮體與打撈船間距越近,浮體受到的波浪力及輻射力也越大,其運動情況也更為劇烈。從操作便利性角度出發(fā),當(dāng)浮體與打撈船間距相對接近時,打撈操作相對便利。因此,從保持浮體在高海況下的穩(wěn)定性與打撈船在高海況下的操作便利性這兩方面考慮,存在一個較為平衡的橫向位置。
針對不同縱向位置處,不同工況間存在一定的相位差,尤其以工況6最為明顯,這應(yīng)該與波浪傳播方向和船舶縱向船長方向一致有關(guān),浮體在不同的縱向位置,其受到波浪力也存在著一定的相位差。除此外,不同工況下的橫向受力大小沒有特別明顯的變化,說明浮體受到的橫向力受浮體與船舶的橫向間距影響較大,受浮體與船舶的縱向間距影響較小??傮w看,工況3處對應(yīng)的浮體與打撈船相對位置較優(yōu),后續(xù)在進(jìn)行網(wǎng)具與浮體分析時,將以工況3位置及浮體對應(yīng)運動響應(yīng)作為設(shè)計輸入。
圖10顯示了網(wǎng)囊中不含浮體與含浮體時網(wǎng)具在波浪中變形響應(yīng)情況。當(dāng)網(wǎng)囊中不包含浮體時,從其運動姿態(tài)中可以看出,網(wǎng)具主要隨波面運動,并被提上綱繩、沉子與浮子共同限制;當(dāng)網(wǎng)囊中包含浮體時,浮體大部分浮出水面且隨波面運動明顯,牽引網(wǎng)囊與其同步運動。另外,網(wǎng)具同時也被提上綱繩、沉子與浮子共同限制。
圖10 網(wǎng)具入水姿態(tài)
圖11給出網(wǎng)具中網(wǎng)囊網(wǎng)線、網(wǎng)身力綱、側(cè)網(wǎng)曳綱以及提上綱繩等構(gòu)件的受力極值時歷曲線,需指出的是,由于模型進(jìn)行了簡化處理,實際每個構(gòu)件均由若干網(wǎng)線組成,這里僅給出出現(xiàn)極值的網(wǎng)線時歷情形。由圖11可知,網(wǎng)囊網(wǎng)身整體受力較為平均,側(cè)網(wǎng)上、下曳綱受力差別較大,且側(cè)網(wǎng)下曳綱為網(wǎng)具各部分中受力最大的構(gòu)件,這是由于側(cè)網(wǎng)上曳綱在浮子的作用下漂浮于水面,而下曳綱在沉子的作用下浸沒入水,波浪載荷及浮體運動作用于下曳綱使其張緊程度高于上曳綱。圖12為網(wǎng)具總受力曲線。
圖11 網(wǎng)具各部位受力時歷曲線
圖12 網(wǎng)具總受力及對應(yīng)波面
由圖12可知,網(wǎng)具總受力為121.11 kN。需指出的是,網(wǎng)具總受力為各個構(gòu)件在時域上的疊加,雖然每個構(gòu)件的極值基本出現(xiàn)在同一個波浪周期內(nèi),但時刻點并非完全重合??梢钥闯?網(wǎng)具受力與波幅有直接關(guān)系,波幅越高其受力越大,而受力極值也出現(xiàn)在最大波幅所在的一個波浪周期內(nèi)。
研究發(fā)現(xiàn),浮體受到的橫向力受浮體與船舶的橫向間距影響較大,受浮體與船舶的縱向間距影響較小;浮體與打撈船存在一個較為平衡的橫向位置,在該位置下浮體受船舶輻射擾流影響較小,便于打撈船高海況下開展打撈作業(yè);網(wǎng)具受力情況與波幅正相關(guān),波幅越高網(wǎng)具受力越大,網(wǎng)具受力極值出現(xiàn)在最大波幅所在的一個波浪周期內(nèi)。本研究將打撈回收分解為打撈船與浮體的耦合計算以及打撈浮體與網(wǎng)具的耦合計算兩部分,傳遞參數(shù)為浮體相對船舶的位置及浮體運動響應(yīng),假設(shè)浮體與尾部網(wǎng)囊無相對位移,這可能導(dǎo)致網(wǎng)囊兜住浮體后的形狀與實際有所出入,若要考慮網(wǎng)囊兜住浮體后的變形,則涉及剛性固體與柔性彈性體在復(fù)雜流體環(huán)境下的多重耦合,后續(xù)將對此做進(jìn)一步研究。