發射筒易碎前蓋開蓋過程的流體與固體耦合仿真方法

張寶振， 王漢平， 竇建中， 程夢文

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

0 引言

易碎蓋技術能極大提升導彈作戰的響應速度，同時兼具質量輕、免維護、可靠性高等優點，是導彈貯運發射箱研制中的關鍵技術。當發動機產生的沖擊波或高壓燃氣沖擊后蓋時，將會產生擾動波[1]，擾動波向前蓋方向傳播，可用于開啟易碎前蓋，從而為導彈的發射出筒打開通道；文獻[2-3]利用計算流體力學(CFD)仿真和實驗手段對易碎蓋開蓋過程進行研究，驗證了上述開蓋機制；黨海燕等[4]探討了易碎后蓋開蓋壓力對易碎前蓋的影響，仿真表明，易碎前蓋的開蓋壓力設計應依據易碎后蓋的開蓋壓力進行合理匹配；牛鈺森等[5]利用動網格技術，將后蓋裂片運動過程與流場耦合，對整個開蓋過程進行仿真分析，結果表明，相比于非耦合工況，流體與固體耦合工況的計算結果與實驗數據更接近；靖建全等[6]采用動網格技術模擬后蓋脫落運動過程，并利用分區邊界的類型變化來模擬后蓋破裂過程，結果顯示，后蓋運動和破碎過程對發射箱內擾動波的產生和傳遞有重要影響；文獻[7-8]利用ABAQUS軟件中耦合的歐拉- 拉格朗日(CEL)算法對地面氣動載荷作用下的大型整流罩分離過程進行模擬，發現基于CEL算法的仿真結果與實驗數據吻合較好；王懿等[9]基于CEL算法建立了拋錨入泥深度的數值模型，仿真研究了土壤不排水抗剪強度、拋錨速度、錨尺寸等參數對船舶拋錨入泥深度的影響；徐文杰[10]則利用CEL算法對滑坡涌浪這種復雜的流體與固體耦合問題進行仿真分析，且仿真結果與實驗數據具有較好的吻合度。由此可見，CEL算法在解決低速流動條件下流體與固體耦合問題具有廣泛的適用性。

國內對易碎蓋開裂過程的研究主要是基于結構的斷裂力學[11-12]，而未考慮流場與易碎蓋碎片運動過程的耦合效應；對易碎蓋開蓋過程中流體與固體耦合問題的研究，則多采用CFD中的動網格技術，卻未細致引入易碎蓋結構的開裂過程，且研究對象也均為易碎后蓋。有鑒于此，本文構建了一種針對易碎前蓋與箱內、箱外流場相耦合的仿真模型，將CFD流場仿真壓強作為輸入條件，利用擴展開發的易碎蓋材料混凝土塑性損傷(CDP)本構模型，對易碎前蓋開蓋過程的流體與固體耦合問題進行了仿真研究。這種將CEL算法應用于耦合面變拓撲的易碎前蓋開裂仿真尚屬先例。仿真結果表明：擴展開發的CDP本構模型可用于易碎蓋的結構設計和強度優化；CEL算法能有效解決易碎前蓋的擾動波開蓋過程的流體與固體耦合問題，且能細致捕捉易碎蓋碎片上的載荷變化特性，可為精準控制易碎蓋的沖破力[13]提供數據支撐。

1 CEL算法

CEL算法最早由Noh[14]提出，其兼具拉格朗日方法和純歐拉方法的優勢，能夠有效克服單元畸變和實現高精度的界面捕捉。該算法所使用的結構網格和流體網格各自獨立，兩域的網格可相互重疊，且對兩套節點的對應與否并無特別要求。CEL算法能在流體域和固體域之間進行載荷、位移、速度等參數信息的實時傳遞，不需要復雜繁瑣的人工干預。因此，CEL算法在模擬流體流動、液面晃動、流體與固體耦合等大變形問題時，具有強大的優勢。

CEL算法的控制方程由連續方程、動量守恒方程和能量守恒方程組成，控制方程需在歐拉坐標下進行時間積分。擾動波開蓋涉及發射箱內外的環境，故需附加理想氣體的狀態方程：

p+pa=ρR(T-Tz)，

(1)

式中：p為氣體壓強；pa為環境壓強；ρ為氣體密度；R為氣體常數；T為氣體溫度；Tz為絕對溫度的0 K值。

ABAQUS軟件中的CEL算法采用罰函數耦合算法[15]來實現流體域和結構域之間的耦合。計算過程中，通過追蹤固體節點對流體歐拉材料界面的貫穿深度δ來確定界面貫穿力的大小。如果發生固體節點對流體材料界面的貫穿，界面力就會分布到流體歐拉材料的錨定點上，如圖1所示。界面力的大小F與貫穿深度δ呈正比：

圖1 罰函數耦合算法Fig.1 Penalty function coupling algorithm

F=kiδ，

(2)

式中：ki為第i個單元的罰剛度系數，其值取決于拉格朗日和歐拉材料特性。

2 易碎前蓋材料的本構模型和校驗

2.1 易碎前蓋材料實驗及本構模型的選擇

在易碎前蓋的擾動波開蓋過程中，前蓋承受的是沖擊載荷，材料的應變強化以及應變率效應影響十分明顯，因此，需要測試不同應變率條件下的擬靜態及動態拉伸、壓縮力學性能[16]。擬靜態拉伸、壓縮實驗試件按國家標準GB 9641—1988《硬質泡沫塑料拉伸性能試驗方法》、GB 8813—1988《硬質泡沫塑料壓縮性能試驗方法》予以制備，實驗在電子萬能實驗機完成；動態拉伸、壓縮實驗試件由航天科工集團第10研究院航天天馬公司提供，實驗在霍普金森拉桿和霍普金森桿壓桿實驗裝置上完成。擬靜態、動態實驗應變率數據如表1所示。

表1 實驗類型及數據Tab.1 Test types and test data

限于篇幅，在此僅列出應變率為1.0的擬靜態拉伸、壓縮實驗的應力- 應變曲線，如圖2所示。從圖2中可以發現：易碎前蓋的材料屬于硬脆性材料，單軸拉伸行為表現準脆性，單軸壓縮行為表現先硬化、后軟化，非常類似于 CDP本構模型。因此，可考慮用擴展開發了斷裂失效的CDP本構模型來表征其材料的力學特性，該模型中的參數可根據實驗數據予以擬合。

圖2 擬靜態拉伸和壓縮的工程應力- 應變曲線Fig.2 Engineering stress-strain curves of quasi-static tension and compressed case

2.2 易碎蓋材料本構模型的擴展開發及驗證

由于CDP本構模型無失效單元的刪除功能，為了附加對失效單元的標記及刪除操作，擴展開發了用戶子程序VUSDFLD.

E=(1-di)E0，

(3)

式中：di中的i分別對應下標t和c，分別表示拉伸、壓縮，dt、dc分別為拉伸和壓縮塑性損傷因子，di的取值范圍為0～1，0表示材料未出現損傷，1表示材料強度完全喪失。

由CDP本構模型受拉、受壓的應力- 應變關系和張勁等[17]基于規范混凝土的應力- 應變提出公式，可建立開裂應變與受拉損傷的聯系、非彈性應變與受壓損傷的聯系，計算公式為

(4)

(5)

CDP本構模型在三維多軸狀態下的應力- 應變關系可通過損傷彈性方程表示：

(6)

通過σ計算靜水壓力以標記積分點的拉伸、壓縮狀態，并將仿真計算的累積損傷因子di與材料實驗擬合出的損傷因子進行比較，若計算的損傷因子大于材料失效的損傷因子，則將材料標記為失效，并完成對單元的刪除。

在ABAQUS軟件中分別進行材料擬靜態拉伸、壓縮狀態和動態拉伸狀態的仿真，將仿真結果和實驗結果對比，以此來對開發的CDP本構模型進行校驗，結果如圖3和圖4所示。結果表明：擴展開發的CDP本構模型能夠較好地吻合試件的失效狀態，且靜態拉壓載荷- 變形量數據也吻合較好。其中，拉伸載荷的差異在于材料的拉伸、壓縮彈性模量不一致，即材料的壓縮彈性模量大于拉伸彈性模量，而仿真是按壓縮彈性模量來賦予材料彈性模量屬性的。壓縮的軟化段仿真載荷值小于實驗值，其原因在于：仿真模型中單元失效后會被刪除，且未設置單元失效后斷裂界面與單元表面接觸，因此仿真模型的承載能力有所降低。

圖3 實驗和仿真破壞形態Fig.3 Experimental and simulated failure modes

圖4 單軸拉伸和壓縮載荷- 變形量的實驗數據與仿真結果對比Fig.4 Comparison of simulated results and experimental data of uniaxial tension/compression loading and deformation

3 基于CFD方法的流場仿真

3.1 仿真模型

考慮到發射筒和導彈的結構特征，將流場仿真模型簡化為二維軸對稱模型，計算域由彈體、噴管及發射筒前、后蓋以及后蓋外的部分區域組成。為監測流場狀態參數的變化情況和擾動波的傳播特性，在發射箱內距噴管出口0～5 m處的筒壁、筒與彈之間以及前蓋上設置壓強、溫度監測點。模型計算域網格均采用結構化網格，網格總數為32萬個。仿真模型的計算域及各監測點的分布如圖5所示。

圖5 模型計算域及監測點分布Fig.5 Computational domain and monitoring point distribution

3.2 條件設定

導彈發動機燃氣參數如表2所示。

表2 燃氣參數表Tab.2 Gas parameters

依據總體要求，破片壓強為3 MPa，發射箱內預先充氣10 kPa，故在初始參數設置時，將導彈、發射筒間的氣體壓強設置為10 kPa. 發動機高壓室以外的計算域初始參數為靜止大氣參數：p0=101 325 Pa，T0=300 K. 發動機燃燒室的總溫Tt=3 000 K，總壓為導彈發動機燃燒室的總壓時間歷程，如圖6所示。

圖6 燃燒室總壓時程曲線Fig.6 Curve of total pressure over time

噴管出口邊界條件為壓力入口。后易碎蓋開啟前，其邊界條件為無滑移絕熱壁面；開啟后，其邊界條件為壓力出口。為獲取后蓋在不同時刻開啟時，前蓋在反射擾動波作用下所能達到的載荷狀態，因此將前蓋設定為壁面。

通過設定后蓋的開啟時間來控制后蓋邊界條件的切換，仿真設定的后蓋開啟時間分別為1 ms、2 ms、3 ms、4 ms、5 ms、6 ms和7 ms，共計7種工況。按前述設定參數，對以上7種工況進行仿真，探究發射過程中的反射擾動波、傳播特性，以及前易碎前蓋的受載特性。

4 基于CEL算法的建模及分析

采用CEL算法來處理擾動波開蓋的流體與固體耦合問題，需要將CFD方法得到的壓強時程用作CEL模型的入口條件。而在建立CEL模型之前，需將CFD模型和ABAQUS歐拉模型的箱壁上1號～5號監測點處的壓強時程曲線進行對比，以驗證兩方法計算出的壓力波在幅值、傳播速度上的一致性，這是進行流體與固體耦合仿真的基礎。

4.1 歐拉單元仿真分析及校驗

4.1.1 基于歐拉方法的流場建模及仿真

由于發射筒的結構特性，可選取結構的1/6進行建模。建立發射箱內距離噴管出口1 m處(即1號監測點截面)至前蓋的歐拉模型，并建立與CFD模型位置相匹配的監測點。所有網格均為EC3D8R單元，節點數250 835個，單元數230 232個。歐拉模型的網格模型及監測點分布如圖7所示。

圖7 歐拉模型的網格及監測點分布Fig.7 Grid and monitoring points distribution of Eulerian model

選用理想氣體狀態方程，歐拉材料的定義與CFD方法中定義的氣體屬性一致。距離噴管出口1 m截面處設置為壓強入口，壓強幅值曲線按CFD仿真提供的監測點1處的壓強時程曲線予以設置；除距離噴管出口1 m截面處，其余外部邊界均設置為絕熱無滑移邊界條件。由于歐拉仿真驗證不需拉格朗日單元，所以初始條件下，歐拉體積分數全部為1.

選擇ABAQUS顯示動力學求解器，仿真時間為0.025 s. 按上述設定參數，對7種工況進行仿真，獲取各監測點處的壓強時程曲線。

4.1.2 歐拉方法與CFD方法的仿真對比

根據CFD方法和歐拉方法的仿真計算，得到各監測點處的壓強時程曲線。限于篇幅，僅對后蓋開啟時間為6 ms這一工況，將兩種方法中筒壁上1號～5號監測點、前蓋3號監測點的壓強時程曲線對比分析，分別如圖8(a)和圖8(b)所示。

圖8 筒壁上各監測點和前蓋3號監測點的壓強時程曲線(CFD方法與歐拉方法對比)Fig.8 Pressures at monitoring points on the canister wall and monitoring point 3 on the front cover (CFD method vs. Eulerian method)

從圖8中可以看出：歐拉方法與CFD方法得到的箱壁、前蓋3號監測點的壓強時程曲線基本重合，可見兩方法得到的壓力擾動波的幅值、傳播速度、反射特性吻合較好，這是進行壓力波開蓋流體與固體耦合分析的有力支撐；各監測點均出現壓強二次峰值，這是因前蓋對壓力波的反射而形成。圖8中同一監測點處的壓力波在傳播速度、壓強幅值方面存在微小差別，這是由歐拉方法和CFD方法確立的監測點位置不完全重合導致的，即CFD方法是根據距離設置的監測點，而歐拉方法的監測點位置是根據網格節點予以定義的。

本設計是基于單一景點深度開發的專門性的游客服務的手機APP。本設計專注于單一景點開發，因此景點的信息覆蓋全面，沒有冗雜信息干擾，目的性明確，垂直服務于景點旅游的游客。本設計的意義在于旅游景區為游客提供更好的服務，為用戶提供更好的指引。同時游客能夠通過APP隨時隨地了解某一景區的最新動態，有利于自己選擇合適的時間出游，并且能夠快速地滿足酒店、旅游社以及旅客三方面的需求[3]。

4.2 易碎前蓋的流體與固體耦合仿真

4.2.1 CEL模型的建立

考慮到發射筒和易碎蓋的結構特征，可用其1/6循環對稱的結構進行建模計算。易碎前蓋的1/6結構示意圖如圖9所示。

圖9 易碎前蓋的1/6結構示意圖Fig.9 1/6 schematic diagram of fragile front cover

CEL模型建模與歐拉模型基本一致，不同之處為：前蓋網格為拉格朗日單元，需將前蓋和發射筒的歐拉域進行裝配。邊界條件除絕熱無滑移和對稱邊界外，需設置歐拉邊界：自由流入和無反射流出；還需對拉格朗日單元附近的歐拉單元進行局部加密，以保證耦合區域的網格密度，同時避免計算量的過分增加。

CEL模型中易碎前蓋的拉格朗日網格(C3D8R)為61 952個，流場域的歐拉網格(EC3D8R)為334 635個，網格模型如圖10所示。

圖10 CEL模型的網格及監測點分布Fig.10 Grid and monitoring points distribution of CEL model

4.2.2 CEL模型仿真結果分析

經過仿真計算和數據處理，得到CEL模型監測點處的壓強時程曲線。限于篇幅，僅對后蓋開啟時間為6 ms這一工況，將CEL算法與CFD方法的結果進行對比分析，如圖11所示。

圖11 箱壁上各監測點和前蓋監測點3的壓強時程對比(CFD方法與CEL算法對比)Fig.11 Pressure curves of monitoring points on the canister wall and No.3 point on the front cover (CFD method vs. CEL algorithm)

CEL模型中易碎前蓋內表面接觸壓力云圖如圖12所示，其中圖12(a)～圖12(j)為選取的特征較為明顯的壓力云圖，對應的時刻依次為：12.875 ms、13.0 ms、13.5 ms、15.0 ms、15.5 ms、16.0 ms、16.5 ms、17.0 ms、17.5 ms、18.0 ms. 從圖12中可以看出：12.875 ms是反射壓力波傳遞到前蓋的第一個輸出步，此時，易碎前蓋內表面并無壓力變化；在13.0 ms時，壓力波已傳遞到前蓋，但只是先傳到前蓋的法蘭處；再經過4個計算步，在13.5 ms時，壓力波完全作用到蓋上，壓力波在不同時刻到達前蓋各個位置，并且向更靠近中軸線的方向匯聚，所以在蓋的中間處會有明顯的壓強增大，這有利于前蓋的開啟；壓力波持續一段時間后，在15.5 ms時，前蓋完全沿預制溝槽開啟，且開裂效果良好。如圖13所示，易碎前蓋已完全碎裂并飛離發射筒。之后，易碎蓋碎片在壓力波沖擊下進一步運動。

圖12 易碎前蓋開裂過程的內表面壓力云圖Fig.12 Pressure nephogram of inner surface of fragile front cover

圖13 前易碎蓋完全破裂Fig.13 Completely broken fragile front cover

5 結論

本文構建了考慮易碎前蓋與壓力波交互作用的CEL流體與固體耦合模型，利用擴展開發的易碎前蓋材料的CDP本構模型，在輔以CFD仿真數據作為壓力邊界的條件下，對易碎前蓋的壓力波開蓋過程進行了流體與固體耦合仿真分析。根據仿真分析，可得出以下主要結論：

1) 擴展開發后的CDP本構模型能有效模擬材料的失效以及失效單元的刪除，也能很好復現材料實驗和產品實驗的結果，適合于易碎前蓋的本構模擬應用。

2) 以CFD仿真結果作為歐拉模型的邊界條件，對壓力波的傳播、反射特性進行了對比仿真驗證。結果表明：CFD方法、歐拉方法的壓力波強度、傳播速度及其反射特性吻合較好，歐拉方法能較好地模擬壓力波的傳播及反射特性，這說明基于有限元的歐拉方法在模擬壓力波的傳播、反射方面具有匹配CFD的模擬精度。

3) 將CFD仿真結果作為CEL流體與固體耦合模型的邊界條件，開展易碎前蓋的開裂過程仿真，在壓力波的自由傳播、反射特性及其與易碎前蓋接觸耦合方面，CEL模型更符合實際狀態。采用CEL算法能夠較好地解決耦合面變拓撲的易碎前蓋開裂的仿真問題，這種將CFD與CEL相結合的仿真方法為易碎前蓋的壓力波開蓋流體與固體耦合分析提供了一種可信的新方法，為工程中易碎蓋的設計和優化提供更優良的解決方案。