CN111368389B - 一种预测复合材料层合板失效强度的方法 - Google Patents

一种预测复合材料层合板失效强度的方法 Download PDF

Info

Publication number
CN111368389B
CN111368389B CN201910961138.0A CN201910961138A CN111368389B CN 111368389 B CN111368389 B CN 111368389B CN 201910961138 A CN201910961138 A CN 201910961138A CN 111368389 B CN111368389 B CN 111368389B
Authority
CN
China
Prior art keywords
plane
damage
stress
fiber
strain
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN201910961138.0A
Other languages
English (en)
Other versions
CN111368389A (zh
Inventor
陈静芬
杨凤祥
陈善富
刘志明
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Jinan University
Original Assignee
Jinan University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Jinan University filed Critical Jinan University
Priority to CN201910961138.0A priority Critical patent/CN111368389B/zh
Publication of CN111368389A publication Critical patent/CN111368389A/zh
Application granted granted Critical
Publication of CN111368389B publication Critical patent/CN111368389B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • YGENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02TCLIMATE CHANGE MITIGATION TECHNOLOGIES RELATED TO TRANSPORTATION
    • Y02T90/00Enabling technologies or technologies with a potential or indirect contribution to GHG emissions mitigation

Landscapes

  • Investigating Strength Of Materials By Application Of Mechanical Stress (AREA)

Abstract

本发明公开的一种预测复合材料层合板失效强度的方法,包含以下顺序的步骤:步骤一、建立复合材料层合板有限元模型;步骤二、建立复合材料损伤本构模型;步骤三、基于ABAQUS‑VUMAT有限元用户动态子程序模块,使用FORTRAN语言编写用户自定义子程序实现提出的损伤本构模型,以求解应力、应变和损伤;步骤四、对有限元模型进行计算,预测复合材料层合板的失效强度。本发明利用ABAQUS‑VUMAT用户自定义子程序来数值实现所建立的三维损伤本构模型,该模型同时考虑了剪切非线性和损伤累积导致材料性能退化的影响,能准则预测复合材料的失效强度。

Description

一种预测复合材料层合板失效强度的方法
技术领域
本发明涉及复合材料力学性能分析领域,特别涉及一种预测复合材料层合板失效强度的方法。
背景技术
纤维增强树脂基复合材料具有高比模量和比强度、优良的能量吸收性能,特别是各向刚度与强度可设计性等特点,广泛应用在航空航天、军事、海洋、土木和机械等工程领域。由于基体材料具有塑性特点,已有大量实验研究表明很多纤维增强树脂基复合材料(如AS4/3501-6、AS4/PEEK和T300/914等)展现出显著的非线性力学行为。Lafarie-Frenot等(LAFARIE-FRENOT M C,TOUCHARD F.Comparative in-plane shear behaviour of long-carbon-fibre composites with thermoset or thermoplastic matrix[J].CompositesScience and Technology,1994,52(3):417-25.)对T300/914碳纤维/环氧树脂及AS4/PEEK碳纤维/热塑性树基[±45°]复合材料层合板进行了加/卸载反复拉伸试验,试验表明,这两种材料展现出显著的残余应变。Carlsson等(CARLSSON L A,ARONSSON C G,B CKLU NDJ.Notch sensitivity of thermoset and thermoplastic laminates loaded intension[J].Journal of Materials Science,1989,24(5):1670-82.)对碳纤维/环氧树脂[±45°]复合材料层合板进行了拉伸试验,试验表明,该复合材料层合板的剪切方向展现出显著的非线性行为。另外,由于复合材料组分纤维和基体存在不同的力学属性导致其力学行为和破坏模式非常复杂,纤维断裂、基体开裂和分层损伤等破坏模式能够独立或同时发生。复合材料层合板内部微裂纹、纤维断裂和基体开裂等会随着外荷载的增加发生损伤初始并演化,由损伤演化导致的材料属性退化是其另一重要力学特性。因此,为了准确预测复合材料的失效强度,在建立描述复合材料力学行为的本构关系时应同时包含非线性和材料属性退化两个特性。
近年来针对复合材料层合板提出的损伤本构模型未能合理考虑复合材料存在的非线性力学行为、面外应力的影响、合适的失效准则及损伤演化,基本集中在:只考虑材料的弹性行为,现有研究表明,如果模型未考虑非线性效应可能会低估复合材料结构的能量吸收能力;使用Hashin准则来判断层内损伤的起始萌生,Hashin准则虽然可以区分纤维与基体的不同损伤模式,但无法说明损伤产生的物理机制(李力,黄哲峰,杨增钦,et al.基于三维Puck失效准则及唯象模量退化的复合材料臂杆屈曲分析[J].复合材料学报,1-11.),对于适当的横向压缩能抑制剪切破坏发生的现象无法给出较为合理的解释,所以现有一小部分技术提出使用Puck准则(PUCK A,SCH RMANN H.Chapter 5.6-Failure analysis ofFRP laminates by means of physically based phenomenological models[M]//HINTONM J,KADDOUR A S,SODEN P D.Failure Criteria in Fibre-Reinforced-PolymerComposites.Oxford;Elsevier.2004:832-76.)来判断,使用该准则需要计算基体的失效断裂面角度,而现有的计算断裂面角度的方法是Puck遍历法(PUCK A,SCH RMANN H.Chapter5.6-Failure analysis of FRP laminates by means of physically basedphenomenological models[M]//HINTON M J,KADDOUR A S,SODEN P D.Failure Criteriain Fibre-Reinforced-Polymer Composites.Oxford;Elsevier.2004:832-76.)、黄金分割法(WIEGAND J,PETRINIC N,ELLIOTT B.An algorithm for determination of thefracture angle for the three-dimensional Puck matrix failure criterion for UDcomposites[J].Composites Science and Technology,2008,68(12):2511-7.)以及分区黄金分割法(SCHIRMAIER F J,WEILAND J,K RGER L,et al.A new efficient andreliable algorithm to determine the fracture angle for Puck’s 3D matrixfailure criterion for UD composites[J].Composites Science and Technology,2014,100(19-25.),这些算法存在计算效率低和计算结果不可靠等问题;损伤导致的材料属性退化使用突降式,这样会使材料提前发生破坏,计算结果偏保守。
ABAQUS软件可以基于二维Hashin失效准则对弹性复合材料的失效强度进行预测,但无法直接运用较有优势的Puck失效准则,更无法直接对有非线性力学行为的复合材料进行失效强度的预测。针对现有模型存在的问题,需要开发一个同时考虑复合材料剪切非线性效应和损伤累积导致材料属性退化的三维损伤本构模型。
发明内容
本发明的目的在于克服现有技术的缺点与不足,提供一种预测复合材料层合板失效强度的方法。
本发明的目的通过以下的技术方案实现:
一种预测复合材料层合板失效强度的方法,包含以下顺序的步骤:
步骤一、建立复合材料层合板有限元模型;
步骤二、建立复合材料损伤本构模型;
步骤三、基于ABAQUS-VUMAT有限元用户动态子程序模块,使用FORTRAN语言编写用户自定义子程序实现提出的损伤本构模型,以求解应力、应变和损伤;
步骤四、对有限元模型进行计算,预测复合材料层合板的失效强度。
所述步骤一,具体如下:
复合材料层合板的铺层角度沿厚度方向中平面对称布置,为减少计算量,有限元分析沿厚度方向采用对称边界条件只模拟1/2厚的复合材料层合板,每一层厚度方向只划分1个单元,复合材料层合板有限元模型包含8层C3D8R实体单元;由于孔边周围存在应力集中现象,且损伤也是由此开始延伸,因此需对孔边周围区域进行网格细化;建立参考点与自由端面之间加载方向位移一致性约束条件:受拉荷载均采用位移加载方式,左加载面施加固支约束,右端自由端面外设置一个参考点,然后把参考点和端面进行绑定,在Abaqus/CAE模块中,采用creat constraint方法建立coupling耦合约束方程,此时,将位移荷载施加在参考点上,同时只要输出参考点上的位移和反力即U1和RF1,就能够获得加载端面上的位移与反力。
所述步骤二具体为:
步骤(1):建立含损伤的复合材料层合板本构关系;
复合材料应力-应变本构方程:σ=C(d):εe,
其中,符号“:”表示对两个张量指标的缩并计算;σ是有效应力张量;是名义应力张量,通常是指Cauchy应力张量;εe是弹性应变张量;e表示弹性;C(d)是含损伤单向复合材料层合板的四阶刚度张量;C是未损伤单向复合材料层合板的四阶线弹性刚度张量;d是一维向量(d1,d2,d3,d23,d13,d12),其中d1、d2、d3分别为纤维方向纤维损伤的损伤变量、平面内垂直于纤维方向基体损伤的损伤变量、层间平面外方向分层损伤的损伤变量;d12、d23、d13分别为12、23、13平面内的剪切损伤变量;定义坐标系x1-x2-x3为单向板的自然坐标系,xl-xn-xt为断裂面的局部坐标系,两个坐标系下的xl轴重合;12、23、13平面分别对应坐标系下的x1x2平面、x2x3平面、x1x3平面;
把损伤变量引入刚度矩阵,使刚度随着损伤的发展而逐渐变弱,即:
C(d)=M-1(d):C:MT,-1(d);
其中,M-1(d)为M(d)的逆矩阵,MT,-1(d)为M(d)转置矩阵的逆矩阵;M(d)为损伤因子张量,其损伤主轴系下矩阵形式表示如下:
复合材料主坐标系中单层板的三维正交各项异性损伤本构模型如下:
所述复合材料主坐标系为单向板的自然坐标系x1-x2-x3
其中:
其中,σ1、σ2和σ3分别为纤维方向、垂直于纤维方向、层间平面外方向的名义正应力(图1中,x1-x2-x3坐标系下x1表示纤维方向、x2表示垂直于纤维方向、x3表示层间平面外方向);τ23、τ12和τ13分别为图1所示的x1-x2-x3坐标系下x1x2平面、x2x3平面、x1x3平面内的剪应力;ε1、ε2和ε3分别为纤维方向、垂直于纤维方向、层间平面外方向的工程正应变;γ23、γ13和γ12为图1所示的x1-x2-x3坐标系下x1x2平面、x2x3平面、x1x3平面内的工程剪应变;E1、E2、E3分别为纤维方向、垂直于纤维方向、层间平面外方向的未损伤单向复合材料单层的弹性模量,G23、G13、G12分别为x1x2平面、x2x3平面、x1x3平面内未损伤单向复合材料单层的剪切模量,ν12、ν13、ν23分别为纤维方向与垂直于纤维方向的、纤维方向与层间平面外方向的、垂直于纤维方向与层间平面外方向的泊松比,ν21、ν31、ν32分别为垂直于纤维方向与纤维方向的、层间平面外方向与纤维方向的、层间平面外方向与垂直于纤维方向的泊松比,满足关系式,/>
步骤(2):建立三维Puck失效准则来判断纤维和基体损伤,三维Hou准则(HOU J P,PETRINIC N,RUIZ C.A delamination criterion for laminated composites underlow-velocity impact[J].2001,61(14):2069-74.)来判断分层损伤,具体建立方式为:
(a)对于纤维拉伸和压缩,损伤初始判据为:
其中,为垂直纤维方向的名义应力张量,Sft、Sfc分别为单向板纵向拉伸强度、压缩强度;φft、φfc分别为纤维拉伸应力危险系数、纤维压缩应力危险系数;
(b)对于基体拉伸和压缩,损伤初始判据为:
图1定义了作用在断裂面上的各种应力,图中坐标系x1-x2-x3为单向板的自然坐标系,xl-xn-xt为断裂面的局部坐标系,两个坐标系下的xl轴重合;θ为σ2到最危险断面的旋转角度,θ的取值区间为[-90°,90°];σn(θ)和εn(θ)分别为潜在断裂面上的法向应力和法向应变(如图1所示断裂面xl-xn-xt局部坐标系下xn方向);τnl(θ)和εnl(θ)分别为潜在断裂面上平行于纤维方向剪应力和剪应变(如图1所示断裂面xl-xn-xt局部坐标系下xl xn平面内的剪应力和剪应变);τnt(θ)和εnt(θ)分别为潜在断裂面上垂直于纤维方向剪应力和剪应变(如图1所示断裂面xl-xn-xt局部坐标系下xt xn平面内的剪应力和剪应变);τnv(θ)为τnl(θ)与τnt(θ)的合力,ψ为τnv(θ)与τnt(θ)的夹角,则断裂面上的有效应力、应变分量计算如下:
分别为x2方向、x3方向的有效应力(图1中x1-x2-x3坐标系下x1表示纤维方向、x2表示垂直于纤维方向、x3表示层间平面外方向),/>分别为图1所示的x1-x2-x3坐标系下x2x3平面、x1x3平面、x1x2平面内的有效剪应力;
式中:φmt和φmc分别为基体拉伸与压缩应力危险系数;和/>为单向板垂直于纤维方向的拉伸与压缩强度,t、c分别表示拉伸和压缩,⊥、P分别表示垂直于纤维方向和平行于纤维方向;/>为单向板的面内剪切强度;/>为只有/>剪切应力作用下的失效应力;A表示断裂面;/>和/>为潜在断裂面上的法向应力对基体失效的促进参数,/>和/>为潜在断裂面上的法向应力对基体失效的抑制参数;对于碳纤维增强复合材料,将分别取值为0.27、0.35、0.27、0.3;
(c)分区二次插值法进行断裂角搜索
基体拉伸与压缩应力危险系数φmt和φmc是断裂面角度θ的一元函数,断裂面角度随应力状态的变化而变化,每一种应力状态下都有其最危险的潜在断裂面,可以通过一维搜索优化算法求得基体应力危险系数最大值的方法来求解断裂面角度θfp
为求解断裂面角度θfp,主要有以下三种方法:
第一种:Puck遍历法,在[-90°,90°]区间内每次以1°的增量逐角度计算基体应力危险系数,最终求得基体最大应力危险系数所在的角度为潜在断裂面角度。该方法用于有限元计算每个增量步内每个积分点上至少需要迭代180次,计算量极大。
第二种:基于黄金分割法和二次插值法提出的一种快速求解基体失效断裂面角度的算法,该方法每个增量步内每个积分点上仅需不超过6次计算即能够求解出基体失效断裂面角度。该方法大大提高了计算效率,但该方法对基体应力危险系数曲线含有多个峰值的情况,有可能得到的是局部最优解而不是全局最优解,因此可靠性较差。
第三种:在Wiegand的基础上提出了分区黄金分割算法,该方法首先分隔区间计算出含有基体应力危险系数极大值的搜索区间,然后在这些有极大值的区间内使用黄金分割法进行区间搜索求解区间内的局部最大值,最后所有局部最大值进行比较得到全局最大值,避免了求解的最优解是局部最优解的问题,但此方法在含有极大值的区间内还需多次采用黄金分割法进行分割点的寻找计算,直到求解得到满足精度的局部基体应力危险系数最大值所在的角度。
需要说明的是,θ为σ2到最危险断面的旋转角度,θ为断裂面角度,表达的是同一个意思:基体发生损伤会有一个损伤断裂面,将这个断裂面与σ2的夹角称为θ,这个角度在不同的应力状态下是不同的,通过求解应力危险系数来搜索求解,只有当φm取得最大值才是最危险的潜在断裂面,这时角度才是潜在断裂面角度,只有当φm取得最大值并且大于等于1,这时的角度才是断裂面角度θfp
为了解决Puck遍历法计算效率低,黄金分割法结合的二次插值法计算效率高但可靠性差,分区黄金分割法在含极大值区间内仍需多次进行黄金分割点的寻找计算等问题,本发明基于Schirmaier的研究结论:即基体损伤阈值的局部极大值不超过3个且极大值对应的角度间距不低于25°,提出了采用分区二次插值法来求解基体断裂面角度,该算法计算效率和计算精度高且可靠性良好。
该方法首先以10°的间隔将区间[-90°,90°]等分为18个子区间,分别计算18个间隔点上的损伤判断阈值,并对每相邻的3个间隔点上的损伤判断阈值进行比较,如果中间点的损伤判断阈值高于两侧间隔点上的数值,则可以确定这个区间内必定存在一个局部极大值,如图2所示:这三个点坐标分别为A(θ1m1))、B(θ2m2))和C(θ3m3)),然后在区间[θ12]内采用二次插值法进行搜索求得局部极大值,以此类推,最后将所有的局部极大值对比求得全局最大值;通过二次插值构造函数为:
对函数F(θ)求最大值,可以得到“潜在”最危险断裂角:
不同应力状态下损伤起始阈值随旋转角度θ的变化规律不同,且峰值点的位置和数量也不一样。为了说明本发明提出的分区二次插值法的优势,本发明以三种典型三维应力状态,即一个峰值(Case1)、两个峰值(Case2)、三个峰值(Case3)应力状态为例,如表1所示,采用Matlab软件进行编程计算,在相同硬件环境下进行测试,对比了Puck遍历法、分区黄金分割法和分区二次插值法在这三种应力状态下断裂面角度的计算精度和计算效率,如表2所示。在1°的精度要求下,采用Puck遍历法需要计算180个状态点,采用分区黄金分割法需要计算24个状态点,而分区二次插值法只需计算19个状态点。对于这三种应力状态,当精度均为1°时,采用三种方法所需的计算时间相接近,而本发明提出的方法具有略微优势;当精度为0.1°时,采用本发明提出的方法所需的计算时间大约为Puck遍历法的1/10,分区黄金分割法的1/2;当精度为0.01°时,采用本发明提出的方法所需计算时间大约为Puck遍历法的1/100,分区黄金分割法的1/2。因每个单元每个增量步都需计算这些状态点数,可见,当模型单元数量和增量步较大时,分区二次插值法能显著提高模型的计算效率且具有良好的计算精度。
表1三种典型的三维应力状态
表2三种典型三维应力状态下断裂角计算精度和计算效率对比
Notes:θ,N and T represent the fracture angle,the number of statepoints and the calculating time by different search methods respectively.
θ、N和T分别表示在不同搜索方法下的断裂角、状态点的个数和计算时间。
(d)分层损伤初始判据
其中,为层间平面外方向的有效正应力,/>和/>分别为图1所示的x1-x2-x3坐标系下x2x3平面、x1x3平面内的有效剪应力;
其中:φzt和φzc分别为拉伸与压缩分层应力危险系数;Szt为单向板的法向拉伸强度;S23和S13分别表示图1所示的x1-x2-x3坐标系下x2x3平面、x1x3平面内单向板的剪切强度;当由有效应力状态确定的φI的值等于损伤初始阈值1时,损伤开始萌生;若由当前有效应力状态得到的φI值超过先前加载的历史损伤阈值时,损伤进一步发展;其中,I={ft,fc,mt,mc,zt,zc},ft、fc、mt、mc、zt和zc分别是指纤维拉伸、纤维压缩、基体拉伸、基体压缩、拉伸分层、和压缩分层;
步骤(3):建立剪切非线性模型,具体建立方式为:
复合材料由于基体的塑性,通常展现出显著的剪切非线性或不可逆的塑性变形行为;本发明采用Soutis等(SHI Y,SWAIT T,SOUTIS C.Modelling damage evolution incomposite laminates subjected to low velocity impact[J].Composite Structures,2012,94(9):2902-13.)提出的半经验表达式来表征材料的剪切非线性,表达式如下所示:
其中,参数A可通过对单向复合材料层合板试件的偏轴拉伸或压缩试验进行线性回归分析获得,γij为图1所示的x1-x2-x3坐标系下x1x2平面、x2x3平面、x1x3平面内的工程剪应变,Gij为图1所示的x1-x2-x3坐标系下x1x2平面、x2x3平面、x1x3平面内的初始剪切模量,Sij为图1所示的x1-x2-x3坐标系下x1x2平面、x2x3平面、x1x3平面内的剪切强度,为图1所示的x1-x2-x3坐标系下x1x2平面、x2x3平面、x1x3平面内的有效剪应力;
剪切非线性损伤初始判断准则如下:当φij=1时,认为剪切损伤开始萌生;式中,ij代表某个平面,φij是判断准则,φij的含义就是这个三个平面内剪切应力是否导致损伤的判断准则;
步骤(4):建立损伤演化模型,具体建立方式为:
纤维损伤演化模型为:
其中,ft、fc分别代表纤维拉伸和纤维压缩,t、c分别指拉伸和压缩,dft(c)是指纤维拉伸/压缩损伤变量;ε1是指纤维方向应变;是指纤维发生初始损伤时(φft(c)=1)的拉伸/压缩应变;/>是指纤维发生最终失效时(dft(c)=1)的拉伸/压缩损伤应变;
基体损伤演化模型为:
其中,“<·>”为Macauley符号,当x∈R时,<x>=(x+|x|)/2;dmt(c)是指基体拉伸/压缩损伤变量;εeq,mt(c)是指损伤断裂面上的等效拉伸/压缩应变;是指断裂面上发生初始损伤时(φmt(c)=1)的等效拉伸/压缩应变;/>是指断裂面上发生最终失效时(dmt(c)=1)的等效拉伸/压缩应变;
分层损伤演化模型为:
其中,dzt(c)是指发生拉伸/压缩分层损伤变量;εeq,zt(c)是指导致拉伸/压缩分层损伤的等效应变;是指发生初始拉伸/压缩分层损伤时(φzt(c)=1)的等效应变;/>是指发生拉伸/压缩分层最终失效时(dzt(c)=1)的等效应变;
剪切损伤演化模型:
其中,dij是指图1所示的x1-x2-x3坐标系下x1x2平面、x2x3平面、x1x3平面内的剪切损伤变量;是指图1所示的x1-x2-x3坐标系下x1x2平面、x2x3平面、x1x3平面内发生剪切初始损伤时(φij=1)的应变;/>是指图1所示的x1-x2-x3坐标系下x1x2平面、x2x3平面、x1x3平面内发生剪切最终失效时(dij=1)的应变。
整个步骤二都在描述本构模型,损伤本构模型包含:损伤判断准则、损伤演化、本构模型。本构模型就是应力应变关系,是abaqus用户动态子程序模块调用子程序进行分析求得。
所述步骤三具体为:
步骤(1):开始当前增量步,读取前一时刻收敛状态量及当前增量步中应变增量,更新应变和有效应力;
步骤(2):根据有效应力代入到步骤二的分步骤(2)、(3)判断是否有损伤出现,如有损伤出现再通过步骤二的分步骤(4)更新损伤变量,再通过有效应力和损伤变量计算名义应力。
所述步骤四具体为:将步骤一建立的复合材料层合板有限元模型文件和步骤三建立的ABAQUS-VUMAT用户子程序联合,完成对复合材料层合板失效强度的预测。
所述步骤四具体为:先在ABAQUS软件里面建立复合材料层合板有限元模型,然后调用编写好的子程序进行应力应变分析,最后得到的应力位移曲线就是该模型的力学行为反应,得到的最大值为失效强度值。
本发明与现有技术相比,具有如下优点和有益效果:
本发明利用ABAQUS-VUMAT用户自定义子程序来数值实现所建立的三维损伤本构模型,该模型同时考虑了剪切非线性和损伤累积导致材料性能退化的影响,能准则预测复合材料的失效强度。对于ABAQUS内嵌的基于Hashin失效准则的二维弹性损伤本构,本发明基于更准确的Puck失效准则的包含剪切非线性效应的三维损伤本构,并且在求解断裂面角度的方法上,本发明提出的分区二次插值法,与目前提出的Puck遍历法、黄金分割法、分区黄金分割法相比,具有更高的计算效率和计算精度。
附图说明
图1为层合板子层应力状态和潜在断裂面的定义图。
图2为实施例的试件几何尺寸图。
图3为[0°/90°/±45°]2S带孔层合板拉伸应力-应变预测曲线。
图4为[±45°]4S带孔层合板拉伸应力-应变预测曲线。
具体实施方式
下面结合实施例及附图对本发明作进一步详细的描述,但本发明的实施方式不限于此。
在ABAQUS/CAE中建立Carlsson等(CARLSSON L A,ARONSSON C G,B CKLU NDJ.Notch sensitivity of thermoset and thermoplastic laminates loaded intension[J].Journal of Materials Science,1989,24(5):1670-82.)所实验的材料含缺口的AS4/3501-6复合材料层合板,试件铺层分别为[0°/90°/±45°]2S及[±45°]4S,单层厚度为0.125mm,承受轴向拉伸荷载。试件总长度C为280mm,其中B试验区长为200mm,宽W=76.2mm,两端A试验夹持端为40mm,中间圆孔直径D分别为6.35、12.7和25.4mm。试件几何尺寸如图2所示。
数值分析采用的AS4/3501-6材料属性及模型参数值均来自于已报道的同种材料的相应数值。材料单向板的材料属性见表3:
表3
利用层合板铺层角度沿厚度方向中平面对称布置,为减少计算量,有限元模型采用对称边界条件只模拟1/2厚层合板,每一层厚度方向只划分1个单元,模型包含8层C3D8R实体单元。由于孔边周围存在应力集中现象,且损伤也是由此开始延伸,因此需对孔边周围区域进行网格细化。受拉荷载均采用位移加载方式,左加载面施加固支约束,右端自由端面外设置一个参考点,将位移荷载施加在参考点上,建立参考点与自由端面之间加载方向位移一致性约束。
利用ABAQUS/EXPLICIT计算模拟复合材料层合板的拉伸失效过程,利用用户子程序VUMAT读取当前应变增量,更新应变和有效应力,根据有效应力判别单元是否进入损伤,进入损伤阶段则根据损伤演化模型计算损伤变量,从而得到名义应力,最后可从ABAQUS-VISUALIZATION模块得到模型的反力-位移曲线。
图3、4为分别采用本发明开发的本构模型和Abaqus v6.14自带的纤维增强复合材料弹性损伤模型对铺层角度为[0°/90°/±45°]2S和[±45°]4S的3种不同孔径的带孔层合板进行了数值分析,得到的拉伸应力-应变预测曲线。曲线1表示Present model D=6.35mm、曲线2表示Present model D=12.7mm、曲线3表示Present model D=25.4mm、曲线4表示Abaqus built-in model D=6.35mm、曲线5表示Abaqus built-in model D=12.7mm、曲线6表示Abaqus built-in model D=25.4mm。表3为预测失效强度与试验值的比较,本发明的预测值比Abaqus v6.14自带的纤维增强复合材料弹性损伤模型的预测值误差要小的多,最大误差为9.8%,可见,本发明提出的损伤本构模型能理想的预测其失效强度。
表4为[0°/90°/±45°]2S及[±45°]4S带孔层合板拉伸破坏强度预测值与试验结果对比。
表4
本发明是在ABAQUS软件的基础上进行的用户子程序的开发,提出的三维损伤本构模型同时考虑了剪切非线性效应和损伤累积导致的材料性能退化的影响,能理想的预测复合材料层合板的失效强度,为深入阐明复合材料结构的损伤失效特性,提升轻量化强度设计水平提供了技术支撑。
上述实施例为本发明较佳的实施方式,但本发明的实施方式并不受上述实施例的限制,其他的任何未背离本发明的精神实质与原理下所作的改变、修饰、替代、组合、简化,均应为等效的置换方式,都包含在本发明的保护范围之内。

Claims (8)

1.一种预测复合材料层合板失效强度的方法,其特征在于,包含以下顺序的步骤:
步骤一、建立复合材料层合板有限元模型;
步骤二、建立复合材料损伤本构模型;具体为:
步骤(1):建立含损伤的复合材料层合板本构关系;
复合材料应力-应变本构方程:σ=C(d):εe
其中,符号“:”表示对两个张量指标的缩并计算;σ是有效应力张量;是名义应力张量;εe是弹性应变张量;e表示弹性;C(d)是含损伤单向复合材料层合板的四阶刚度张量;C是未损伤单向复合材料层合板的四阶线弹性刚度张量;d是一维向量(d1,d2,d3,d23,d13,d12),其中d1、d2、d3分别为纤维方向纤维损伤的损伤变量、平面内垂直于纤维方向基体损伤的损伤变量、层间平面外方向分层损伤的损伤变量;d12、d23、d13分别为12、23、13平面内的剪切损伤变量;定义坐标系x1-x2-x3为单向板的自然坐标系,xl-xn-xt为断裂面的局部坐标系,两个坐标系下的xl轴重合;12、23、13平面分别对应坐标系下的x1x2平面、x2x3平面、x1x3平面;
把损伤变量引入刚度矩阵,使刚度随着损伤的发展而逐渐变弱,即:
C(d)=M-1(d):C:MT,-1(d);
其中,M-1(d)为M(d)的逆矩阵,MT,-1(d)为M(d)转置矩阵的逆矩阵;M(d)为损伤因子张量,其损伤主轴系下矩阵形式表示如下:
复合材料主坐标系中单层板的三维正交各项异性损伤本构模型如下:
所述复合材料主坐标系为单向板的自然坐标系x1-x2-x3
其中:
其中,σ1、σ2和σ3分别为纤维方向、垂直于纤维方向、层间平面外方向的名义正应力;τ23、τ12和τ13分别为x1-x2-x3坐标系下x1x2平面、x2x3平面、x1x3平面内的剪应力;ε1、ε2和ε3分别为纤维方向、垂直于纤维方向、层间平面外方向的工程正应变;γ23、γ13和γ12为x1-x2-x3坐标系下x1x2平面、x2x3平面、x1x3平面内的工程剪应变;E1、E2、E3分别为纤维方向、垂直于纤维方向、层间平面外方向的未损伤单向复合材料单层的弹性模量,G23、G13、G12分别为x1x2平面、x2x3平面、x1x3平面内未损伤单向复合材料单层的剪切模量,ν12、ν13、ν23分别为纤维方向与垂直于纤维方向的、纤维方向与层间平面外方向的、垂直于纤维方向与层间平面外方向的泊松比,ν21、ν31、ν32分别为垂直于纤维方向与纤维方向的、层间平面外方向与纤维方向的、层间平面外方向与垂直于纤维方向的泊松比,满足关系式,/>
步骤(2):建立三维Puck失效准则来判断纤维和基体损伤,三维Hou准则来判断分层损伤,具体建立方式为:
(a)对于纤维拉伸和压缩,损伤初始判据为:
其中,为垂直纤维方向的名义应力张量,Sft、Sfc分别为单向板纵向拉伸强度、压缩强度;φft、φfc分别为纤维拉伸应力危险系数、纤维压缩应力危险系数;
(b)对于基体拉伸和压缩,损伤初始判据为:
坐标系x1-x2-x3为单向板的自然坐标系,xl-xn-xt为断裂面的局部坐标系,两个坐标系下的xl轴重合;θ为σ2到最危险断面的旋转角度,θ的取值区间为[-90°,90°];σn(θ)和εn(θ)分别为潜在断裂面上的法向应力和法向应变;τnl(θ)和εnl(θ)分别为潜在断裂面上平行于纤维方向剪应力和剪应变;τnt(θ)和εnt(θ)分别为潜在断裂面上垂直于纤维方向剪应力和剪应变;τnv(θ)为τnl(θ)与τnt(θ)的合力,ψ为τnv(θ)与τnt(θ)的夹角,则断裂面上的有效应力、应变分量计算如下:
分别为x2方向、x3方向的有效应力,/>分别为x1-x2-x3坐标系下x2x3平面、x1x3平面、x1x2平面内的有效剪应力;
式中:φmt和φmc分别为基体拉伸与压缩应力危险系数;和/>为单向板垂直于纤维方向的拉伸与压缩强度,t、c分别表示拉伸和压缩,⊥、P分别表示垂直于纤维方向和平行于纤维方向;/>为单向板的面内剪切强度;/>为只有/>剪切应力作用下的失效应力;A表示断裂面;/>和/>为潜在断裂面上的法向应力对基体失效的促进参数,/>和/>为潜在断裂面上的法向应力对基体失效的抑制参数;对于碳纤维增强复合材料,将/>分别取值为0.27、0.35、0.27、0.3;
(c)分区二次插值法进行断裂角搜索
基体拉伸与压缩应力危险系数φmt和φmc是断裂面角度θ的一元函数,断裂面角度随应力状态的变化而变化,每一种应力状态下都有其最危险的潜在断裂面,通过一维搜索优化算法求得基体应力危险系数最大值的方法来求解断裂面角度θfp
(d)分层损伤初始判据
其中,为层间平面外方向的有效正应力,/>和/>分别为x1-x2-x3坐标系下x2x3平面、x1x3平面内的有效剪应力;
其中:φzt和φzc分别为拉伸与压缩分层应力危险系数;Szt为单向板的法向拉伸强度;S23和S13分别表示x1-x2-x3坐标系下x2x3平面、x1x3平面内单向板的剪切强度;当由有效应力状态确定的φI的值等于损伤初始阈值1时,损伤开始萌生;若由当前有效应力状态得到的φI值超过先前加载的历史损伤阈值时,损伤进一步发展;其中,I={ft,fc,mt,mc,zt,zc},ft、fc、mt、mc、zt和zc分别是指纤维拉伸、纤维压缩、基体拉伸、基体压缩、拉伸分层、和压缩分层;
步骤(3):建立剪切非线性模型,具体建立方式为:
复合材料由于基体的塑性,通常展现出显著的剪切非线性或不可逆的塑性变形行为;采用半经验表达式来表征材料的剪切非线性,表达式如下所示:
其中,参数A可通过对单向复合材料层合板试件的偏轴拉伸或压缩试验进行线性回归分析获得,γij为x1-x2-x3坐标系下x1x2平面、x2x3平面、x1x3平面内的工程剪应变,Gij为x1-x2-x3坐标系下x1x2平面、x2x3平面、x1x3平面内的初始剪切模量,Sij为x1-x2-x3坐标系下x1x2平面、x2x3平面、x1x3平面内的剪切强度,为x1-x2-x3坐标系下x1x2平面、x2x3平面、x1x3平面内的有效剪应力;
剪切非线性损伤初始判断准则如下:当φij=1时,认为剪切损伤开始萌生;式中,ij代表某个平面,φij是判断准则,φij的含义就是这个三个平面内剪切应力是否导致损伤的判断准则;
步骤(4):建立损伤演化模型,具体建立方式为:
纤维损伤演化模型为:
其中,ft、fc分别代表纤维拉伸和纤维压缩,t、c分别指拉伸和压缩,dft(c)是指纤维拉伸/压缩损伤变量;ε1是指纤维方向应变;是指纤维发生初始损伤时的拉伸/压缩应变;是指纤维发生最终失效时的拉伸/压缩损伤应变;
基体损伤演化模型为:
其中,“<·>”为Macauley符号,当x∈R时,<x>=(x+|x|)/2;dmt(c)是指基体拉伸/压缩损伤变量;εeq,mt(c)是指损伤断裂面上的等效拉伸/压缩应变;是指断裂面上发生初始损伤时的等效拉伸/压缩应变;/>是指断裂面上发生最终失效时的等效拉伸/压缩应变;
分层损伤演化模型为:
其中,dzt(c)是指发生拉伸/压缩分层损伤变量;εeq,zt(c)是指导致拉伸/压缩分层损伤的等效应变;是指发生初始拉伸/压缩分层损伤时的等效应变;/>是指发生拉伸/压缩分层最终失效时的等效应变;
剪切损伤演化模型:
其中,dij是指x1-x2-x3坐标系下x1x2平面、x2x3平面、x1x3平面内的剪切损伤变量;是指x1-x2-x3坐标系下x1x2平面、x2x3平面、x1x3平面内发生剪切初始损伤时的应变;/>是指x1-x2-x3坐标系下x1x2平面、x2x3平面、x1x3平面内发生剪切最终失效时的应变;
步骤三、基于ABAQUS-VUMAT有限元用户动态子程序模块,使用FORTRAN语言编写用户自定义子程序实现提出的损伤本构模型,以求解应力、应变和损伤;
步骤四、对有限元模型进行计算,预测复合材料层合板的失效强度。
2.根据权利要求1所述预测复合材料层合板失效强度的方法,其特征在于,所述步骤一,具体如下:
复合材料层合板的铺层角度沿厚度方向中平面对称布置,有限元分析沿厚度方向采用对称边界条件只模拟1/2厚的复合材料层合板,每一层厚度方向只划分1个单元,复合材料层合板有限元模型包含8层C3D8R实体单元;对孔边周围区域进行网格细化;建立参考点与自由端面之间加载方向位移一致性约束条件:受拉荷载均采用位移加载方式,左加载面施加固支约束,右端自由端面外设置一个参考点,然后把参考点和端面进行绑定,在Abaqus/CAE模块中,采用creat constraint方法建立coupling耦合约束方程,此时,将位移荷载施加在参考点上,同时只要输出参考点上的位移和反力即U1和RF1,就能够获得加载端面上的位移与反力。
3.根据权利要求1所述预测复合材料层合板失效强度的方法,其特征在于,所述求解断裂面角度θfp,具体为:
在[-90°,90°]区间内每次以1°的增量逐角度计算基体应力危险系数,最终求得基体最大应力危险系数所在的角度为潜在断裂面角度。
4.根据权利要求1所述预测复合材料层合板失效强度的方法,其特征在于,所述求解断裂面角度θfp,具体为:
每个增量步内每个积分点上仅需不超过6次计算即能够求解出基体失效断裂面角度。
5.根据权利要求1所述预测复合材料层合板失效强度的方法,其特征在于,所述求解断裂面角度θfp,具体为:
首先分隔区间计算出含有基体应力危险系数极大值的搜索区间,然后在这些有极大值的区间内使用黄金分割法进行区间搜索求解区间内的局部最大值,最后所有局部最大值进行比较得到全局最大值,直到求解得到满足精度的局部基体应力危险系数最大值所在的角度;
该方法首先以10°的间隔将区间[-90°,90°]等分为18个子区间,分别计算18个间隔点上的损伤判断阈值,并对每相邻的3个间隔点上的损伤判断阈值进行比较,如果中间点的损伤判断阈值高于两侧间隔点上的数值,则可以确定这个区间内必定存在一个局部极大值:这三个点坐标分别为A(θ1m1))、B(θ2m2))和C(θ3m3)),然后在区间[θ12]内采用二次插值法进行搜索求得局部极大值,以此类推,最后将所有的局部极大值对比求得全局最大值;通过二次插值构造函数为:
对函数F(θ)求最大值,可以得到“潜在”最危险断裂角:
不同应力状态下损伤起始阈值随旋转角度θ的变化规律不同,且峰值点的位置和数量也不一样。
6.根据权利要求1所述预测复合材料层合板失效强度的方法,其特征在于,所述步骤三具体为:
步骤(1):开始当前增量步,读取前一时刻收敛状态量及当前增量步中应变增量,更新应变和有效应力;
步骤(2):根据有效应力代入到步骤二的分步骤(2)、(3)判断是否有损伤出现,如有损伤出现再通过步骤二的分步骤(4)更新损伤变量,再通过有效应力和损伤变量计算名义应力。
7.根据权利要求1所述预测复合材料层合板失效强度的方法,其特征在于,所述步骤四具体为:将步骤一建立的复合材料层合板有限元模型文件和步骤三建立的ABAQUS-VUMAT用户子程序联合,完成对复合材料层合板失效强度的预测。
8.根据权利要求7所述预测复合材料层合板失效强度的方法,其特征在于,所述步骤四具体为:先在ABAQUS软件里面建立复合材料层合板有限元模型,然后调用编写好的子程序进行应力应变分析,最后得到的应力位移曲线就是该模型的力学行为反应,得到的最大值为失效强度值。
CN201910961138.0A 2019-10-11 2019-10-11 一种预测复合材料层合板失效强度的方法 Active CN111368389B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201910961138.0A CN111368389B (zh) 2019-10-11 2019-10-11 一种预测复合材料层合板失效强度的方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201910961138.0A CN111368389B (zh) 2019-10-11 2019-10-11 一种预测复合材料层合板失效强度的方法

Publications (2)

Publication Number Publication Date
CN111368389A CN111368389A (zh) 2020-07-03
CN111368389B true CN111368389B (zh) 2023-11-07

Family

ID=71210393

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201910961138.0A Active CN111368389B (zh) 2019-10-11 2019-10-11 一种预测复合材料层合板失效强度的方法

Country Status (1)

Country Link
CN (1) CN111368389B (zh)

Families Citing this family (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112052614A (zh) * 2020-08-31 2020-12-08 暨南大学 山棕纤维床垫的横观各向同性材料有限元模型建立方法
CN112100744B (zh) * 2020-09-03 2023-06-06 上海交通大学 叠层板材集群孔试验系统
CN112329205B (zh) * 2020-10-12 2022-04-29 湖北航天技术研究院总体设计所 一种复合材料结构低速冲击损伤确定方法及装置
CN112487645B (zh) * 2020-12-01 2022-05-10 中国科学院软件研究所 一种统一的各向同性、各向异性虚拟材料的能量建模方法和装置
CN112906263B (zh) * 2021-01-28 2022-06-28 天津大学 含制孔分层损伤的复合材料层合板强度预测方法
CN113609688B (zh) * 2021-08-09 2024-03-29 大连理工大学 复合材料细观切削仿真中纤维剪断、弯断失效的精确判定方法
CN113642169B (zh) * 2021-08-09 2024-03-15 大连理工大学 一种适用于横观各向同性碳纤维的多模式失效判定方法
CN114611389A (zh) * 2022-03-04 2022-06-10 北京航空航天大学 一种基于人工智能的复合材料失效高效模拟方法

Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN105740566A (zh) * 2016-02-18 2016-07-06 浙江大学 一种预测层状复合材料层内损伤和层间分层的有限元方法
CN106503292A (zh) * 2016-09-20 2017-03-15 浙江大学 预测低速冲击下复合材料层合板渐进失效的有限元方法

Patent Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN105740566A (zh) * 2016-02-18 2016-07-06 浙江大学 一种预测层状复合材料层内损伤和层间分层的有限元方法
CN106503292A (zh) * 2016-09-20 2017-03-15 浙江大学 预测低速冲击下复合材料层合板渐进失效的有限元方法

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
Glass/Epoxy复合材料叠层板面内剪切连续损伤模型;庞宝君等;《复合材料学报》;20130514(第05期);第222-207页 *
基于弹塑性损伤本构模型的复合材料层合板破坏荷载预测;陈静芬;《复合材料学报》;20170430;第547-552页 *
考虑剪切非线性影响的复合材料连续损伤模型及损伤参数识别;刘伟先等;《复合材料学报》;20130514(第06期);第222-224页 *

Also Published As

Publication number Publication date
CN111368389A (zh) 2020-07-03

Similar Documents

Publication Publication Date Title
CN111368389B (zh) 一种预测复合材料层合板失效强度的方法
Kam et al. Predictions of deflection and first-ply failure load of thin laminated composite plates via the finite element approach
Qu et al. Three-dimensional free and transient vibration analysis of composite laminated and sandwich rectangular parallelepipeds: Beams, plates and solids
Cazzani et al. A refined assumed strain finite element model for statics and dynamics of laminated plates
Thai et al. A layerwise C0-type higher order shear deformation theory for laminated composite and sandwich plates
Beklemysheva et al. Numerical simulation of the failure of composite materials by using the grid-characteristic method
Biswas et al. Comparative study on transient response analysis of hybrid laminated composite plates with experimental verification
Wu et al. Effects of higher-order global–local shear deformations on bending, vibration and buckling of multilayered plates
Noh et al. Failure analysis of glass/epoxy and graphite/epoxy laminates due to the effect of variation in lamination scheme and angle of fibre orientation
Li et al. A novel hybrid auxetic honeycomb with enhanced load-bearing and energy absorption properties
Jun et al. Free vibration analysis of third-order shear deformable composite beams using dynamic stiffness method
CN112926244A (zh) 一种复合材料层合板开孔件极限载荷确定方法
Zhang et al. Applications of VAM-based homogenization model in free and forced vibrations of sandwich plates with bowtie-shaped auxetic core
Karkkainen et al. A direct micromechanical approach toward the development of quadratic stress gradient failure criteria for textile composites
Nayak et al. Transient response of initially stressed composite sandwich plates
Semmani et al. Analysis and optimization of composite kagome grid panels subjected to the low velocity impact
Al-Masri et al. Buckling solutions of clamped-pinned anisotropic laminated composite columns under axial compression using bifurcation approach and finite elements
Nader et al. Probabilistic finite element analysis of modified ASTM D3039 tension test for marine grade polymer matrix composites
Li et al. Application of the Burgers model and elasticity theory to obtain viscoelastic solutions for flexural creep of layered beams
Liu et al. Failure progression and mesh sensitivity analyses by the plate element-failure method
Diveyev et al. Dynamic properties identification for laminated plates
Naghipour et al. Numerical simulation of composite plates to be used for optimization of mobile bridge deck
Lestari et al. Dynamic characteristics and effective stiffness properties of honeycomb composite sandwich structures for highway bridge applications
Sadr et al. Fundamental frequency optimization of angle-ply laminated plates using Elitist-Genetic algorithm and finite strip method
Mahdi et al. Influence of Fiber Orientation on the Natural Frequencies of Laminated Composite Beams

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant