楊博睿，張桂勇，2，王雙強(qiáng)，嚴(yán)博騫，王 鵬

（1.大連理工大學(xué)工業(yè)裝備結(jié)構(gòu)分析國家重點實驗室；船舶工程學(xué)院，遼寧大連 116024；2.高新船舶與深海開發(fā)裝備協(xié)同創(chuàng)新中心，上海 200240）

0 引 言

流固耦合是船舶與海洋工程領(lǐng)域中的一種常見現(xiàn)象，它涵蓋了非定常流體流動和固體大變形分析，具有十分重要的研究價值。1972 年，Peskin 提出浸沒邊界法（IBM）[1]，并首次應(yīng)用了非貼體網(wǎng)格。該方法無需進(jìn)行網(wǎng)格重構(gòu)操作，極大地節(jié)省了計算資源，因此十分適合應(yīng)用于求解固體域發(fā)生較大變形或移動的流固耦合問題中。在此基礎(chǔ)上，浸沒類方法的思想也與有限元法相結(jié)合，應(yīng)用到流固耦合問題的求解當(dāng)中，取得了較好的模擬效果[2]。

光滑點插值法（S-PIM）是由Liu和Zhang[3]基于廣義梯度光滑度技術(shù)[4]提出的一種插值方法。與傳統(tǒng)的有限元方法（FEM）[5]相比，S-PIM 可以在使用簡單的低階單元時進(jìn)行高階插值，軟化了模型剛度，提高了計算精度并降低了對網(wǎng)格質(zhì)量的依賴性。本文采用的浸沒光滑點插值法（IS-PIM），便是在浸沒類方法的框架下，采用光滑點插值法作為固體求解器的一種流固耦合計算方法[6]。

在傳統(tǒng)的ISPIM 求解過程中，流固耦合力是通過虛擬流體拉格朗日網(wǎng)格計算的[7]，并沒有考慮固體邊界的速度梯度。所以這種計算方式會導(dǎo)致邊界處粘性力被忽略。本文提出一種基于真實流體歐拉網(wǎng)格求解，并通過插值傳遞到固體節(jié)點的流固耦合力施加方法。該方法直接通過流體節(jié)點的額外動量變化求解流固耦合力，可規(guī)避需單獨進(jìn)行粘性力計算的缺點，程序流程簡單。本文對新方法與傳統(tǒng)求解方式進(jìn)行對比介紹，并通過數(shù)值算例展示新方法在提高計算精度方面的作用，同時驗證流固耦合力條件和速度條件在節(jié)點動量改變方面的等效性，強(qiáng)化新方法的理論基礎(chǔ)。

1 浸沒光滑點插值方法介紹

1.1 流固耦合控制方程

浸沒光滑點插值法采用與直接力法相似的思路，將邊界處原本復(fù)雜的流固耦合作用轉(zhuǎn)化為動量方程中的體力，施加到不可壓縮粘性流體控制方程，即Navior-Stokes 方程中，以表示耦合作用對流場的影響：

文中統(tǒng)一采用右上角標(biāo)f、s、fc 分別表示所屬區(qū)域為流體域、固體域和虛擬流體域，F(xiàn)SI表示流固耦合。上式中，vf表示流體速度，Pf表示流體壓力，F(xiàn)f,FSI表示流體受到的流固耦合力，g為重力加速度。整個流固耦合系統(tǒng)控制方程可以分成三部分。

（1）不可壓縮粘性流體的控制方程：

（2）固體運(yùn)動方程：

式中，ρs、us和Ps分別為固體的密度、位移和第一類皮奧拉-基爾霍夫應(yīng)力。

（3）流固耦合系統(tǒng)的速度條件和力條件：

1.2 流體控制方程的求解

本文在流體求解時采用了基于特征線算子分離的半隱式CBS方法[8]。該方法基于特征線理論，以流體守恒方程為基礎(chǔ)，可以得到計算簡單的顯式格式，同時還具有較好的收斂性，消除了穩(wěn)定性條件的限制[9]。通過使用三節(jié)點三角形單元對流體域進(jìn)行離散，單元的速度和壓力值可以表示為

文中統(tǒng)一采用右下角標(biāo)表示變量指代的單元和方向。其中大寫字母代表節(jié)點序號，小寫字母表示方向。左上角標(biāo)則表示時刻。如式（5）中，ΦfJ為流體單元第J號節(jié)點的形函數(shù)，nPfJ、nvfJi分別表示t=n時刻流體節(jié)點J的壓力值和i方向的速度（i=1,2分別表示x方向和y方向）。

完成空間和時間離散之后，CBS方法對流體從n時刻到n+1時刻的求解可簡單概括為三個步驟。

步驟二：利用中間時刻速度和壓力邊界條件求解第（n+1）時刻的流體壓力n+1PfJ，

步驟三：利用當(dāng)前壓力值和速度邊界條件求解第（n+1）時刻處的流體速度，

1.3 固體控制方程的求解

非線性固體的求解采用光滑點插值方法S-PIM。光滑域構(gòu)造格式如圖1所示。

固體單元采用的光滑函數(shù)W(Xs)為

圖1 邊基光滑域構(gòu)造格式[10]Fig.1 Construction of edge-based smoothing domain[10]

式中：Ji表示光滑域內(nèi)的第i個節(jié)點；Sˉ為光滑皮奧拉-基爾霍夫第二應(yīng)力，為偏應(yīng)力Sˉdev和體積應(yīng)力Sˉvol之和。

2 改進(jìn)的流固耦合力求解方式

2.1 施加流固耦合速度條件

虛擬流體節(jié)點需要施加耦合速度條件，以此來保證與固體的同步運(yùn)動。由于采用非貼體網(wǎng)格，節(jié)點速度信息的傳遞需要通過形函數(shù)插值實現(xiàn)：

2.2 施加流固耦合力條件

在求解固體域時，需要將流體求解所得到的力條件傳遞給虛擬流體域，作為流固耦合力條件。虛擬流體的動量方程可寫為

式中，力Ffc保證了虛擬流體節(jié)點速度與固體節(jié)點速度相同的虛擬流體假設(shè)，其本質(zhì)可看作為流固耦合力。在原始IS-PIM 的力條件施加過程中，粘性力成分是在保證虛擬流體與固體擁有相同的運(yùn)動響應(yīng)的基礎(chǔ)上，通過速度梯度求得的：

2.3 基于歐拉網(wǎng)格的流固耦合力求解

原始IS-PIM 方法是基于虛擬流體拉格朗日網(wǎng)格求解流固耦合力，進(jìn)而施加力條件的。由于真實流體網(wǎng)格與固體網(wǎng)格不完全重合，因此會導(dǎo)致流固邊界處的粘性力被忽略，產(chǎn)生數(shù)值誤差，在模擬大剛度或低雷諾數(shù)流動問題時尤為明顯。因此，本文提出一種基于歐拉網(wǎng)格的求解思路。不同于傳統(tǒng)的求解方式，該方法直接通過流體控制方程求解各個流體單元所受到的流固耦合力，然后作為耦合力條件施加到固體節(jié)點上。

根據(jù)CBS計算過程，由式（6）可得到基于歐拉網(wǎng)格求解時的動量方程為

再考慮流固耦合速度條件的作用，將動量方程改寫為

聯(lián)立以上二式，可以得到歐拉網(wǎng)格下計算出的流體節(jié)點的流固耦合力：

該公式形式上與牛頓第二定律公式類似，體現(xiàn)了流固耦合力與流體節(jié)點動量變化的關(guān)系，并將在后續(xù)的數(shù)值算例部分進(jìn)行驗證。接下來只需通過形函數(shù)進(jìn)行單元應(yīng)力的插值，將流體節(jié)點流固耦合力傳遞給固體節(jié)點，作為求解固體域時的耦合力條件即可，過程如圖2所示。

圖2 應(yīng)力插值過程示意圖Fig.2 Schematic diagram of stress interpolation process

3 數(shù)值算例

3.1 圓柱繞流

首先對低雷諾數(shù)下圓柱繞流問題進(jìn)行模擬，驗證本方法對流固耦合力計算的精度及其有效性。設(shè)置圓柱直徑D=1 m，流域長度L=32 m，寬度B=16 m。幾何域信息如圖3 所示。物理屬性方面，流體為不可壓縮粘性流體，密度設(shè)為ρf=1000 kg/m3；固體圓柱為固定的剛體。流場左側(cè)邊界有沿X方向、速度為Vf=1 m/s 的來流。分別采用固體節(jié)點數(shù)為380、575 和743，以及流體節(jié)點數(shù)為22 141、31 161 和55 940的網(wǎng)格對雷諾數(shù)為100時圓柱繞流的阻力系數(shù)Cd進(jìn)行了對比。網(wǎng)格劃分如圖4所示，最終模擬結(jié)果如表1 所示。可以看到隨著固體節(jié)點或流體節(jié)點數(shù)的增加，相同雷諾數(shù)下圓柱的阻力系數(shù)變化并不大，結(jié)果已收斂。

圖3 圓柱繞流的幾何域信息Fig.3 Geometric model of flow past a cylinder

表1 Re=100時，不同方法在各網(wǎng)格下求得的CdTab.1 Summary of Cd obtained by 2 methods under different mesh sizes when Re=100

驗證網(wǎng)格的收斂性后，選用網(wǎng)格5 對雷諾數(shù)為20、40、100 和200 時的圓柱繞流現(xiàn)象進(jìn)行了模擬。將圓柱繞流的阻力系數(shù)Cd與其他研究學(xué)者的結(jié)果以及基于拉格朗日網(wǎng)格進(jìn)行額外修正粘性力的ISPIM算法[12]進(jìn)行了對比，結(jié)果如表2所示。

表2 兩種方法在不同雷諾數(shù)下求得的CdTab.2 Summary of Cd obtained by 2 methods under different Re

結(jié)果顯示，新提出的基于歐拉網(wǎng)格的求解方式無需進(jìn)行邊界粘性力的修正，即可較好地模擬低雷諾數(shù)下的圓柱繞流現(xiàn)象，而原始ISPIM 由于無法計算流固邊界的粘性力，對阻力系數(shù)的模擬偏小很多。圖5 為各雷諾數(shù)下圓柱繞流穩(wěn)定后阻力系數(shù)的對比。

接下來對新方法模擬的流場速度、壓力云圖進(jìn)行分析。圖6～7 分別為不同雷諾數(shù)下流場穩(wěn)定后速度云圖及流線?？梢钥闯?，在低雷諾數(shù)流動下，速度云圖在圓柱上下兩側(cè)呈對稱分布，且圓柱兩側(cè)的流線較為平滑。隨著雷諾數(shù)的增大，流場波動變得劇烈，圓柱后形成的旋渦結(jié)構(gòu)也越發(fā)明顯。

圖5 各雷諾數(shù)下求得的Cd與其他學(xué)者結(jié)果對比Fig.5 Comparison of Cd under different Re with those of other researches

圖6 基于歐拉網(wǎng)格模擬速度云圖及流線（左：Re=20，右：Re=40）Fig.6 Velocity contour and streamline of fluid for Re=20(left)and Re=40(right)

圖7 基于歐拉網(wǎng)格模擬速度云圖及流線（左：Re=100，右：Re=200）Fig.7 Velocity contour and streamline of fluid for Re=100(left)and Re=200(right)

圖8 基于歐拉網(wǎng)格模擬Re=100時圓柱繞流的壓力云圖Fig.8 Pressure contour of fluid for Re=100

圖9 基于歐拉網(wǎng)格模擬Re=200時圓柱繞流的壓力云圖Fig.9 Pressure contour of fluid for Re=200

圖8～9 分別展示了Re=100 和Re=200 的壓力云圖，圓柱后側(cè)能觀察到明顯的卡門渦街現(xiàn)象。在圖9 中，圓柱尾流上下交替出現(xiàn)的低壓區(qū)顏色較圖8 中的更深，體現(xiàn)了雷諾數(shù)增大后旋渦整體強(qiáng)度的提高。上述模擬現(xiàn)象均與實際物理現(xiàn)象相符合。

圖10 為對Re=40 時圓柱后尾渦長度的測量，表3 為測量結(jié)果與其他學(xué)者結(jié)果的對比。結(jié)果顯示基于歐拉網(wǎng)格求解的IS-PIM能夠較準(zhǔn)確地模擬圓柱繞流的尾流和旋渦結(jié)構(gòu)。

圖10 尾渦長度測量示意圖Fig.10 Measurement of vorticity length

表3 Re=40時流場穩(wěn)定后尾渦長度大小比較Tab.3 Summary of the length of the vorticity for Re=40

3.2 圓盤沉降

圓盤受重力作用在流場中自由沉降是一個經(jīng)典的流固耦合算例。設(shè)置流域?qū)? cm，高5 cm。圓盤直徑D=0.25 cm。物理屬性方面，設(shè)置固體密度為2 g/cm3，泊松比為0.3，彈性系數(shù)Es=1.0×104g/(cm·s2)。流體密度為1 g/cm3。初始時刻流體和固體速度均為0。首先以粘性系數(shù)μ=1.25 g/(cm·s)為例，選用不同節(jié)點數(shù)的網(wǎng)格對圓盤在下落過程中所達(dá)到的穩(wěn)定速度進(jìn)行了模擬計算，結(jié)果如表4所示。圖11以基于歐拉網(wǎng)格的計算為例，展示了不同網(wǎng)格下的模擬結(jié)果，可以看出隨著節(jié)點數(shù)的增加，模擬結(jié)果已收斂。

表4 各網(wǎng)格的節(jié)點數(shù)量及不同方法模擬的下落穩(wěn)定速度Tab.4 Node number of meshes and steady velocities in falling simulated by different methods

經(jīng)收斂性驗證后，選擇網(wǎng)格5 進(jìn)行后續(xù)的數(shù)值模擬。采用新提出的基于歐拉網(wǎng)格求解方式和原始IS-PIM 對粘性系數(shù)為1 g/(cm·s)、1.25 g/(cm·s)和1.5 g/(cm·s)三種情況進(jìn)行模擬，幾何域信息及網(wǎng)格劃分情況如圖12所示。

圖11 基于歐拉網(wǎng)格在粘性系數(shù)1.25 g/(cm·s)時圓盤速度時間曲線Fig.11 Velocity-time curve with μ=1.25 g/(cm·s)solved based on Euler grid

圖12 圓盤沉降的幾何域及圓盤附近的網(wǎng)格情況Fig.12 Geometric domain of disk falling and the grid near the disk

不同粘性系數(shù)下圓盤從開始下落到達(dá)到穩(wěn)定下落速度時的速度時間歷程曲線如圖13～15 所示，表5為最終結(jié)果及誤差對比。從速度時間曲線中可以看出，0.15 s過后，圓盤的下落速度基本穩(wěn)定，而原始IS-PIM 的計算方式誤差非常大，這是由于流體的粘性在圓盤下落中起到非常重要的作用，而基于歐拉網(wǎng)格的計算方式在很大程度上彌補(bǔ)了這一缺陷。

圖13 粘性系數(shù)為1時圓盤速度時間曲線Fig.13 Velocity-time curve with μ=1 g/(cm·s)

圖14 粘性系數(shù)為1.25時圓盤速度時間曲線Fig.14 Velocity-time curve with μ=1.25 g/(cm·s)

圖15 粘性系數(shù)為1.5時圓盤速度時間曲線Fig.15 Velocity-time curve with μ=1.5 g/(cm·s)

為了說明耦合速度條件在改變流體節(jié)點動量方面與耦合力具有等效性，驗證公式（18）中對流固耦合力作用的解釋。在基于歐拉網(wǎng)格的計算模擬中，每隔0.01 s 對粘性系數(shù)為1.25 g/(cm·s)的圓盤落水的固體節(jié)點動量進(jìn)行了監(jiān)測，并與同時刻所求得的流固耦合力進(jìn)行對比，繪制成時間歷程曲線，如圖16 所示。二者吻合較好，說明了該方法中所采用的施加耦合邊界條件的方式來考慮原本較為復(fù)雜的流固耦合力作用的正確性。

4 結(jié) 語

表5 不同粘性系數(shù)下求得的下落中穩(wěn)定速度（cm/s）Tab.5 Summary of steady velocities in falling with different viscous coefficients(cm/s)

圖16 基于歐拉網(wǎng)格模擬下流體節(jié)點在y方向上的額外動量變化與流固耦合力對比圖Fig.16 Comparison of extra momentum change and FSI force of fluid nodes in y direction simulated by Euler grid

本文針對原始IS-PIM 無法計算流固邊界粘性力的問題，提出了一種基于真實流體歐拉網(wǎng)格求解并施加流固耦合力的新方法。該方法從計算流體域全域出發(fā)求解動量方程，方法理論和求解步驟更加簡單，從根本上避免了流固邊界處粘性力被忽略的問題。通過數(shù)值算例的驗證結(jié)果可以說明，基于歐拉網(wǎng)格計算流固耦合力的方法在模擬固定圓柱的繞流問題時，可以準(zhǔn)確地體現(xiàn)粘性力作用，因此阻力系數(shù)的計算結(jié)果與原始IS-PIM 相比更加精確。在模擬固體運(yùn)動的問題時，如圓盤沉降，精度提升的效果也十分明顯，有效解決了原始IS-PIM 在模擬中忽略邊界粘性力而導(dǎo)致圓盤下落過快的問題。在后續(xù)的研究中，可以基于該求解方式對流固邊界的信息交換方式加以優(yōu)化，進(jìn)一步提高模擬精度。