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

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

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

0 引言

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

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

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

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

1 數學模型

1.1 柔性飄帶穩定彈丸原理

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

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

1.2 飄帶材料本構模型

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

(1)

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

1.3 控制方程

1.3.1 流場控制方程

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

(2)

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

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

(3)

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

1.3.2 拉格朗日結構控制方程

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

(4)

式中:ρs為結構材料密度;yi為結構域中的節點位移。

1.4 流體與結構的耦合

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

(5)

(6)

式中:vr為流體節點與結構節點的相對速度;vs為結構節點速度;vf為流體節點速度。

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

Fn=kidn

(7)

式中:ki為流體節點與結構節點質量模型特征的剛度系數。

2 數值方法驗證

2.1 數值方法驗證模型

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

圖3 飄帶模型及網格劃分Fig.3 Ribbon model and grid division

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

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

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

2.2 數值方法驗證結果

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

表2 仿真結果與實驗數據對比Table 2 Comparison between simulation results and experimental data

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

3 仿真模型與結果

3.1 仿真模型

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

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

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

圖7 流場網格劃分Fig.7 Grid division of flow field

3.2 仿真結果

3.2.1 單根飄帶阻力特性

為方便氣動分析,將阻力轉化為阻力系數,阻力系數cd的計算公式為:

(8)

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

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

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

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

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

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

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

3.2.2 雙根飄帶以及環狀飄帶的阻力特性

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

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

圖11 雙根飄帶與環狀飄帶阻力系數對比Fig.11 Comparison of drag coefficients between double ribbons and annular ribbon

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

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

3.3 結果分析

3.3.1 單根飄帶結果分析

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

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

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

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

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

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

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

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

3.3.2 雙根飄帶以及環狀飄帶結果分析

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

圖15 環狀飄帶環形部壓差阻力與總阻力隨速度的變化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 結論

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

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

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

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

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

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