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

底瑛棠,趙蘭浩,毛 佳

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

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

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

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

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

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

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

1 數(shù)值計算方法

1.1 流體控制方程

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

流體運動由不可壓縮Navier-Stokes方程描述:

(1)

?·u=0

(2)

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

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

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

(3)

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

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

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

(4)

從而可以將附加體力項寫作下式:

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

(5)

1.1.2 ICLS界面捕捉方法

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

(6)

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

CLS方法的對流輸運及重新初始化方程分別為:

(7)

(8)

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

(9)

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

LS方法中φ的對流輸運及重新初始化方程為:

(10)

(11)

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

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

(12)

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

1.2 固體控制方程

剛性塊體的運動符合牛頓第二定律:

(13)

(14)

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

(15)

(16)

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

(Fct)max=Fcn·tanφμ

(17)

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

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

(18)

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

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

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

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

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

(Δ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分別為流體方程與固體方程的右端項。當耦合系統(tǒng)迭代過程滿足收斂條件時,跳出循環(huán),否則執(zhí)行第n+1步第k+1次迭代。

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

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í)行下一迭代步的計算。

2 數(shù)值仿真結果

2.1 圓柱常速入水

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

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

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

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

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

2.2 楔形體對稱垂直入水

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

圖5 楔形體垂直入水初始時刻示意

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

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

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

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

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

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

圖8 楔形體對稱入水垂向砰擊力時程圖

圖9 楔形體對稱入水垂向速度時程圖

2.3 楔形體非對稱斜入水

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

圖10 楔形體斜入水初始時刻示意

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

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

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

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

圖13 斜入水楔形體角加速度時程圖

2.4 多體入水

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

圖14 多圓柱入水初始時刻示意

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

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

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

圖16 入水圓柱加速度時程圖

3 結 論

1)高解析度CFD-DEM-IBM方法能夠合理描述結構移動邊界,準確計算流固耦合作用力,獲得高解析度流場,實現(xiàn)對自由面演化過程的精細刻畫。

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

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

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