考慮彈-筒柔性的潛載導彈彈射多體動力學分析方法

張子豪,王漢平

(北京理工大學 宇航學院,北京 100081)

0 引言

潛射彈道導彈具有隱蔽性好、生存力強、突防性佳、射程遠、破壞性大等特點,是大國博弈的鎮山之寶,一直深受青睞,其發射方式一般采用燃氣-蒸汽式垂直彈射。在垂直彈射過程中,其彈道過程通常分為入水段、水彈道、出水段、空中段,其中前3 個階段為潛射導彈所獨有,也是關系其發射成敗的關鍵。導彈和發射筒之間借助多圈適配器、密封環支撐,當導彈彈射入水時,導彈慣性較大,彈射過程中所受到的彈射力、重力、附連水慣性力以及水動阻力將導致較大的軸向載荷,對彈體剛性以及橫向振動會產生不可忽視的影響。另外,由于載體運動以及導彈相對載體的運動,導彈還受橫向流、波浪、適配支撐等條件影響,最終彈體和發射筒之間會產生相對旋轉運動,因此水下彈射發射過程的動力學問題是一個復雜多場耦合且高度非線性的問題。

目前國內針對該問題的研究有兩種思路:其一是對流體作用進行簡化處理,列出導彈運動過程的顯式方程,并結合計算流體動力學(CFD)軟件的用戶自定義函數(UDF)開發進行仿真分析;其二是對彈體與發射筒的接觸作用進行簡化處理,基于有限元技術開展流體與固體耦合的仿真研究。文獻[5 -14]采用多相流模型,在考慮空化的條件下運用UDF 和動網格技術對彈體運動過程進行了數值仿真,并討論了波浪、橫向流等因素對彈體姿態的影響,但這些研究都只將發射筒簡化為一種約束,沒有考慮實際彈體和發射筒間復雜的相互作用;文獻[15-17]利用牛頓定律列出了導彈的運動學和動力學方程組,并編程進行了仿真計算,但對水動力載荷及彈體與適配器間作用力的簡化過于簡單,與實際物理過程相差較大。文獻[18 -19]對水動力和附連水慣性力進行了較為詳細的表達,但是沒有考慮發射筒和彈體之間的相互作用,只是簡單認為發射筒對彈體施加了支撐力和摩擦力。程載斌等使用任意拉格朗日-歐拉(ALE)模型對導彈無攻角潛射過程進行了仿真,該模型無法適用于有攻角的情形,而實際的潛射狀態攻角是一個無法回避的因素。劉傳龍等使用切片法在考慮適配器作用的條件下進行了流體與固體耦合計算。裴譞等考慮了氣密環的柔性與剛性彈體的耦合作用,研究了氣密環的抗壓剛度和結構尺寸對彈體出筒過程中橫向振動的影響。尚書聰等面向不同艇速下適配器和氣密環-減震墊兩種橫向支撐方式開展了對比仿真,研究了支撐方式對彈體姿態的影響,但模型中只考慮了適配器的抗壓剛度,忽略了適配器抗彎剛度。

綜上所述可知,已有研究均沒有考慮發射筒的柔性影響。而彈體的柔性變形也將導致彈體的應力剛化,且發射筒柔性和彈體柔性將會產生交互耦合作用,這也是上述研究尚未涉及的問題。關于彈體與定向器間的相互作用,文獻[24 -26]基于Euler-Bernouli 理論建立了柔性梁的偏微分方程,并自編程實現了移動物體作用下柔性梁振動的仿真研究,但是模型過于簡單,與實際應用相去甚遠。近年來,多體動力學仿真軟件雖然具備了處理剛體與剛體、柔性體與剛體以及雙柔性體之間接觸建模的能力,但當復雜柔性體參與接觸時,收斂性、穩定性較差,而且計算量巨大,直接制約了其應用范圍。針對此問題,比較通用的方法是借助特定開發以簡化接觸建模,在保證精度和穩定性的同時大幅度降低計算量。Ambrogi 等基于ADAMS 平臺開發了一種柔性點線高副,在卡車通過柔性橋梁的耦合模型中進行了驗證,但是并未公開算法細節。王林鵬等基于三次樣條曲線開發了一種柔性點線高副,并對移動導彈激勵下的彈體發射架耦合振動模型進行了仿真分析,但該算法在由型值點坐標換算控制點坐標的過程中需要求解大型代數方程組,計算量較大。竇建中等開發了一種基于Hermite 插值函數的柔性點線高副,該算法省去了樣條曲線繁瑣的計算,但是應用對象是剛性彈體定向鈕與柔性梁間的點線約束,而非雙柔性結構。

本文面向潛射導彈彈體和發射筒的雙柔性結構,提出一種高效、實用的動力學建模分析方法,采用修正的Craige-Bampton 法構建兩結構的柔性,并基于Hermite 插值函數構建柔性點線高副,以實現柔性彈體與柔性發射筒間基于適配器的接觸約束關系。基于Hermite 形狀函數實現了適配器對發射筒的載荷等效作用、水動力以及附連水慣性力對彈體的等效簡化處理,構建了同時考慮彈體和發射筒柔性的潛載導彈垂直彈射發射動力學模型,并以某算例進行了仿真校驗,驗證了該方法的有效性和優越性。

1 彈射發射動力學建模

某型導彈彈體通過4 圈適配器與發射筒相接觸來實現適配導向約束,在發射起始時刻,導彈、發射筒和潛艇均以相同速度勻速航行。假設潛艇航行速度與洋流速度在同一平面,于是潛射導彈的彈射發射運動也可看作是一個平面運動。基于發射筒固連坐標系進行彈射發射動力學建模:原點位于發射筒筒底中心,導彈初始軸向為發射筒固連坐標系的軸方向,法向作為發射筒固連坐標系的軸方向,同時軸方向也是坐標系的運動方向。全局坐標系的方向與發射筒固連坐標系方向一致。彈體和發射筒的參數都按實際參數設置,導彈的速度按照實際測得的數據進行加載驅動。彈體和發射筒均按柔性體處理,彈體和發射筒之間的接觸適配導向關系則通過柔性點線高副實現,其中柔性線以發射筒節點為型值點構造,彈體與適配器接觸部位設置啞物體,啞物體中心就是點線高副中的點(本文僅闡述建模仿真的思路,在此忽略洋流的影響)。圖1 所示為坐標系示意圖,圖2 所示為導彈結構及適配器分布圖。圖2 中,為彈體長度,為支撐面到彈尾距離,、、、分別為后適配器、第3 圈適配器、中適配器、上適配器到支撐面的距離。

圖1 坐標系示意圖Fig.1 Schematic diagram of coordinate systems

圖2 導彈結構及適配器分布圖Fig.2 Missile structure and adapter arrangement

1.1 發射筒和彈體的柔性化

彈體按離散質量點和等效梁相連接的方式進行簡化,劃分為41 段,梁參數由彈體截面半徑、當量厚度、截面當量彈性模量以及截面當量剪切模量予以等效;發射筒按梁單元模型予以簡化,其參數按發射筒的材料屬性及截面參數予以賦值。將兩部件用修正的Craige-Bampton 法進行柔性化處理,保留前20 階固定界面的正交模態,并記錄發射筒節點的初始坐標和對應的初始斜率,作為柔性線的初始型值點。其中,彈體的保留節點在4 圈適配器的掛接位置,發射筒的保留節點為發射筒在潛艇艙段中的懸掛處。

1.2 適配器作用的簡化

整圈適配器的力學性能表現為抗壓和抗彎剛度。適配器的變形情況可用其與彈體、發射筒相接觸的兩點與間相對位置關系來表示:dr為適配器的壓縮量,在坐標系中適配器的初始厚度方向與全局坐標系軸方向重合,因此就是整圈適配器的彎曲角度,按線性處理,則有

式中:F為沿方向的抗壓力;M為垂直平面與反向的抗彎力矩;為抗壓剛度,k為抗彎剛度,可通過實驗數據擬合得到。圖3 所示為適配器作用示意圖。圖3 中,x、y分別為適配器與彈體接觸截面的軸向和法向,x、y分別為適配器與發射筒接觸截面的軸向和法向。表1 所示為各適配器剛度。

表1 各圈適配器抗壓剛度和抗彎剛度Tab.1 Compressional and flexural stiffness es of adapters

圖3 適配器作用示意圖Fig.3 Adapter force between missile and canister

1.3 發射筒和彈體基于4 圈適配器的接觸約束

發射筒和彈體之間的真實狀態是借助4 圈適配器的接觸適配實現支撐和導向約束,在建模時按GForce 力元予以等效處理。簡單起見,僅以某一適配器為例對模型進行闡述(見圖4):由于動力學仿真ADAMS 軟件中基于模態中性文件的柔性體不支持浮動Marker,按柔性體處理的彈體和發射筒均不能直接用作力的反作用物體,必須先在彈體與適配器的接口位置(即彈體的保留節點)處構建一個體積、慣性均非常小的啞物體1,并將其固定在彈體上,從而借助啞物體就可將適配器的載荷傳遞到彈體上,并且也可借助該啞物體來檢測彈體適配器位置(即點)處的絕對位移和相對偏角。發射筒的節點坐標和對應位置處偏角可以實時提取,用以刷新柔性線,柔性線的動態變化就表征了發射筒的動態特性。建立啞物體2,啞物體2 與柔性線之間是點線高副約束(PTCV Joint),啞物體2 對應于適配器與發射筒的接觸適配點,即點(初始狀態與啞物體1 重合)。在啞物體1 與啞物體2 之間是在面約束(Inplane Joint),將啞物體2 約束在過啞物體1 中心且垂直于彈體軸線的平面上。提取啞物體1 與啞物體2 間的相對位移dr和相對轉角,即可利用自定義的Gforce 力元將適配器的作用力等效加載到發射筒和啞物體1 上。

圖4 彈體、適配器及發射筒間的力學模型拓撲結構Fig.4 Topological structure of mechanical model of missile,adapter and canister

1.3.1 柔性線的實現

柔性線的實現參考文獻[28 -29],通過CBKSUB 子程序實時提取發射筒節點的坐標及對應點的偏角數據,利用公共變量將坐標及偏角數據傳遞給CURSUB 子程序,在CURSUB 中用Hermite 插值函數定義柔性線。圖5 所示為彈體、發射筒間適配環節的數據流。

圖5 彈體、發射筒間適配環節的數據流Fig.5 Data flow diagram of contact between missile and canister

1.3.2 載荷的加載

有關載荷加載,以前述適配器為例,需要在柔性發射筒的每個節點處都建立一個按GFOSUB 定義的GForce 力元。在GFOSUB 子程序中,首先提取dr與,計算適配器的抗壓力和抗彎力矩。之后提取該節點坐標,根據坐標來判定該節點是否受力,如果受力,則從能量角度利用Hermite 形狀函數將適配器作用力等效分配在發射筒跨越啞物體2 的鄰近兩節點上(見圖6)。圖6 中,F為施加在節點上沿全局軸方向的作用力,F為沿全局軸方向的作用力,為垂直平面的力矩,F為力F在全局坐標系下軸方向投影,F為力F在全局坐標系下軸方向投影。形函數表達式為

圖6 載荷加載邏輯圖Fig.6 Logic diagram of adapter loading

式中:為無量綱參數,

為有限段的長度。

1.4 彈體水動力、附連水慣性力的定義及加載

對彈體作用力的加載是在單元局部坐標系Oxy中開展的,表示第個節點,x軸為單元軸,y軸為單元法向。使用Gforce 力元進行加載,受力物體為彈體,施力物體為大地,將水動力及附連水慣性力的分布等效簡化為單元兩節點的集中力。

1.4.1 水動力的定義

1.4.1.1 整彈的水動力

水動力和水動力矩是與導彈姿態密切相關的。其中入水過程導彈所受橫向阻力為

式中:為入(出)水深度占全長的比例;為導彈質心速度與軸向的夾角(攻角);為彈質心速度;ξ()為橫向水動力入水系數;C(,)為橫向水動力系數;為海水密度;為導彈的參考面積。

入水過程中,由橫向阻力產生的水動力矩為

式中:()為俯仰水動力矩入水系數;(,)為俯仰水動力矩系數。

入水過程軸向水動阻力為

式中:ξ()為軸向水動力入水系數;C(,)為軸向水動阻力系數。

在出水過程導彈所受橫向阻力為

出水過程中,軸向水動阻力為

1.4.1.2 水動力的分布

在入水段,橫向水動力在軸向方向上的分布為

水動力矩的產生主要是由于水動力載荷的軸向分布引起,因此在考慮橫向水動力載荷的軸向分布之后自然就包含了水動力矩載荷。

入水段軸向阻力在軸向的分布為

在出水段階段,橫向水動力在軸向的分布為

出水段軸向阻力在軸向方向上的分布為

考慮到艇體的橫向速度,彈體的攻角和速度模的平方表達為

式中:v、v、θ分別為發射筒坐標系下單元的橫向速度、軸向速度、角位移。

1.4.2 附連水質量慣性力的定義

入水過程中軸向附連水慣性力為

入水過程中橫向附連水慣性力為

式中:λ(,)、λ(,)分別為軸向和橫向附連水質量的入水系數;a、a分別為當地節點在其局部坐標系中的軸向、橫向加速度。

入水過程中軸向附連水質量慣性力在軸向的分布為

入水過程中橫向附連水質量慣性力在軸向方向上的分布為

出水過程中軸向附連水質量慣性力為

出水過程中橫向附連水質量慣性力為

出水過程中軸向附連水質量慣性力在軸向的分布為

出水過程中橫向附連水質量慣性力在軸向的分布為

1.4.3 分布力加載

彈體有限段的劃分利用1.1 節柔性化處理時進行的分段。有限元變形的描述仍采用Hermite 插值函數,但是坐標系使用有限段局部坐標系,形函數同(3)式。

假設單元12 段和23 段上節點1、節點2 和節點3處的分布力大小分別為、、,而節點處的分布力f的表達式在入水段為

而在出水段則為

式中:h為對應于節點的入(出)水高度比。

將分布力按線性分布予以近似,則在12 段上單元分布力表達式為

在23 段上,單元分布力表達式為

圖7 為有限元示意圖,其中,、分別為23 單元、12 單元實際入水部分占有限段長度的無量綱參數。

圖7 有限段示意圖Fig.7 Schematic diagram of finite element

把發射筒筒口坐標和作為研究對象的該有限段節點坐標、、進行比較,根據比較結果得出該有限段所處的位置。圖8 展示了有限段位置的4 種情況,其中(1)、(2)、(3)、(4)分別表示兩個單元都入水、右單元入水、左單元部分入水、右單元部分入水;兩個單元都在筒內。

圖8 不同位置有限段示意圖Fig.8 Schematic diagram of finite elements in different positions

根據所處位置,、的計算方式如下:

式中:、分別為12 段和23 段的長度。

根據1.3.2 節中利用Hermite 形狀函數對有限段變形的描述,可計算出入水段節點2 處等效的集中力與集中力矩:

使用GFOSUB 子程序定義GForce 力元,即可完成載荷的加載。

2 仿真校驗

對所有數據都進行無量綱化處理,處理方式如下。

無量綱時間為

式中:為彈體出筒時刻速度。

無量綱力為

式中:為彈體質量;為重力加速度。

無量綱力矩為

無量綱位移為

無量綱速度為

無量綱角度為

式中:為彈體俯仰角;為彈體出筒時刻質心的角速度。

無量綱角速度為

式中:為俯仰角速度。

初始時刻彈體與發射筒重合,計算時采用實測的速度-時間歷程作為導彈的驅動,艇速為3.5 kn,仿真計算導彈入水、水彈道、出水的全過程,結果如圖9～ 圖12 所示,關鍵時間節點導彈姿態如表2所示。

由圖9 可知,彈體在推力作用下加速離筒。由圖11、圖12 結合表2 可以看出,隨著彈體入水體積的增加,彈體受到的水動力載荷會不斷增加,導致姿態變化,隨著出筒過程中適配器的脫離,彈體運動更加自由,姿態變化會進一步加劇;出水后,隨著出水體積的增加,彈體的姿態變化會變緩。

表2 關鍵時間節點的導彈姿態Tab.2 Missile attitude at critical time

圖9 全過程仿真彈體質心速度-時間歷程Fig.9 Velocity-time curve of the center of mass of missile (whole trajectory)

圖10 全過程仿真彈體質心位移-時間歷程Fig.10 Displacement-time curve of the center of mass of missile(whole trajectory)

圖11 全過程仿真彈體質心俯仰角速度-時間歷程Fig.11 Pitch angular velocity-time curve of the center of mass of missile(whole trajectory)

圖12 全過程仿真彈體質心俯仰角-時間歷程Fig.12 Pitch angle-time curve of the center of mass of missile(whole trajectory)

將仿真數據與實驗結果進行對比可知,實驗測得的數據只有的前2.26,包括入水段和水彈道的一部分,截取柔體仿真數據的前2.26,并附上剛體仿真數據。對比結果如圖13、圖14 及表3 所示。

表3 出筒時刻彈體運動信息對比Tab.3 Comparision of attitudes of missiles out of cainster

圖13 俯仰角速度-時間歷程對比Fig.13 Simulated and experimental pitch angular velocity-time curves

由圖13、圖14 可以看出:柔體仿真曲線與實驗曲線吻合良好,從而證明了仿真模型的有效性。對比柔體仿真與剛體仿真,柔體仿真平均誤差和最大誤差是剛體仿真的0.2 倍。相比剛體模型,柔體模型的結果與實測數據更加接近。值得一提的是,柔體仿真占用計算資源并不多,本算例用筆記本電腦I5 六核處理器10 min 完成計算,與剛體模型的計算量接近。

圖14 俯仰角-時間歷程對比Fig.14 Simulated and experimental pitch angle-time curves

仿真給出了四圈適配器的受力(見圖15),中適配器受力狀況最為惡劣,受到的最大無量綱載荷為0.078,后適配器的受載方向與前3 個相反,符合實際彈體產生俯仰角后適配器的變形情況。圖16 和圖17 分別為全過程俯仰角、俯仰角速度隨時間變化曲線。

圖15 全過程適配器載荷-時間歷程Fig.15 Curves of load on adapters (whole trajectory)

圖16 全過程俯仰角-時間歷程對比Fig.16 Comparison of pitch angles at the center of mass,head and tail of missile (whole trajectory)

圖17 全過程俯仰角速度-時間歷程對比Fig.17 Comparison of pitch angular velocities at the center of mass,head and tail of missile (whole trajectory)

由圖16、圖17 可以看出,彈體不同位置處的運動情況是有明顯差異的,證明柔體仿真是有必要的。結合圖11、圖13、圖16、圖17 可以看出,柔體仿真相比傳統剛體仿真不僅可以捕捉彈體不同位置的運動狀態,還可以提供豐富的振動特性。彈體柔性是由固定界面模態綜合法計算出的前20 階模態表達的。

為了驗證本文方法的適用性,分別對航速為1.5 kn、2.5 kn、3.5 kn 進行仿真,研究潛艇航行速度對俯仰角和俯仰角速度的影響,仿真結果如圖18、圖19 所示。

由圖18 可以看出,航行速度對彈體出筒姿態有較大的影響,隨著航行速度的增加,彈體的俯仰角速度和俯仰角會顯著增加,但是曲線的變化趨勢基本相同。

圖18 全過程俯仰角-時間歷程對比Fig.18 Comparison of pitch angles in different simulation cases (whole trajectory)

由圖11、圖13、圖19 可知:前2.5 彈體有明顯的振動,這是由適配器的脫離引起載荷劇變所致;后半段彈體已完全出筒,部分出水,隨著水動力載荷的減小,響應幅值減小,表現為俯仰角速度曲線的相對平滑。從圖19 中可以看出,質心的振動不如其他位置振動明顯,可見彈體振型疊加后在質心位置有相互削弱甚至抵消的效果。

圖19 全過程俯仰角速度-時間歷程對比Fig.19 Comparison of pitch angular velocities in different simulation cases(whole trajectory)

3 種航速下適配器的受載極值都發生在適配器中(見表4)。從表4 中不難看出,隨著航行速度的增加,適配器的載荷會顯著增加,但各圈適配器受載時間歷程的趨勢相近,且后適配器的受載方向與前3 圈適配器相反。

表4 各工況適配器載荷極值Tab.4 Max load on adapters in different conditions

3 結論

本文針對潛載導彈彈射發射提出了一種可以同時考慮發射筒、彈體柔性的動力學快速分析方法。將柔性點線高副模型拓展應用于潛射導彈的垂直彈射發射動力學建模,模型中同時考慮彈體和發射筒的柔性,基于Hermite 形狀函數對適配器的抗壓力和抗彎力矩進行等效分配,并編程加載到發射筒中跨越適配點的鄰近兩節點。基于Hermite 形狀函數及彈體所受水動力、附連水慣性力的線性分布,假設對水動力、附連水慣性力進行等效加載,實現了潛射導彈垂直彈射發射動力學分析。針對某潛射彈的動力學仿真結果表明仿真姿態與實驗數據吻合較好,驗證了建模方法的實用可信,且模型簡潔、穩定、高效,為潛載彈的彈射發射建模提供了一種新思路,為拓展柔體之間接觸的簡化建模提供了新途徑。得出如下主要結論:

1) 在水動力、附連水慣性力的單元線性分布假設基礎上,采用Hermite 形狀函數和有限元理論將分布的水動力和附連水慣性力向彈體節點等效簡化的處理方法,具有足夠的精度,滿足使用要求。

2) 將彈體在適配器位置處相對于發射筒的運動位移及姿態變化所計算得到的適配器力和力矩,用Hermite 形狀函數向發射筒跨越適配器點的最鄰近兩節點予以等效簡化,算法可信,結果合理。

3) 彈體不同位置的運動角度差異明顯,表明彈體在運動過程中存在高頻振動,也進一步表明柔體仿真十分必要。

4) 本文仿真方法相比剛體仿真計算結果更接近實驗數據,相對于傳統有限元方法僅需很少的計算資源。

5) 潛艇行速度研究結果表明:潛艇航行速度會顯著影響導彈姿態以及適配器的受載大小,但導彈的姿態歷程趨勢相同,僅影響具體的量值,速度越大,姿態偏差越大。