推進劑“脫濕”損傷研究的內聚力單元方法

2022-06-08 09:10:58崔輝如許玉榮

國防科技大學學報 2022年3期

崔輝如，呂軒，許玉榮

(1. 陸軍工程大學國防工程學院, 江蘇南京 210007； 2. 湖北航天技術研究院總體設計所，湖北武漢 430040)

近年來，細觀力學的發展為研究推進劑非線性力學性能提供了新的研究手段。在此基礎上，可以結合數值損傷模型，開展推進劑細觀力學性能研究，揭示推進劑“脫濕”損傷機理，對推進劑非線性力學性能的預示具有較大的推動意義[1-7]。

通過顆粒填充手段模擬生成的推進劑細觀結構為推進劑松弛模量的預測[4, 8]、推進劑細觀損傷的數值模擬研究[5]提供了便利條件。推進劑細觀結構中，顆粒和基體之間黏接界面的破壞是“脫濕”產生的主要原因。Zhang[9]以及Han[10]等建立了固體推進劑的細觀力學模型，對大變形下的推進劑內部界面“脫濕”現象進行模擬，采用內聚單元對黏接界面進行建模。為了模擬不同固體含量的推進劑中細觀損傷的發展歷程，職世君等[6]通過分子動力學的方法得到了推進劑的顆粒裝填模型，其基于內聚接觸手段結合雙線性內聚力本構關系，建立了AP顆粒和基體之間的損傷模型；基于有限元方法，對比分析了不同固體含量對細觀損傷形態的影響并指出高體分比狀態下顆粒排列非常致密，界面損傷的相互影響較大是加快固體推進劑產生損傷至形成宏觀裂紋的過程的主要因素；此外，研究還表明，內聚強度越小，復合固體推進劑的損傷起始時間越早，損傷程度越大，相應的拉伸應力越小。韓龍等[2]開展了復合固體推進劑細觀界面率相關參數的反演分析并揭示了黏接界面性能隨加載速率的變化規律。封濤等[1]采用雙線性內聚力模型來描述顆粒與基體黏接單元的力學響應，深入研究了細觀結構對端羥基聚丁二烯(hydroxyl terminated polybutadiene)推進劑力學性能的影響。周紅梅等[11]同樣采用了內聚力單元的方式，開展了推進劑細觀結構在常溫以及低溫下的損傷機理研究，研究表明低溫下的推進劑細觀結構破壞形式更為復雜。此外，還有研究[12]指出，在細觀結構下，推進劑一般首先在極區產生脫黏，顆粒的大小對推進劑的力學性能有較大影響。內聚力單元方法是推進劑細觀“脫濕”損傷研究的重要途徑之一[13]。在推進劑細觀幾何的處理、內聚力單元模型的構建以及邊界條件的設計方面，上述文獻中的處理方式還不是特別全面和完善。通過內聚力單元方式進行推進劑“脫濕”損傷機理的研究還有待進一步深入。

為了進一步研究推進劑細觀“脫濕”損傷機理，本文在合理的代表體積單元尺寸的基礎上，采用內聚力單元方法對顆粒和基體之間的黏接界面進行建模。結合周期幾何和周期邊界條件，采用單軸拉伸和剪切試驗的手段，開展推進劑細觀“脫濕”損傷的機理研究，并對不同推進劑體分比、應變速率以及內聚強度進行“脫濕”損傷力學行為影響分析。

1 推進劑細觀有限元模型

1.1 細觀幾何模型的構建

推進劑是典型的高填充型復合材料，顆粒的體積分數往往會達到60%，甚至更高。目前分子動力學方法是行之有效的方法，Kochevets[14]、申柳雷[15]的博士學位論文均給出了利用分子動力學進行推進劑顆粒裝填模型建立的具體細節。

為了描述一個真實的細觀結構的力學響應，細觀模型必須要滿足周期幾何以及周期邊界條件的要求。在進行顆粒裝填時，如果發生顆粒的“壓線”事件，那么對應幾何邊(或者三維模型中的面)上就要有相應的顆粒，才能滿足周期幾何的限制條件。此外，在進行顆粒裝填之前，顆粒的數量以及相應的幾何尺寸是按照材料的配方計算出來的固定的值，并不能因為裝填的過程而改變。為了解決上述的兩個問題，以模型的中心構建直角坐標系，第i個顆粒的坐標(顆粒邊界頂點的坐標)為ηi(xi,yi,zi)，其中無量綱的坐標范圍為[-1,1]。上述問題可以轉化為

(1)

通過這樣的手段，可以使得所有“壓線”的顆粒都在模型的對應幾何邊界上有對稱的顆粒，盡管顆粒的數量比初始計算的數量增加，但這并不影響總的裝填分數，如圖1所示。

圖1 顆粒發生“壓線”事件以及周期幾何的處理Fig.1 "Line pressing" event of particles and periodic geometry processing

1.2 推進劑細觀“脫濕”行為仿真

推進劑細觀“脫濕”損傷仿真的核心是顆粒和基體界面分離行為的模擬。應用比較成功的主要有內聚接觸模型和內聚力單元模型兩種。內聚接觸是ABAQUS自帶的功能，其本質是對接觸行為的模擬。當然，這種接觸是一種廣義的“接觸”，既可以描述界面之間的相互擠壓行為，又可以描述界面之間的相互分離行為。內聚力單元會將單元上分布的內聚力等效到單元的節點上，代入總體剛度陣的計算。然而，內聚接觸是分布在界面上的一連串的離散的彈簧，被黏接的物體是通過一系列離散的彈簧相連接，至于沒有彈簧連接的地方，是不存在任何相互作用的；內聚力單元則是界面上分布的連續的彈簧，被黏接物之間沒有任何的間隙，黏接界面之間任何微小的分離都將被捕捉，并被代入單元節點力的計算中。正因為兩種不同的建模理念，內聚接觸方法會比內聚力單元方法更容易收斂，但是精度會大大降低，甚至會導致錯誤的結果，畢竟分段式的連接并不符合黏接的實際情況。因此，本文將采用更符合實際的內聚力單元進行界面損傷的建模。圖2(b)給出了圖2(a)中顆粒的有限元模型以及界面上分布的零厚度內聚力單元，可以發現，界面上分布的內聚力單元和顆粒單元的邊界是完全吻合的。

圖2 顆粒網格以及零厚度內聚力單元分布Fig.2 Distribution of mesh of particles and zero thickness cohesive elements

PPR(Park-Paulino-Roesler)內聚力模型是考慮物理斷裂參數以及連續斷裂邊界的勢相關內聚力模型[16]，模型的勢函數表達式為

(2)

其中：Δn和Δt分別為法向和切向分離位移；Γn和Γt為能量常數；m和n為無量綱參數；〈·〉為MacAuley算子；αk和βk為模型的形狀參數，分別控制法向和切向的內聚力和分離位移曲線下降段的凹凸特性；δn和δt為模型的法向和切向的臨界位移；φt和φn分別代表切向和法向內聚能，其表達式為

(3)

本文中，黏接界面初始內聚強度σmax和τmax為0.2 MPa、初始臨界位移δn和δt為50 μm、初始無量綱參數m和n為0.01、形狀參數αk和βk為2.10。此外，顆粒采用線彈性材料屬性，基體采用線黏彈性材料屬性。

1.3 單元平均應力應變的后處理方法

商業有限元軟件并不能直接給出模型的等效應力應變關系。為了便于進一步分析模型的力學響應特性，需要對計算結果進行處理。對于三角形網格模型，其平均應力應變關系有如下表達式

(4)

(5)

考慮到網格數量達到一定的規模時，網格對計算結果的影響差異不大，進行有限元的分析過程中，節點的位移是計算出來的直接的結果。但是考慮到部分節點既屬于基體又屬于顆粒夾雜，節點處應力的處理方法并不統一。因此，本文直接采用積分點出的應力和應變結果代入式(4)～(5)進行計算，即

(6)

(7)

2 推進劑“脫濕”損傷機理

2.1 單軸拉伸下的推進劑細觀“脫濕”損傷分析

從單軸拉伸試驗出發，著重分析考慮和不考慮損傷狀態下，推進劑拉伸階段的應力應變曲線。“脫濕”現象和應變速率具有較大的關聯性，因為基體和黏接界面的力學性能對應變速率較為敏感。本文暫時不考慮基體和黏接界面的率相關力學性能。在細觀推進劑結構拉伸過程中，結構整體的應變速率為0.2%/s。

圖3(a)首先給出了結構在定速拉伸階段，縱向應力應變曲線隨時間的關系。由圖可以發現，顆粒和基體黏接界面的損傷效應會明顯地降低推進劑的剛度和承載能力。圖3(b)為拉伸過程中縱向“脫濕”因子的變化曲線。“脫濕”因子具體的定義為

(8)

其中，D代表“脫濕”因子，σwithout和σwith分別表示不考慮損傷和考慮損傷狀態下的應力。

在加載開始階段，結構的“脫濕”因子在2 s的時間內就增加到了0.11左右。這說明，在加載的開始階段，界面的損傷效應就已經存在。值得注意的是，在“脫濕”因子隨加載時間穩定增長之前，曲線存在一個“凹陷”。本文認為“凹陷”的存在是局部黏接界面處界面剛度符號的改變導致的：在加載的初始階段，黏接界面處的分離位移都比較小，界面剛度基本處于上升階段。隨著時間的推移，局部黏接界面處的分離位移加速增長，界面剛度開始出現下降，這一局部剛度的突然改變導致了“脫濕”因子曲線上“凹陷”的存在。局部黏接界面處的損傷效應向外傳遞時，大部分界面處的剛度還處于增長階段，局部的負剛度效應會被整體的剛度抵消，因此，在很短的時間內，“脫濕”因子曲線繼續隨時間不斷增長。

(a) 應力應變隨時間分布(a) Stress and strain versus time

(b) “脫濕”因子隨時間分布(b) "Dewetting" factor versus time圖3 定速拉伸階段縱向應力應變和“脫濕”因子變化曲線Fig.3 Longitudinal stress/strain and "dewetting" factor versus time in stretch period with constant speed

圖4分別給出了四種典型時刻推進劑細觀結構的變形。在t=20 s時刻，整個結構并沒有發生相對于結構尺寸而言較為明顯的裂紋。結合“脫濕”因子曲線，可以判斷這一時刻，結構內部細微的裂紋早就已經讓結構的承載能力產生了較大幅度的削減。在t=40 s時刻，結構內部已經開始出現明顯的裂紋。如圖4(b)所示，在圖中的實線框，裂紋呈現“成對”出現的現象，排列順序和拉伸方向一致。也就是說，某個界面處的裂紋并不是“孤立”存在的，在裂紋的附近，沿著結構受力的方向，總能找到另一個明顯的裂紋。圖4(c)中，t=60 s時刻，裂紋“成對”出現的現象更加明顯。這樣的“裂紋對”可以在兩個顆粒與基體之間的黏接界面發生，也可以在三個顆粒與基體之間的黏接界面發生。此外，這樣的“裂紋對”的產生還需要界面與界面之間有足夠近的距離，圖中粗虛線框的部分，由于在界面附近沒有距離足夠近的另一個界面，“裂紋對”現象并沒有出現。圖4(d)中，t=70 s時刻，“裂紋對”已經布滿整個結構表面了。部分“裂紋對”的參與界面甚至從兩個發展到了三個。此時，黏接界面處的裂紋已經有擴展到基體中的趨勢。

(a) t=20 s (b) t=40 s

2.2 純剪試驗下的推進劑細觀“脫濕”損傷分析

本小節將進一步利用純剪試驗對推進劑細觀結構的“脫濕”機理進行研究。圖5分別給出了兩種剪切方向(模式A和模式B)下的變形示意，切應變的應變速率均為0.02%/s。可以明顯地看到，“裂紋對”在沿著剪切的方向(雙箭頭線)上有序分布。在垂直于剪切方向的位置上，幾乎看不到明顯的裂紋。此外，部分界面由于在沿著剪切方向上的黏接界面數目較少，距離較遠，同樣沒有“裂紋對”的出現。

(a) 模式A(a) Model A (b) 模式B(b) Model B圖5 t=40 s時刻模式A和B的變形Fig.5 Diagram of deformed model A and B at t=40 s

圖6(a)所示為剪切階段切應力和應變隨時間的變化。以沿著二四象限的剪切模式為例，在考慮和不考慮界面損傷狀態下，結構在xy以及yx方向上的應力應變曲線是一致的。這就從細觀層面上驗證了切應力互等原理，進一步說明了細觀結構和仿真手段的準確性。另外，考慮界面損傷效應的應力應變曲線明顯低于無損傷狀態下的應力應變曲線。這說明界面結構的損傷會削弱推進劑的抗剪能力。圖6(b)給出了兩種加載方向下的“脫濕”因子變化曲線，同一模式(模式A或模式B)下不同方向上的“脫濕”因子變化曲線幾乎是重合的。盡管不同模式下(模式A和模式B)的“脫濕”因子變化曲線不一致，但是誤差保持在1%左右。此外，類似于單軸拉伸，“凹陷”現象同樣存在，這里不再重復解釋。

(a) 應力應變隨時間分布(a) Stress and strain versus time

(b) “脫濕”因子隨時間分布(b) "Dewetting" factor versus time圖6 純剪作用下推進劑的應力應變和“脫濕”因子變化曲線Fig.6 Stress/strain and "dewetting" factor versus time under pure shear effect

3 考慮“脫濕”損傷的推進劑力學行為分析

3.1 推進劑體分比對“脫濕”損傷的影響規律

推進劑中顆粒的體積與推進劑總體積的比值稱為體分比。按照同一級配(粒徑分布)生成四種不同體分比的推進劑細觀有限元模型，對應的體分比Vf分別為60%、65.3%、70%以及75%。采用同樣的單軸拉伸試驗，縱向應變速率固定為0.2%/s。圖7分別給出了這些有限元模型在考慮和不考慮界面“脫濕”損傷效應下的應力應變曲線。可以發現，隨著體分比的增加，結構中固體填料增多，整體剛度變大，應力應變曲線整體水平獲得了提升。但是，隨著體分比的增大，考慮界面“脫濕”損傷效應下的應力應變水平比不考慮損傷下的應力應變曲線之間差異愈發明顯。

(a) Vf=60%

(b) Vf=65.3%

(d) Vf=75%圖7 不同體分比下的應力應變隨時間分布Fig.7 Stress/strain versus time under different volume factors

圖8 不同體分比下的推進劑“脫濕”因子隨時間變化曲線Fig.8 "Dewetting" factor versus time under different volume factors

圖8分別給出了這些不同體分比下的推進劑細觀結構“脫濕”因子隨時間的變化關系。可以看到，“凹陷”現象在初始階段已經變成異常的上下波動，但是這樣的波動之后，“脫濕”因子曲線會繼續平穩發展。不同體分比下的“脫濕”因子曲線進入穩定發展階段的時刻并沒有特別明顯的差異。此外，盡管體分比越高，“脫濕”因子曲線提升現象越明顯，但是曲線的斜率之間并沒有較大的差異。

3.2 應變速率對“脫濕”損傷的影響規律

推進劑在不同的應變速率下往往會表現出不同的力學性能。圖9分別給出了體分比為65.3%、內聚強度為0.2 MPa時，不同應變速率v下的推進劑細觀結構應力應變曲線。可以看到，隨著應變速率的提升，應力應變曲線不斷提升，相應的“脫濕”損傷效應愈明顯。

圖10給出了不同應變速率下的“脫濕”損傷因子變化曲線。這里，高應變速率下的“脫濕”因子曲線明顯大于低應變速率下的“脫濕”因子曲

(a) v=0.1%/s

(b) v=0.2%/s

(d) v=0.4%/s圖9 不同應變速率下的應力應變曲線Fig.9 Stress/strain versus time under different tensile rates

線。此外，高應變速率下的“脫濕”因子曲線斜率明顯大于低應變率下的“脫濕”因子曲線。這是因為，應變速率提升的狀態下，結構中基體和顆粒間的“脫濕”來不及“愈合”，那么相應的損傷效應就會加劇。

圖10 應變速率對推進劑“脫濕”因子的影響Fig.10 Effects of tensile rates on "dewetting" factor

3.3 內聚強度對“脫濕”損傷的影響規律

(a) σmax=τmax=0.1 MPa

(b) σmax=τmax=0.2 MPa

(d) σmax=τmax=0.4 MPa圖11 不同內聚強度下應力應變隨時間分布Fig.11 Stress/strain versus time under different cohesive strength

為了探究不同的內聚強度對應力應變曲線的影響，圖11給出了不同內聚強度下推進劑細觀結構在0.2%/s的應變速率下的應力應變曲線。在考慮界面“脫濕”損傷狀態時，內聚強度越小，細觀結構的“脫濕”損傷效應越明顯。如果界面的內聚強度接近無窮大，結構的計算結果應該和不考慮“脫濕”損傷的狀態一致。這一經驗理論和仿真預示是吻合的。

圖12給出了拉伸階段，不同內聚強度下的“脫濕”因子變化曲線。可以看到，當內聚強度達到0.4 MPa時，“脫濕”因子曲線在“凹陷”點之后一直呈現下降的趨勢。在PPR內聚關系中，內聚強度的提升會使得內聚剛度增大，因此在內聚強度很高時，隨著分離位移的提高，結構整體的剛度加大，界面的“脫濕”損傷效應就會得到抑制。因此，在內聚強度為0.3 MPa和0.4 MPa時，結構整體的“脫濕”因子曲線會出現持續下降的趨勢。也就是說，在考慮界面“脫濕”損傷狀態時，應變速率越大、內聚強度越小，細觀結構的“脫濕”損傷效應越明顯。

圖12 內聚強度對推進劑“脫濕”因子的影響Fig.12 Effects of cohesive strength on “dewetting” factor

4 結論

從推進劑細觀分析模型出發，基于分子動力學的推進劑顆粒裝填方法、周期幾何以及周期邊界的處理手段，在單軸拉伸和純剪試驗下，進行了基于內聚力模型的推進劑“脫濕”損傷機理研究，開展了推進劑體分比、應變速率、內聚強度參數對推進劑細觀結構仿真結果的影響分析。主要結論如下：

1)在外載荷作用下，沿著細觀結構的受力方向，“裂紋對”現象的出現是推進劑發生“脫濕”損傷的內因。

2)隨著體分比的增大，推進劑細觀結構“裂紋對”現象增多，“脫濕”因子曲線提升現象明顯。體分比由60%提升至75%的過程中，“脫濕”因子增加近一倍。

3)應變速率的減小和內聚強度的提升會使得“脫濕”因子減小，相應地削弱推進劑細觀結構的“脫濕”損傷效應。在內聚強度為0.3 MPa和0.4 MPa時，結構整體的“脫濕”因子曲線會出現持續下降的趨勢。