一種基于改進流固耦合力求解的浸沒光滑點插值法

楊博睿，張桂勇，2，王雙強，嚴博騫，王 鵬

（1.大連理工大學工業裝備結構分析國家重點實驗室；船舶工程學院，遼寧大連 116024；2.高新船舶與深海開發裝備協同創新中心，上海 200240）

0 引 言

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

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

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

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

1.1 流固耦合控制方程

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

文中統一采用右上角標f、s、fc 分別表示所屬區域為流體域、固體域和虛擬流體域，FSI表示流固耦合。上式中，vf表示流體速度，Pf表示流體壓力，Ff,FSI表示流體受到的流固耦合力，g為重力加速度。整個流固耦合系統控制方程可以分成三部分。

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

（2）固體運動方程：

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

（3）流固耦合系統的速度條件和力條件：

1.2 流體控制方程的求解

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

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

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

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

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

1.3 固體控制方程的求解

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

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

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

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

2 改進的流固耦合力求解方式

2.1 施加流固耦合速度條件

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

2.2 施加流固耦合力條件

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

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

2.3 基于歐拉網格的流固耦合力求解

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

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

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

聯立以上二式，可以得到歐拉網格下計算出的流體節點的流固耦合力：

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

圖2 應力插值過程示意圖Fig.2 Schematic diagram of stress interpolation process

3 數值算例

3.1 圓柱繞流

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

3.2 圓盤沉降

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

表4 各網格的節點數量及不同方法模擬的下落穩定速度Tab.4 Node number of meshes and steady velocities in falling simulated by different methods

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

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

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

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

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

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

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

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

4 結 語

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

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

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