王雙強(qiáng)，回 達(dá)，張桂勇，3，張智謙，宗 智，3

（1.貝德福德海洋研究所漁業(yè)與海洋部，加拿大達(dá)特茅斯B2Y 4A2；2. 大連理工大學(xué)運(yùn)載工程與力學(xué)學(xué)部船舶工程學(xué)院，遼寧大連116024；3.大連理工大學(xué)工業(yè)裝備結(jié)構(gòu)分析國家重點(diǎn)實(shí)驗(yàn)室，遼寧大連116024；4. 新加坡高性能計(jì)算研究院，新加坡117576）

0 引 言

流場中距離較近的兩個(gè)或多個(gè)物體相對運(yùn)動(dòng)時(shí)會(huì)產(chǎn)生相互干擾，這種干擾現(xiàn)象廣泛存在于船舶與海洋工程領(lǐng)域，例如多船開采、多船裝運(yùn)以及海上作業(yè)所依靠的多柱形結(jié)構(gòu)海上平臺(tái)等。這些問題不但涉及到流體與固體間的相互作用，而且還涉及固體與固體之間的相互影響，是一種復(fù)雜的流固耦合現(xiàn)象，對這些現(xiàn)象所蘊(yùn)含的科學(xué)問題進(jìn)行研究極具實(shí)際意義[1-2]。

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

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

1 控制方程

1.1 控制方程分解

流固耦合系統(tǒng)控制方程可表示為如下：

浸沒光滑有限元方法引入虛擬流體，并假設(shè)虛擬流體有兩大特性：與真實(shí)流體相同的密度和粘性以及與固體粒子相同的運(yùn)動(dòng)特性[10]?；谶@兩個(gè)假設(shè)可將流體域與固體域分離，固體所占區(qū)域Ωs用虛擬流體域Ωfc代替，滿足Ωs=Ωfc和Γfc=Γs，這樣整個(gè)流體域Ωf*=Ωf?Ωfc可采用歐拉網(wǎng)格描述，固體域Ωs采用拉格朗日網(wǎng)格描述，其基本思路如圖1所示?？刂品匠桃卜纸鉃槿糠郑翰豢蓧嚎s粘性流體的N-S方程、固體控制方程和流固耦合邊界條件。以下公式均采用Ωf代替Ωf*表示整個(gè)流體域。

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

式中，Ωfc代表虛擬流體區(qū)域。

1.2 控制方程離散

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

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

1.3 梯度光滑技術(shù)

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

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

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

2.1 施加流固耦合力條件

2.2 施加流固耦合速度條件

2.3 流固耦合力計(jì)算

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

2.4 程序流程

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

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

3 算 例

3.1 流域中自由下落的圓盤

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

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

圖4 為流體域不同時(shí)刻的垂向速度分布。從圖4 可以看出，隨著運(yùn)動(dòng)物體速度的增加，流域內(nèi)固體尾部交錯(cuò)形成渦。圖5 為不同時(shí)刻固體域內(nèi)垂向速度Vy分布。從圖5 可以看出，固體運(yùn)動(dòng)至t=0.04 s 時(shí)速度分布表現(xiàn)出一定的不均勻性，這是因?yàn)殡S著時(shí)間的增加，各個(gè)節(jié)點(diǎn)受力不同及加速度有差異引起的；但隨著時(shí)間推移，在t=0.06 s至t=0.08 s時(shí)這種不均勻性逐漸減弱。

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

圖4 不同時(shí)刻流域內(nèi)垂向速度Vy分布Fig.4 Vertical velocity distribution of fluid domain at different times

圖5 不同時(shí)刻固體域內(nèi)垂向速度Vy分布Fig.5 Vertical velocity distribution at different times

圖6 計(jì)算結(jié)果與理論解對比Fig.6 Calculation results in comparison with theoretical prediction

3.2 流域內(nèi)雙圓柱體干擾

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

圖7 流域內(nèi)雙圓柱運(yùn)動(dòng)計(jì)算模型Fig.7 Calculation model of two cylinders

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

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

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

4 結(jié) 論

本文采用浸沒光滑有限元法模擬二維流域內(nèi)自由下落的圓盤以及雙圓柱相對運(yùn)動(dòng)時(shí)引起的流場變化特性，并通過與文獻(xiàn)對比得出以下結(jié)論:

（1）引用虛擬流體，將流固分離后分別求解可以避免網(wǎng)格重構(gòu)，同時(shí)采用簡單的三角形單元離散計(jì)算域，可以簡化前處理過程。

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

（3）浸沒光滑有限元法能很好地模擬雙圓柱運(yùn)動(dòng)時(shí)的干擾作用。從流線圖可以看出，虛擬流體速度可通過插值得到有效傳遞；從壓力云圖可以看出，圓柱相對運(yùn)動(dòng)時(shí)區(qū)間壓力變大；距離較近時(shí)正壓區(qū)域減少，負(fù)壓值得以加強(qiáng)；逐漸遠(yuǎn)離時(shí)尾部壓力受干擾作用影響變?nèi)?，符合物理本質(zhì)特性。