平面應變各向異性本構關系及在應力波傳播模擬中的應用*

黃 霞,湯文輝,蔣邦海

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

1 引 言

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

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

2 平面應變正交各向異性彈塑性本構關系

2.1 平面應變正交異性彈性本構關系

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

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

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

偏應力s 各分量也可根據相應的公式求出。

以沖擊絕熱線為參考的Grü neisen 物態方程為

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

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

2.2 平面應變正交異性塑性本構關系

在塑性變形時,應力狀態與變形路徑或歷史有關,應力、應變間沒有一一對應關系,但應力增量與彈性應變增量之間滿足H ooke 定律,因此本構關系用增量形式表達。應變增量可分解為彈性應變增量和塑性應變增量,即,有

同樣引入物態方程修正平均正應力,可得到修正后塑性變形時的平均應力增量

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

至此已經給出平面應變條件下,既考慮了高壓下體積壓縮非線性,在低壓下又能反映各向異性力學性能的正交異性彈塑性本構關系。

2.3 Tsai-Hill 屈服準則

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

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

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

3 材料主軸的旋轉及客觀應力率

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

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

式中:R 為正交旋轉矩陣,U 為正定拉伸矩陣,根據U=RTF 的正定性,可得R 的表達式。

獲得旋轉矩陣R 后,就可以將系統坐標系中求得的應變、變形率等物理量轉換到材料坐標系中,使得材料本構關系能夠正確使用。2 個坐標系之間的應變、變形率、柯西應力張量轉換有如下關系式

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

該式并不滿足式(11),因此柯西應力率張量﹒σ不是客觀的,不能代表反映物體應力狀態變化的應力變化率。而彈塑性本構關系是以率的形式給出的,為此引入客觀應力率。通常采用的客觀應力率有Jaumann 率、Truesdell 率和Green-Naghdi 率,由于簡單剪切的大變形彈塑性計算中Jaumann 率造成不正確的響應,因此彈塑性材料使用Green-Naghdi率[9]。基本形式為

式中:Ω=﹒RRT, σ為系統坐標系中的柯西應力。

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

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

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

4 模擬平面應變條件下各向異性材料中的應力波傳播的數值算例

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

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

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

表1 材料參數Table 1 Material constants

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

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

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

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

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

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

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

5 結 論

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

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

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

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

(4)應力波在傳播過程中表現出彈塑性傳播特點。

圖7 沿TF 材料經向碰撞,使用彈塑性模型和彈性模型時對稱軸線y=3 cm 上的應力空間分布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] 李永池,王紅五,袁福平,等.碳酚醛各向異性本構關系和波傳播[J].寧波大學學報(理工版),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] 李永池,譚福利,姚磊,等.含損傷材料的熱粘塑性本構關系及其應用[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] 蔣邦海,張若棋.一種碳纖維織物增強復合材料應變率相關的各向異性強度準則[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] 蔣邦海.正交織物復合材料的動態本構模型及熱激波研究[D].長沙:國防科學技術大學,2006.