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

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

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

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

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

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

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

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

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

1 數學模型

1.1 FVM方法

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

1.1.1 控制方程

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

(1)

(2)

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

1.1.2 VOF方法

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

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

(3)

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

(4)

(5)

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

1.2 FEM方法

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

(6)

C=α1M+α2K

(7)

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

(8)

(9)

(10)

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

(11)

(12)

1.3 FVM-FEM耦合方法

圖1 雙向流固耦合模型

2 算例分析

2.1 潰壩水流沖擊剛性結構

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

圖2 潰壩水流沖擊剛性結構

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

圖3 水流沖擊剛性結構物的自由表面分布

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

圖4 計算域中液體質量的變化

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

圖5 左側壁處的液面變化

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

圖6 A點處壓力變化

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

2.2 潰壩水流沖擊彈性結構

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

圖7 潰壩水流沖擊彈性結構

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

圖8 彈性結構左上角的位移變化

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

圖9 不同時刻的彈性板變形和水流形態

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

圖10 不同時刻流體的壓力和彈性結構的應力分布

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

圖11 剛性結構和彈性結構左側中點處的壓力

3 結論

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

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

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

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