結(jié)構(gòu)物砰擊入水過程的CFD-DEM-IBM數(shù)值仿真

底瑛棠,趙蘭浩,毛 佳

(河海大學(xué) 水利水電學(xué)院,南京 210098)

入水問題廣泛存在于自然界和工程的眾多領(lǐng)域,如導(dǎo)彈入水[1]、船舶砰擊[2]、滑坡涌浪[3]甚至打水漂[4]等,具有重要的實(shí)際應(yīng)用價(jià)值。然而,入水問題的精確高效數(shù)值仿真仍面臨艱巨的技術(shù)挑戰(zhàn)。當(dāng)具有復(fù)雜幾何外形的結(jié)構(gòu)撞擊水面時(shí),固體大運(yùn)動(dòng)邊界的處理成為首要難點(diǎn);同時(shí),運(yùn)動(dòng)物體對(duì)液面擠壓與排開,必然伴隨一系列的自由液面劇烈演化過程,如液面飛濺、空泡和射流[5]等現(xiàn)象,給數(shù)值模擬帶來額外的困難;另外,結(jié)構(gòu)的運(yùn)動(dòng)決定流場(chǎng)的演化,流場(chǎng)的演化又反過來影響結(jié)構(gòu)的運(yùn)動(dòng)特性,呈現(xiàn)流固之間強(qiáng)烈的耦合效應(yīng)。歸納起來,結(jié)構(gòu)移動(dòng)邊界的合理描述、劇烈變形自由液面的準(zhǔn)確捕捉及各相之間耦合作用的精細(xì)刻畫是入水現(xiàn)象數(shù)值仿真面臨的核心難題。

為解決結(jié)構(gòu)入水過程中的移動(dòng)邊界問題,王永虎等[6]運(yùn)用動(dòng)網(wǎng)格技術(shù),嚴(yán)格按照結(jié)構(gòu)邊界構(gòu)造貼體網(wǎng)格,并隨邊界移動(dòng)進(jìn)行網(wǎng)格變形或重建,構(gòu)建了三維圓柱自由入水模型;然而,當(dāng)結(jié)構(gòu)發(fā)生大位移時(shí),動(dòng)網(wǎng)格逐時(shí)間步的網(wǎng)格更新消耗大,還可能產(chǎn)生網(wǎng)格扭曲等問題。Ma等[7]、謝路毅等[8]采用嵌套網(wǎng)格方法,采用固定的背景網(wǎng)格離散計(jì)算域,無需實(shí)時(shí)更新,但需采取合理的挖洞措施,且難以考慮多體情況下子網(wǎng)格間的相互作用。相比之下,采用固定笛卡爾網(wǎng)格,并利用拉格朗日點(diǎn)來表征復(fù)雜結(jié)構(gòu)邊界的浸入邊界法(IBM)[9]更加簡(jiǎn)潔:通過在流體動(dòng)量方程中附加體力項(xiàng)來施加固體邊界條件,IBM物理概念明確,網(wǎng)格生成簡(jiǎn)單,對(duì)于結(jié)構(gòu)復(fù)雜外形與大幅位移問題處理優(yōu)勢(shì)顯著,目前已被廣泛應(yīng)用于包含復(fù)雜邊界的流固耦合問題[10]。

針對(duì)強(qiáng)非線性的自由面問題,Kleefsman等[11]采用單相流簡(jiǎn)化,僅考慮水相演化并構(gòu)造自由面邊界條件以描述氣相的影響。然而,單相流簡(jiǎn)化無法反映實(shí)際入水過程中氣液兩相劇烈的相互摻雜,也不能解釋氣液耦合流動(dòng)的潛在機(jī)制。目前常用的多相流界面捕捉方法為VOF方法[12]和LS(level set)方法[13]。VOF方法通過構(gòu)造流體體積函數(shù)捕捉界面,以其高度守恒性得到廣泛應(yīng)用,但需要復(fù)雜的界面重構(gòu)技術(shù)。LS方法概念簡(jiǎn)單,通過距離指示函數(shù)隱式追蹤界面,但具有嚴(yán)重的質(zhì)量損失問題。Olsson等[14]提出守恒式LS(CLS)方法,采用雙曲正切函數(shù),以保守的方式對(duì)流,顯著提高了質(zhì)量守恒性,然而,該方法界面邊緣法向?qū)_動(dòng)極為敏感,可能使法向計(jì)算失準(zhǔn)。據(jù)此,Zhao等[15]提出改進(jìn)的守恒式LS方法(ICLS),以CLS方法捕獲界面,并以LS方法獲得界面信息,在保證方法簡(jiǎn)單性的同時(shí)兼顧了守恒性及界面參數(shù)的準(zhǔn)確性要求。

結(jié)構(gòu)入水沖擊過程中的流固耦合機(jī)制十分復(fù)雜。何春濤等[16]采用常速入水模型,只考慮結(jié)構(gòu)位移對(duì)流體的作用,實(shí)際上只求解了包含動(dòng)邊界的純流體力學(xué)問題。而工程中以自由入水為主,結(jié)構(gòu)運(yùn)動(dòng)軌跡及流體的響應(yīng)事先未知,有學(xué)者考慮流固耦合效應(yīng)開展了對(duì)自由入水問題的探索。如Calderer等[10]結(jié)合曲線IBM及LS方法提出流固耦合模型,預(yù)測(cè)浮體與自由面的相互作用;Sun等[17]采用光滑粒子動(dòng)力學(xué)方法處理自由面與剛體間的耦合作用。然而,現(xiàn)有研究普遍局限于單體入水,王聰?shù)萚18]雖對(duì)雙體入水作出嘗試,但均不考慮固體間的碰撞。對(duì)于滑坡體入水或拋石等問題,散體材料自身的相互作用描述也是不可或缺的。

為合理考慮包含離散介質(zhì)相互作用的流固耦合問題,有學(xué)者采用流體動(dòng)力學(xué)-離散單元法(CFD-DEM)[19]。CFD-DEM方法雖引入了接觸模型考慮離散介質(zhì)的物理碰撞,但流固耦合作用力通過Di Felice[20]提出的經(jīng)驗(yàn)公式估算,只適用于圓形顆粒,無法準(zhǔn)確滿足流固耦合邊界條件,流場(chǎng)也是非解析的。底瑛棠等[21]曾提出高解析度CFD-DEM-IBM方法解決了上述問題,但缺乏對(duì)自由面的合理描述,不適用于入水現(xiàn)象仿真。

為實(shí)現(xiàn)具有普遍工程適用價(jià)值的入水問題數(shù)值仿真,本文提出一種解析的CFD-DEM-IBM方法,通過Navier-Stokes方程描述流體,采用固定的笛卡爾網(wǎng)格離散計(jì)算域,依據(jù)距離勢(shì)離散單元法分析固體運(yùn)動(dòng)特性并有效處理固體間的相互作用,引入隱式直接力IBM追蹤固體移動(dòng)邊界并展開相間作用力的精確求解,確保同時(shí)滿足固體表面無滑移邊界條件及連續(xù)性條件,并采用ICLS方法捕捉自由液面,在保證守恒性的前提下獲得準(zhǔn)確的界面信息。為實(shí)現(xiàn)流固間的強(qiáng)耦合,將變量通過插值函數(shù)與分布函數(shù)在流固耦合面上傳遞,采用交錯(cuò)迭代格式求解流固耦合系統(tǒng),在每個(gè)時(shí)間步內(nèi)反復(fù)迭代多次以精確滿足收斂條件。從而,本文所提出的CFD-DEM-IBM方法不僅能夠獲得準(zhǔn)確的耦合作用力及高解析度的流場(chǎng)信息,克服入水問題數(shù)值仿真固有的困難,還能夠合理反映離散塊體間的碰撞,深入描述多體入水現(xiàn)象的復(fù)雜行為特征,為相關(guān)工程解答提供可靠的技術(shù)支撐。

1 數(shù)值計(jì)算方法

1.1 流體控制方程

1.1.1 Navier-Stokes方程與浸入邊界法

流體運(yùn)動(dòng)由不可壓縮Navier-Stokes方程描述:

(1)

?·u=0

(2)

式中:u為流速,t為時(shí)間,ρ為密度,p為壓力,τ=μ(?u+?Tu)為黏性應(yīng)力張量,μ為動(dòng)力黏度,fb為體力項(xiàng),f為附加體力項(xiàng),反映固體邊界對(duì)流體的作用。

根據(jù)CBS(Characteristic-based split scheme)算法[22],時(shí)間離散后式(1)可以寫作下式:

Δu=Gn+Hn+1+fn+1Δt

(3)

如圖1所示,固體由分布在其邊界上的浸入邊界(IB)點(diǎn)來表征,流體域則采用固定的笛卡爾網(wǎng)格離散,流體網(wǎng)格無需隨固體運(yùn)動(dòng)實(shí)時(shí)更新。根據(jù)直接力浸入邊界法[23],流固耦合面上需滿足無滑移邊界條件Vn+1=Un+1,其中Vn+1是IB點(diǎn)的期望速度,Un+1為同一點(diǎn)由周圍流體網(wǎng)格節(jié)點(diǎn)插值得到的速度。定義分布函數(shù)D將IB點(diǎn)上的物理量分布到網(wǎng)格節(jié)點(diǎn)及插值函數(shù)I將網(wǎng)格節(jié)點(diǎn)上的物理量插值到IB點(diǎn)[24],結(jié)合式可得

圖1 入水問題的CFD-DEM-IBM方法示意

Vn+1=Un+1=I(un+Gn+Hn+1+fn+1Δt)

(4)

從而可以將附加體力項(xiàng)寫作下式:

fn+1Δt=D[Vn+1-I(un+Gn+Hn+1)]

(5)

1.1.2 ICLS界面捕捉方法

ICLS方法[15]融合LS及CLS方法的優(yōu)勢(shì),采用守恒性好的CLS方法對(duì)流運(yùn)輸和捕捉界面,LS方法計(jì)算準(zhǔn)確的界面法向。CLS方法中雙曲正切函數(shù)H為

(6)

式中:φ為符號(hào)距離函數(shù),ε為H的彌散寬度。

CLS方法的對(duì)流輸運(yùn)及重新初始化方程分別為:

(7)

(8)

式中:τ為虛擬時(shí)間步,n為重新初始化開始時(shí)刻的界面法向。n可由下式計(jì)算:

(9)

式中δ為給定距離,與界面附近網(wǎng)格尺寸Δh有關(guān)。

LS方法中φ的對(duì)流輸運(yùn)及重新初始化方程為:

(10)

(11)

界面附近φ值根據(jù)式(6)由H反算:

φ=2εtanh-1(2H-1)

(12)

流體密度ρ=ρ1+(ρ2-ρ1)H(φ)和黏性系數(shù)μ=μ1+(μ2-μ1)H(φ)等材料性質(zhì)即可由H計(jì)算,其中,下標(biāo)1和2分別為第1種和第2種流體。

1.2 固體控制方程

剛性塊體的運(yùn)動(dòng)符合牛頓第二定律:

(13)

(14)

根據(jù)勢(shì)函數(shù)離散單元法[25],接觸力Fc由法向接觸力Fcn和切向接觸力Fct構(gòu)成:

(15)

(16)

根據(jù)庫(kù)侖摩擦定律,切向力需要滿足下式:

(Fct)max=Fcn·tanφμ

(17)

式中:(Fct)max為最大允許靜摩擦,φμ為最大靜摩擦角。

注意流體施加在固體上的耦合作用力Ff是基于拉格朗日描述作用在IB點(diǎn)上的,與式中歐拉描述下施加于網(wǎng)格節(jié)點(diǎn)上的附加體力f不同。由文獻(xiàn)[26]有

(18)

式中:Ω″為邊界Γs所圍固體域,Ω為計(jì)算域邊界Γout以內(nèi)、固體邊界Γs以外的流體域,Ω′為由流體域Ω及固體域Ω″共同組成的整個(gè)計(jì)算域,如圖1所示。

1.3 流固耦合系統(tǒng)迭代求解方案

由式(5)可見,引入浸入邊界法后,流體速度、壓力及體力三者相耦合,需要進(jìn)行迭代計(jì)算。本文采用分區(qū)強(qiáng)耦合算法求解流固耦合系統(tǒng),在一個(gè)時(shí)間步內(nèi)迭代數(shù)次以滿足收斂條件,提高解的收斂性與穩(wěn)定性。現(xiàn)假設(shè)第n步的計(jì)算已經(jīng)完成,第n+1步的第k次迭代步驟如圖2所示,圖中‖·‖為任一范數(shù),ζ可以是u、p或δ,εζ為對(duì)應(yīng)項(xiàng)給定的小值。

圖2 CFD-DEM-IBM方法迭代求解步驟

依據(jù)標(biāo)準(zhǔn)伽遼金有限單元法[27]對(duì)流體控制方程進(jìn)行空間離散、速度Verlet算法[28]對(duì)固體控制方程離散,流固控制方程的迭代格式分別為:

(Δun,k,Δpn,k)=Rf(un,pn,δn+1,k-1)

(19)

Δδn,k=Rs(δn,un+1,k,pn+1,k)

(20)

式中:Rf、Rs分別為流體方程與固體方程的右端項(xiàng)。當(dāng)耦合系統(tǒng)迭代過程滿足收斂條件時(shí),跳出循環(huán),否則執(zhí)行第n+1步第k+1次迭代。

采用投影法[29]求解流體相,將速度增量拆分為中間速度場(chǎng)增量Δu*=Gn+fn+1和速度增量的修正量Δu**=Hn+1,如圖2所示,中間速度場(chǎng)u*,k不含壓力和附加體力項(xiàng)。內(nèi)循環(huán)中的附加體力fn+1,k,i由下式計(jì)算:

fn+1,k,iΔt=fn+1,k,i-1Δt+D[Vn+1-I(un+Gn+Hn+θ2+fn+1,k,i-1Δt)]

(21)

內(nèi)、外循環(huán)的收斂判斷條件分別由下式給出:

‖Vn+1-I(un+Gn+Hn+θ2+fn+1,k,i-1Δt)‖<ε

(22)

‖I(un+Kn+Rn+θ2,k)-I(un+Kn+Rn+θ2,k-1)‖<ε

(23)

式中ε為容差。若滿足收斂性要求則退出迭代,否則繼續(xù)執(zhí)行下一迭代步的計(jì)算。

2 數(shù)值仿真結(jié)果

2.1 圓柱常速入水

為驗(yàn)證CFD-DEM-IBM方法對(duì)自由面的準(zhǔn)確捕捉及高解析度流場(chǎng)的反映,開展文獻(xiàn)[10,30]中圓柱常速入水現(xiàn)象仿真。半徑1 m圓柱以v=1 m/s的速度入水,初始時(shí)刻圓心距離水面h=1.25 m,水深20 m,水體密度與黏度分別為1 kg/m3和1.0×10-3Pa·s,對(duì)應(yīng)氣體取值為1×10-3kg/m3和1.8×10-5Pa·s,重力加速度設(shè)為1 m/s2。計(jì)算域?qū)?0 m,高30 m,依據(jù)文獻(xiàn)[10],將計(jì)算域剖分為600×480個(gè)單元,圓柱運(yùn)動(dòng)路徑附近采用Δx=Δy=0.025 m的均勻網(wǎng)格,圓柱表面布置502個(gè)IB點(diǎn),計(jì)算域邊界均為無滑移固壁。

圖3給出了圓柱入水的自由液面位置及渦量云圖,T=vt/h為量綱一的時(shí)間。當(dāng)圓柱接觸靜止水面時(shí),水面由于固體邊界的擠壓而變形,產(chǎn)生沿圓柱邊界的兩股射流;隨著圓柱浸入,射流逐漸變陡并與圓柱表面分離,直至圓柱完全浸沒在水面以下時(shí),液面呈現(xiàn)凹陷狀態(tài);隨后,不穩(wěn)定的液面持續(xù)變形并最終合并,在中間形成向上的飛濺。圖3成功再現(xiàn)了圓柱入水全過程,與文獻(xiàn)[10,30]中的數(shù)值結(jié)果符合,渦量在自由液面曲率發(fā)展的區(qū)域及圓柱附近產(chǎn)生,與文獻(xiàn)[10]吻合,驗(yàn)證了本文方法能夠獲得高解析度流場(chǎng),精細(xì)描述流場(chǎng)結(jié)構(gòu)。

圖3 圓柱入水自由液面位置及渦量云圖

為詳細(xì)論證本文方法捕獲自由液面的準(zhǔn)確性,圖4對(duì)比了文獻(xiàn)[10,30]中的數(shù)值結(jié)果與本文得到的不同時(shí)刻界面輪廓。注意到T=3.0時(shí)刻的差異以及T=4.0時(shí)刻空腔關(guān)閉時(shí)產(chǎn)生的氣泡可能與接觸角邊界條件的設(shè)置有關(guān),在文獻(xiàn)[10,30]中均有所提及。考慮到自由液面的網(wǎng)格敏感性,網(wǎng)格分辨率也對(duì)計(jì)算結(jié)果存在影響。總的來說,文中結(jié)果與以往研究符合較好,展示了本文方法準(zhǔn)確刻畫瞬態(tài)演化自由液面的能力。

圖4 圓柱入水自由液面形態(tài)對(duì)比

2.2 楔形體對(duì)稱垂直入水

圓柱常速入水實(shí)際只是包含移動(dòng)邊界的純流體力學(xué)問題。本文開展對(duì)稱楔形體垂直自由入水?dāng)?shù)值仿真,旨在揭示本文方法處理入水問題時(shí)對(duì)強(qiáng)烈的流固耦合作用的準(zhǔn)確反映及自由液面與移動(dòng)邊界的妥善處理能力。算例取自文獻(xiàn)[11,31-34],計(jì)算參數(shù)如圖5所示,水體密度與黏度分別為1 000 kg/m3和1.0×10-3Pa·s,對(duì)應(yīng)空氣取值為1 kg/m3及1.8×10-5Pa·s。三維楔形體長(zhǎng)1 m,其中測(cè)量截面長(zhǎng)度為0.2 m,本文依據(jù)文獻(xiàn)[31-32]將該三維入水問題簡(jiǎn)化為二維考慮,僅考慮豎向自由度。參照文獻(xiàn)[32],將計(jì)算域剖分為520×280個(gè)網(wǎng)格,楔形體邊界上布置862個(gè)IB點(diǎn),邊界條件設(shè)置同圓柱常速入水。

圖5 楔形體垂直入水初始時(shí)刻示意

楔形體起初由靜止距離水面2.0 m處自由下落,接觸水面時(shí)具有6.15 m/s的初速度。為了驗(yàn)證數(shù)值方法的網(wǎng)格收斂性,保持其他條件不變,對(duì)時(shí)間步長(zhǎng)為Δt=0.000 1 s,最小網(wǎng)格尺寸Δh分別為0.010 00、0.005 00、0.002 50 、0.001 25 m 4種情況下的楔形體入水砰擊力進(jìn)行計(jì)算,結(jié)果如圖6(a)所示。由圖6(a)可見,本文方法具有良好的網(wǎng)格收斂性,4種不同網(wǎng)格尺寸得到的計(jì)算結(jié)果基本一致,后續(xù)計(jì)算選取Δh=0.002 50 m。

圖6 數(shù)值方法時(shí)空離散收斂性驗(yàn)證

同樣地,為驗(yàn)證本文方法的時(shí)間步長(zhǎng)收斂性,保持其他條件不變,采用最小網(wǎng)格尺寸Δh=0.002 50 m,對(duì)比Δt分別為0.001 00、0.000 50、0.000 10、0.000 05 s 4種時(shí)間步長(zhǎng)下的楔形體入水砰擊力。如圖6(b)所示,本文所提出的高解析度CFD-DEM-IBM方法具有良好的時(shí)間步長(zhǎng)收斂性,分別采用4種時(shí)間步長(zhǎng)計(jì)算得到的結(jié)果吻合度高。后續(xù)數(shù)值計(jì)算的時(shí)間步長(zhǎng)統(tǒng)一選為Δt=0.000 10 s。

本文計(jì)算得到的自由液面及壓力云圖如圖7所示,可以看出,楔形體入水后,在楔形體兩側(cè)激起兩股射流,隨后射流沿楔形體邊界上升并在重力作用下發(fā)生彎曲,與邊界表面發(fā)生分離;壓力峰值在入水初期位于楔形體底部尖端附近,隨后移動(dòng)到射流區(qū)根部并逐漸降低。總的來說,壓力分布與Oger等[31]的數(shù)值結(jié)果具有較高的一致性。

圖7 楔形體對(duì)稱入水自由液面及壓力云圖

圖8、9分別給出了垂向砰擊力和楔形體下落速度的定量對(duì)比。由圖8、9可見,楔形體入水過程經(jīng)歷了兩個(gè)階段:第1階段(0 s

圖8 楔形體對(duì)稱入水垂向砰擊力時(shí)程圖

圖9 楔形體對(duì)稱入水垂向速度時(shí)程圖

2.3 楔形體非對(duì)稱斜入水

本文對(duì)三自由度傾斜楔形體自由入水問題展開仿真。計(jì)算參數(shù)與文獻(xiàn)[35]一致,如圖10所示。水體和空氣的密度分別為1 000 kg/m3和1.2 kg/m3,黏度分別為1.0×10-3Pa·s和1.8×10-5Pa·s。楔形體表面分布1 511個(gè)IB點(diǎn),計(jì)算域剖分為796×464個(gè)單元,邊界均為無滑固壁。

圖10 楔形體斜入水初始時(shí)刻示意

楔形體初始由靜止在距離水面0.610 0 m處受重力自由下落,自由液面及渦量云圖如圖11所示。當(dāng)楔形體在空氣中下落時(shí),其左、右尖端處形成渦結(jié)構(gòu)并逐漸脫落,至楔形體撞擊水面時(shí),大的渦結(jié)構(gòu)與液面相互作用破碎成小渦,與文獻(xiàn)[35]的數(shù)值結(jié)果符合,表明本文方法能夠精細(xì)刻畫流體結(jié)構(gòu),獲得高解析度流場(chǎng)信息。非對(duì)稱楔形體入水造成非對(duì)稱的液面飛濺,及楔形體由傾斜“回正”的現(xiàn)象也與文獻(xiàn)[31]的數(shù)值結(jié)果描述相符。

圖11 楔形體斜入水自由液面與渦量云圖

為定量對(duì)比計(jì)算結(jié)果,圖12、13分別將量綱一的垂向加速度與角加速度與已有文獻(xiàn)結(jié)果展開比較。其中,量綱一的垂向加速度/g中y為楔形體垂向位移,g=9.81 m/s2為重力加速度。由圖12、13可見,本文曲線與已有文獻(xiàn)結(jié)果整體趨勢(shì)類似,幅值與文獻(xiàn)[35]中的數(shù)值結(jié)果尤為一致,而與文獻(xiàn)[36]的試驗(yàn)結(jié)果和文獻(xiàn)[37]中的數(shù)值結(jié)果存在一定差異,這是由于本文和文獻(xiàn)[35]中都未考慮壓縮性,且文獻(xiàn)[36-37]也未定義詳細(xì)的計(jì)算域尺寸及水深參數(shù)。總體而言,本文準(zhǔn)確再現(xiàn)了文獻(xiàn)[35]中的斜入水問題,進(jìn)一步驗(yàn)證了CFD-DEM-IBM方法的正確性,也說明了方法對(duì)于復(fù)雜單體入水問題的普適性。

圖12 斜入水楔形體量綱一的加速度(/g)時(shí)程圖

圖13 斜入水楔形體角加速度時(shí)程圖

2.4 多體入水

涉及塊體間相互作用的多體入水問題數(shù)值仿真尤為復(fù)雜,在以往的研究中鮮見報(bào)道。本文針對(duì)5個(gè)相同剛性圓柱的入水問題展開模擬,以揭示本文方法對(duì)復(fù)雜多體入水問題乃至更一般的工程問題的適用性。計(jì)算參數(shù)如圖14所示,固體密度為2 500 kg/m3,各圓柱表面布置471個(gè)浸入邊界點(diǎn),計(jì)算域由500×400個(gè)單元離散,流體特性與邊界條件設(shè)置同楔形體非對(duì)稱斜入水。圓柱的最大靜摩擦角取為0,法向與切向剛度系數(shù)均取為2×107N/m。

圖14 多圓柱入水初始時(shí)刻示意

圖15給出了多圓柱入水過程中的自由水面和壓力云圖。初始時(shí)刻,5個(gè)分散的圓柱在一定排列方式下由靜止受重力作用下落,如圖15(a)所示;隨后在t=0.638 s時(shí)刻,底部3個(gè)圓柱撞擊自由水面,產(chǎn)生巨大的砰擊力,如圖15(b)中圓柱底部呈現(xiàn)高壓區(qū);從而,底部圓柱速度驟降,而上方圓柱仍持續(xù)下落,造成原來分散的5個(gè)圓柱發(fā)生碰撞,如圖15(c)所示,此時(shí)水面發(fā)生變形,出現(xiàn)水花飛濺現(xiàn)象;緊接著,不穩(wěn)定的自由面持續(xù)變形,如圖15(d)所示,由于復(fù)雜的流固耦合作用及固體間的碰撞,壓力的分布與單體入水時(shí)全然不同,下方圓柱間飛濺的液面也由于上方圓柱的阻礙發(fā)生變形彎曲,并推動(dòng)上方圓柱,使5個(gè)圓柱重新變得分散。總之,本文對(duì)多體入水時(shí)劇烈的自由面演化、塊體間的碰撞及強(qiáng)烈的流固耦合效應(yīng)都得到了很好的反映,驗(yàn)證了本文CFD-DEM-IBM方法處理多體入水問題的靈活性與普適性。

圖15 多圓柱入水自由液面與壓力云圖

圖16給出了圓柱的加速度時(shí)程圖。由圖16可見,初始5個(gè)圓柱在重力作用下勻加速規(guī)則下落;t=0.620 s時(shí)下方圓柱砰擊水面急劇減速,而上方圓柱仍保持以重力加速度向下運(yùn)動(dòng);t=0.730 s時(shí)刻,下方3個(gè)圓柱在中間射流作用下分散,從而上方圓柱與圓柱4發(fā)生碰撞,使得上方圓柱速度驟降而圓柱4受到高速撞擊加速下落;在t=0.780 s時(shí),圓柱間發(fā)生第2次碰撞,下方兩側(cè)圓柱受到上方圓柱的撞擊后開始加速下落。結(jié)合圖15可見,多體入水現(xiàn)象呈現(xiàn)出高度對(duì)稱性,且在不同入水階段存在不同的力的主導(dǎo)趨勢(shì),而本文所提出的CFD-DEM-IBM方法能夠兼顧離散體間的碰撞處理、自由液面的瞬態(tài)演化及復(fù)雜的流固耦合效應(yīng)描述,在涉及固體間相互作用的多體入水問題求解方面優(yōu)勢(shì)顯著。

圖16 入水圓柱加速度時(shí)程圖

3 結(jié) 論

1)高解析度CFD-DEM-IBM方法能夠合理描述結(jié)構(gòu)移動(dòng)邊界,準(zhǔn)確計(jì)算流固耦合作用力,獲得高解析度流場(chǎng),實(shí)現(xiàn)對(duì)自由面演化過程的精細(xì)刻畫。

2)通過4種不同網(wǎng)格尺寸與時(shí)間步長(zhǎng)的敏感性分析,驗(yàn)證了CFD-DEM-IBM方法在處理入水問題方面具有良好的收斂性。

3)本文再現(xiàn)了圓柱1 m/s常速入水、底升角30°的楔形體對(duì)稱自由入水及底升角20°的楔形體非對(duì)稱自由入水過程,與以往研究符合較好,揭示了本文方法對(duì)入水問題數(shù)值仿真的準(zhǔn)確性與可靠性。

4)多體入水與單體入水現(xiàn)象迥異,通過設(shè)計(jì)5個(gè)圓柱入水過程仿真,證明了CFD-DEM-IBM方法能夠合理考慮固體間的碰撞,在處理復(fù)雜多體入水問題方面具有獨(dú)到的優(yōu)越性,能夠推廣應(yīng)用于更加一般的實(shí)際工程問題。