基于浸沒光滑有限元法的雙圓柱體運動問題研究

王雙強，回 達，張桂勇，3，張智謙，宗 智，3

（1.貝德福德海洋研究所漁業與海洋部，加拿大達特茅斯B2Y 4A2；2. 大連理工大學運載工程與力學學部船舶工程學院，遼寧大連116024；3.大連理工大學工業裝備結構分析國家重點實驗室，遼寧大連116024；4. 新加坡高性能計算研究院，新加坡117576）

0 引 言

流場中距離較近的兩個或多個物體相對運動時會產生相互干擾，這種干擾現象廣泛存在于船舶與海洋工程領域，例如多船開采、多船裝運以及海上作業所依靠的多柱形結構海上平臺等。這些問題不但涉及到流體與固體間的相互作用，而且還涉及固體與固體之間的相互影響，是一種復雜的流固耦合現象，對這些現象所蘊含的科學問題進行研究極具實際意義[1-2]。

求解流固耦合問題的數值方法可分為動網格法和固定網格算法。耦合信息傳遞時涉及到固體移動邊界，拉格朗日移動網格方法[3]需要基于物體邊界不斷的網格重構，會給計算帶來很大困難。浸沒邊界法[4]的提出有效地避免了這個問題，為求解流固耦合問題提供了新的思路，并由此發展了一系列的浸沒類方法，如改進的浸沒邊界法[5]、虛擬流體方法[6]、浸沒有限元法[7]等。而原始浸沒邊界法的不足之處在于基于固體簡化的模型難以處理復雜的本構關系。浸沒有限元法將流體與固體統一為一個整體，固體采用有限元方法能夠反映真實的本構關系。在歐拉坐標系下基于Galerkin思想的不可壓縮粘性流體N-S方程的求解也有諸多方法，其中基于CBS算法的Galerkin方法求解流體流動時速度和壓力采用相同的插值函數，可以求解層流、湍流、可壓和不可壓縮等一系列流體問題[8]。

在浸沒有限元方法的基礎上，通過與最新推出的基于梯度光滑技術的光滑有限元方法[9]相結合，Zhang 等[10]提出了新穎高效的浸沒光滑有限元法，在求解典型的流固耦合問題中得到了驗證[11-12]。該方法采用浸沒光滑有限元法建立固體模型，用基于半隱式CBS 特征線法的Galerkin 方法模擬流體運動，為解決復雜流固耦合問題提供了有效途徑。基于浸沒光滑有限元法，本文以二維流域內雙圓柱為例，研究多體運動時的相互干擾作用及相應流場變化特性。

1 控制方程

1.1 控制方程分解

流固耦合系統控制方程可表示為如下：

浸沒光滑有限元方法引入虛擬流體，并假設虛擬流體有兩大特性：與真實流體相同的密度和粘性以及與固體粒子相同的運動特性[10]。基于這兩個假設可將流體域與固體域分離，固體所占區域Ωs用虛擬流體域Ωfc代替，滿足Ωs=Ωfc和Γfc=Γs，這樣整個流體域Ωf*=Ωf?Ωfc可采用歐拉網格描述，固體域Ωs采用拉格朗日網格描述，其基本思路如圖1所示。控制方程也分解為三部分：不可壓縮粘性流體的N-S方程、固體控制方程和流固耦合邊界條件。以下公式均采用Ωf代替Ωf*表示整個流體域。

圖1 浸沒光滑有限元基本思路Fig.1 Basic idea of immersed smoothed FEM

式中，Ωfc代表虛擬流體區域。

1.2 控制方程離散

1.2.1 求解N-S方程的CBS算法

本文采用半隱式的CBS算法離散不可壓縮粘性流體N-S方程，將速度和壓力場分離求解，并根據壓力對速度場進行修正，分為三步迭代求解：動量計算、壓力計算以及動量修正[13-14]。

1.3 梯度光滑技術

對于二維問題光滑域的構造主要有基于邊的光滑有限元和基于點的光滑有限元，本文固體求解采用基于邊的光滑有限元。固體模型位移梯度變量光滑后可表示為[10]

2 浸沒光滑有限元方法求解

考慮流體和固體相互作用時，需要將流固耦合作用補充到上述求解方程中。浸沒光滑有限元處理交界面的基本思想即是將耦合力等效為體力施加到固體上，對固體所覆蓋的虛擬流體區域施加流固耦合速度條件[10]，這樣可以分別求解流體和固體控制方程。

2.1 施加流固耦合力條件

2.2 施加流固耦合速度條件

2.3 流固耦合力計算

對流體考慮流固耦合力作用時，式（8）可表示為

2.4 程序流程

浸沒光滑有限元方法求解多體問題流程如圖2所示。

圖2 浸沒光滑有限元計算流程圖Fig.2 Calculation procedure of the immersed smoothed FEM

3 算 例

3.1 流域中自由下落的圓盤

如圖3所示，圓盤在重力作用下在流域內自由下落，所有邊界均滿足不可滑移邊界條件。流域尺寸為W=2.0 cm，H=5.0 cm。圓盤直徑D=0.25 cm，ρf=1 g/cm3，楊氏模量Es=104g/（cm?s2），泊松比ν=0.3。圓心初始位置距計算域頂端和左右邊界為L=1.0 cm。文獻[11]采用浸沒光滑有限元方法模擬圓盤在水中自由運動現象，固體達到穩定時的速度與經驗公式吻合。本算例計算流域介質為空氣ρf=10-3g/cm3，以驗證固液密度比較大時方法的適用性。重力加速度為g=980 cm/s2。網格采用任意生成的三角形單元，固體網格采用301個節點離散，流域網格采用25 351個節點離散。

圖3 流域內自由下落圓盤問題描述及流固離散模型Fig.3 Problem setting of a falling disk in the fluid domain and the corresponding calculation model

圖4 為流體域不同時刻的垂向速度分布。從圖4 可以看出，隨著運動物體速度的增加，流域內固體尾部交錯形成渦。圖5 為不同時刻固體域內垂向速度Vy分布。從圖5 可以看出，固體運動至t=0.04 s 時速度分布表現出一定的不均勻性，這是因為隨著時間的增加，各個節點受力不同及加速度有差異引起的；但隨著時間推移，在t=0.06 s至t=0.08 s時這種不均勻性逐漸減弱。

圖6為固體質點平均速度和位移的時間歷程圖，與理論解進行對比，結果吻合很好。由于采用浸沒方法的思想，計算此問題時我們無需進行網格重構，而穩定準確的計算結果也證明了浸沒光滑有限元法方法在計算固液大密度比問題時的可靠性。

圖4 不同時刻流域內垂向速度Vy分布Fig.4 Vertical velocity distribution of fluid domain at different times

圖5 不同時刻固體域內垂向速度Vy分布Fig.5 Vertical velocity distribution at different times

圖6 計算結果與理論解對比Fig.6 Calculation results in comparison with theoretical prediction

3.2 流域內雙圓柱體干擾

如圖7 所示，兩圓柱以1 m/s 的速度在流域內做相對運動，雷諾數Re=40。其中速度邊界條件為頂部和底部滿足不可滑移條件，壓力邊界條件為進口與出口處壓力為0。網格采用三角形單元，圓柱網格采用396 個節點；流域網格采用237 429 個節點。流固平均網格尺寸比分別為1:1。將t=8 s時兩圓柱上下相對的計算結果與文獻[19-20]對比，結果如圖8所示。從圖8可以看出本文計算結果與文獻吻合較好。同時本文計算結果顯示圓柱內部流線分布較為均勻，由于流域內固體覆蓋處的速度從固體域插值而來，說明速度可以通過插值算法得到有效傳遞。

圖7 流域內雙圓柱運動計算模型Fig.7 Calculation model of two cylinders

圖8 t=8 s時兩圓柱橫坐標相同時流線圖Fig.8 Streamline at t=8 s when two cylinders have the same x-coordinate value

圖9 不同時刻速度和壓力分布云圖Fig.9 Contours of velocities and pressures at different times

圖9顯示了三個不同時刻時的速度和壓力分布云圖，其中左側為速度云圖，右側為壓力云圖。從圖中可以看出，t=3.3 s時兩圓柱已有干擾作用但不明顯，兩圓柱間區域壓力值增大；t=7.9 s時兩圓柱距離較近，干擾作用較強，壓力云圖表現明顯，物體附近高壓區域減小而低壓區得以加強；t=14 s 時干擾作用主要存在于尾部區域且已經比較微弱。

4 結 論

本文采用浸沒光滑有限元法模擬二維流域內自由下落的圓盤以及雙圓柱相對運動時引起的流場變化特性，并通過與文獻對比得出以下結論:

（1）引用虛擬流體，將流固分離后分別求解可以避免網格重構，同時采用簡單的三角形單元離散計算域，可以簡化前處理過程。

（2）浸沒光滑有限元法能夠模擬固液大密度比的流固耦合現象。

（3）浸沒光滑有限元法能很好地模擬雙圓柱運動時的干擾作用。從流線圖可以看出，虛擬流體速度可通過插值得到有效傳遞；從壓力云圖可以看出，圓柱相對運動時區間壓力變大；距離較近時正壓區域減少，負壓值得以加強；逐漸遠離時尾部壓力受干擾作用影響變弱，符合物理本質特性。