平面应变各向异性本构关系及在应力波传播模拟中的应用*

黄 霞,汤文辉,蒋邦海

(国防科技大学理学院技术物理研究所,湖南 长沙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.