基于雙向流固耦合模型的潰壩水流數(shù)值模擬

王春正, 孫亮*, 黃玲玲, 陳麗芬

(1. 武漢理工大學船海與能源動力工程學院, 武漢 430063; 2. 大連理工大學港口海岸及近海工程國家重點實驗室, 大連 116024; 3. 大連市海洋可再生能源重點實驗室, 大連 116024)

流固耦合是一類多學科交叉融合的問題,廣泛存在于工程、科學和自然界中。例如,液化天然氣船在運輸過程中遭遇的液艙晃蕩問題;海洋立管在復雜海洋環(huán)境下發(fā)生的渦激振動問題;風電機組葉片在運轉(zhuǎn)過程中面臨的氣動彈性問題;超大型浮式結(jié)構(gòu)物在復雜波浪條件下產(chǎn)生的水彈性響應問題以及盾構(gòu)隧道施工[1]等場景。這些問題的共性是流體作用于結(jié)構(gòu)而引起彈性結(jié)構(gòu)發(fā)生位移或者變形,同時結(jié)構(gòu)的動力響應又對流體的運動產(chǎn)生影響。

流固耦合問題較為復雜且涉及結(jié)構(gòu)變形破壞、流體劇烈變化等因素,往往難以通過嚴格的公式推導獲得解析解,因此,模型試驗和數(shù)值模擬成為研究流固耦合問題的兩種常用方法。然而模型試驗受試驗成本、儀器精度、環(huán)境因素干擾的制約,常難以獲得較為精確的結(jié)果。與模型試驗相比,數(shù)值模擬具有成本低、周期短、物理過程清晰等優(yōu)勢,從而成為解決流固耦合問題的主流。

對于流固耦合問題的數(shù)值模擬,根據(jù)是否在同一組控制方程中求解流體域和固體域,可以將研究方法分為強耦合和弱耦合兩種。其中強耦合方法是在同一組控制方程中求解,又被稱為直接耦合法。Walhorn等[2]提出了時空有限元法來分析具有自由液面流動的流固耦合問題,并對潰壩水流沖擊彈性結(jié)構(gòu)問題進行了數(shù)值模擬;Idelsohn等[3]提出了一個以統(tǒng)一形式處理流體子域和固體子域的拉格朗日公式,并基于PFEM(particle finite element method,粒子有限元法)對若干流固耦合問題進行了數(shù)值模擬;Rafiee等[4]提出了一種基于光滑流體動力學(smoothed particle hydrodynamics,SPH)的數(shù)值方法,并對彈性板在水壓下變形和潰壩水流沖擊彈性結(jié)構(gòu)等問題進行了模擬;Sun等[5]基于材料點法(material point method,MPM)同樣提出了一種模擬流固耦合的方法,并對液倉晃蕩、彈性水壩溢流等問題進行了數(shù)值模擬。

弱耦合方法是在不同組控制方程中分別求解流體和固體域,又被稱為間接耦合法。Cerquaglia等[6]基于PFEM和有限元法(finite element method,FEM)相耦合,以分區(qū)的方法對鳥類撞擊、潰壩等流固耦合問題進行了數(shù)值模擬;張友林等[7]基于移動粒子半隱式方法(moving particle semi-implicit method,MPS)和FEM,開發(fā)出了MPS-FEM耦合求解器,并使用該求解器數(shù)值研究了潰壩流動對彈性結(jié)構(gòu)的沖擊問題;張智瑯等[8]基于解耦有限粒子法(decoupled finite particle method,DFPM)和光滑有限元(smoothed finite element method,SFEM)相耦合,以分區(qū)的方法對結(jié)構(gòu)入水,潰壩水流沖擊彈性結(jié)構(gòu)等流固耦合問題進行了數(shù)值模擬;趙冉等[9]基于格子玻爾茲曼法(lattice Boltzmann method,LBM)和浸入邊界法(immersed boundary method,IBM)相耦合的方法,對重力場內(nèi)單根柔性纖維與均勻流場的耦合行為進行了數(shù)值模擬;姚學昊等[10]提出了一種基于SPH和虛粒子和排斥力的近場動力學(peridyna-mics,PD)相耦合的方法,利用流體粒子-虛粒子接觸算法處理流-固界面,并使用該方法對彈性結(jié)構(gòu)在水壓下變形和潰壩沖擊彈性結(jié)構(gòu)等問題進行分區(qū)求解。

強耦合方法因流固控制方程組同時收斂較為困難,故難以針對復雜的物理問題建立統(tǒng)一的控制方程組且消耗計算資源巨大,而弱耦合方法使用不同的控制方程分離求解流體域和固體域,因此,可以組合使用各種成熟的流固求解器,靈活分析復雜的流固耦合問題,對計算機性能的需求也大幅降低。故而,本文研究采用弱耦合的方法對流固耦合問題進行分析。

長期以來,有限體積法(finite volume method,FVM)在解決帶有自由表面的兩相流問題中占據(jù)了主導地位,FEM可以精確、穩(wěn)定、高效地處理結(jié)構(gòu)的變形問題。現(xiàn)使用基于FVM的開源計算流體動力學(computational fluid dynamics,CFD)軟件OpenFOAM和基于FEM的開源計算固體力學(computational solid mechanics,CSM)軟件CalculiX,通過開源耦合庫preCICE對雙向流固耦合問題展開分析。首先對代數(shù)流體體積法(volume of fluid,VOF)和幾何VOF兩種自由液面捕捉的方法進行對比,隨后使用幾何VOF方法對潰壩水流沖擊彈性結(jié)構(gòu)的雙向流固耦合問題進行模擬,通過與已有的數(shù)據(jù)進行對比,驗證本文所使用分析方法的可靠性,為今后研究相關(guān)問題提供參考。

1 數(shù)學模型

1.1 FVM方法

在本文所采用的方法中,將氣體和液體視為單一可變密度的連續(xù)流體,并考慮黏性,其控制方程和自由液面捕捉的VOF方法如下所述。

1.1.1 控制方程

開源CFD軟件OpenFOAM采用有限體積法求解連續(xù)方程[式(1)]和N-S方程[式(2)][11]。

(1)

(2)

式中:U為速度矢量;ρ為流體密度;μ為動力黏性系數(shù);g為重力加速度;p為流體壓力;fσ為表面張力;t為時間。

1.1.2 VOF方法

在使用OpenFOAM中的兩相流求解器分析潰壩水流沖擊結(jié)構(gòu)物時,需要追蹤流體區(qū)域內(nèi)液體和氣體的交界面,即自由液面。Hirt等[12]提出的VOF方法根據(jù)每個時刻網(wǎng)格中流體體積的變化來追蹤自由液面,是目前使用最廣泛的自由液面追蹤方法。在分析中引入流體體積分數(shù)α表示流體計算域內(nèi)氣體和液體的分布:當α=1時,表示流體區(qū)域內(nèi)全部為液體;當α=0時,表示流體區(qū)域內(nèi)全部為氣體;當0<α<1時,表示流體區(qū)域內(nèi)既有氣體也有液體。根據(jù)求解方法的不同,VOF方法可以分為代數(shù)VOF方法和幾何VOF方法。

(1)代數(shù)VOF方法。代數(shù)VOF方法為獲取更小的界面厚度,引入了人工壓縮項來建立體積分數(shù)運輸方程[13],即

(3)

(2)幾何VOF方法。幾何VOF方法(isoAdvector方法)是Roenby等[14]提出的一種新型幾何VOF方法。在每個分析步內(nèi),該方法首先會在氣體和液體交界的網(wǎng)格單元中,找到Isoface(最適合劃分氣體和液體的平面)[14],并基于這些Isoface對自由液面進行幾何重構(gòu)。之后,會使用這些重構(gòu)的自由液面對流體體積分數(shù)場進行更新,下一分析步的流體體積分數(shù)場由式(4)[14]給出。

(4)

(5)

式中:αi為第i網(wǎng)格處的流體體積分數(shù);Vi為第i個網(wǎng)格處的體積;Bi為組成網(wǎng)格單元i的所有面集合;Sij為確定法向量方向的輔助因子;Sj為網(wǎng)格表面的面積法向量;Aj(τ)為τ時刻網(wǎng)格表面Fj處的瞬時淹沒面積;φj(t)為通過Fj處的容積通量;x為笛卡爾坐標;u(x)為網(wǎng)格表面Fj中x處的流速矢量。

1.2 FEM方法

使用基于FEM的開源CSM軟件CalculiX[15]對結(jié)構(gòu)場分析,其控制結(jié)構(gòu)單元變形的動力學方程由式(6)給出。

(6)

C=α1M+α2K

(7)

式中:M為結(jié)構(gòu)的質(zhì)量矩陣;C為瑞麗阻尼矩陣;K為結(jié)構(gòu)剛度矩陣;F(t)為施加在結(jié)構(gòu)上隨時間t變化的外力;y為結(jié)構(gòu)單元的位移矢量;系數(shù)α1和α2的取值與結(jié)構(gòu)的阻尼比和固有頻率有關(guān)。

(8)

(9)

(10)

式中:參數(shù)α、β、γ控制α-Method求解的精度和穩(wěn)定性。在CalculiX中參數(shù)β、γ由參數(shù)α表示。

(11)

(12)

1.3 FVM-FEM耦合方法

圖1 雙向流固耦合模型

2 算例分析

2.1 潰壩水流沖擊剛性結(jié)構(gòu)

為了比較代數(shù)VOF方法和幾何VOF方法捕捉自由液面在實際應用中的效果,首先分別使用基于代數(shù)VOF法的interFoam求解器和基于幾何VOF法的interIsoFoam求解器對該潰壩水流沖擊剛性結(jié)構(gòu)的過程[18]進行研究。如圖2所示,初始時刻左側(cè)區(qū)域有寬0.146 m、高0.292 m的水柱,在下壁面的中間位置有一個寬0.024 m、高0.048 m的剛性結(jié)構(gòu),A點位于剛性結(jié)構(gòu)的左側(cè)中點(距離底邊界0.024 m)。在該算例中,空氣的密度為1 kg/m3,液體的密度為1 000 kg/m3。計算過程中經(jīng)過收斂性分析最終使用的網(wǎng)格大小為0.002 m,時間的步長設置為自適應時間步長,該時間步長由最大庫朗數(shù)控制,本文采用的最大庫朗數(shù)為0.5。

圖2 潰壩水流沖擊剛性結(jié)構(gòu)

圖3給出了分別使用基于代數(shù)VOF法和幾何VOF法的求解器進行數(shù)值模擬,得到的不同時刻(t=0.2、0.3、0.4、0.5 s)液體的界面形態(tài)與試驗[18]的對比。從圖3中可以看出,隨著坍塌的水流抨擊剛性結(jié)構(gòu),使用代數(shù)VOF方法捕捉的自由液面可以觀察到非常明顯的液體體積損失,并且由于碰撞本應該出現(xiàn)的飛濺液珠和復雜尖銳的自由液面沒有被精確捕捉,反而因為數(shù)值模擬的誤差部分水體消失,而使用幾何VOF法可以將自由液面很好地捕捉并且沒有觀察到有明顯的水體損失。因此,幾何VOF方法對于劇烈變化的自由液面具有更好的捕捉能力。

圖3 水流沖擊剛性結(jié)構(gòu)物的自由表面分布

圖4為基于代數(shù)VOF方法和幾何VOF方法得到的計算域內(nèi)水體質(zhì)量的變化情況。初始時,積分得到水體的質(zhì)量為0.085 264 kg,與計算得到的理論值相同,即初始誤差為0。在3.0 s時,使用幾何VOF法得到的水體質(zhì)量穩(wěn)定在0.084 71 kg,誤差為0.649 7%,而使用代數(shù)VOF法得到的水體質(zhì)量穩(wěn)定在0.083 48 kg,誤差為2.092 3%。因此,基于幾何VOF方法得到的結(jié)果具有更好的質(zhì)量守恒性。

圖4 計算域中液體質(zhì)量的變化

圖5給出了使用兩種VOF方法模擬左邊壁上的液面變化。從圖5可以看出,與試驗值[18]和Issakhov等[19]的數(shù)值模擬值對比,使用兩種VOF方法對左邊壁上液面的模擬結(jié)果沒有較大差別,原因在于左邊壁處液體的運動比較平緩,并沒有出現(xiàn)碰撞、飛濺等強非線性現(xiàn)象,這也說明了對于比較緩慢變化的沖擊過程幾何VOF方法和代數(shù)VOF方法均可以給出高精度的預測值。

圖5 左側(cè)壁處的液面變化

圖6給出了使用兩種VOF方法模擬圖2中A點處的壓力隨時間變化的對比。從圖6可以看出當水柱突然坍塌沖擊剛性結(jié)構(gòu)的A點時,該處的壓力突然增大并出現(xiàn)最大峰值,之后剛性障礙物左側(cè)的水體會繼續(xù)沖擊A點,并在0.6 s內(nèi)逐漸形成較為明顯的三個峰值。使用幾何VOF方法所得的結(jié)果與Issakhov等[19]的結(jié)果基本吻合,而使用代數(shù)VOF方法所得的結(jié)果與Issakhov等[19]的結(jié)果差別較大。這在一定程度上與圖3中觀察到的自由液面破碎現(xiàn)象以及圖4中使用基于代數(shù)VOF法模擬會出現(xiàn)較大的液體質(zhì)量損失有關(guān)。

圖6 A點處壓力變化

以上的分析表明:基于代數(shù)VOF方法的流體求解器在模擬潰壩水流沖擊等強非線性問題時的精度要低于基于幾何VOF方法所得到的相應結(jié)果(自由表面分布、壓力和質(zhì)量守恒)。在分析雙向流固耦合問題時,對流體區(qū)域壓力的較大計算誤差會影響結(jié)構(gòu)的響應,而結(jié)構(gòu)響應的分析誤差又會影響流體的運動,這些誤差會在耦合過程中不斷累積,從而導致最終的模擬結(jié)果產(chǎn)生較大的差異。為保證雙向流固耦合的分析精度,進一步的研究中將使用具有更強液面捕捉能力、質(zhì)量守恒性更好的幾何VOF方法來模擬雙向流固耦合問題。

2.2 潰壩水流沖擊彈性結(jié)構(gòu)

為驗證本文基于FVM-FEM方法在模擬帶有自由液面的雙向流固耦合問題時的準確性,同時進一步揭示此類流固耦合問題的內(nèi)在機理,對潰壩水流沖擊彈性結(jié)構(gòu)問題[2]進行研究。如圖7所示,水箱寬0.584 m,水箱左下腳存有寬0.146 m,高0.292 m的水柱,在下壁面的中間位置有一個寬0.012 m、高0.08 m的彈性結(jié)構(gòu)。在該算例中,空氣的密度為1 kg/m3,液體的密度為1 000 kg/m3,彈性結(jié)構(gòu)的密度2 500 kg/m3,彈性模量為1 MPa,泊松比為0。

圖7 潰壩水流沖擊彈性結(jié)構(gòu)

圖8給出了本文方法(FVM-FEM)與Walhorn等[2]、Idelsohn等[3]、Rafiee等[4]所計算得出的彈性結(jié)構(gòu)左上角水平位移的對比。從圖8可以看出,使用本文的方法FVM-FEM與其他3種方法所得的結(jié)果曲線,在總體趨勢上具有較好的一致性,且在關(guān)鍵點處的數(shù)值基本吻合。在0.14 s左右,4種方法的結(jié)果中均出現(xiàn)向左的水平位移(出現(xiàn)負值),本文結(jié)果為0.068 cm;在0.24 s左右時,4種方法所模擬的彈性障礙物的水平位移均達到向右的最大值,本文結(jié)果為4.88 cm,與Idelsohn等[3]的結(jié)果幾乎完全一致;在0.65 s左右,彈性障礙物的水平位移達到第二個峰值,本文結(jié)果為1.68 cm,與Walhorn等[2]和Idelsohn等[3]的結(jié)果吻合較好,優(yōu)于Rafiee等[4]的結(jié)果,驗證了本文方法的可靠性;在第二個峰值之后,4種數(shù)值模擬方法得到的水平位移曲線具有較為明顯的差別,這是因為潰壩水流沖擊彈性結(jié)構(gòu)是一個流體和固體相互耦合的過程,同時液體區(qū)域具有更強的非線性,求解難度更大,對求解的精確度和穩(wěn)定性要求高,但可以看到的是本文的方法與其他3種方法所得的結(jié)果在總體趨勢上保持一致,再次驗證了本文所采用的方法在模擬具有自由液面的雙向流固耦合問題時具有較好的可靠性。

圖8 彈性結(jié)構(gòu)左上角的位移變化

圖9給出了不同方法對于彈性結(jié)構(gòu)變形和流體分布在不同時刻的比較?？梢钥闯?在t=0.14 s時,彈性結(jié)構(gòu)下端受到水流的沖擊作用,結(jié)構(gòu)下部向右彎曲,頂部出現(xiàn)向左的偏轉(zhuǎn),因此結(jié)構(gòu)左上角的位移為負值(圖8);在t=0.26 s時,水流完全覆蓋彈性結(jié)構(gòu)左側(cè),出現(xiàn)明顯的“射流”現(xiàn)象,此時彈性結(jié)構(gòu)具有最大向右的彎曲幅度,同時在計算域內(nèi)可以觀察到飛濺的液珠和尖銳的自由液面形狀,這與參考文獻[3- 4]中的結(jié)果非常接近,本文方法在捕捉劇烈變化的自由液面時具有較高的精度;t=0.34 s時,水體越過彈性結(jié)構(gòu)的抨擊右側(cè)壁面,彈性結(jié)構(gòu)開始出現(xiàn)向左的回彈;t=0.42 s時,水流越過彈性結(jié)構(gòu)繼續(xù)砰擊右側(cè)壁面,形成較為明顯的“分流”現(xiàn)象,部分液體沿右側(cè)壁面開始回流,彈性結(jié)構(gòu)繼續(xù)向左回彈。

圖9 不同時刻的彈性板變形和水流形態(tài)

圖10給出了不同時刻(對應圖8中位移的峰值和谷值,即t=0.14、0.23、0.53、0.65、0.84、0.94 s)液體壓力和彈性板應力的分布。從圖10可以看出,t=0.14 s時,潰壩水流沿著彈性結(jié)構(gòu)左側(cè)向上流動,水流受到彈性結(jié)構(gòu)的阻礙作用,匯聚在彈性結(jié)構(gòu)的左下角,水體在該處產(chǎn)生壓力的最大值(約為4.2 kPa),彈性結(jié)構(gòu)在中部偏下向右側(cè)彎曲,而彈性結(jié)構(gòu)的頂端受到的壓力很小,因而發(fā)生了向左側(cè)的較小偏移(圖8);t=0.23 s時,潰壩水流越過彈性結(jié)構(gòu)頂端沖出去,形成“射流”現(xiàn)象,此時彈性結(jié)構(gòu)的頂端出現(xiàn)向右水平位移的最大值,同時彈性結(jié)構(gòu)左下角的壓力因為液體的流動而變小,彈性結(jié)構(gòu)左右兩側(cè)的應力均出現(xiàn)了整個過程中的最大值,約為221 kPa和-221 kPa,此時也是彈性結(jié)構(gòu)最可能發(fā)生破壞的時刻;t=0.53 s時,由于水流與右側(cè)墻壁碰撞出現(xiàn)“分流”現(xiàn)象,一部分水體沿右側(cè)的水槽底邊壁出現(xiàn)回流,在彈性結(jié)構(gòu)右下角出現(xiàn)液體壓強的較大區(qū)域,此時彈性結(jié)構(gòu)的頂端出現(xiàn)向左回彈;t=0.65 s時,右側(cè)的水體沿彈性結(jié)構(gòu)右側(cè)邊壁爬升,與結(jié)構(gòu)左側(cè)流動的水體產(chǎn)生“對流現(xiàn)象”,在彈性結(jié)構(gòu)頂端形成了一個漩渦,彈性結(jié)構(gòu)物開始晃動,解釋了圖8中本文方法得到的位移曲線在t=0.65 s附近出現(xiàn)輕微震蕩的原因,而其他3種方法沒有捕捉到該現(xiàn)象,同時可以觀察到在彈性結(jié)構(gòu)的頂端形成局部壓力較大的區(qū)域;t=0.84 s和t=0.94 s時,受彈性結(jié)構(gòu)運動的影響,左側(cè)的水體被不斷擠壓,右側(cè)墜落的水體和回流的水體開始發(fā)生沖撞,并且這種不規(guī)則非線性作用不斷加強,導致圖8中4種方法所得的位移曲線在t=0.65 s以后存在較大的差異。

圖10 不同時刻流體的壓力和彈性結(jié)構(gòu)的應力分布

圖11給出了具有相同布置的彈性結(jié)構(gòu)(圖7)和剛性結(jié)構(gòu)在結(jié)構(gòu)左側(cè)中點處壓力值的對比。從圖11可以看到,在水流第一次沖擊結(jié)構(gòu)時,由于彈性結(jié)構(gòu)產(chǎn)生彎曲,相比于剛性結(jié)構(gòu)減小了水流的沖擊作用,所以彈性結(jié)構(gòu)左側(cè)中點處所受到的壓力更小。圖11中彈性結(jié)構(gòu)左側(cè)中點處壓力值震蕩較大,最終趨向于剛性結(jié)構(gòu)左側(cè)中點處的壓力值。

圖11 剛性結(jié)構(gòu)和彈性結(jié)構(gòu)左側(cè)中點處的壓力

3 結(jié)論

基于雙向流固耦合模型對潰壩水流進行了數(shù)值分析,研究中使用了開源CFD求解器OpenFOAM和開源CSM求解器CalculiX,得出如下結(jié)論。

(1)分別采用基于代數(shù)VOF方法的interFoam求解器和基于幾何VOF方法的interIsoFoam求解器模擬了潰壩水流沖擊剛性結(jié)構(gòu)問題,通過比較液體的流態(tài)、質(zhì)量守恒性、指定位置處的壓力和液面高度,表明幾何VOF方法在模擬潰壩沖擊等強非線性、具有復雜劇烈變化的自由液面問題時具有更好的精度。

(2)基于FVM-FEM的弱耦合方法,并使用幾何VOF方法對潰壩水流沖擊彈性結(jié)構(gòu)問題進行了數(shù)值模擬,所得的結(jié)構(gòu)響應和流體運動過程與已有文獻結(jié)果吻合良好,驗證了本文所采用的雙向耦合方法在模擬帶自由液面的復雜流固耦合問題時具有較好的可靠性。

(3)對于潰壩水流沖擊彈性結(jié)構(gòu)問題中流體的壓力場和結(jié)構(gòu)的應力場進行了進一步分析,同時對比了潰壩水流沖擊剛性結(jié)構(gòu)和彈性結(jié)構(gòu)在相同位置處的壓力。研究較為深入地揭示了流固雙向耦合中的機理問題,為后續(xù)相關(guān)研究打下了良好的基礎。