正交各向異性CFRP材料的本構(gòu)關(guān)系及其在平板撞擊模擬中的應(yīng)用

張 昆， 湯文輝， 冉憲文

(國防科技大學(xué) 文理學(xué)院，長沙 410073)

碳纖維增強(qiáng)樹脂基類復(fù)合材料(Carbon-Fiber Reinforced Polymer, CFRP)以低密度、高比模量以及在高溫高壓等極端條件下良好的力學(xué)性能，在航空航天及軍工領(lǐng)域獲得了廣泛的應(yīng)用[1-2]，如用于制造航天飛機(jī)機(jī)翼，尾舵及導(dǎo)彈殼體等。在空天環(huán)境中，上述部件可能隨時受到飛鳥，空間碎片等的高速撞擊，相對撞擊速度從幾百米每秒到幾千米每秒不等，因此研究該類材料在高速撞擊下的動態(tài)力學(xué)響應(yīng)行為，對于提高飛行器整體的安全系數(shù)具有重要意義[3-6]。通常這一類復(fù)合材料具有力學(xué)性能的各向異性，這會對材料的強(qiáng)度，屈服，及應(yīng)力波傳播等產(chǎn)生顯著影響。與傳統(tǒng)金屬類材料相比較，數(shù)值計算上該類材料在本構(gòu)模型及物態(tài)方程(Equation-of State, EOS)的處理上也較為復(fù)雜。對于準(zhǔn)靜態(tài)及低壓情況，應(yīng)力應(yīng)變計算采用各向異性本構(gòu)關(guān)系就能達(dá)到很好的效果。而對于動高壓情況，需要引入物態(tài)方程來計算。值得注意的是，對于各向異性材料，在高壓段材料的體積形變會對壓力產(chǎn)生貢獻(xiàn)，即材料的體變和形變發(fā)生耦合，因此對于這一類材料，EOS需要作出一定修正。對于這一問題，Anderson等[7-8]給出了考慮正交各向異性修正的靜水壓力及偏應(yīng)力，并得到了多項式形式的Grüneisen物態(tài)方程，使得壓力計算可以同時考慮高壓下的體積壓縮非線性及低壓下不同方向的強(qiáng)度差異性。Lukyanov等[9-10]提出了一種對應(yīng)力應(yīng)變張量的廣義分解方法，并同樣采用Grüneisen物態(tài)方程進(jìn)行了平板撞擊問題的模擬計算。美國Sandia實驗室的沖擊動力學(xué)軟件CTH[11]在正交各向異性復(fù)合材料的處理上收錄了上述兩種方法，但目前研究大多數(shù)采用了第一種修正模型。

對于CFRP類材料，在強(qiáng)沖擊載荷作用下，壓力的急劇升高會導(dǎo)致材料內(nèi)部層裂產(chǎn)生及裂紋的擴(kuò)展，這一損傷在宏觀上表現(xiàn)為與金屬材料屈服類似的不可逆形變。文獻(xiàn)[12-13]對這一問題進(jìn)行了研究，嘗試采用Tsai-Hill屈服準(zhǔn)則來描述碳纖維增強(qiáng)復(fù)合材料的不可逆形變問題。事實上屈服準(zhǔn)則的選取對于材料塑性狀態(tài)的判定及后續(xù)塑性應(yīng)變計算具有重要意義，文獻(xiàn)[14]對于各向異性材料的屈服準(zhǔn)則選取進(jìn)行了較全面的介紹。在塑性屈服準(zhǔn)則在這一類材料的使用上，Anderson等采用了一種廣義von-Mises屈服準(zhǔn)則，并成功運(yùn)用于沖擊動力學(xué)程序EPIC92。蔣邦海等[15-16]完成了一系列從準(zhǔn)靜態(tài)到中高應(yīng)變率加載下的某型碳纖維增強(qiáng)酚醛樹脂材料(C/PF)的力學(xué)實驗，得到了該材料的Tsai-Hill屈服準(zhǔn)則參數(shù)并成功運(yùn)用在一維下的高速撞擊模擬中?；诖它S霞等[17-19]構(gòu)建了該材料的二維各向異性動態(tài)本構(gòu)模型，通過自行編寫有限元程序?qū)崿F(xiàn)了矩形，圓柱殼形靶板在平板撞擊，X射線熱擊波輻照等外載荷加載下C/PF材料的熱力學(xué)響應(yīng)模擬。目前關(guān)于CFRP材料的平板撞擊研究大多局限于一維或二維研究，而三維正交各向異性材料本構(gòu)關(guān)系，物態(tài)方程的構(gòu)建及其在動高壓條件下的數(shù)值模擬計算中的應(yīng)用鮮有文獻(xiàn)報道。

基于上述情況，本文以文獻(xiàn)[20]中C/PF材料為研究對象，建立了一套完整的三維各向異性材料本構(gòu)模型，在此基礎(chǔ)上編寫了一顯式拉格朗日有限元程序，對一有限尺寸平板撞擊問題進(jìn)行了數(shù)值計算。

1 理論模型

1.1 彈性段本構(gòu)關(guān)系

對于CFRP類材料，碳纖維作為增強(qiáng)基承擔(dān)了絕大部分的力學(xué)載荷，因此碳纖維的力學(xué)性能及其編織方式?jīng)Q定了材料整體的力學(xué)性能。本文所研究C/PF材料是由碳纖維二維正交平紋機(jī)織成碳布后層疊，再加入酚醛樹脂熱壓固化生成。在主軸坐標(biāo)系下材料存在三個力學(xué)主軸方向，分別為碳纖維經(jīng)向，碳纖維緯向，碳布厚度方向，如圖1所示。

圖1 C/PF材料示意圖Fig.1 A schematic of C/PF material

三維應(yīng)變條件下， 彈性段應(yīng)力分量σij與應(yīng)變分量εij的本構(gòu)關(guān)系采用廣義Hook定律描述

(1)

式中：Cij為剛度矩陣系數(shù)，其分量由材料在三個主軸方向上的楊氏模量、泊松比及剪切模量確定。對應(yīng)力及應(yīng)變進(jìn)行球-偏分解，正交各向異性材料彈性段的靜水壓p表示為

(2)

(3)

(4)

1.2 塑性本構(gòu)模型及Tsai-Hill屈服準(zhǔn)則

當(dāng)材料應(yīng)力狀態(tài)在應(yīng)力空間超越屈服面時，材料進(jìn)入塑性變形階段。根據(jù)塑性增量理論，應(yīng)力是過程量與應(yīng)變歷史相關(guān)，而應(yīng)力增量[dσ]與彈性應(yīng)變增量[dεe]之間仍滿足廣義Hook定律

[dσ]=[C][dεe]=[C]([dε]-[dεp])

(5)

式中： [C]為剛度矩陣； [dε]為應(yīng)變增量張量； [dεp]為塑性應(yīng)變增量張量。同樣對[dσ]進(jìn)行分解，得到塑性段內(nèi)壓力增量為

(6)

可以看到，在塑性段壓力不但與體應(yīng)變相關(guān)，同時也會受到偏應(yīng)變增量及塑性應(yīng)變增量的影響。對于塑性應(yīng)變的計算，材料的應(yīng)力狀態(tài)始終在屈服面上，根據(jù)正交性法則和一致性法則，塑性應(yīng)變增量可表示為

(7)

式中： [?F/?σ]為塑性勢函數(shù)即屈服函數(shù)對應(yīng)力的偏導(dǎo)張量。對于正交各向異性復(fù)合材料。本文中對于C/PF材料采用一9參數(shù)平方形式的Tsai-Hill準(zhǔn)則

(8)

(9)

式中：Yij分別為主軸方向上的單軸壓縮強(qiáng)度和相對應(yīng)方向上的純剪強(qiáng)度。λ為塑性流動因子

(10)

1.3 Grüneisen物態(tài)方程(EOS)的修正

在式(6)中， dpp與θ之間的線性關(guān)系只有在低壓段成立。在高壓下，二者之間呈非線性關(guān)系，因此需要額外引入物態(tài)方程來計算壓力。固體材料壓縮常采用Grüneisen物態(tài)方程

p=pH(v)+ρΓ(e-eH)

(11)

式中：ρ為密度；v為比體積；e為比內(nèi)能；Γ為Grüneisen系數(shù)；pH,eH組成沖擊壓縮線。

(12)

(13)

式中：c0和s為Hugoniot參數(shù)；ρ0為初始密度；v，v0分別為當(dāng)前及初始比體積。定義壓縮比

(14)

在小應(yīng)變條件下忽略應(yīng)變高次項有

(15)

對式(11)做關(guān)于μ的多項式展開，三階截斷并代入式(15)，得到多項式形式的物態(tài)方程為

(16)

其中，各項參數(shù)可表示為

(17)

對式(16)求全微分，得到增量形式的壓力dp為

(18)

式(16)和式(18)目前僅僅適用于各向同性材料，并沒有反映材料的各向異性特點。對于沖擊動力學(xué)問題，壓力由物態(tài)方程計算，在低壓彈性段，由本構(gòu)關(guān)系得到的pe的表達(dá)式(6)應(yīng)當(dāng)是物態(tài)方程式(18)在低壓下的極限情況，同時也要考慮到高壓段偏應(yīng)變及塑性應(yīng)變增量對壓力的影響，因此對式(15)可做如下修正：①用廣義體積模量A′替換物態(tài)方程中體應(yīng)變一階項系數(shù)A1(事實上A1即為傳統(tǒng)體積模量)；②在方程中添加偏應(yīng)變對壓力貢獻(xiàn)項；③方程中添加塑性應(yīng)變對壓力貢獻(xiàn)項。基于上述修正可得到正交各向異性修正下的壓力為

(19)

目前涉及到正交各向異性CFRP材料的本構(gòu)方程及物態(tài)方程的研究，在工程上多采用與式(16)相似的體應(yīng)變?nèi)康亩囗検叫问?，或只考慮了偏應(yīng)變項而忽略塑性應(yīng)變項的影響[21]。而本文給出的三維應(yīng)變條件下正交各向異性的本構(gòu)模型能充分體現(xiàn)壓力隨體應(yīng)變的非線性關(guān)系，同時材料正交各向異性力學(xué)性能對壓力的各個影響因素也已考慮在內(nèi)。在后續(xù)計算中本文將對各個修正項的影響進(jìn)行定量分析。

2 算例及分析

為了對本文提出的模型進(jìn)行驗證，本文對該材料的平板撞擊(Planar Plate Impact, PPI)問題進(jìn)行了模擬。初始設(shè)定靶板及飛片中碳纖維經(jīng)向，緯向及碳布厚度三個主軸方向依次與實驗室坐標(biāo)系下xyz坐標(biāo)軸方向重合，靶板與飛片尺寸相同均為0.25 cm×1.0 cm×1.0 cm。飛片沿x坐標(biāo)軸方向，以速度w撞擊靜止放置的靶板，如圖2所示。

圖2 平板撞擊實驗示意圖Fig.2 A schematic of planar plate impact experiment

有限元前處理采用Ansys生成標(biāo)準(zhǔn).k格式文件進(jìn)行讀入，網(wǎng)格尺寸為0.012 5 cm×0.025 cm×0.025 cm。計算采用自行編寫的顯式動態(tài)拉格朗日描述有限元程序[22]。材料物態(tài)方程參數(shù)為：ρ0=1.38 g/cm3，c0=2.35 km/s，s=1.66，Γ=2.32。材料在其主軸方向的力學(xué)參數(shù)及屈服函數(shù)參數(shù)見表1和表2。

表1 C/PR材料在主軸方向上的力學(xué)參數(shù)Tab.1 Mechanical parameters of C/PR material in principal axis directions

表2 C/PR材料在主軸方向上的Tsai-Hill屈服準(zhǔn)則參數(shù)Tab.2 Tsai-Hill yield parameters of C/PR material in principal axis directions GPa

以飛片速度w=1 500 m/s為例，圖3和圖4分別給出飛片與靶板整個模型在t=0.3 μs時刻的三維壓力分布以及z=0.5 cm處xy平面內(nèi)壓力分布。

圖3 飛片及靶板三維壓力云圖Fig.3 3D contour of pressure in flyer and target

圖4 xy截面壓力云圖Fig.4 Contour of pressure in xy section

從圖3和圖4可以看出，碰撞發(fā)生后對稱的壓力波分別向靶板及飛片的前后自由面處傳播，同時在碰撞面附近產(chǎn)生了材料的擠壓變形和側(cè)向膨脹，在碰撞面周邊的邊側(cè)自由面上產(chǎn)生了側(cè)向稀疏波，波的傳播呈現(xiàn)出明顯的三維效應(yīng)。

圖5給出了該時刻，飛片-靶板中軸線上的三個主應(yīng)力分布情況?？梢钥吹?，中軸線上三個主應(yīng)力的分布并不相同，x方向即碰撞方向應(yīng)力值最高，y方向應(yīng)力高于z方向應(yīng)力。這是因為碳纖維緯度方向與碳布厚度方向這兩個主軸方向雖都與速度方向垂直，但不同主軸方向上的模量，泊松比等存在明顯差異，經(jīng)計算得到的主應(yīng)力也完全不同。這一應(yīng)力分布顯著不同于各向同性材料，反映了各向異性材料的力學(xué)特點。

為了對本文所建立本構(gòu)模型做進(jìn)一步分析，在上述工況下本文比較了采用原始各向同性物態(tài)方程式(18)及修正后的物態(tài)方程式(19)計算的結(jié)果。t=0.3 μs時刻靶板中軸線上的壓力分布如圖6所示。

取靶板內(nèi)中軸線上距碰撞面0.068 cm處的單元作為觀測點，記錄得到兩種模型的壓力歷史曲線如圖7所示。

圖5 中軸線應(yīng)力分布Fig.5 The distribution of stresses alone the axis

圖6 壓力沿軸線分布Fig.6 The distribution of pressure alone the axis of two models

圖7 固定點處壓力歷史曲線Fig.7 The history curve of pressure at the fixed point

從圖7可以看到，無論從某一時刻壓力的空間分布還是固定點處的壓力歷史曲線都顯示采用原始物態(tài)方程計算得到的壓力較高。為了更進(jìn)一步分析這一差異特性與飛片加載速度的關(guān)系，本文設(shè)定飛片速度從100～3 000 m/s，對這一平板撞擊實驗行了多次數(shù)值計算，計算結(jié)果見表3。

表3 不同飛片速度下原始物態(tài)方程及修正物態(tài)方程計算結(jié)果Tab.3 Results of the original EOS and modified EOS with different flyer velocity

當(dāng)w=100 m/s，壓力相對偏差最大達(dá)到27.2%，隨著飛片速度增加逐步降低。當(dāng)飛片速度達(dá)到3 000 m/s時，壓力的相對偏差小于3%，兩種模型計算結(jié)果趨于一致。本文分析認(rèn)為，對于平板撞擊問題，低速低壓段，材料的體應(yīng)變較小，此時p與θ近似成線性關(guān)系，比例系數(shù)為體模量。式(18)和式(19)得到的體積模量A1與廣義體積模量A′分別為7.62 GPa和4.76 GPa差異較大，因此在低速狀態(tài)下，兩種模型壓力相差較大。

圖8 壓力相對偏差隨飛片速度的變化規(guī)律Fig.8 Variation of p relative error with the velocity of flyer

隨著飛片速度w的上升，材料體應(yīng)變增加，θ的高階項對壓力的影響隨之增大， 壓力p與體應(yīng)變θ成非線性關(guān)系。而兩種模型中表征高壓下材料特性的物態(tài)方程參數(shù)ρ0，c0，s，Γ是相同的，即p與θ的二階三階關(guān)系是一致的。因此兩種模型在低壓段差異較大而高壓區(qū)趨于一致的計算結(jié)果是合理的。另一方面，為了定量的研究偏應(yīng)變修正項及塑性應(yīng)變修正項對壓力的影響，因此本文分別定義

(20)

(21)

其中積分時間為0時刻至壓力穩(wěn)定為平臺壓力的時刻ts。從表3數(shù)據(jù)計算可得，在w=100 m/s時，偏應(yīng)變修正項p_εd對壓力p的貢獻(xiàn)絕對值最高， 接近-10 %。但隨著速度的上升，p_εd的絕對值增高但在壓力中所占百分比下降，材料在高壓下呈現(xiàn)出流體特性，各向異性特性趨于消失。而塑性應(yīng)變項p_εp對壓力的貢獻(xiàn)在計算中均未超過-2%，所占百分比非常有限，因此在各個階段的計算中均可以忽略不計。

3 結(jié) 論

(1) 基于彈塑性理論，從廣義Hook定律出發(fā)，導(dǎo)出了正交各向異性材料在彈性和塑性階段應(yīng)力增量與應(yīng)變增量之間的關(guān)系式；采用Tsai-Hill準(zhǔn)則描寫了各向異性材料的塑性應(yīng)變增量；采用Grüneisen物態(tài)方程描述了壓力與體應(yīng)變之間非線性關(guān)系。

(2) 將正交各向異性材料的三維本構(gòu)模型與三維有限元計算結(jié)合起來，對C/PR正交各向異性材料飛片與靶板的碰撞過程進(jìn)行了數(shù)值模擬。結(jié)果表明，所建立的本構(gòu)模型能夠反映出材料力學(xué)性能的各向異性特性。

(3) 數(shù)值模擬結(jié)果表明，飛片速度較低時，正交各向異性模型得到的壓力顯著小于各向同性模型。隨著飛片速度升高，平臺壓力上升，兩種模型計算結(jié)果趨于一致。

(4) 對于正交各向異性模型中引入的偏應(yīng)變修正項，在低壓條件下會對壓力產(chǎn)生一定影響。而塑性應(yīng)變增量項無論飛片速度高低，均可忽略而不造成太大誤差。