基于浸入邊界法的雙層柔性蒙皮減阻性能研究

任洵濤，趙丹，陳輝，崔進

1. 中國船舶集團第七〇四研究所，上海 200031

2. 哈爾濱工程大學 機電工程學院，黑龍江 哈爾濱 150001

水下航行器所受流體阻力的大小直接影響其航速以及能源消耗速度。目前主要的流體減阻技術有溝槽表面減阻技術[1]、低表面能減阻技術[2]、微氣泡減阻技術[3]和外加射流減阻技術[4]等。但是從水下航行器的實際工程應用角度來看，很多減阻技術因為減阻機理不明確、使用維護困難等制約因素而不適合用于航行器的流體減阻上。因此提出一種新型而有效的減阻結構具有重要意義。

柔性表面減阻技術起源于研究人員對海豚皮膚的仿生研究[5]，指的是在水下航行器外殼上敷設一層柔性蒙皮，柔性蒙皮會隨著來流產生相應的形貌變化，進而產生減阻效果。值得注意的是，柔性表面減阻技術雖然能起到一定的減阻效果，但由于其表面減阻的結構參數通常來源于自然界中動植物表面特定的連續性幾何形貌，導致減阻表面抽象度低、在實際工程中應用難度大以及結構參數設定具有一定盲目性等問題，不便于對研究對象進行減阻特性的深入研究。

近年來，雙層隔振系統廣泛應用于艙室內的減振研究中[6]，雙層隔振系統是在單層隔振系統的基礎上增加一個中間質量塊及一層彈性支撐形成的，即雙層隔振系統由2 層彈簧阻尼隔振器和隔振器之間的中間質量塊組成。相比于單層隔振系統，雙層隔振系統能在較寬的頻率范圍內發揮更好的減振效果，現已被廣泛應用在汽車、航天和船舶等領域，是振動控制系統研究的重點領域。

本文針對減小水下航行器阻力問題，考慮到當前表面減阻技術的局限性以及雙層隔振系統的廣泛應用，將柔性表面減阻技術和微小化、陣列化的雙層隔振系統相結合，提出一種新型的雙層柔性蒙皮結構。針對此結構復雜且具有大變形的特點，本文通過基于浸入邊界法的程序框架IB2d[7-8]建立仿真模型，并對其進行數值計算，以研究雙層柔性蒙皮各結構參數對減阻性能的影響。

1 雙層柔性蒙皮模型及控制方程

1.1 雙層柔性蒙皮物理模型

本文設計的柔性蒙皮結構由柔性壁面和彈性支撐組成，彈性支撐則是由中間質量塊和上下層彈簧及阻尼構成，其物理模型如圖1 所示。

圖1 柔性蒙皮物理模型

雙層柔性蒙皮外層的柔性壁面一方面作為傳統柔性蒙皮發揮流體減阻作用，另一方面發揮約束作用，用以固定彈性支撐。

1.2 浸入邊界法的控制方程

雙層柔性蒙皮結構復雜且存在大變形問題，常規的數值方法在處理這類問題時，需要生成復雜的貼體網格，且需要不斷在新舊網格上交換各種數據。這兩方面不但加大了計算量，而且降低了仿真計算的精度和穩定性。浸入邊界法通過在N-S 方程上附加一個體積力項來表示固體對流體的影響，無需生成復雜的貼體網格，大大簡化了計算過程，適合對雙層柔性蒙皮進行仿真計算。

浸入邊界法控制方程以N-S 方程為出發點，描述了具有慣性的牛頓流體的運動。包括雙層柔性蒙皮和流體在內的整個物理區域，都假定為不可壓縮牛頓流體的黏性流動，其連續性和動量方程可分別表示為

式中：ρ為流體的密度；u(x,t)=(u(x,t),v(x,t))為歐拉描述下的流體速度矢量；p(x,t)為流體壓力；f(x,t)為雙層柔性蒙皮的單位力浸入邊界施加到流體上的體積力；自變量是位置x=(x,y)和時間t。

式(1)等價于流體的動量守恒，而式(2)則是要求流體不可壓縮的條件。

流體與雙層柔性蒙皮浸入邊界之間的相互作用方程為

式中：X(r,t)為拉格朗日參數r標記的材料點在時間t的笛卡爾坐標；F(r,t)為雙層柔性蒙皮浸入邊界中的彈性變形對函數施加到流體上的每單位面積的力的拉格朗日位置r和時間t的乘積。式(3)中的歐拉力f(x,t)定義在整個流體域 Ω上，而拉格朗日力F(r,t)定義在雙層柔性蒙皮浸入邊界上。

正則化的增量函數為

基于浸入邊界法控制方程的離散空間上采用的是有限差分法，時間上采用的是基于中點規則的二階龍格–庫塔方法，具體離散過程參見文獻[9]。

1.3 雙層柔性蒙皮仿真模型

本文利用基于浸入邊界法的IB2d 軟件包，通過修改相應的代碼程序來建立仿真模型、設置邊界約束和仿真參數等。

1)柔性壁面仿真模型

柔性蒙皮外層柔性壁面的材料為聚氨酯彈性體的柔韌薄膜，其尺寸為1 m×1 m×4 mm，使用建模單元中的拉格朗日點、彈簧和質量單元建立柔性壁面仿真模型，如圖2 所示。

圖2 柔性壁面仿真模型示意

由圖2 可知，柔性壁面的浸入邊界由若干拉格朗日點組成，其上、下兩面為浸入邊界，左、右兩側面為邊界約束。模型中每2 個相鄰點的距離為2 mm，柔性壁面x方向長為1 m，沿x方向的上、下2 個浸入邊界各由501 個拉格朗日點組成；默認z方向為單位長度；y方向厚4 mm，左右兩側設置3 個拉格朗日點作為固定約束，在每個拉格朗日點上附加相同質量單元，模擬柔性壁面本身的質量。每2 個相鄰拉格朗日點之間連接彈簧單元，模擬柔性壁面x和y方向上的剛度。

2)彈性支撐仿真模型

柔性蒙皮內層的彈性支撐是分布式雙層隔振系統。使用建模單元中的拉格朗日點、阻尼彈簧和質量單元建立彈性支撐，如圖3 所示。

圖3 彈性支撐仿真建模示意

根據上述建模過程，建立的雙層柔性蒙皮仿真模型如圖4 所示。

圖4 雙層柔性蒙皮仿真模型

2 算例驗證

2.1 圓柱繞流

圓柱繞流是流體力學中的經典問題，對其有許多的理論以及實驗研究[10-12]，方便對仿真結果進行驗證。因此，本文按照文獻[10-12]中的圓柱繞流試驗參數進行了圓柱繞流仿真計算，在仿真與試驗參數一致的情況下，驗證了IB2d 數值計算結果的準確性。

設置計算域尺寸為1 m×0.25 m，選擇直徑為D=0.08 m 的圓柱，以計算域左下角為原點，圓柱中心位于(0.1 m，0.125 m)。流體密度 ρ設定為100 kg/m3，流場速度v=1 m/s。取網格為512×128，時間步長為5×10-5s。確定上述的參數后，設定計算時間為27 s，即無量綱時間T=180，測量渦脫落時間并計算渦脫落頻率，最后求得斯特勞哈爾數、阻力系數、升力系數。按上述的參數設置并完成仿真計算，得到雷諾數分別為100、150 和200 時的斯特勞哈爾數，并與其他研究人員實驗測量數據[10-12]進行對比，如表1 所示。由表可知，基于浸入邊界法數值計算得到的斯特勞哈爾數與實驗測量數據具有很好的一致性。

表1 斯特勞哈爾數比較

Re=150 時圓柱繞流的升力系數和阻力系數曲線分別如圖5 和圖 6 所示。

圖5 升力系數曲線

圖6 阻力系數曲線

將本文得到的平均升力和阻力系數與Lai 等[10]的結果進行對比，結果如表2 所示。由表可知，本文計算結果與其計算結果吻合良好。

表2 阻力系數對比

本文在Re=150、計算時間T=180 時，根據數值計算結果繪制渦脫落的瞬時渦度等值線，得到圓柱尾渦狀態如圖7 所示。由圖可知，沿圓柱表面分布的邊界層、分離剪切層的卷起以及圓柱后側的分離區域都清晰可見。圓柱尾渦出現明顯正、負交替向下游發展的卡門渦街現象。

圖7 圓柱繞流尾渦狀態

綜上，根據圓柱繞流算例中的升力系數、阻力系數、斯特勞哈爾數以及圓柱繞流尾渦狀態，驗證了IB2d 數值計算結果的準確性。

2.2 柔性絲線

為進一步說明IB2d 軟件包解決柔性大變形流固耦合問題的準確性和有效性，本文利用柔性絲線案例對其進行了驗證。

本算例來源于Zhu[13]的數值仿真?？紤]一種不可壓縮黏性流體在重力驅動下的二維垂直流動，彈性纖維最初水平放置在實驗通道中心，中點固定在頂部邊界附近。通道的2 個側壁是剛性的，側壁上的速度采用無滑移邊界條件。其計算模型示意圖如圖8 所示。仿真參數設置如表3所示。

圖8 柔性絲線計算模型示意

表3 仿真參數設置

按照表3 的參數進行仿真計算，所得到的計算結果如圖9 所示。由圖可知，仿真結果與文獻[13]中的實驗結果相吻合，進一步說明了IB2d 軟件包能夠處理大變形問題，且具有較好的精度。

圖9 柔性絲線仿真計算與實驗數據對比

3 雙層柔性蒙皮仿真結果分析

3.1 仿真參數設置

仿真模型邊界約束為柔性壁面左、右兩側面以及各彈性支撐底端，邊界約束的虛擬彈簧剛度取2×109N/m。在流體域尺寸為1.6 m×0.6 m，網格密度為400×150 時，滿足了剛性邊界條件。同理，質量單元虛擬彈簧剛度也設置為2×109N/m。流體密度ρ 設置為100 kg/m3，動力黏度為0.001 N·s/m2，背景流速u分別設置為2、5、10、20 m/s，則蒙皮中點處的雷諾數分別為1×105、 2.5×105、 5×105、1×106。

3.2 仿真參數可行性驗證

為了驗證本文中的網格密度的可行性，取時間步長為5×10-5s，流速為5 m/s，網格密度分別取400×150、400×300、800×300 以及800×600 時進行仿真計算，并以網格密度為800×600 時對應的平均阻力系數為基準，計算了另外3 種網格密度相對于800×600 時平均阻力系數的誤差，結果如表4 所示。

表4 不同網格密度下的仿真結果

由表4 可知，網格密度為800×600 時的平均阻力系數相對于800×300 時不發生變化，并且平均阻力系數在網格密度為400×150 時相對于網格密度為800×600 時的誤差僅為1.8%，滿足仿真計算要求。為了提高計算效率，選擇400×150 的網格密度。

為了驗證本文中時間步長的取值合理，在網格密度為400×150 的情況下，取流速為5 m/s，時間步長分別取5×10-5、2.5×10-5、1×10-5s 進行仿真計算，并以1×10-5s 對應的平均阻力系數為基準，計算了另外2 種時間步長相對于1×10-5s 時平均阻力系數的誤差，結果如表5 所示。

表5 不同時間步長下的仿真結果

由表5 可知，隨著時間步長的增加，誤差不發生變化，并且隨著時間步長的增加平均阻力系數相差較小，即3 種時間步長都滿足仿真計算要求。為提高計算效率，本文選擇時間步長為5×10-5s 對雙層柔性蒙皮進行仿真計算。

3.3 不同參數下雙層柔性蒙皮減阻性能分析

本文以阻力系數衡量雙層柔性蒙皮的減阻性能。柔性蒙皮在流體作用下所受的阻力來自2 個方面，即剪切應力和分布在柔性蒙皮浸入邊界上的壓力。根據蒙皮浸入邊界對流體網格施加的力的x分量，即仿真結果中x方向的歐拉力f1，計算柔性蒙皮的阻力。柔性蒙皮阻力FR的計算公式為

3.3.1 彈性支撐間距對蒙皮減阻性能的影響

設置蒙皮無量綱結構參數為：m1=0.4、m2=0.4、k1=0.4、k2=0.4，對雙層柔性蒙皮模型進行了數值計算，得到不同間距l下蒙皮的阻力系數CD如圖10所示。同時，為更加直觀地反映蒙皮的減阻效果，在相同流場環境下對剛性平板的阻力系數進行了數值計算。由圖可知，l=1～4 時，雙層柔性蒙皮的阻力系數CD隨間距l的增大而增大，且l=7 時的阻力系數最小，l=9 時的阻力系數最大。值得注意的是，l=1～8 時，雙層柔性蒙皮的阻力系數CD均小于剛性平板，表明蒙皮在此間距范圍內具備減阻效果。

圖10 不同間距l 對應的阻力系數

為探究間距l對阻力系數CD的影響機理，給出雙層柔性蒙皮近壁面速度分布曲線如圖11 所示。由圖可知，l=4 時的近壁面流速比l=1 時的要大，且l=7 時近壁面流速最小，l=9 時近壁面流速最大。而近壁面流速降低會使黏性阻力減小，進而印證了圖10 所示的蒙皮阻力系數CD隨間距l變化的規律。

圖11 雙層柔性蒙皮近壁面流速

下面從結構變形的角度探究蒙皮阻力系數CD隨間距l變化的規律。給出雙層柔性蒙皮表面流速矢量圖以及速度流線如圖12 和圖13 所示。由圖12 可知，l=1、4 時，雙層柔性蒙皮幾乎不產生明顯的變形。l=7 時，雙層柔性蒙皮發生類似于波浪型的變形，且在其下凹處存在著穩定、低速的漩渦(圖13)。根據胡海豹等[14]的研究結果表明，穩定、低速的漩渦可以起到類似于滾動軸承的作用，將流體之間的滑動摩擦轉變為滾動摩擦，從而減小了蒙皮受到的摩擦阻力。l=9 時，由于支撐間距較大，雙層柔性蒙皮結構在流體的作用下被破壞，因而使其阻力系數急劇增大。

圖12 不同間距下雙層柔性蒙皮表面流速矢量

圖13 雙層柔性蒙皮表面速度流線

3.3.2 彈性支撐剛度比對蒙皮減阻性能的影響

設置蒙皮無量綱結構參數為l=7、m1=0.4、m2=0.4、k2=0.4，對雙層柔性蒙皮模型進行了數值計算。得到不同剛度比γ下蒙皮的阻力系數CD如圖14 所示。由圖可知，γ=0.25～1.5 時，阻力系數CD隨著剛度比γ的增大逐漸減小，而繼續增大剛度比γ=2 時，雙層柔性蒙皮的阻力系數幾乎未發生變化。為探究剛度比γ對阻力系數CD影響的內在機理，給出圖15 所示的雙層柔性蒙皮近壁面流速圖。由圖可知，γ=0.25～1.5 時，雙層柔性蒙皮近壁面流速隨剛度比的增大而減??；而繼續增大剛度比γ=2 時，蒙皮近壁面流速幾乎不變。這與圖14 所示的蒙皮阻力系數CD隨剛度比γ的變化規律一致。進而說明了在一定范圍內增大剛度比可以減小蒙皮近壁面流速，從而使阻力系數減小。

圖14 不同剛度比γ 對應的阻力系數

圖15 蒙皮近壁面流速

3.3.3 流速對蒙皮減阻性能的影響

設置蒙皮無量綱結構參數為l=7、m1=0.4、m2=0.4、k1=0.4，k2=0.4，分別在流速為2、5、10、20 m/s 情況下，對雙層柔性蒙皮進行建模以及仿真計算。得到相對減阻率η隨流速u的變化情況如圖16 所示。由圖可知，u=2、5 m/s 時，蒙皮的相對減阻率隨流速的增大而增大，且最大相對減阻率可達14.48%。當流速u增大到10 m/s 時，蒙皮的相對減阻率降低至5.61%，而繼續增大流速u至20 m/s 時，相對減阻率為-13.36%，表明蒙皮在此流速下不具備減阻性能。

圖16 相對減阻率隨流速的變化

為探究流速對蒙皮相對減阻率影響的內在機理，給出不同流速下的近壁面流速曲線如圖17 所示。由圖可知，u 為2、5 m/s 時，雙層柔性蒙皮近壁面流速均值小于剛板近壁面流速均值，表明雙層柔性蒙皮在此流速下具備減阻效果，且蒙皮近壁面平均流速減小率隨外加流場流速的增大而增大；外加流場流速增大至u=10 m/s 時，平均流速減小率降低；而增大外加流場流速至u=20 m/s時，剛板近壁面流速均值小于雙層柔性蒙皮流速均值，說明此時雙層柔性蒙皮不具備減阻效果。這與圖16 所示的相對減阻率隨流速的變化規律一致。進而說明了在一定范圍內增大流速會使雙層柔性蒙皮的減阻性能先提高后降低。

圖17 不同流速下蒙皮和剛板近壁面流速分布曲線

4 結論

本文從水下航行器流體減阻的設計需求出發，將雙層隔振系統微小化、陣列化，結合柔性壁面減阻技術構建出了一種雙層柔性蒙皮結構，并基于浸入邊界法對其減阻性能進行了仿真計算，結論如下：

1）間距l=1～8 時，雙層柔性蒙皮的阻力系數小于剛板的阻力系數，表明雙層柔性蒙皮結構在特定的間距下具備減阻效果。

2）在一定范圍內增大剛度比會使近壁面流速降低，從而起到減小黏性阻力、提高蒙皮減阻性能的作用。

3）特定結構參數的雙層柔性蒙皮在13 m/s 的流速范圍內具備減阻效果，對于更高流速下的減阻性能還有待對其結構參數優化后進一步的探究。