水中高壓脈動氣泡與浮體流固耦合特性研究1)

胡振宇 曹卓爾 李 帥 張阿漫

(哈爾濱工程大學船舶工程學院,哈爾濱 150001)

引言

脈動氣泡廣泛地存在于自然界中,并且被應用于水下爆炸[1-2]、空化空蝕[3-4]、超聲清洗[5]、海底資源勘探[6-7]和醫(yī)療[8-9]等領域.迄今為止,前人對氣泡動力學的研究主要集中于自由場[10]、壁面[11-12]和自由液面[13-14]等簡單邊界條件下的氣泡;對復雜邊界下的氣泡尤其是自由液面?浮體?氣泡三者非線性耦合作用的研究相對匱乏.氣泡與近水面結構物的非線性耦合作用是海洋工程和軍事研究領域的重要科學問題.浮體在氣泡脈動壓力的影響下運動并伴隨著自由液面大變形,這種運動響應反過來又對氣泡動力學特性產生影響.目前氣泡?浮體耦合作用中所蘊含的力學機理仍未完全揭示,因此本文針對該問題進行數(shù)值研究,同時開展放電氣泡實驗驗證數(shù)值模型,旨在增強對氣泡與近水面結構物非線性流固耦合效應的認識,為相關領域的研究和應用提供基礎性支撐.

Rayleigh[10]最早基于球狀假設提出著名的Rayleigh-Plesset(RP) 方程,隨后Herring[15]和Keller等[16]在RP 方程的基礎上進一步考慮流體的黏性和可壓縮性的影響,得到各種不同形式的改進RP 方程.近邊界(自由液面、剛/彈性壁面、其他氣泡界面)脈動氣泡往往無法保持球狀,這種非球型脈動極大地增加了理論研究的難度.在實驗研究方面,研究者使用水下爆炸實驗[17-19]、放電氣泡實驗[20-21]和激光氣泡實驗[22-23]來研究氣泡與各種邊界的耦合作用,高速攝影技術也被用于捕捉氣泡的動力學行為.Lauterborn 和Bolle[24]使用激光聚焦產生單個氣泡以研究氣泡與剛性壁面的耦合作用,發(fā)現(xiàn)在氣泡坍塌后期產生指向壁面的高速射流,其速度可達120 m/s.高速水射流與結構表面的砰擊被認為是空蝕的重要成因[11,25].然而由于實驗攝影中時空分辨率的限制,精細地捕捉射流內部形態(tài)仍然是一個挑戰(zhàn).通過數(shù)值模擬,能得到更多的流場信息以分析氣泡與流場邊界的非線性耦合作用.

Chahine 等[26]最早針對放電氣泡與柱形浮體的流固耦合作用進行了實驗和數(shù)值研究,分別使用剛固浮體與可動浮體進行研究發(fā)現(xiàn)浮體的運動響應對氣泡坍塌模式、射流形成和氣泡脈動周期具有重要影響.Klaseboer 等[27-28]對水下爆炸氣泡與漂浮船體的耦合作用進行了數(shù)值研究,發(fā)現(xiàn)氣泡射流方向受到自由液面、船體表面曲率和重力的影響.在數(shù)值計算中使用了鏡像法來模擬自由液面,該方法對于γf>1.4(γf為無量綱氣泡與自由液面距離,取氣泡最大半徑為特征長度) 的工況,精度較好,但對于近距離的強耦合工況(γf<1.4) 需要考慮全非線性自由液面條件.Zong 等[29]使用邊界積分法(BIM)聯(lián)合有限元法建立了氣泡與剛性/可變形浮體耦合模型,考慮全非線性自由表面條件以及結構的運動和彈性變形,并采用傳統(tǒng)的差分法求解結構表面壓力(伯努利方程中??/?t項),然而這種方法可能不適用于結構運動很快的強耦合工況[30-31].Li 等[30]首次將輔助函數(shù)法應用于氣泡與水中懸浮球體的非線性耦合效應研究中,相比傳統(tǒng)差分法,計算更為準確,本文將該方法拓展至氣泡與浮體的流固耦合特性研究中.近年來,有限體積法[32]、有限元法[33]等方法也被成功的應用于氣泡動力學的數(shù)值研究中,相比于這些域離散方法,BIM 僅需對流場邊界進行離散,這種降維處理使BIM 計算效率具有極大的優(yōu)勢.

本文聯(lián)合BIM 與輔助函數(shù)法建立脈動氣泡與浮體全非線性耦合模型.首先開展了一系列放電氣泡與柱形浮體耦合作用機理實驗,然后對比數(shù)值與實驗結果,驗證計算模型與方法的有效性.采用本文的數(shù)值模型系統(tǒng)地研究無量綱氣泡與浮體距離參數(shù)γs對氣泡動力學特性和結構響應的影響規(guī)律.本文的研究內容安排如下;第二章給出本文數(shù)值模型的基本理論;第三章首先驗證數(shù)值模型然后分析討論實驗與數(shù)值結果,主要關注氣泡射流模式、射流速度、氣泡遷移以及浮體運動響應;第四章給出本文研究的結論.

1 基本理論和數(shù)值模型

1.1 邊界元法

本文采用基于勢流理論的軸對稱邊界元法,建立脈動氣泡與浮體非線性流固耦合作用數(shù)值模型.基于流體無黏、無旋、不可壓縮假設,流場(水域)可由拉普拉斯方程描述

上式可轉化為積分形式的邊界積分方程[34-35]

式中,p為控制點,q為源點,A為p點觀察流場(水域)的立體角,S為包括氣泡表面Sb、自由液面Sf、浮體濕表面Sw以及無窮大球面S∞的流體邊界,n表示邊界法向,以指向流場外為正.圖1 所示為氣泡與浮體非線性流固耦合作用數(shù)值模型示意圖,引入柱坐標系O(r,θ,z),其坐標原點位于初始氣泡中心,z軸指向重力的反方向且為對稱軸,定義浮體的入水深度為d,初始氣泡中心至初始自由液面的距離為h;本文考察圓柱型浮體與氣泡的耦合作用,圓柱浮體半徑為Rs,氣泡最大等效半徑為Rm.

圖1 氣泡與浮體非線性流固耦合作用數(shù)值模型示意圖Fig.1 Sketch of the numerical model for investigating the nonlinear fluid-structure interaction between a bubble and a floating structure

離散氣泡、浮體和自由液面邊界,邊界積分方程轉化為一個線性方程組,求解此線性方程組后可獲得氣泡和自由液面節(jié)點的法向速度與浮體濕表面節(jié)點速度勢,通過有限差分可獲得氣泡與自由液面的切向速度從而合成總速度[36].

對浮體附近脈動氣泡動力學特性的時域模擬需顯式更新邊界表面節(jié)點位置和速度勢.給出氣泡表面與自由液面滿足的運動學邊界條件如下

其中x為節(jié)點位置矢量.

氣泡表面和自由液面滿足的動力學邊界條件可由非定常伯努利方程給出,分別為

式中,p∞為z=0 平面內流場中無窮遠處靜水壓力,pb為氣泡內壓,g表示重力加速度,ρ 為流體密度.氣泡與自由液面滿足的動力學邊界條件用于更新節(jié)點速度勢.由于氣泡內部不可冷凝氣體在脈動過程中滿足理想氣體絕熱定律[37],忽略表面張力作用,pb可由下式計算

式中,pv為氣體的飽和蒸汽壓,本文取2337 Pa;p0為初始氣泡內部不可冷凝氣體壓力,V0為初始氣泡體積,V為瞬時氣泡體積;κ 為氣泡內部不可冷凝氣體的比熱比,本文取κ=1.4.

速度勢? 需滿足浮體濕表面不可穿透條件如下

式中,n表示邊界單位法向量,Vs表示浮體運動速度.

1.2 輔助函數(shù)法

通過牛頓第二定律

可求得浮體加速度a.式中m為浮體質量,G為浮體重力,空氣壓力Fa通過在浮體干表面上對壓力(大氣壓patm) 積分即可獲得,水動力Fh通過伯努利方程獲得[29]

傳統(tǒng)方法通過對相鄰時間步的節(jié)點速度勢采取向前差分近似求解??/?t項,但是時間步過小或浮體運動劇烈時,精度不足.理論上這種方法得到的是d?/dt,可通過物質導數(shù)與局部導數(shù)的關系得到修正的差分法即??/?t=d?/dt?Vs·??.本文采用的輔助函數(shù)法可精確求解結構響應,后面將對比修正后的傳統(tǒng)差分法、未修正的差分法和輔助函數(shù)法的計算精度.

本文中浮體僅存在z向運動,作用在浮體上的外力沿z軸,故可將式(8)改寫為

上式左側第一項代表慣性力項(其中積分項代表附加慣性力),左側第二項代表重力,χ 與ξ 為兩個輔助函數(shù),Sd為浮體與大氣接觸的表面.

在氣泡表面與自由液面,χ 滿足邊界條件[30]

在結構濕表面,χ 滿足邊界條件如下[38]

ξ 在氣泡表面與自由液面滿足的邊界條件可分別表示為[30]

ξ 在浮體濕表面滿足邊界條件[31]

輔助函數(shù)χ 與ξ 滿足拉普拉斯方程,并已知其在氣泡和自由液面滿足的狄利克萊邊界條件與在浮體濕表面上滿足的諾伊曼邊界條件,采用類似于求解邊界積分方程的方式可以得到在浮體濕表面上的值,代入式(10)中可求解出浮體加速度.

1.3 環(huán)狀氣泡模型

射流砰擊時,原來的單連通氣泡轉變?yōu)殡p連通氣泡,即變?yōu)榄h(huán)狀氣泡.此時氣泡射流砰擊處的速度勢不連續(xù).為保證計算的穩(wěn)定性,需要解決砰擊處節(jié)點的速度勢突跳問題.本文采用Wang 等[35]提出的渦環(huán)模型,在氣泡內部布置一個渦環(huán),其環(huán)量為砰擊處的速度勢突跳.將氣泡表面速度勢分解為[34]

其中,?vor為渦環(huán)誘導速度勢,?res為殘余速度勢.殘余速度勢?res在環(huán)狀氣泡表面連續(xù),采用邊界積分法進行計算,渦環(huán)誘導速度勢采用一種半解析法求解.詳細的計算過程可參考文獻[39],本文不再贅述,這里僅對渦環(huán)模型的速度勢分解思想進行介紹.

1.4 三相接觸模型

對于水面漂浮結構存在氣?液?固三相接觸線,如圖2 所示,接觸線需要特殊處理以同時滿足自由液面狄利克萊邊界條件與浮體表面諾伊曼邊界條件[40].本文采用雙節(jié)點法追蹤交界線,這種方法同樣適用于氣泡與結構接觸問題.將接觸點k分解為k1和k2兩點,分別屬于自由液面和柱體;使用下標fk1和sk2分別表示屬于自由液面、浮體的交界點.下標b,f 和s分別表示除了交界點外氣泡、自由液面和浮體上的節(jié)點.下標n表示速度勢的法向導數(shù).

圖2 三相接觸模型示意圖Fig.2 Sketch of the solid-liquid-gas contact model

離散邊界積分方程(2) 后可得到如下線性方程組

由于k1和k2兩點在空間位置上完全相同,影響系數(shù)矩陣H和G第3 和第4 行完全相同,整理可得

求解線性方程組后可得到k1點的法向速度.

當氣泡產生位置距浮體表面很近時,由于氣泡膨脹或者氣泡朝壁面遷移,導致氣泡附著到浮體表面.若氣泡節(jié)點與浮體節(jié)點距離小于單元尺度,傳統(tǒng)BIM 產生嚴重的數(shù)值不穩(wěn)定性;如圖3(a) 所示,部分氣泡節(jié)點穿過結構表面,產生了計算錯誤.因此本文將氣泡與浮體間的薄水膜去除,連接氣泡與浮體節(jié)點,進而保證數(shù)值計算的穩(wěn)定性,如圖3(b) 所示.Wang 等[41]的研究表明,這層薄水膜壓力趨近氣泡內壓,將其去除未對氣泡的動力學特性產生可見的影響.節(jié)點去除的范圍和接觸判定時刻由氣泡最小單元長度?s控制,當氣泡節(jié)點與浮體節(jié)點最小距離小于?s時,即認為氣泡即將附著在浮體表面,并將距離小于2.5?s[38]的氣泡與浮體節(jié)點刪除并連接氣泡與浮體,采用雙節(jié)點法處理產生的氣?液?固三相交界線.氣泡附著在浮體后,氣泡內壓直接作用在浮體表面上,此時浮體運動方程如下

式中,Sg為浮體與氣泡接觸表面.

圖3 傳統(tǒng)BIM 計算結果(a)與采用三相接觸模型計算結果(b)對比Fig.3 Comparison between the traditional BIM(a)and the gas-liquid-solid contact model(b)

1.5 無量綱化

本文將所有物理量進行無量綱化,并使用無量綱量進行計算與討論.取Rm,?p=p∞?pv和ρ 分別作為長度、壓力和密度的特征量,其他物理量如時間、速度、速度勢和加速度的無量綱特征量分別為Rm(ρ/?p)1/2,(?p/ρ)1/2,Rm(?p/ρ)1/2和?p/(ρRm)[34-35,37].本文后續(xù)出現(xiàn)的物理量如無特殊說明均為無量綱形式.

定義浮力參數(shù)δ

強度參數(shù)ε

初始氣泡中心位置與自由液面距離參數(shù)γf

初始氣泡中心位置與浮體距離參數(shù)γs

尺度比ζ

由于放電氣泡實驗中很難得到氣泡內壓,并且Turangan 等[42]的研究發(fā)現(xiàn)強度參數(shù)ε 在100 ～500 范圍內變化,氣泡射流速度變化很小,氣泡的動力學特性變化也很小.經過計算發(fā)現(xiàn),本文選取ε=300,實驗與數(shù)值結果吻合良好.計算中氣泡初始半徑與強度參數(shù)相關[43],本文選取R0=0.113.計算中所用時間步長由下式計算得到[37]

其中C為單步最大的速度勢改變量.

2 結果分析與討論

2.1 收斂性分析

為了驗證本文數(shù)值模型的正確性和有效性,首先進行收斂性分析,包括網格與時間步的收斂性分析.分別使用氣泡節(jié)點數(shù)量N等于80,120,160,200和240 的模型進行計算.圖4(a) 展示了不同網格密度下射流砰擊瞬時氣泡形態(tài),節(jié)點數(shù)為80,120 和160 的結果中氣泡形態(tài)具有一定的差異,當節(jié)點數(shù)N200 時氣泡輪廓基本重合,證明當氣泡網格節(jié)點取200 時能獲得收斂的結果.

本文采用四階Runge-Kutta 法對脈動氣泡與浮體相互耦合的時域運動過程進行模擬,顯式算法的收斂性依賴于時間步長.取C分別為0.04,0.02,0.01和0.005 分別進行計算.圖4(b) 描述了浮體加速度隨時間的變化,4 種方案取得了一致的結果,進一步證明了本文提出的數(shù)值模型的有效性.考慮程序計算的精度以及所需的計算成本,本文后續(xù)計算選取N=200,C=0.01.

圖4 模型的收斂性分析(a)計算網格的收斂性(b)計算時間步的收斂性Fig.4 Convergence study with respect to(a)the mesh and(b)the time step

2.2 典型工況對比分析

本節(jié)將分別對3 個典型距離參數(shù)γs下氣泡與浮體非線性耦合作用進行研究和討論,展示實驗中的物理現(xiàn)象,利用數(shù)值模擬分析現(xiàn)象背后的力學機理.本文實驗方法詳見已發(fā)表文獻[31,44].

2.2.1 弱耦合工況

第一個實驗中初始氣泡深度h=74.7 mm,初始浮體入水深度d=45.0 mm,浮體半徑Rs=16.0 mm,氣泡最大等效半徑Rm=16.8 mm.根據(jù)實驗設定計算初始參數(shù)為:ε=300,R0=0.113,κ=1.4,γf=4.45,γs=1.77,ζ=0.95,δ=0.04.圖5 所示為本工況數(shù)值與實驗結果對比圖,在此對本文的數(shù)值與實驗對比圖進行說明;圖中左側為高速攝影捕捉的氣泡圖像,右側為對應時刻氣泡運動模擬結果.每一張圖片左上方顯示了實驗中氣泡脈動時間,右側數(shù)字為數(shù)值計算結果對應的氣泡脈動時間,圖5(a) 左下方給出了比例尺,圖5(b) 中箭頭指出了實驗圖像中自由液面位置,為了更好地分析浮體與氣泡的非線性耦合效應,數(shù)值結果中還給出無量綱壓力與速度場.

圖5(a)中,氣泡在內部高壓的推動下猛烈膨脹.圖5(b) 所示為氣泡膨脹達到其最大體積時刻,此時氣泡近似球狀,自由液面和浮體未對氣泡脈動產生明顯的影響.由于過度膨脹氣泡周圍壓力比靜水壓力小,在接下來的脈動過程中在流體慣性的作用下氣泡開始坍塌.浮體對氣泡界面收縮具有延時作用,在氣泡坍塌前期,氣泡頂部略微凸起,而底部具有較大坍塌速度因此表面較為平坦,如圖5(c)所示.這種非球狀坍塌導致氣泡頂部在坍塌中期獲得更大的曲率,如圖5(d)所示.在氣泡坍塌后期,氣泡周圍流體壓力迅速增大,但氣泡頂部壓力梯度更大,使頂部界面向下(逆壓力梯度方向)快速收縮,如圖5(e)所示.第一個周期結束,并未形成射流,而在回彈階段氣泡頂部產生遠離浮體的射流并穿透另一側氣泡表面,如圖5(f) 所示.氣泡頂部射流現(xiàn)象與Lauterborn[45]的預測一致,即氣泡具有高曲率的部分更容易產生射流.此工況中氣泡產生背離浮體的射流(稱為“反射流”),說明浮體對氣泡的吸引力較小,無法誘導指向壁面的射流,故將本工況定義為弱耦合工況.總體上,數(shù)值計算得到的氣泡運動特性與實驗吻合良好.

圖5 弱耦合工況數(shù)值與實驗結果對比Fig.5 Comparison between the numerical results and the experimental images for the weakly coupling case

2.2.2 中度耦合工況

第二個實驗中氣泡產生位置距離浮體更近,耦合作用更強,氣泡初始深度h=60.0 mm,初始浮體入水深度d=42.6 mm,浮體半徑Rs=16.0 mm,氣泡最大等效半徑Rm=16.4 mm.依據(jù)實驗設定計算初始參數(shù)如下:ε=300,R0=0.113,κ=1.4,γf=3.66,γs=1.06,ζ=0.98,δ=0.04.

圖6 給出了本工況中幾個典型時刻數(shù)值與實驗結果對比.銅絲燃燒生成一個高壓氣核,在內部高壓氣體作用下迅速膨脹,如圖6(a) 所示.圖6(b)表示氣泡達到最大體積時刻,氣泡上方流體運動受到浮體的阻礙,氣泡上表面被壓平,此時氣泡上方流體繼續(xù)向上運動推動浮體遠離氣泡.在氣泡坍塌前期,遠離浮體的氣泡底部近似球形坍塌,如圖6(c)所示.柱體阻礙了流體從頂部流向氣泡中心,而對側方的流體的運動影響較小,因此形成了“梨”形氣泡.圖6(d)～圖6(f)描述了氣泡坍塌中后期的運動,氣泡兩側以及底部收縮的更快,頂部運動速度很小,“梨”形氣泡逐漸轉變?yōu)椤八蟆毙螝馀?氣泡在垂向被拉長.在坍塌后期氣泡底部壓力梯度變大,產生高壓區(qū),使底部氣泡表面快速收縮形成射流.圖6(f)的實驗圖像中氣泡底部留下銅絲燃燒產物,表明形成了一束向上的射流.此工況中浮體對氣泡的動力學行為影響更大,并能誘導氣泡產生指向壁面的射流(稱為“正射流”),因此將此工況命名為中度耦合工況,從實驗中也能觀察到浮體在膨脹階段被氣泡排斥而在收縮階段被吸引的現(xiàn)象.

圖6 中度耦合工況數(shù)值與實驗結果對比Fig.6 Comparison between the numerical results and the experimental images for the intermediately coupling case

2.2.3 強耦合工況

為研究浮體與脈動氣泡強烈的非線性耦合作用,第三個實驗進一步減小了浮體與氣泡間的距離.此時氣泡在膨脹階段附著到浮體底面,具有明顯的流固耦合效應,固將此工況命名為“強耦合工況”.實驗中參數(shù)如下;初始氣泡深度h=42.6 mm,浮體入水深度d=36.0 mm,浮體半徑Rs=14.0 mm,氣泡最大等效半徑Rm=13.1 mm.依據(jù)實驗設定計算初始參數(shù)如下:ε=300,R0=0.113,κ=1.4,γf=3.25,γs=0.5,ζ=1.07,δ=0.04.

圖7 給出了強耦合工況中幾個典型時刻實驗與數(shù)值計算結果對比.如圖7(a)所示,在氣泡膨脹初期,由于柱形浮體的阻礙作用,氣泡頂部變平.在膨脹階段氣泡上表面逐漸靠近浮體底面,需啟動三相接觸模型以保證數(shù)值穩(wěn)定性.圖7(b)所示為氣泡與浮體表面接觸時刻,氣泡即將附著到浮體表面.圖7(c)表示氣泡最大體積時刻,此時氣泡頂部仍有一定的膨脹速度,而底部與側面在收縮,進一步說明浮體對氣泡運動具有延時的作用.側面流體在氣泡坍塌時能自由的流向氣泡中心,頂部流體受到阻礙,如圖7(d)所示.在接下來的坍塌過程中,氣泡底部產生高壓區(qū),驅動射流直接砰擊在浮體上.實驗中還觀察到射流前進階段氣泡撕裂為頂部大氣泡和底部微小氣泡,計算中并未產生撕裂現(xiàn)象,如圖7(e)和圖7(f)所示.圖中最后兩幀數(shù)值與實驗結果具有較大偏差,這可能是由于氣泡脈動過程中浮體底面附著的空化氣泡(見圖7(c))與放電氣泡產生耦合作用所致.并且三相接觸線附近氣泡表面產生翻卷和破碎的現(xiàn)象,這可能與固體材料的潤濕性能相關,而本文計算中未考慮接觸角的影響.總體上氣泡形態(tài)與周期與實驗結果吻合良好.值得注意的是,柱形浮體在此強耦合工況中產生了顯著的運動.

圖7 強耦合工況數(shù)值與實驗結果對比Fig.7 Comparison between the numerical results and the experimental images for the strongly coupling case

圖8 給出了強耦合工況中浮體位移和浮體速度時歷曲線,分別對比了本文輔助函數(shù)法模型、傳統(tǒng)修正后的差分模型以及實驗結果.實驗結果中誤差棒的中心代表4 次測量的平均值,誤差棒的范圍代表標準差;菱形標記表示氣泡最大體積時刻.由圖8(a)可知,浮體在氣泡膨脹階段具有向上的運動,由于慣性作用氣泡達到最大體積后浮體仍繼續(xù)向上運動.在坍塌階段浮體向下加速運動再逐漸減速.可以看出,相比于傳統(tǒng)方法本文采用的輔助函數(shù)方法與實驗結果更加一致.在氣泡膨脹早期,數(shù)值結果中浮體運動較快,這是由于數(shù)值計算中忽略了銅絲燃燒放熱的過程;而在后期隨著時間步減小,傳統(tǒng)差分法結果逐漸偏離實驗結果.從圖8(b) 可以看出在氣泡膨脹前期,兩種方法求得的結構速度基本一致,而在后期傳統(tǒng)差分法精度降低,在氣泡附著于浮體表面后(對氣泡拓撲結構進行“數(shù)值切割”可能進一步增大了不穩(wěn)定性),速度時歷曲線出現(xiàn)劇烈震蕩(圖中虛線表示氣泡與浮體接觸時刻),而本文方法能得到穩(wěn)定且準確的結果.故對于研究近距離的脈動氣泡與浮體非線性耦合作用,本文提出的數(shù)值模型更為有效.

圖8 柱形浮體(a)位移與(b)速度時歷曲線(菱形標記代表最大體積時刻,虛線表示氣泡與浮體接觸時刻)Fig.8 Comparison of(a)the displacement and(b)the velocity histories of the cylindrical floating structure(the diamonds denote the moment of reaching the maximum bubble volume,the dashed line denotes the moment of the bubble-structure attachment)

2.3 無量綱參數(shù)討論

2.3.1 γs對氣泡射流模式的影響

影響氣泡與浮體非線性耦合效應的參數(shù)眾多,主要有無量綱氣泡與浮體距離參數(shù)γs、無量綱氣泡與自由面距離參數(shù)γf與尺度比ζ.眾多的影響參數(shù)導致了這一問題的復雜性.限于篇幅,本文僅研究γs的影響,其余參數(shù)與實驗參數(shù)一致,本文提出的數(shù)值模型也可研究其余參數(shù)的影響.依據(jù)實驗設定后續(xù)數(shù)值研究的初始參數(shù):ε=300,R0=0.113,γ=1.4,d=2.6,ζ=1,δ=0.04,γs=0.2 ～2.

圖9 展示了不同距離參數(shù)γs下,射流砰擊時刻的氣泡形態(tài),明顯氣泡射流模式分為5 種類型,分別為頸縮型環(huán)狀射流(0.2γs0.3)、接觸射流(0.4γs0.6)、非接觸射流(0.7γs1)、對射流(1.1γs1.3)和反射流(1.4γs2),圖中虛線代表浮體.

圖9 不同無量綱氣泡與浮體距離γs 下氣泡射流模式Fig.9 Jetting patterns for different dimensionless distances γs

圖10 給出了γs=0.2(頸縮型環(huán)狀射流)工況中典型時刻的壓力和速度場.氣泡收縮時牽引浮體側面流體向下運動,同時水平方向氣泡能自由收縮,如圖10(a)所示.水平方向流體與垂向流體在柱體底部最大半徑處匯聚,在此處形成高壓區(qū)(見圖10(b) 和圖10(c)).在壓力梯度作用下流體向斜下方運動,形成頸縮現(xiàn)象.類似于Cui 等[44]的分析,浮體的向下運動可能進一步改變流體運動,使其偏離形成環(huán)狀向外的射流導致氣泡撕裂為一個大單連通氣泡與小環(huán)狀氣泡.

圖10 γs=0.2 工況中典型時刻壓力與速度場(分別對應最大體積時刻t=0.910,氣泡坍塌時刻t=1.132,環(huán)向射流砰擊時刻t=1.403)Fig.10 Pressure and velocity fields at some typical moments for γs=0.2(corresponding times are t=0.910 at the moment of reaching the maximum bubble volume,t=1.132 during the bubble collapse stage and t=1.403 at the moment of the annular jet impact respectively)

圖9 第3 ～5 幀中氣泡附著在浮體表面,此時產生的射流稱為接觸射流.圖11 給出了γs=0.6(接觸射流) 工況中典型時刻壓力和速度場.在氣泡坍塌后期,底部附近形成了局部高壓區(qū)域,驅動氣泡界面快速收縮形成高速細射流直接砰擊在浮體上,如圖11(a)和圖11(b)所示.射流砰擊到浮體底部后,在浮體底面產生很高的壓力.流體受到壁面的作用變?yōu)閺较蛄鲃?最終導致氣泡脫離浮體表面,如圖11(c)所示.在氣泡脫離瞬時,浮體底面中心附近壓力比射流砰擊時刻更高,這是射流砰擊產生的局部高壓與氣泡收縮產生的脈動壓力疊加導致的結果.

圖11 γs=0.6 工況中典型時刻壓力與速度場(分別對應氣泡坍塌時刻t=1.792，射流砰擊時刻t=1.810 和氣泡與浮體脫離時刻t=1.838)Fig.11 Pressure and velocity fields at some typical moments for γs=0.6(corresponding times are t=1.792 during the bubble collapse stage,t=1.810 at the moment of the jet impact and t=1.838 at the moment of the bubble-floating structure detachment respectively)

圖9 第6 ～9 幀氣泡射流模式稱為非接觸射流.氣泡底部形成的高速射流首先砰擊在氣泡上表面,穿過水膜后砰擊在浮體底面.圖12 給出了γs=0.9(非接觸射流)工況中典型時刻壓力和速度場.在坍塌階段氣泡下方壓力梯度最大,導致正射流的形成.射流砰擊在氣泡上表面后,在氣泡與浮體間形成局部高壓區(qū)域,如圖12(a)和圖12(b)所示.環(huán)狀氣泡繼續(xù)收縮同時向浮體運動,在頂部形成突出物最終與浮體底面接觸.值得注意的是射流砰擊后浮體底面中心附近壓力比上文接觸射流工況更大.環(huán)狀氣泡繼續(xù)坍塌,壁面附近壓力超過140 MPa,如圖12(c)所示.

圖12 γs=0.9 工況中典型時刻壓力與速度場(分別對應氣泡坍塌時刻t=1.806,射流砰擊時刻t=1.807 和氣泡與浮體接觸時刻t=1.809)Fig.12 Pressure and velocity fields at some typical moments for γs=0.9(corresponding times are t=1.806 during the bubble collapse stage,t=1.807 at the moment of the jet impact and t=1.809 at the moment of the toroidal bubble-floating structure attachment respectively)

圖9 第10 ～12 幀展示了反射流與正射流相撞的現(xiàn)象(稱為“對射流”).如圖13(a)所示,在氣泡坍塌末期,氣泡頂部首先產生一個寬度較大的凹陷.氣泡底部形成很高的曲率,在其下方產生一個高壓區(qū)域.高壓區(qū)使氣泡底部迅速收縮形成一束指向壁面的正射流.最終正射流與寬而慢的反射流相撞(圖13(b)和圖13(c)).射流相互撞擊后在撞擊區(qū)域附近產生了很高的壓力.反射流內部未與正射流砰擊的部分繼續(xù)向下運動,導致砰擊區(qū)域環(huán)狀凹陷的形成,如圖13(d)所示.這種由射流相互撞擊產生的環(huán)狀凹陷無法發(fā)展成射流,在回彈階段逐漸消失,如圖13(e) 和圖13(f) 所示.在坍塌和回彈階段,氣泡和高壓區(qū)逐漸向浮體遷移.

圖13 γs=1.3 工況中典型時刻壓力與速度場(分別對應氣泡坍塌時刻t=1.806,t=1.807,射流砰擊時刻t=1.809,氣泡坍塌時刻t=1.810 和氣泡回彈時刻t=1.820,t=1.853)Fig.13 Pressure and velocity fields at some typical moments for γs=1.3(corresponding times are t=1.806,t=1.807 during the bubble collapse stage;t=1.809 at the moment of the jet impact;t=1.810 during the bubble collapse stage;t=1.820,t=1.853 during the bubble rebounding stage respectively)

圖9 第13 幀中氣泡產生了反向射流與徑向射流,徑向射流砰擊將氣泡撕裂為兩個單連通的氣泡.圖9 第14 ～15 幀中氣泡只產生徑向凹陷,無法穿透氣泡.圖14 給出了γs=1.5(反射流) 工況中典型時刻壓力和速度場.在氣泡坍塌后期,氣泡頂部壓力比周圍壓力更高導致氣泡形成向下的凹陷.然而第一個周期內,向下的凹陷無法形成射流穿透氣泡表面,如圖14(a)所示.隨后氣泡兩側收縮較快形成徑向的氣泡凹陷.在氣泡回彈階段,上方凹陷發(fā)展為一束向下的反射流.由于氣泡膨脹射流變得越來越細而徑向凹陷開始消失,如圖14(b)所示.在反射流尖端形成了非常奇特的壓力分布,即在射流尖端具有較高的壓力,而在兩側卻形成較低的環(huán)形壓力區(qū)域.反射流與下方氣泡界面砰擊后,在氣泡底部形成了局部高壓,而其余方向壓力由于氣泡過度膨脹降至較低水平.此時氣泡呈現(xiàn)為“蝴蝶”形,如圖14(c)所示.

圖9 第16 ～19 幀中氣泡僅產生反射流.圖15 給出了γs=1.9 工況中典型時刻的壓力與速度場.與圖14 工況相同,在第一個脈動周期內射流無法穿透氣泡表面,而在回彈階段射流變細并加速向下運動,如圖15(a)和圖15(b)所示.射流砰擊產生后氣泡下方形成局部高壓,并且形成了向下的突出物,如圖15(c)和圖15(d)所示.

圖14 γs=1.5 工況中典型時刻壓力與速度場(分別對應氣泡坍塌時刻t=1.803,氣泡回彈時刻t=1.833,t=1.860)Fig.14 Pressure and velocity fields at some typical moments for γs=1.5(corresponding times are t=1.803 during the bubble collapse stage;t=1.833 and t=1.860 during the bubble rebounding stage respectively)

圖15 γs=1.9 工況中典型時刻壓力與速度場(分別對應氣泡最小體積時刻t=1.816,氣泡回彈時刻t=1.825,t=1.837,t=1.891)Fig.15 Pressure and velocity fields at some typical moments for γs=1.9(corresponding times are t=1.816 at the moment of reaching the minimum bubble volume;t=1.825,t=1.837,and t=1.891 during the bubble rebounding phase respectively)

圖15 γs=1.9 工況中典型時刻壓力與速度場(分別對應氣泡最小體積時刻t=1.816,氣泡回彈時刻t=1.825,t=1.837,t=1.891)(續(xù))Fig.15 Pressure and velocity fields at some typical moments for γs=1.9(corresponding times are t=1.816 at the moment of reaching the minimum bubble volume;t=1.825,t=1.837,and t=1.891 during the bubble rebounding phase respectively)(continued)

當γs很小時(γs0.3),浮體底面兩側產生的局部高壓使氣泡產生頸縮,最后導致氣泡撕裂.γs增大至(0.4γs1.3),氣泡與浮體具有較強的流固耦合效應.由于Bjerknes 力[46-47]作用,在坍塌階段氣泡產生指向浮體的正射流,而當氣泡與浮體較遠時(1.1γs2.0),在坍塌階段氣泡頂部受到浮體的延時作用具有較高曲率并且在氣泡上方壓力梯度更大,使氣泡上表面產生向下的反射流.自由液面對氣泡具有排斥作用,進一步促進了反射流的形成[48-49].

2.3.2 γs對射流速度的影響

由于本文實驗裝置和條件的限制,很難精確測量射流速度,為驗證本文模型能預測射流速度,首先對比本文模型計算結果與Cui 等[44]放電氣泡與不同厚度浮冰耦合作用的實驗.圖16 所示為最大射流速度Vjet(正射流)與浮冰厚度的關系(γs=0.93).由圖可知,當浮冰厚度大于28.1 mm 時,計算所得最大射流速度基本在實驗測量范圍內,而隨著浮冰厚度減小或者射流速度增大,兩者偏差增大,最大誤差為16.8%(定義為實驗與計算偏差/實驗值).實驗中由于氣泡壁面作為一個“凹透鏡”折射光線導致測量的射流速度可能偏小[43,50].并且當射流速度很大時,需要極高的拍攝幀率才能精確測量最大射流速度值,Cui等[44]的實驗中拍攝幀率可能不足.這兩種因素導致了數(shù)值計算結果相比實驗結果偏大.浮冰厚度增大,射流速度減小,并且數(shù)值與實驗結果吻合程度更高,一定程度上驗證了本文計算模型能較準確的預測射流速度.本文后續(xù)進一步研究γs對脈動氣泡與浮體非線性耦合效應的影響.

圖16 最大射流速度數(shù)值與Cui 等[44] 實驗結果對比Fig.16 Comparison of the maximum jet velocity between the present model and the experiment[44]

圖17 給出了3.3.1 節(jié)中工況(γs=0.4 ～2)最大射流速度Vjet和最大反射流速度Vcounter與γs的關系.正射流速度遠高于反射流且隨γs先增大后減小再增大,在γs=0.8 時達到最大值136.8(1368 m/s);反射流速度隨γs增大而增大,最大值為27.2(272 m/s).圖12 給出的γs=0.9 工況中,高速細射流砰擊在另一側氣泡表面,導致浮體底面附近具有很高的壓力(超過140 MPa).而當0.7γs0.9 時,射流速度均大于100,可能在浮體底部產生巨大的射流砰擊載荷,造成結構局部損傷.與Lechner 等[51]對超近壁面附近脈動氣泡動力學特性的數(shù)值研究結果相似,在氣泡坍塌階段與浮體底面平行的流動比垂直于底面的流動更快,導致氣泡在遠端形成很高的曲率,從而形成了速度約1000 m/s 的高速細射流(Lechner 等[51]發(fā)現(xiàn)的射流速度超過2000 m/s).本文數(shù)值模型中最大射流速度遠大于Philipp 和Lauterborn[50]通過激光氣泡與固定剛性板耦合作用實驗測量的射流速度(150 m/s),可見浮體與脈動氣泡的非線性流固耦合效應能顯著影響氣泡的動力學特性,在實際應用中需要考慮這種耦合作用的影響.但實際中如此高速的細射流可能存在不穩(wěn)定性[51],并且流體可壓縮性具有重要影響,這不在本文的討論范圍內.

圖17 射流速度與無量綱氣泡與浮體距離γs 的關系Fig.17 Maximum jet velocity versus the dimensionless bubble-floating structure distance γs

2.3.3 γs對氣泡遷移和結構運動的影響

圖18 展示了不同γs下氣泡形心的位移時歷曲線(計算至射流穿透前).當γs<1.5 時,在氣泡膨脹階段,由于浮體對氣泡上方流體的運動具有很強的阻礙,導致氣泡下表面相比上表面膨脹得更快,類似于氣泡被浮體推開,形心遠離結構.而在坍塌階段氣泡受到浮體Bjerknes 吸引[46]作用強于自由液面的Bjerknes 排斥[35,48]作用,氣泡又被浮體吸回.在氣泡坍塌階段近浮體與遠離浮體側不平衡壓力使氣泡整體受到一個朝向壁面的力,從而導致了氣泡朝浮體遷移,如圖6(d)～圖6(f)所示.由圖18 可知隨著γs的增大,這種相互排斥、吸引作用減小.當γs1.5 時(此時γf4.1),氣泡在膨脹階段基本保持其形心位置不變,而在坍塌和回彈階段氣泡逐漸遠離浮體.這可能是因為自由液面對氣泡的排斥力比浮體的吸引力更強.

圖18 不同γs 下氣泡形心的位移時歷曲線Fig.18 Time histories of the displacement of the bubble centroid for different γs

圖19 浮體最大速度與無量綱氣泡與浮體距離γs 的關系Fig.19 Maximum velocity of the floating body versus the dimensionless bubble-floating structure distance γs

3 結論

本文采用邊界積分法聯(lián)合輔助函數(shù)法,建立了高壓脈動氣泡與浮體非線性流固耦合數(shù)值模型,同時開展了柱形浮體與放電氣泡流固耦合作用實驗,對比數(shù)值與實驗結果,二者吻合良好.采用該數(shù)值模型,本文對氣泡與浮體無量綱距離參數(shù)γs進行了系統(tǒng)研究,得到的主要結論如下;

(1)將本文流固耦合模型、傳統(tǒng)流固耦合模型(差分法計算速度勢偏導數(shù))與實驗結果進行對比分析,發(fā)現(xiàn)傳統(tǒng)方法在處理氣泡與浮體接觸工況時,精度明顯不足,而本文的輔助函數(shù)法大大提高了計算穩(wěn)定性與精度.

(3) 正射流速度隨著γs增大先增大后減小再增大,而反射流速度隨γs增大而增大.正射流速度能達到約1000 m/s 的量級,說明脈動氣泡與浮體的流固耦合效應在實際應用中影響重大,其形成機理為氣泡底部高曲率與浮體的Bjerknes 吸引力聯(lián)合作用.該現(xiàn)象需要未來更精密的實驗來驗證.

(4) 總體而言,在氣泡膨脹階段,浮體排斥氣泡,氣泡形心向下運動;在坍塌階段,當0.2γs1.4 時,浮體對氣泡的Bjerknes 吸引力強于自由液面Bjerknes 排斥力,導致氣泡向浮體遷移;當1.5γs2,浮體對氣泡的Bjerknes 吸引力比自由液面Bjerknes 排斥力弱,導致氣泡遠離自由液面.