粘彈性應變能釋放率的虛擬網格計算方法①

崔輝如,姜強強,沈 棟,呂 軒,任孝宇

(1.陸軍工程大學 國防工程學院,南京 210007;2.陸軍工程大學 土木工程博士后科研流動站,南京 210007;3.中國人民解放軍91112部隊,寧波 315046;4.航天動力先進技術湖北省重點實驗室,武漢 430040)

0 引言

當前,粘彈性材料已經被廣泛地用于各種工程結構中。伴隨著結構尺寸以及載荷環境的日益極端化,粘彈性結構的斷裂問題顯得尤為重要。以固體火箭發動機為例,其藥柱屬于典型的粘彈性材料。在日益頻繁的戰備值班以及實彈演練過程中,藥柱結構內部很容易出現微裂紋,進而直接或間接導致藥柱結構完整性的失效[1-2]。不同于彈性材料,粘彈性斷裂問題尤其是粘彈性斷裂參數的數值計算方法還有待深入發掘和研究。

斷裂力學參數的計算是裂紋擴展判斷的重要依據。應力強度因子以及應變能釋放率等斷裂力學參數的計算大部分依賴于原始的網格或者構建虛擬的區域進行求解。在原始網格的基礎上,最早最直接的方法是基于裂紋尖端前沿的單元應力或裂紋尖端后的節點位移的外推法[3]。王珊[4]利用有限元節點位移確定含裂紋Kirchhoff板辛本征解中的待定系數,計算出了裂紋尖端附近應力場的顯式表達式。段靜波等[5]將加料有限元法擴展應用于粘彈性介質中裂紋應變能釋放率的計算。祁勇峰等[6]在裂紋尖端引入Williams解析解級數,應用解析覆蓋與周邊數值覆蓋聯合求解三維穿透直裂紋的應力強度因子。RAJARAM等[7]采用積分點處的攝動梯度與直接解析微分之間的聯系,DAIMON等[8]引入一個使得輔助解滿足積分平衡關系式的修正項來進行強度因子的求解。JIN等[9]利用彈性-粘彈性對應原理研究了粘彈性梯度材料Ⅰ型及Ⅱ型裂紋的應力強度因子解析解。另外,還有研究采用虛擬網格來進行斷裂參數的求解。例如,采用沿裂紋前緣的虛擬圓柱單元,利用極坐標下的高斯點在圓域中進行J積分的計算[10-11]。為避免網格中應力強度因子沿裂紋前沿的振蕩,可以利用被積函數的移動最小二乘近似[12]。類似地,還有研究采用盤形域積分法逐點計算J積分和強度因子[17]。GROSSMAN-PONEMON等[14]在虛擬網格中展示了一種基于獨立網格的多網格交互作用積分方法。CHOI等[15]利用虛擬網格插值計算虛擬網格裂紋尖端位移,通過計算虛擬網格J積分,進而獲取了低質量網格下二維以及三維彈性材料應力強度因子。

裂紋尖端斷裂力學參數計算已經有不少的實現辦法和實例支撐,但是目前還有以下不足:一是研究對象大多針對線彈性材料,對粘彈性介質斷裂問題缺乏深入研究;二是伴隨著裂紋擴展過程中非結構且低質量網格的使用,上述方法實現流程相對復雜,并且難以確保仿真精度。本文在前期研究的基礎上[15],在裂紋尖端構建了一種新的虛擬網格形式,結合虛擬裂紋閉合方法,發展了一種針對粘彈性材料裂紋尖端應變能釋放率的高精度數值計算方法。利用兩種典型裂紋案例與文獻驗證了本文方法的有效性。

1 線粘彈性本構模型

在泊松比為常數的條件下,推進劑三維線粘彈性本構關系表達式為[16]

(1)

其中,

(2)

(3)

式中ν為泊松比;σij、Sij及σkk分別為應力張量、偏應力張量及球應力張量;eij及εkk分別為偏應變張量以及球應變張量;t為加載時間;E(t)為松弛模量。

E(t)Prony級數表達式為

(4)

2 粘彈性裂紋尖端虛擬網格

如圖1,以原始網格裂紋尖端P點為中心構建虛擬網格。虛擬網格由2×2個正方形單元構成(D2PKJ、PHNK、EFPD1以及FGHP),虛擬單元的邊長為原始網格裂紋尖端平均單元邊長的2倍[15]。與此同時,虛擬網格需要在原始裂紋路徑上設置同樣的裂縫,因此虛擬網格一共有10個虛擬節點。

圖1 原始網格(黑色三角形網格)與虛擬網格(藍色四邊形網格)示意圖

虛擬節點的位移計算遵循以下原則。假設,原始網格某單元的組成節點為A、B和C,那么虛擬節點N剛好位于該單元區域的判斷依據是

SABN+SBCN+SCAN=SABC

(5)

式中SABN、SBCN、,SCAN、SABC分別為三角形ABN、BCN、CAN及ABC的面積。

假設原始網格上單元A、B和C三個節點位移分別是(uA,vA),(uB,vB)以及(uC,vC)。N點、A點、B點以及C點的坐標為分別為(xN,yN)、(xA,yA)、(xB,yB)以及(xC,yC)。

利用面積法進行虛擬節點N位移(uN,vN)的插值計算

(6)

其中,

(7)

以此類推,不難計算出任意時間步上虛擬網格的位移場。

考慮到粘彈性材料具有典型的時間相關特性,為此在進行虛擬網格應力場計算的時候,不僅需要代入上一步計算中的應力、應變增量等數據還要將當前步產生的應力、應變增量等數據傳遞到下一計算步中。有關粘彈性本構模型的數值化計算方法可以參考相關文獻,這里不再詳細展開[17-18]。

虛擬網格應力計算以及原始網格的應力計算都需要進行相關狀態變量的數據傳遞,上述操作會給主程序的設計和編程帶來很大的不便。考慮虛擬網格中節點和單元總數較少,可以在每一步計算中利用相同的位移插值計算方案,計算出虛擬網格在當前步以及前序時間步中的位移場,進而直接計算出當前載荷步上虛擬網格的應力以及應變場。上述操作可以有效降低時間步之間狀態變量傳遞的工作量,簡化編程任務。為了方便表示,將該虛擬網格應力重構的方法記作(Stress Recovery Technique,SRT)。圖2給出了本文中粘彈性應變能釋放率數值計算的流程示意。

圖2 應變能釋放率數值計算流程示意圖

本文中的虛單元均采用常規一階四節點四邊形單元。虛單元的單元積分方案有兩種,分別是2×2高斯積分點的完全積分以及1×1高斯積分點的縮減積分。如無特殊交代,虛單元采用完全積分方案。此外,本文中的原始網格均采用三節點三角形單元,積分點數為1。

3 虛擬裂紋閉合方法

Irwin在假設“勢能改變與裂紋閉合一個擴展增量所需的功等效”基礎上,利用裂紋閉合積分計算裂紋尖端的能量釋放率[3]:

(8)

式中(見圖3)a和Δa分別為裂紋長度和裂紋擴展增量;σyy和τxy分別為沿著閉合裂紋面上的法向和切向應力;T為裂紋體厚度;Δv和Δu分別為當閉合裂紋張開時閉合裂紋面上的法向和切向位移。

圖3 閉合的虛擬裂紋示意圖

上式在實際的計算過程中需要開展兩次分析計算,這給裂紋擴展問題的研究帶來很大的不便。為此,RYBICKI等[19]首次提出了應變能釋放率計算的虛擬裂紋閉合方法。該方法假設虛擬裂紋尖端后面的張開位移和實際裂紋尖端后面的張開位移近似相等。此外,為了避免積分號中出現不精確的應力值,采用節點力代替應力的積分,即

(9)

式中(見圖4)Fx1和Fy1分別為節點1處內力沿著水平和垂直方向的分力;Δv2,3和Δu2,3分別為節點2和3之間垂直和水平方向的位移差。

圖4 虛擬裂紋閉合方法計算應變能釋放率示意圖

在上節虛擬網格的基礎上,利用編程很容易實現節點1處支反力以及節點2和3處位移的計算,進而計算出任意時刻裂紋尖端的應變能釋放率。

4 驗證算例

圖5 三參數線粘彈性模型

對于上述粘彈性結構受階躍載荷時,張淳源[21]給出了結構裂紋尖端應變能釋放率的解析表達式:

(10)

式中C(t)為材料的蠕變柔量;平面應變時,χ=1.0-v2,平面應力時,χ=1.0;KⅠ和KⅡ分別為結構所受階躍載荷變成定常載荷時裂紋尖端的Ⅰ型和Ⅱ型裂紋強度因子。

對于經典的三參數線粘彈性模型,其蠕變柔量的形式為

(11)

4.1 單邊裂紋受拉測試

圖6(a)為一單邊含裂紋缺陷的粘彈性平板,板的長度和寬度分別2h=20和b=10。平板的厚度為1,板中間位置設置長度為a=1的水平裂紋。平板兩端受均布拉應力σ(t)=σ0H(t)作用,H(t)代表單位階躍函數,σ0=1,作用時間為200 s。對于定常均布拉應力σ0作用,裂紋尖端的Ⅰ型強度因子滿足[22]:

(a) Geometric model of tensile specimen with single edge crack (b) Semi-structured low-quality grid model

(12)

其中,

(13)

為方便,定義相對誤差用以表示計算結果的準確性。各個時刻相對誤差表達式如下

(14)

圖6(b)給出了平板的非結構網格劃分情況,網格劃分均采用一階三角形單元。在遠離裂紋尖端的平板邊緣,單元尺寸為2.0,在裂紋尖端附近對網格進行了加密,單元平均尺寸為0.05。

結合裂紋尖端的虛擬網格技術(虛擬單元尺寸為0.1)與虛擬裂紋閉合方法,分別開展泊松比為0.3時,平面應力以及平面應變條件下的Ⅰ型應變能釋放率的計算。從圖7可以明顯地看到,本文提出的虛擬網格技術比文獻[20]中的加料單元方法相比精度大幅提高,相對誤差均在5%以內。加料單元方法的本質是將斷裂參數作為單元的附加自由度,即所謂的“加料”。該方法需要對裂紋尖端相鄰單元的位移模式進行精細化設計。含裂紋結構有限元分析結束后,主程序將會同時獲得斷裂參數和位移場的結果。但是,加料后的單元并沒有改善裂紋尖端應力奇異的現狀,斷裂參數的計算還是在裂紋尖端奇異應力的基礎上獲得的。這是本文提出的虛擬網格技術比加料單元方法精度更高的主要原因。以上算例說明,本文提出的方法在應對低泊松比條件下的Ⅰ型應變能釋放率計算是準確可行的。

(a)State of plane stress (b)State of plane strain

對于推進劑這一類典型的粘彈性材料,其泊松比一般大于0.495,甚至接近0.5。在利用虛擬網格技術時,分別采用兩種數值積分手段,進行積分點上應力應變的數值計算。一種是采用全積分形式,即采用4個積分點進行計算。考慮到虛擬單元方法的本質是對裂紋尖端應力場的重構,為此將全積分形式的虛擬網格方法記作(Stress Recovery Technique-Fully Integration,SRT-F)。另一種方案采用縮減積分的形式,即采用一個高斯積分點,將該方法記作(Stress Recovery Technique-Reduced Integration, SRT-R)。

圖8給出了平面應變狀態時兩種積分方案下,應變能釋放率的計算結果與文獻解的對比情況。由圖8可以看到,虛擬網格應力的全積分計算方案效果較差,相對誤差接近100%。與同樣采用縮減積分方案的文獻解相比,本文的虛擬網格技術獲得的應變能釋放率精度更高,在加載后期,相對誤差低于1%。由此可見,虛擬網格縮減積分技術可以應對高泊松比條件下的Ⅰ型應變能釋放率的數值計算。

圖8 平面應變狀態階躍加載下應變能釋放率相對誤差隨時間變化(v=0.495)

4.2 單邊裂紋受剪測試

圖9(a)為一單邊含裂紋缺陷的粘彈性平板,板的長度和寬度分別2h=20和b=10。平板的厚度為1,板中間位置設置長度為a=3的水平裂紋。平板左側邊緣中間位置分別受一對相互平行集中力Q(t)=Q0H(t)作用,Q0=1,作用時間為200 s。對于定常集中力Q0作用,裂紋尖端的Ⅱ型強度因子滿足[23]:

(a) Geometric model of shear specimen with single edge crack (b) Semi-structured low-quality grid model

(15)

其中,fⅡ=2.5935。

圖9(b)給出了平板的非結構網格劃分情況,網格劃分同樣采用一階三角形單元。在遠離裂紋尖端的平板邊緣,單元尺寸為2.0,在裂紋尖端附近對網格進行了加密,單元平均尺寸同樣為0.05。

圖10給出了平面應變狀態時兩種積分方案下,應變能釋放率的計算結果與文獻解的對比情況。虛擬網格應力的全積分計算方案效果較差,相對誤差在90%左右。本文的虛擬網格技術獲得的應變能釋放率相對誤差低于7%。另外一方面,文獻解的精度較差,加載后期一度超過40%。可以認為,本文提出的虛擬網格技術能夠應對高泊松比條件下的Ⅱ型應變能釋放率的高精度數值計算。

圖10 平面應變狀態階躍加載下應變能釋放率

進一步地,為了檢驗復雜載荷歷程下,虛擬網格技術的計算能力,按照文獻[20],設計了如圖11兩種載荷歷程形式,其中相互平行集中力Q(t)=Q0f(t)。

(a)Load history I (b)Load history II

(a)Load history I (b)Load history II

因此,文獻解與本文提出的虛擬網格全積分求解方案相對誤差較大。第50 s時刻,本文提出的虛擬網格縮減積分求解方案計算出的應變能釋放率為0.499 2,可以認為這是三種方案中,精度最高的求解手段。類似地,在載荷歷程二中,應變能釋放率的最大值不會超過0.538 8,文獻解與虛擬網格全積分求解方案存在較大偏差。

5 結論

本文提出了一種粘彈性應變能釋放率的虛擬網格求解方案,通過兩種典型含裂紋缺陷,結合解析表達式,驗證了求解方案的準確性和有效性,主要結論如下:

(1)虛擬網格技術結合虛擬裂紋閉合方法可以有效地利用低質量網格實現粘彈性裂紋尖端應變能釋放率的高精度數值計算。與文獻提出的加料單元方法相比,實現思路簡潔、計算精度更高并且容易編程實現。此外,由于不需要單獨設計網格,本文提出的方案更加適用于粘彈性裂紋擴展過程中的應變能釋放率計算。

(2)對于泊松比接近于0.5的粘彈性材料,仿真結果表明,本文提出的虛擬網格縮減積分方案同樣可以高效地應對。

(3)本文提出的二維粘彈性應變能釋放率求解方案可以移植和擴展到任意材料二維裂紋前緣的應變能釋放率計算問題中。對于三維應變能釋放率的求解,可以參考本文方法,在裂紋面前緣構建三維虛擬網格并結合虛擬裂紋閉合方法,實現三維裂紋面上斷裂參數的數值計算。此外,本文方法將應用到粘彈性材料的裂紋擴展仿真中,為任意低質量網格條件下的裂紋擴展提供更加精確和科學的判斷。

(4)與文獻[5]求解方案相比,本文提出的粘彈性應變能釋放率虛擬網格計算方法實現方便、精度高并且易于在后期的裂紋擴展仿真研究中實現。