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

任洵濤，趙丹，陳輝，崔進(jìn)

1. 中國船舶集團(tuán)第七〇四研究所，上海 200031

2. 哈爾濱工程大學(xué) 機(jī)電工程學(xué)院，黑龍江 哈爾濱 150001

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

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

近年來，雙層隔振系統(tǒng)廣泛應(yīng)用于艙室內(nèi)的減振研究中[6]，雙層隔振系統(tǒng)是在單層隔振系統(tǒng)的基礎(chǔ)上增加一個(gè)中間質(zhì)量塊及一層彈性支撐形成的，即雙層隔振系統(tǒng)由2 層彈簧阻尼隔振器和隔振器之間的中間質(zhì)量塊組成。相比于單層隔振系統(tǒng)，雙層隔振系統(tǒng)能在較寬的頻率范圍內(nèi)發(fā)揮更好的減振效果，現(xiàn)已被廣泛應(yīng)用在汽車、航天和船舶等領(lǐng)域，是振動(dòng)控制系統(tǒng)研究的重點(diǎn)領(lǐng)域。

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

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

1.1 雙層柔性蒙皮物理模型

本文設(shè)計(jì)的柔性蒙皮結(jié)構(gòu)由柔性壁面和彈性支撐組成，彈性支撐則是由中間質(zhì)量塊和上下層彈簧及阻尼構(gòu)成，其物理模型如圖1 所示。

圖1 柔性蒙皮物理模型

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

1.2 浸入邊界法的控制方程

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

浸入邊界法控制方程以N-S 方程為出發(fā)點(diǎn)，描述了具有慣性的牛頓流體的運(yùn)動(dòng)。包括雙層柔性蒙皮和流體在內(nèi)的整個(gè)物理區(qū)域，都假定為不可壓縮牛頓流體的黏性流動(dòng)，其連續(xù)性和動(dòng)量方程可分別表示為

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

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

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

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

正則化的增量函數(shù)為

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

1.3 雙層柔性蒙皮仿真模型

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

1)柔性壁面仿真模型

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

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

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

2)彈性支撐仿真模型

柔性蒙皮內(nèi)層的彈性支撐是分布式雙層隔振系統(tǒng)。使用建模單元中的拉格朗日點(diǎn)、阻尼彈簧和質(zhì)量單元建立彈性支撐，如圖3 所示。

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

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

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

2 算例驗(yàn)證

2.1 圓柱繞流

圓柱繞流是流體力學(xué)中的經(jīng)典問題，對其有許多的理論以及實(shí)驗(yàn)研究[10-12]，方便對仿真結(jié)果進(jìn)行驗(yàn)證。因此，本文按照文獻(xiàn)[10-12]中的圓柱繞流試驗(yàn)參數(shù)進(jìn)行了圓柱繞流仿真計(jì)算，在仿真與試驗(yàn)參數(shù)一致的情況下，驗(yàn)證了IB2d 數(shù)值計(jì)算結(jié)果的準(zhǔn)確性。

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

表1 斯特勞哈爾數(shù)比較

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

圖5 升力系數(shù)曲線

圖6 阻力系數(shù)曲線

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

表2 阻力系數(shù)對比

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

圖7 圓柱繞流尾渦狀態(tài)

綜上，根據(jù)圓柱繞流算例中的升力系數(shù)、阻力系數(shù)、斯特勞哈爾數(shù)以及圓柱繞流尾渦狀態(tài)，驗(yàn)證了IB2d 數(shù)值計(jì)算結(jié)果的準(zhǔn)確性。

2.2 柔性絲線

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

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

圖8 柔性絲線計(jì)算模型示意

表3 仿真參數(shù)設(shè)置

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

圖9 柔性絲線仿真計(jì)算與實(shí)驗(yàn)數(shù)據(jù)對比

3 雙層柔性蒙皮仿真結(jié)果分析

3.1 仿真參數(shù)設(shè)置

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

3.2 仿真參數(shù)可行性驗(yàn)證

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

表4 不同網(wǎng)格密度下的仿真結(jié)果

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

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

表5 不同時(shí)間步長下的仿真結(jié)果

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

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

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

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

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

圖10 不同間距l(xiāng) 對應(yīng)的阻力系數(shù)

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

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

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

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

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

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

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

圖14 不同剛度比γ 對應(yīng)的阻力系數(shù)

圖15 蒙皮近壁面流速

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

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

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

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

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

4 結(jié)論

本文從水下航行器流體減阻的設(shè)計(jì)需求出發(fā)，將雙層隔振系統(tǒng)微小化、陣列化，結(jié)合柔性壁面減阻技術(shù)構(gòu)建出了一種雙層柔性蒙皮結(jié)構(gòu)，并基于浸入邊界法對其減阻性能進(jìn)行了仿真計(jì)算，結(jié)論如下：

1）間距l(xiāng)=1～8 時(shí)，雙層柔性蒙皮的阻力系數(shù)小于剛板的阻力系數(shù)，表明雙層柔性蒙皮結(jié)構(gòu)在特定的間距下具備減阻效果。

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

3）特定結(jié)構(gòu)參數(shù)的雙層柔性蒙皮在13 m/s 的流速范圍內(nèi)具備減阻效果，對于更高流速下的減阻性能還有待對其結(jié)構(gòu)參數(shù)優(yōu)化后進(jìn)一步的探究。