平面應(yīng)變各向異性本構(gòu)關(guān)系及在應(yīng)力波傳播模擬中的應(yīng)用*

黃 霞,湯文輝,蔣邦海

(國防科技大學(xué)理學(xué)院技術(shù)物理研究所,湖南 長沙410073)

1 引 言

近年來,復(fù)合材料在國防領(lǐng)域得到了越來越廣泛的應(yīng)用,以高強度、高剛度、低密度等特點,已成為航空、航天等國防工業(yè)部門的一種重要工程材料。在航空、航天等領(lǐng)域中,材料的外在環(huán)境非常復(fù)雜,可能面臨高速撞擊、射線輻射等動載荷環(huán)境,因此對復(fù)合材料動態(tài)響應(yīng)特征的研究,對提高材料性能、加強航天器的安全性有著非常重要的作用。在研究復(fù)合材料力學(xué)性能的過程中,必須考慮各向異性力學(xué)特征,它會對強度、應(yīng)力波傳播等帶來影響,為了分析復(fù)合材料的各向異性響應(yīng)特征,必須使用各向異性本構(gòu)模型。為了處理各向異性材料本構(gòu)模型中容變和畸變的耦合效應(yīng),C.E.Anderson 等[1]、P.E.O'Donoghue 等[2]將各向異性條件下的靜水壓及應(yīng)力偏量表達式進行了修正;另外,他們將物態(tài)方程引入到各向異性本構(gòu)模型中,使得修改后的Grüneisen 物態(tài)方程既能反映高壓下的體積壓縮非線性,又能考慮低壓下材料的各向異性強度性能。各向異性強度準(zhǔn)則是各向異性本構(gòu)模型研究中的一個重要問題,從各向同性強度準(zhǔn)則基礎(chǔ)上發(fā)展起來的適用于復(fù)合材料的強度準(zhǔn)則已有十幾種[3],最常用的是T sai-Hill 屈服準(zhǔn)則、Tsai-Wu 屈服準(zhǔn)則等。C.E.Anderson 等[4]將Tsai-Hill 強度準(zhǔn)則作為各向異性理想塑性屈服準(zhǔn)則,給出了各向異性塑性變形的計算方法,并將這個方法應(yīng)用到了大型沖擊動力學(xué)程序EPIC 當(dāng)中。李永池等[5-6]發(fā)展了以Tsai-Hill 屈服準(zhǔn)則和Johnson-Cook 模型為基礎(chǔ)的橫觀各向同性粘塑性本構(gòu)模型,并用于應(yīng)力波傳播的數(shù)值模擬。蔣邦海等[7]研究了碳纖維增強型復(fù)合材料率相關(guān)的Tsai-Hill 屈服準(zhǔn)則,并用于一維應(yīng)力波傳播的數(shù)值模擬。此外,由于復(fù)合材料在宏觀上表現(xiàn)出類似金屬的彈塑性特征,可把用于金屬動態(tài)性能方面的彈塑性理論方法用于復(fù)合材料上[8]。自20 世紀(jì)90 年代,就開展了關(guān)于各向異性本構(gòu)模型在數(shù)值模擬中應(yīng)用的研究,其中部分研究成果已用在了較新版本的有限元大型沖擊動力學(xué)數(shù)值模擬軟件當(dāng)中,例如LSDYNA950(1999)中就嵌入了橫觀各向同性與正交各向異性彈塑性本構(gòu)模型,使得這些軟件具有分析纖維增強復(fù)合材料動態(tài)響應(yīng)的能力。但是,復(fù)合材料力學(xué)性能除了具有各向異性特性,還具有應(yīng)變率效應(yīng),同時纖維與基體的界面效應(yīng)、損傷、溫度等都將對力學(xué)性能產(chǎn)生影響。因此建立計及多種因素的復(fù)合材料動態(tài)本構(gòu)模型,使得復(fù)合材料動態(tài)響應(yīng)數(shù)值模擬結(jié)果更加可靠,將是今后研究工作的努力方向。

本文中,利用Hooke 定律、Grü neisen 物態(tài)方程及Tsai-Hill 屈服準(zhǔn)則,建立計及體積壓縮非線性,能處理彈性、塑性變形的二維應(yīng)變正交異性彈塑性本構(gòu)模型。討論在該本構(gòu)模型中容變率和畸變率耦合處理的問題、物態(tài)方程引入的問題以及坐標(biāo)軸旋轉(zhuǎn)、應(yīng)力修正等方面的內(nèi)容。最后,結(jié)合這些理論采用顯式有限元方法,自行編寫程序模擬某纖維增強型復(fù)合材料碰撞過程中平面應(yīng)力波的傳播。

2 平面應(yīng)變正交各向異性彈塑性本構(gòu)關(guān)系

2.1 平面應(yīng)變正交異性彈性本構(gòu)關(guān)系

對于正交各向異性材料,彈性階段的應(yīng)力應(yīng)變關(guān)系可用Hooke 定律描述,考慮到在平面應(yīng)變(假設(shè)為(1,2)平面)條件下,ε33=ε13=ε23=0,σ13=σ23=0,材料參數(shù)具有對稱性

式中:cij為與材料彈性模量、泊松比和剪切模量相關(guān)的剛度矩陣系數(shù)。

對于各向同性材料,容變率和畸變率解耦處理已應(yīng)用于大多數(shù)沖擊動力學(xué)程序,反映在算法上就是靜水壓和偏應(yīng)力分別利用物態(tài)方程和本構(gòu)關(guān)系計算。但是對于各向異性材料他們不能簡單解耦,因為靜水壓可能導(dǎo)致形狀改變,而應(yīng)力偏量也會造成體積變化,因此采用平均正應(yīng)力代替各向同性材料中靜水壓的概念,使得正交異性材料的平均應(yīng)力和應(yīng)力偏量在形式上可以解耦,便于在計算程序中應(yīng)用[2]。于是定義拉為正,壓為負(fù),將應(yīng)力σ分解為平均正應(yīng)力p=-(σ11+σ22+σ33)/3 和偏應(yīng)力s,同時將應(yīng)變ε分解為體應(yīng)變θ=ε11+ε22和偏應(yīng)變e,則根據(jù)式(1)可得到

偏應(yīng)力s 各分量也可根據(jù)相應(yīng)的公式求出。

以沖擊絕熱線為參考的Grü neisen 物態(tài)方程為

將上式表示成下面的多項式形式

為了使得所計算的平均正應(yīng)力既能反映高壓下的體積壓縮非線性效應(yīng),又能體現(xiàn)出低壓下材料的各向異性特征[2],結(jié)合式(2)和(4),得到

2.2 平面應(yīng)變正交異性塑性本構(gòu)關(guān)系

在塑性變形時,應(yīng)力狀態(tài)與變形路徑或歷史有關(guān),應(yīng)力、應(yīng)變間沒有一一對應(yīng)關(guān)系,但應(yīng)力增量與彈性應(yīng)變增量之間滿足H ooke 定律,因此本構(gòu)關(guān)系用增量形式表達。應(yīng)變增量可分解為彈性應(yīng)變增量和塑性應(yīng)變增量,即,有

同樣引入物態(tài)方程修正平均正應(yīng)力,可得到修正后塑性變形時的平均應(yīng)力增量

同理可以得到相應(yīng)的應(yīng)力偏量增量表達式,進而得到修正后的應(yīng)力值。

至此已經(jīng)給出平面應(yīng)變條件下,既考慮了高壓下體積壓縮非線性,在低壓下又能反映各向異性力學(xué)性能的正交異性彈塑性本構(gòu)關(guān)系。

2.3 Tsai-Hill 屈服準(zhǔn)則

判斷材料是處于彈性變形階段還是塑性變形階段需要用到屈服準(zhǔn)則,各向異性材料常用的屈服判據(jù)是Tsai-Hill 屈服準(zhǔn)則。在平面應(yīng)變條件下,Tsai-Hill 屈服準(zhǔn)則的基本形式為

當(dāng)材料進入塑性變形階段后,由式(6)計算應(yīng)力增量時需要計算塑性應(yīng)變增量,于是從Tsai-Hill 屈服準(zhǔn)則出發(fā),利用經(jīng)典塑性流動理論計算塑性應(yīng)變增量[4]。根據(jù)正交性法則和一致性法則得到

式(9)等價于給出了塑性應(yīng)變增量,因此當(dāng)給定某一時刻的應(yīng)力狀態(tài)時,即可由屈服準(zhǔn)則判斷是否處于塑性變形階段。如果發(fā)生塑性變形,則由式(9)求解塑性應(yīng)變增量,再根據(jù)各向異性塑性本構(gòu)關(guān)系得到的偏應(yīng)力增量和修正后的平均應(yīng)力增量獲得應(yīng)力增量,從而求得下一時刻的應(yīng)力狀態(tài),如此循環(huán)便可獲得整個時間域上的解。

3 材料主軸的旋轉(zhuǎn)及客觀應(yīng)力率

各向異性材料的物理性能及本構(gòu)關(guān)系都是基于材料主軸方向成立的,對于二維問題,材料變形會導(dǎo)致材料主軸的偏轉(zhuǎn),因此必須考慮材料主軸與系統(tǒng)坐標(biāo)之間的轉(zhuǎn)換問題。在二維問題中,用(x,y)表示系統(tǒng)坐標(biāo)系,用(1,2)表示材料主軸坐標(biāo)系,并且假定初始時刻系統(tǒng)坐標(biāo)系與材料主軸坐標(biāo)系重合,任一時刻材料某處的主軸坐標(biāo)系相對整體坐標(biāo)系旋轉(zhuǎn)α角度(以逆時針為正),關(guān)于二維問題中α的求解可采取下面方法。

按照連續(xù)介質(zhì)力學(xué)的極分解定理,任一時刻材料某處變形梯度矩陣可分解為如下形式

式中:R 為正交旋轉(zhuǎn)矩陣,U 為正定拉伸矩陣,根據(jù)U=RTF 的正定性,可得R 的表達式。

獲得旋轉(zhuǎn)矩陣R 后,就可以將系統(tǒng)坐標(biāo)系中求得的應(yīng)變、變形率等物理量轉(zhuǎn)換到材料坐標(biāo)系中,使得材料本構(gòu)關(guān)系能夠正確使用。2 個坐標(biāo)系之間的應(yīng)變、變形率、柯西應(yīng)力張量轉(zhuǎn)換有如下關(guān)系式

從式(11)可以看出應(yīng)變、變形率、應(yīng)力張量均為標(biāo)價無關(guān)張量,即他們是客觀的。但是對于柯西應(yīng)力率張量﹒σ有

該式并不滿足式(11),因此柯西應(yīng)力率張量﹒σ不是客觀的,不能代表反映物體應(yīng)力狀態(tài)變化的應(yīng)力變化率。而彈塑性本構(gòu)關(guān)系是以率的形式給出的,為此引入客觀應(yīng)力率。通常采用的客觀應(yīng)力率有Jaumann 率、Truesdell 率和Green-Naghdi 率,由于簡單剪切的大變形彈塑性計算中Jaumann 率造成不正確的響應(yīng),因此彈塑性材料使用Green-Naghdi率[9]?；拘问綖?/p>

式中:Ω=﹒RRT, σ為系統(tǒng)坐標(biāo)系中的柯西應(yīng)力。

為了便于描述,將主軸坐標(biāo)系中的σ1-2記為,系統(tǒng)坐標(biāo)系中的σx-y記為σ,其他量也采取同樣記法。事實上材料主軸坐標(biāo)系相當(dāng)于嵌在材料中的一個旋轉(zhuǎn)坐標(biāo)系,σ1-2相當(dāng)于旋轉(zhuǎn)柯西應(yīng)力,根據(jù)式(12)～(13)容易得到,從而有

由于Green-Naghdi 應(yīng)力率與旋轉(zhuǎn)柯西應(yīng)力率之間的聯(lián)系,在應(yīng)用Green-Naghdi 應(yīng)力率時,材料特性在非線性旋轉(zhuǎn)構(gòu)形中處理。因此從n 時刻的應(yīng)力σn修正到n+1 時刻的應(yīng)力σn+1的計算步驟為:

(1)已知n 時刻的σn、Rn,計算旋轉(zhuǎn)柯西應(yīng)力

4 模擬平面應(yīng)變條件下各向異性材料中的應(yīng)力波傳播的數(shù)值算例

用一種纖維增強型復(fù)合材料(TF 材料)作各向異性材料,3 個主方向分別為材料厚度方向、纖維布經(jīng)向和緯向,如圖1 所示。TF 材料參數(shù)分別為:ρ=1.38 g/cm3, c0=2.35 km/s,s=1.66, γ=2.32。假設(shè)有2 個形狀和結(jié)構(gòu)完全相同的TF 材料A、B,x、y 方向尺寸為2 cm×6 cm,z 方向尺寸無限大,B 初始靜止,A 以300 m/s 的初始速度沿x 方向與B 發(fā)生正碰撞。該問題可簡化為平面應(yīng)變碰撞問題,模型如圖2 所示。為了驗證模型及程序的正確性及更好地說明材料各向異性彈塑性力學(xué)性能,作了3 次不同的模擬。第1 次將碰撞方向即x 方向取為TF 材料厚度方向,y 方向取為材料經(jīng)向,z 方向取為材料緯向;第2 次將x 方向取為材料經(jīng)向,y 方向取為材料厚度方向,z 方向取為材料緯向;第3 次模擬僅考慮彈性本構(gòu)模型,其他條件與前2 次模擬相同。3 次模擬相應(yīng)的材料參數(shù)列于表1 中[10]。

圖1 復(fù)合材料主方向及其鋪層Fig.1 Principal directions and lamina of the composite

圖2 用于數(shù)值模擬的簡化模型Fig.2 A simplified model for numerical simulation

表1 材料參數(shù)Table 1 Material constants

使用Tecplot 圖形處理軟件給出了沿TF 材料厚度方向和經(jīng)向碰撞時σxx的等值云圖,如圖3 ～4 所示?？梢钥闯?應(yīng)力波在傳播過程中表現(xiàn)出明顯的二維效應(yīng)和各向異性特征,材料中有正向沖擊波和側(cè)向稀疏波的傳播,沿經(jīng)向碰撞時產(chǎn)生的正向沖擊波速度比沿厚度方向碰撞時大,而側(cè)向稀疏波傳播速度比沿厚度方向時小,這正是材料經(jīng)向彈性模量比厚度方向高的結(jié)果。為了更好地定量分析TF 材料中應(yīng)力波的傳播特征,考慮圖1 中的對稱軸線y=3 cm 上的點沒有y方向的位移,在上下兩側(cè)稀疏波到達以前應(yīng)力波沿該軸線相當(dāng)于在一維應(yīng)變條件下傳播。該線上的應(yīng)力空間分布如圖5 ～7 所示。

由圖5～6 可以看出,在t=4.8 μs 時上下兩側(cè)向稀疏波尚未到達中心軸線,軸線上表現(xiàn)為一維應(yīng)變狀態(tài)。在t=11.6 μs 時,由于應(yīng)力波已到達左右自由邊界,在左右兩端分別向右向左傳播稀疏波,使波寬逐漸減小。圖5(b)中上下兩側(cè)向稀疏波已到達中心軸線,使得各方向主應(yīng)力值均增加,并且沿y 方向的拉伸作用最強,體現(xiàn)出二維特點,而圖6(b)中則沒有。從理論上估算,圖5 和圖6 所對應(yīng)材料y 方向的彈性波速分別為2.607 和2.136 km/s,側(cè)向稀疏波最先到達軸線的時間分別為約11 和14 μs,因此在11.6 μs 時圖5(b)已經(jīng)出現(xiàn)明顯的拉伸作用,而圖6(b)中沒有。此外注意到沿TF 材料經(jīng)向碰撞時,彈性前驅(qū)波比較明顯,材料發(fā)生塑性變形更早,且正向應(yīng)力峰值(絕對值)略大,σyy與σzz差異更大,與沿TF 材料經(jīng)向碰撞時有所不同,這樣的差異性正是纖維鋪層厚度方向與纖維增強方向力學(xué)性質(zhì)不同的體現(xiàn)。

圖3 沿TF 材料厚度方向碰撞,使用彈塑性本構(gòu)模型時σxx 的等值云圖Fig.3 σxx contour for elastic-plastic constitutive model while compacting along TF thickness direction

圖4 沿TF 材料經(jīng)向碰撞,使用彈塑性本構(gòu)模型時σxx 的等值云圖Fig.4 σxx contour for elastic-plastic constitutive model w hile com pacting along TF w arp direction

圖5 沿TF 材料厚度方向碰撞,使用彈塑性本構(gòu)模型時對稱軸線y=3 cm 上的應(yīng)力空間分布Fig.5 Stress waves along y=3 cm for elastic-plastic model w hile compacting along thickness direction

圖6 沿TF 材料經(jīng)向方向碰撞,使用彈塑性本構(gòu)模型時對稱軸線y=3 cm 上的應(yīng)力空間分布Fig.6 Stress waves along y=3 cm for elastic-plastic model w hile compacting along warp direction

從圖7 中看出,當(dāng)各向異性材料按照彈塑性本構(gòu)模型計算時,應(yīng)力波在傳播過程中出現(xiàn)明顯的彈性前驅(qū)波,正向應(yīng)力峰值(絕對值)比按照純彈性本構(gòu)模型的計算值小,表現(xiàn)出彈塑性傳播特點。以上數(shù)值模擬結(jié)果與理論分析的一致性驗證了模型的正確性及程序的可靠性。

5 結(jié) 論

給出了平面應(yīng)變條件下正交各向異性材料彈塑性本構(gòu)模型,討論了在該正交異性本構(gòu)模型中容變律和畸變律耦合處理、物態(tài)方程引入以及坐標(biāo)軸旋轉(zhuǎn)、應(yīng)力修正等方面的問題。并以TF 材料碰撞問題為例,將該本構(gòu)模型嵌入自行編制的動態(tài)顯式有限元程序中,模擬平面應(yīng)變條件下應(yīng)力波傳播規(guī)律。通過對數(shù)值模擬結(jié)果的分析表明:

(1)數(shù)值模擬結(jié)果與理論結(jié)果符合良好,驗證了本構(gòu)模型的正確性及程序的可靠性。

(2)材料中除正向沖擊波傳播還有側(cè)向稀疏波傳播,應(yīng)力波在傳播過程中具有二維特點。

(3)應(yīng)力波在TF 材料中傳播時表現(xiàn)出各向異性特點。沿TF 的厚度方向和經(jīng)向正碰撞時所激發(fā)的沖擊波具有不同的動力學(xué)參量,包括正應(yīng)力平臺峰值、碰撞方向的正向應(yīng)力與側(cè)向應(yīng)力、沖擊波速度,側(cè)向稀疏波速度等。

(4)應(yīng)力波在傳播過程中表現(xiàn)出彈塑性傳播特點。

圖7 沿TF 材料經(jīng)向碰撞,使用彈塑性模型和彈性模型時對稱軸線y=3 cm 上的應(yīng)力空間分布Fig.7 S tress w aves along y=3 cm for elastic-plastic model and elastic model respectively while compacting along w arp direction

[1] Anderson C E,O'Donoghue P E,Skerhut D.A mixture theory approach for the shock response of composite materials[J].Journal of Composite Materials,1990,24:1159-1178.

[2] O'Donoghue P E,Anderson C E,Friesenhahn G J,et al.A constitutive formulation for anisotropic materials suitable for wave propagation computer programs[J].Journal of Composite Materials, 1992,26(13):1860-1884.

[3] Rowlands E R.Strength(failure)theories and their experimental correlation[C]∥Sih G C, Skudra A M.Handbook of Composite(3):Failure Mechanics of Composites.Elsevier Science Publishers, 1985:71-125.

[4] Anderson C E,Cox P A,Johnson G R,et al.A constitutive formulation for anisotropic materials suitable for wave propagation computer programs Ⅱ[J].Computational Mechanics,1994,15:201-223.

[5] 李永池,王紅五,袁福平,等.碳酚醛各向異性本構(gòu)關(guān)系和波傳播[J].寧波大學(xué)學(xué)報(理工版),2000,13(B12):71-76.LI Yong-chi,WANG Hong-wu,YUAN Fu-ping,et al.The anisotropic constitutive relations and stress wave propagation of carbon-reinforced bakelite[J].Journal of Ningbo University(NS EE),2000,13(B12):71-76.

[6] 李永池,譚福利,姚磊,等.含損傷材料的熱粘塑性本構(gòu)關(guān)系及其應(yīng)用[J].爆炸與沖擊,2004,24(4):289-298.LI Yong-chi, TAN Fu-li, YAO Lei,et al.Thermo-viscoplastic constitutive relation of damaged materials with application[J].Explosion and Shock Waves,2004,24(4):289-298.

[7] 蔣邦海,張若棋.一種碳纖維織物增強復(fù)合材料應(yīng)變率相關(guān)的各向異性強度準(zhǔn)則[J].爆炸與沖擊,2006,26(4):333-338.JIANG Bang-hai, ZH ANG Ruo-qi.Strain rate-dependent Tsai-Hill st rength criteria for a carbon fiber woven reinforced composite[J].Explosion and Shock Waves,2006,26(4):333-338.

[8] Yoon K J,Sun C T.Characterization of elastic-viscoplastic properties of an AS4/PEEK thermoplastic composite[J].Journal of Composite Materials, 1991,25:1277-1298.

[9] Belytschko T, Liu W K, Moran B.Nonlinear finite element for continua and structures[M].Beijing:Tsinghua Press, 2002:108-119.

[10] 蔣邦海.正交織物復(fù)合材料的動態(tài)本構(gòu)模型及熱激波研究[D].長沙:國防科學(xué)技術(shù)大學(xué),2006.