基于ALE方法的柔性飄帶阻力特性研究

譚涵林，阮文俊，步鵬飛，孫鵬飛

(南京理工大學能源與動力工程學院，江蘇 南京 210094)

0 引言

隨著空氣動力學以及外彈道測試技術(shù)的不斷發(fā)展,彈丸受力和飛行的研究日趨深入,武器在提高射程的同時,還需不斷提高射擊精度,這對彈丸飛行穩(wěn)定性有極高的要求。在彈丸尾部安裝柔性飄帶是保證彈丸穩(wěn)定飛行的一種方式,圖1為典型的柔性飄帶穩(wěn)定彈丸-M42子彈[1]。柔性飄帶在提供穩(wěn)定力矩的同時,具有便于大批量集裝、存放、運輸和成本低等優(yōu)點[2],是一種極其實用的彈丸飛行穩(wěn)定裝置。

對柔性飄帶穩(wěn)定彈丸進行研究的關(guān)鍵在于對飄帶流固耦合過程的研究,由于柔性飄帶極易變形,其在氣流中的擺動是復(fù)雜的流固耦合過程,因此對于柔性飄帶穩(wěn)定彈丸的數(shù)值研究較少,大多為理論分析和實驗研究。

Dahlke等[1]通過風洞試驗和自由飛行試驗,研究了飄帶的長度、寬度、厚度和數(shù)量等參數(shù)對彈丸飛行的影響。Auman等[2]通過水平風洞試驗證明了飄帶阻力是影響彈丸穩(wěn)定飛行的關(guān)鍵,而飄帶的阻力與飄帶的寬度有很大關(guān)系。唐良銳[3]進行了柔性飄帶彈的氣動特性試驗,結(jié)果表明增加飄帶長度、寬度會增強穩(wěn)定效果,但飄帶過長會導(dǎo)致飄帶扭轉(zhuǎn)、擺動加劇,不能獲得可靠的穩(wěn)定力矩。張維全[4-6]在試驗與理論分析的基礎(chǔ)上,建立了飄帶式子彈氣動力的工程計算公式,提出了柔性飄帶穩(wěn)定器的設(shè)計原則。孫傳寶等[7-9]通過風洞實驗研究了柔性旗幟的流固耦合過程。

先前大量的實驗研究了柔性飄帶對彈丸穩(wěn)定的影響,但實驗方法成本極高,要耗費大量人力、物力及財力。而數(shù)值計算由于其經(jīng)濟性和可重復(fù)性已經(jīng)成為一種重要的研究方法,因此需要找到一種合適的數(shù)值方法對柔性飄帶穩(wěn)定彈丸進行研究,以降低研發(fā)成本。文中基于任意拉格朗日-歐拉(ALE)方法,對柔性飄帶的流固耦合過程進行計算。ALE方法兼有拉格朗日和歐拉方法的優(yōu)點,能夠有效地防止網(wǎng)格畸變,在處理柔性大變形物體具有很大的優(yōu)勢,ALE方法已經(jīng)廣泛應(yīng)用于柔性降落傘充氣過程的研究[10-15]。文中利用ALE方法對柔性飄帶流固耦合過程進行求解,將求解結(jié)果與風洞實驗結(jié)果進行對比,驗證了方法的有效性,可為柔性飄帶彈丸的研發(fā)提供參考。

1 數(shù)學模型

1.1 柔性飄帶穩(wěn)定彈丸原理

柔性飄帶穩(wěn)定彈丸的原理如圖2所示,由于飄帶質(zhì)地柔軟,其對彈丸的阻力大體沿著來流方向,且作用于彈丸尾部與飄帶的連接處,當彈丸有攻角α時,飄帶阻力F在彈軸法線方向的分量Fy對質(zhì)心產(chǎn)生使攻角減小的穩(wěn)定力矩Mz,從而起到穩(wěn)定彈丸的作用。

圖2 柔性飄帶提供穩(wěn)定力矩的原理Fig.2 The principle of flexible ribbon to provide stabilizing moment

1.2 飄帶材料本構(gòu)模型

柔性飄帶采用FABRIC織物材料模型,該模型可用于具有非線性動力學特性的柔性織物,常用于模擬降落傘以及安全氣囊,其本構(gòu)關(guān)系[16]為:

(1)

式中:σ1、σ2為材料縱向應(yīng)力和橫向應(yīng)力;v1、v2為材料縱向泊松比和橫向泊松比;E1、E2為材料縱向彈性模量和橫向彈性模量;τ12為剪切應(yīng)力;G12為剪切模量;α為非線性剪切應(yīng)力系數(shù)。當把材料近似為各向同性彈性材料時,則只需定義一個楊氏模量和泊松比。

1.3 控制方程

1.3.1 流場控制方程

在ALE方法描述中,拉格朗日和歐拉坐標間的關(guān)系[17]為:

(2)

式中:Xi為拉格朗日坐標;xi為歐拉坐標;ui為流體速度;wi為網(wǎng)格點移動速度。

則在ALE算法中,不可壓縮流體的N-S方程[17]為:

(3)

式中:ρ為流體密度;e為物質(zhì)內(nèi)能;σij為j方向的應(yīng)力作用于垂直于i軸的平面上;fi為單位體積力。

1.3.2 拉格朗日結(jié)構(gòu)控制方程

飄帶是具有典型非線性動力學特性的柔性大變形體,在固體邊界上,其動力學控制方程[18]為:

(4)

式中:ρs為結(jié)構(gòu)材料密度;yi為結(jié)構(gòu)域中的節(jié)點位移。

1.4 流體與結(jié)構(gòu)的耦合

對上述控制方程進行耦合計算,采用罰函數(shù)法對飄帶與流體間節(jié)點力信息進行傳遞,采用顯式時間法,對于任意時間步tn,流體節(jié)點與結(jié)構(gòu)節(jié)點的相對位移為dn,則dn的迭代方程[16]為:

(5)

(6)

式中:vr為流體節(jié)點與結(jié)構(gòu)節(jié)點的相對速度;vs為結(jié)構(gòu)節(jié)點速度;vf為流體節(jié)點速度。

罰函數(shù)法追蹤結(jié)構(gòu)節(jié)點和流體物質(zhì)位置間的相對位移d。檢查每一個結(jié)構(gòu)節(jié)點對流體物質(zhì)表面的貫穿,如果貫穿,界面力F就會分布到流體單元的節(jié)點上:

Fn=kidn

(7)

式中:ki為流體節(jié)點與結(jié)構(gòu)節(jié)點質(zhì)量模型特征的剛度系數(shù)。

2 數(shù)值方法驗證

2.1 數(shù)值方法驗證模型

孫傳寶等[7-9]對柔性飄帶進行了風洞測力實驗,且詳細描述了飄帶的材料參數(shù),實驗中所用柔性飄帶厚度為0.15 mm,密度為1 380 kg/m3,彈性模量約為1.8 GPa。根據(jù)實驗建立柔性飄帶以及外流場仿真模型,利用Hypermesh進行網(wǎng)格劃分。飄帶使用薄殼單元和FABRIC材料,采用四邊形網(wǎng)格進行劃分,為防止飄帶變形過大而發(fā)生網(wǎng)格穿透,需對飄帶設(shè)置自接觸;流場使用實體單元和NULL材料,采用六面體網(wǎng)格進行劃分,對靠近飄帶的流場網(wǎng)格進行細化以保證計算結(jié)果準確,流場邊界設(shè)置為無反射邊界條件。模型及網(wǎng)格如圖3、圖4所示。

圖3 飄帶模型及網(wǎng)格劃分Fig.3 Ribbon model and grid division

圖4 流場網(wǎng)格劃分以及飄帶(紅色)在流場中的位置Fig.4 Grid division of flow field and position of ribbon (red part) in flow field

實驗研究了長寬比以及來流風速對柔性飄帶阻力的影響。根據(jù)實驗結(jié)果,分別進行了3組柔性飄帶的仿真計算,以驗證數(shù)值方法的有效性。3組仿真計算的工況如表1所示。仿真計算的氣壓為1 atm,溫度為20 ℃,空氣密度為1.205 kg/m3。對每組飄帶設(shè)置阻力監(jiān)測以及位移監(jiān)測,每計算0.005 s儲存一次數(shù)據(jù),對每組飄帶進行時長為0.8 s的仿真計算。

表1 仿真計算工況Table 1 Simulation calculation conditions

2.2 數(shù)值方法驗證結(jié)果

計算結(jié)果與文獻實驗結(jié)果對比如表2所示。數(shù)值計算所得時均阻力Fx以及振幅A與實驗結(jié)果吻合較好;而飄帶擺動頻率相較實驗結(jié)果偏小,但整體變化趨勢相同,造成這種現(xiàn)象的原因可能是數(shù)值計算中沒有考慮重力等因素的影響,此外模擬中溫度、氣壓等大氣條件與實驗可能有所偏差。

表2 仿真結(jié)果與實驗數(shù)據(jù)對比Table 2 Comparison between simulation results and experimental data

圖5為飄帶尾部位移與時間的關(guān)系,可以看出數(shù)值計算能夠一定程度上模擬出飄帶擺動的周期性;其中圖5(a)為數(shù)值模擬中對工況2飄帶尾部節(jié)點進行位移監(jiān)測所得位移-時間曲線,圖5(b)為實驗工況下飄帶尾部位移-時間圖像,圖中x1為實驗中的位移監(jiān)測點。由表2及圖5可知,利用數(shù)值方法對飄帶阻力及飄帶擺動進行模擬具有一定的可靠性。

3 仿真模型與結(jié)果

3.1 仿真模型

使用常見的軍用織物材料為基準進行仿真研究,其密度為555 kg/m3,彈性模量為0.43 GPa,厚度為0.1 mm。通過修改其材料參數(shù)以及邊界條件,研究來流速度、材料密度、長寬比、材料剛度以及飄帶形狀對飄帶阻力的影響。仿真計算中所有工況的大氣條件均相同,氣壓為1 atm,溫度為20 ℃,空氣密度為1.205 kg/m3。

仿真模型與2.1節(jié)所述類似,飄帶使用薄殼單元和FABRIC材料,采用四邊形網(wǎng)格進行劃分,設(shè)置自接觸防止網(wǎng)格穿透;流場使用實體單元和NULL材料,采用六面體網(wǎng)格進行劃分,由于飄帶彈所用飄帶長寬比較大,需加大飄帶擺動方向(z方向)的流場,設(shè)置為飄帶長度的4倍,對靠近飄帶的流場網(wǎng)格進行細化,流場邊界設(shè)置為無反射邊界條件。彈體尾部飄帶一般有3種形狀—單根飄帶、雙根飄帶以及環(huán)狀飄帶。不同形狀的飄帶模型以及流場網(wǎng)格劃分如圖6、圖7所示。

圖6 飄帶形狀Fig.6 Shape of ribbon

圖7 流場網(wǎng)格劃分Fig.7 Grid division of flow field

3.2 仿真結(jié)果

3.2.1 單根飄帶阻力特性

為方便氣動分析,將阻力轉(zhuǎn)化為阻力系數(shù),阻力系數(shù)cd的計算公式為:

(8)

式中:Fx為阻力值;ρ為空氣密度;v為來流速度;S為飄帶面積。

首先在恒定單根飄帶的彈性模量以及長寬比情況下,研究來流速度以及飄帶密度對單根飄帶阻力特性的影響。飄帶的尺寸恒定為0.2 m×0.08 m,厚度為0.1 mm,彈性模量恒定為0.43 GPa,采用與2.1節(jié)相同的數(shù)據(jù)采集方式,分別得到了不同來流速度以及不同密度下單根飄帶的阻力系數(shù),計算結(jié)果如圖8所示。由圖可知,阻力系數(shù)總體上隨著飄帶密度的增大而增大,隨著來流速度的增大而減小。

圖8 來流速度和材料密度對阻力系數(shù)的影響Fig.8 Effect of inflow velocity and material density on drag coefficient

隨后將來流速度恒定為30 m/s,飄帶密度恒定為555 kg/m3,飄帶寬度恒定為0.08 m,計算不同長寬比以及彈性模量下單根飄帶的阻力系數(shù),計算結(jié)果如圖9所示。由圖可知,阻力系數(shù)總體上隨彈性模量的增大而增大,隨飄帶長寬比的增大而減小。

圖9 長寬比和彈性模量對阻力系數(shù)的影響Fig.9 Effect of length-width ratio and elastic modulus on drag coefficient

圖中有部分工況與總體規(guī)律不同,這是因為發(fā)生了如圖10所示的扭轉(zhuǎn)現(xiàn)象,在多篇文獻的實驗中也會發(fā)生飄帶扭轉(zhuǎn)現(xiàn)象[1,2,7],實驗以及計算結(jié)果均顯示扭轉(zhuǎn)現(xiàn)象的發(fā)生會降低阻力。

圖10 飄帶擺動圖像Fig.10 Image of ribbon swing

3.2.2 雙根飄帶以及環(huán)狀飄帶的阻力特性

對不同來流速度下雙根飄帶以及環(huán)狀飄帶進行仿真計算,飄帶間的間距均為0.04 m,二者單邊飄帶的尺寸均為0.2 m×0.08 m,厚度為0.1 mm,材料密度為555 kg/m3,彈性模量為0.43 GPa。

在不同來流速度下兩種飄帶的阻力系數(shù)如圖11所示,雙根飄帶的阻力系數(shù)隨速度的增大而減小,速度每增大10 m/s,雙根飄帶阻力系數(shù)平均減少8.2%;環(huán)狀飄帶的阻力系數(shù)隨速度的增大而略微增大,速度每增大10 m/s,環(huán)狀飄帶阻力系數(shù)平均增大2.6%。

圖11 雙根飄帶與環(huán)狀飄帶阻力系數(shù)對比Fig.11 Comparison of drag coefficients between double ribbons and annular ribbon

將二者阻力系數(shù)與相同條件下的單根飄帶阻力系數(shù)進行比較,結(jié)果如圖12所示,其中環(huán)狀飄帶阻力系數(shù)相較于單根飄帶阻力系數(shù)更大,雙根飄帶阻力系數(shù)相較于單根飄帶阻力系數(shù)更小。

圖12 不同形狀飄帶阻力系數(shù)對比Fig.12 Comparison between the drag coefficient of ribbon with different shapes

3.3 結(jié)果分析

3.3.1 單根飄帶結(jié)果分析

根據(jù)前人的實驗研究以及理論分析,飄帶的阻力由表面摩擦阻力以及飄帶波動阻力構(gòu)成,其中飄帶波動阻力占主體部分[4]。

對于單根飄帶的情況,來流速度的增大會同時增加飄帶表面摩擦以及飄帶波動,使得飄帶阻力增大。圖13為飄帶尺寸為0.2 m×0.08 m,厚度為0.1 mm,材料密度為555 kg/m3,彈性模量為0.43 GPa時,飄帶表面摩擦阻力及總阻力隨來流速度的變化曲線,可以發(fā)現(xiàn),來流速度的增大會增大飄帶表面摩擦阻力以及總阻力值,飄帶表面摩擦阻力約為總阻力的10.4%,這與文獻中所述的飄帶波動阻力占總阻力的主體部分相符。根據(jù)阻力計算結(jié)果,飄帶阻力隨來流速度的增大近似呈線性關(guān)系,但根據(jù)式(8),無量綱阻力系數(shù)基于速度的二次方,因此阻力系數(shù)隨來流速度的增大而減小。速度每增大10 m/s,飄帶阻力系數(shù)平均降低9.7%。

圖13 飄帶表面摩擦阻力與總阻力隨速度的變化Fig.13 Changes in surface friction drag and total drag of ribbon with Inflow velocity

由于飄帶波動由上游傳輸?shù)较掠?飄帶尾部波動最為劇烈,尾部的波動阻力占整體波動阻力的主體部分,彈性模量較大的飄帶,參與尾部波動的面積更大,會產(chǎn)生更大的波動阻力[19];在長寬比恒定為2.5,來流速度恒定為30 m/s,飄帶密度恒定為555 kg/m3時,根據(jù)飄帶擺動仿真圖像以及節(jié)點位移測量,彈性模量為0.23 GPa的飄帶,振幅超過0.03 m的區(qū)域約占飄帶總面積的30%,而彈性模量為0.63 GPa的飄帶,振幅超過0.03 m的區(qū)域約占飄帶總面積的42%,其阻力系數(shù)相較于彈性模量為0.23 GPa的飄帶增大了17%;因此飄帶阻力系數(shù)隨彈性模量的增大而增大。彈性模量每增加0.1 GPa,阻力系數(shù)平均增大3.9%。

此外,長寬比越大的飄帶,其尾部波動的面積占整體飄帶面積的比重越小,在彈性模量恒定為0.43 GPa,來流速度恒定為30 m/s,飄帶密度恒定為555 kg/m3時,根據(jù)飄帶擺動仿真圖像以及節(jié)點位移測量,長寬比為1.5的飄帶,其振幅超過0.025 m(最大振幅的一半)的區(qū)域約占飄帶總面積的47%,而長寬比為3的飄帶,其振幅超過0.032 m(最大振幅的一半)的區(qū)域約占飄帶總面積的26%,其阻力系數(shù)相較于長寬比為1.5的飄帶下降了40%;因此長寬比越大的飄帶,基于飄帶面積的阻力系數(shù)越小。長寬比每增加0.5,阻力系數(shù)平均減小17.6%。

當飄帶彈性模量、長寬比以及來流速度相同時,飄帶在氣流中的波動程度相近,但密度更高的飄帶,參與波動的質(zhì)量更大。如圖14所示,根據(jù)飄帶擺動仿真圖像以及節(jié)點位移測量,在彈性模量恒定為0.43 GPa,來流速度恒定為30 m/s,長寬比恒定為2.5時,密度為755 kg/m3的飄帶振幅約為0.055 m,密度為455 kg/m3的飄帶振幅約為0.053 m,兩種飄帶的振幅較為相近,但密度更大的飄帶參與尾部波動的質(zhì)量也越大,因此密度更大的飄帶阻力系數(shù)越大。密度每增加100 kg/m3,飄帶阻力系數(shù)平均增大16.3%。

圖14 不同密度的飄帶尾部位移曲線Fig.14 Tail displacement curves of ribbons with different densities

而當發(fā)生飄帶扭轉(zhuǎn)現(xiàn)象時,會使飄帶波動嚴重減緩,圖9中長寬比為3,彈性模量為0.63 GPa的飄帶發(fā)生了扭轉(zhuǎn)現(xiàn)象,根據(jù)飄帶擺動仿真圖像以及節(jié)點位移測量,發(fā)生扭轉(zhuǎn)時飄帶擺動的最大振幅約下降39%,使其阻力系數(shù)相較彈性模量為0.53 GPa的飄帶下降了24%,由此可見,飄帶的扭轉(zhuǎn)會使阻力明顯降低。

3.3.2 雙根飄帶以及環(huán)狀飄帶結(jié)果分析

來流速度為30 m/s時,3種飄帶氣流速度云圖如圖16所示。對于雙根飄帶,由于氣流速度在雙根飄帶之間減速較大,由30 m/s下降至17 m/s左右,導(dǎo)致飄帶波動減緩以及表面摩擦減少,因此雙根飄帶的阻力系數(shù)小于相同條件下的單根飄帶阻力系數(shù);如圖15所示,對于環(huán)狀飄帶,其環(huán)形部迎風面和背風面具有壓差,會形成額外的環(huán)形部壓差阻力[2],且速度越大,產(chǎn)生的環(huán)形部壓差阻力越大,因此環(huán)狀飄帶的阻力系數(shù)大于相同條件下的單根飄帶阻力系數(shù),且隨著來流速度的增加,環(huán)形部壓差阻力對飄帶整體阻力的影響越大,環(huán)狀飄帶阻力增長的斜率也隨之增大,使得環(huán)狀飄帶阻力系數(shù)隨來流速度的增大而略微增大。

圖15 環(huán)狀飄帶環(huán)形部壓差阻力與總阻力隨速度的變化Fig.15 Changes in pressure difference drag of the annular part and total drag of the annular ribbon with speed

圖16 飄帶附近流場速度云圖(xz平面)Fig.16 Velocity nephogram of flow field near ribbon (xz plane)

4 結(jié)論

利用ALE方法,對柔性飄帶的阻力特性進行仿真模擬,計算了不同來流速度、不同材料屬性以及不同形狀下飄帶的阻力特性,得出以下結(jié)論:

1)通過與文獻風洞實驗對比,驗證了ALE方法能夠有效預(yù)測飄帶阻力值及飄帶擺動的周期性。

2)單根飄帶的阻力系數(shù)隨來流速度以及飄帶長寬比的增大而減小;速度每增大10 m/s,飄帶阻力系數(shù)平均降低9.7%,長寬比每增加0.5,阻力系數(shù)平均減小17.6%。單根飄帶的阻力系數(shù)隨飄帶材料密度以及彈性模量的增大而增大,彈性模量每增加0.1 GPa,阻力系數(shù)平均增大3.9%,密度每增加100 kg/m3,飄帶阻力系數(shù)平均增大16.3%。這些因素的改變會影響飄帶表面摩擦以及波動面積,從而影響飄帶阻力系數(shù)。

3)雙根飄帶的阻力系數(shù)隨速度的增大而減小,速度每增大10 m/s,雙根飄帶阻力系數(shù)平均減少8.2%;環(huán)狀飄帶的阻力系數(shù)隨速度的增大而略微增大,速度每增大10 m/s,環(huán)狀飄帶阻力系數(shù)平均增大2.6%。

4)由于存在環(huán)形部壓差阻力,環(huán)狀飄帶的阻力系數(shù)大于單根飄帶阻力系數(shù);由于雙根飄帶間的氣流減速,雙根飄帶的阻力系數(shù)小于單根飄帶阻力系數(shù)。

5)飄帶的扭轉(zhuǎn)會使飄帶波動嚴重減緩,飄帶波動面積大幅減小,從而降低飄帶阻力。