复合材料失效模拟，利用abaqus用户子程序模拟复合材料层合板失效
1136 Z.-M. HUANG Unfortunately, the nonlinear mechanical behavior of a composite is not well understood in the current literature. Recently, the present author proposed a general constitutive theory, the bridging model [3, 4], for composites. By combining the bridging model with the classical lamination theory, a general constitutive relationship was established for any fibrous laminate, and was summarized in References [1, 2]. a recently completed worldwide failure exercise [5-8] has indicated that the bridging model constitutive theory had two unique features The first feature was that it was the only model in the exercise that could be used to calculate the thermal stresses in the fiber and matrix materials due to a thermal load(temperature variation) applied on the composite (see Reference [6],p. 450). The second feature lay in the fact that the model is consistent in that the laminate nonlinear constitutive equations automatically deteriorate to those of the isotropic matrix material when the fiber content becomes zero or when the fiber becomes the same as that of the matrix. it is noted that the other nonlinear constitutive models used in the exercise [9-15] described the composite constitutive equations in a way somewhat similar to Hooke's law. Namely, a shear stress would not cause an extensional strain whereas a normal stress had no contribution to a shear strain according to a classical plasticity theory such as Prandtl-reus theory for isotropic materials, however, a normal (or shear)stress will generate a shear (or extensional) plastic strain if the material is subjected to a combination of the shear and normal stresses The finite element method has been realized to be one of the most powerful tools in improving design quality and reducing development time for engineering structures Several well-known commercial finite element codes, such as ABAQUs, are capable of analyzing the nonlinear response and strength behavior of a very complicated structure However, their applicability to the structures made of laminated composites is restricted due to a shortage of an efficient inelastic constitutive module for composites in the material libraries of those codes. The purpose of this paper is to incorporate the bridging model based nonlinear constitutive relationship for multidirectional tape laminates into ABAQUS through programming a user subroutine UGENS, as though a new module is added into the aBaQUs material library. In this way, any structure with tape laminates involved, either a sole laminated structure or a structure consisting of laminates and other materials for which constitutive modules are provided in the ABAQUS library, which is subjected to arbitrary load condition can be Fe analyzed readily. In the following sections, the routine development, original subroutine code instructions for providing input data, and an illustration example are presented The original computer code and an input data file are provided in Appendices a and b f the paper SUMMARY OF CONSTITUTIVE RELATIONS For convenience of expression, necessary equations are summarized below more details can refer to References [1-4]. Within the scope of the classical laminate theory, the nonzero stress and strain increments, (dolk and (delk, of the kth lamina of a laminate in the laminate coordinate system(Figure 1)are given by dok=(T])k(S)(71) Blk dT =ICOkldajk-iBjk dT, Analysis of laminate Structures by ABAQUS 1137 where Dark 0 deu+ zk t 2k-l dr 0 少+-4+之 0 dk0,2d。0,+(k+2k-1)dk0/ (2.1) B)k=(B1),(B2),(B3)=(T1)k(Sk){l)k (2.2) yy,and de and dko whry, and dkv are the laminate in-plane strain and curvature Xx2 increments. " L''refers to the laminate coordinate system shown in Figure 1. Zk and zk are the z coordinates of the top and the bottom surfaces of the kth lamina in the laminate system. [Tl is a coordinate transformation matrix between the laminate coordinate system and the material principal coordinate system of the lamina, i.e,(x1, x2)in Figure 2(b). dT is a temperature increment sustained by the laminate. By the bridging Model, the current compliance matrix and thermal expansion coefficient of the lamina in the material principal coordinate system are derived as [3] [S]=(S勹+Vm[S"IAB 小=va+m")+Vm(S"]-S勹)b / dN.dM ← Vector of a moment Figure 1. Laminate coordinate system together with applied force and moment increments 6=01<0 →xd← 6=62=0 y▲ 6=02>0 X1 11 (a) (b) Figure 2.(a) Laminate system with ply-angle,(b) lamina in laminate system, (c)lamina in material principal coordinate system 1138 Z.-M. HUANG where [b]=(v/n+VmlAd, V and Vm denote the volume fractions of the fiber and matrix materials, [S]and [s are the current compliance matrices of the fiber and matrix materials,a and a are the current thermal expansion coefficients of the fiber and matrix materials respectively, and [n is a unit matrix Internal stress increments in the matrix and fiber are correlated with the stress increments applied to the lamina, do=do 1, do22, do123, in the material principal coordinates via d}=[B]{do}+仂b/}dT, (5.1) Ido=[[BIdo)+b"dT, (5.2 fb/+Vm{b"}={0 (5.3) {b"}=([AI[B)([S-[S")(mn-{ (54 where the stress increments in the lamina principal system are obtained from the ones in the laminate system, given by Equation(1), through a coordinate transformation dokk=(Toldo [TIs is another coordinate transformation matrix different from [c. The bridging matrix [a together with [ B are expressed as 11 12a16 b1 b12 b 0 2226 and B 0 622 b (7 00b Explicit expressions of aii and bii can be found in, e.g., Reference 3] FAILURE CRITERIA The total stresses in the fiber and matrix materials as well as on the lamina are updated in the following way Lo"k= o k +ldo"lk, (8.1) JJ=(o/k+(do/j (8.2 an Lokk= lok+] 8.3) As both the internal stresses and the lamina stresses are already known, either a micromechanical or a macromechanical failure criterion can be employed to detect the Analysis of laminate Structures by ABAQUS 1139 lamina failure. In this paper, three failure criteria are incorporated. The first is the modified maximum normal stress criterion [3, 4], the second is the Tsai-Wu criterion, and the third is the hashin -Rotem criterion With the first criterion used, the lamina is considered to have failed as long as the fiber or matrix attains failure, i.e. if the fiber or matrix stresses satisfy either of inequalities(9.1)and(9.2) O (9.1,9.2) where when g(2)<0 )9+((2) 1/q when q O11+ +y(o11-02)2+4(o12) (10.1 O11+σ (011-022)2+4(012) (10.2) In the above, q is a power index, which is introduced to account for a reduction in the material load carrying capacity due to a bi-axial tension in comparison with a uniaxial tension. ou and ou.c are the ultimate tensile and compressive strengths of the fiber or matrix material under uniaxial tension. For the fibers, ou and ou. c are measured along the fiber axial direction. When the power-index q= oo, ( 9. 1) together with(9.2) is equivalent to the classical maximum normal stress criterion In contrast to the maximum normal stress criterion in which the fiber and matrix ultimate strengths are required, both the Hashin-Rotem and the Tsai-Wu criteria rely on five individual strengths of the lamina. Hence, failure experiments on the lamina must be performed in order to obtain those strength parameters The Hashin-Rotem criterion is expressed as [16 2 max where X, X, y, r, and S are longitudinal tension, longitudinal compression, transverse tension, transverse compression, and in-plane shear strengths of the unidirectional lamina, respectively. Further, the Tsai-Wu criterion is given by F1(O11)2+F2(022)2+F301022+F4(012)2+F501+F6022≥1, (12) where [17 Ⅹ-X rr, F2= YY F1F2,F4=,F5 an XX 1140 Z.-M. HUANG It should be pointed out once more that the stresses used in Equations(11) and(12) are those on the lamina, i.e., the stresses determined by Equation(8.3), whereas the stresses used in Equations(9) and (10) are determined either by Equation(8.1)or by Equation(8.2) FE SIMULATION Let us assume that the composite laminate be discretized by three-dimensional (3D) shell elements, and the quantities at every integration point of the element be treated incrementally. Suppose that the external forces and moments of a unit length are denoted by dNxx, dNyy, dNxy, dMxx, dMyy, and dMxy, as shown in Figure 1. They are correlated with the middle surface strain and curvature increments and equivalent thermal load increments through dNxx +dsi 21 212 213 211 212 2 d dNu tds ol2 022 023 012 022 023 d dNv+dss 213 023033 233 deo rv (13 dMx+d2|ag"q"?ng巛|dk dMv+dQ2 212 222023 2120229 dMx+ ds3 Q13Q23Q3Q1923Q lI 2d =2>(C)(k-1)2=2∑(c)(-3).9=3∑(C)(-n (14 d2=∑(B 1)dT,dg2=∑(B)(-=1)d (15) Equation(13)can be more conveniently expressed as dF}=[Q]{d8}, (16) where IdF)=(dPi+ldQ) (17) The generalized sectional applied forces and strains of the laminate at a given point are updated through {P}={P}+{dP},{8}={8}+{d8} (18.1),(18.2 In Abaqus, a 3D shell element such as 4-nodal s4R has 3 translational and 3 rotational degrees of freedom (DOFs). When used to discretize a composite structure, the shell element serves as a local laminate. at any integration point of the element the best way to Analysis of laminate Structures by ABAQUS 1141 incorporate the bridging model based laminate nonlinear constitutive relationship into ABAQUS is through programming a user subroutine UGENS. In the UGENS, arrays FORCE(6 and stran( are passed in as the generalized forces and strains, and will be updated at the end of the routine. The FOrCe array should exclude the contribution from a temperature variation, and hence should be updated through Equation(18.1) On the other hand, the Stran array should contain the contribution of thermal strains and thus updated by(18.2). In addition, the internal stresses in the fiber and matrix of each ply of the element, om and o, are passed in through a state variable array state (NSTATV), where NSTaTV is the length of the STATEV, and will be updated at the end of the routine. An array dStRan6, denoting the generalized strain increments (d8) is passed in as given information The main purpose of the UGENS is to form the elemental current material stiffness matrix [@] in Equation(13), and is stored in an array ddNddE(6, 6)at the end of the routine. This array will be used by abaqus to define the elemental stiffness matrix in the FE simulation. Another purpose is to update the internal stresses in the fiber and matrix materials of all the plies of the element. s the middle strain and curvature increments have already been passed in through DSTRaN(6), the incremental stresses in the fiber and matrix of each ply are calculated straightforwardly using Equations (1),(6), and (5) The third purpose of the UGENS is to check how many plies of the element have failed at this load increment. If the koth ply has failed, the material stiffness elements in Equation(13) should be defined using the following formulas ∑(C)(zk-z kg kol (19) ∑(c) ∑(Ch)( ∑(B1)(k dt. d 21=2∑()(a zk_DdT (20) instead of using Equations(14)and(15), where ko refers to all the failed plies Only when an element attains an ultimate failure, will the ABaQUs evaluation be stopped. In this work, the ultimate failure is defined as a pre-specified ply failure (in a input date file). In general, if the laminate(element) is subjected to an in-plane load the last-ply failure is the case of ultimate failure. If, however, the laminate is subjected to a bending load, the ultimate failure corresponds to an intermediate-ply failure(e.g. the third ply out of a 12-ply laminate)[18, 19] COORDINATE TRANSFORMATION In reality, a composite structure may consist of many different laminates Each laminate, e. g, that in Figure 1, has its own coordinate system and can be discretized into an arbitrary number of 3D shell elements. By convention, the local coordinate system 1142 Z.-M. HUANG of an element is specified as the corresponding laminate system (x,v, z), which is not necessarily coincident with the global one(X,Y, z of the structure. Thus, a coordinate transformation between the elemental local and the structural global coordinate systems is necessary. The local coordinate system of each element can be defined using "ORIENTATION option provided in ABAQUS if necessary. Then, the lamination angle of each ply of the element can be specified and the material stiffness matrix of the element (i.e, laminate)is obtained However, with the UGENS subroutine, the quantities passed in such as dstRan(6) FORCE(6, and STRaN(6) are all expressed in the global coordinate system, whereas ddNdDE(6, 6)to be returned back should also be given in the global system. This means that at the beginning of the UGEns the quantities passed in should be transformed into those in the local (i.e, laminate) system, and before returning the routine all the relevant quantities should be transformed back into the global system. Fortunately, ABAQUS has already provided an array basis(3, 3) into the UGENS routine to accomplish the required coordinate transformations According to a coordinate transformation rule [20], one has {8}=[]168}°,{P}=[L]P (21) where the superscriptsL''and"G''refer to the local and global coordinate systems, respectively n is the directional cosine matrix of the local coordinates in relation to the global one iven by the array basis(3, 3). Because [L=[L], one further has g=[L][gy[L,{8°=[L]{8},{P=[L]{P}2 (22) USER SUBROUTINE UGENS The main subroutine UGENS was programmed in FORTRAN77 language, consistin of 20 other subroutines and a total number of 1273 lines including remark and illustration ines. In the following, the function of each routine together with some major variables is briefly described. In the program, the characteristic of a variable is implicitly defined Namely, those starting with“I,J”,K”,L”’,M”,andN"' are integers or integer arrays, whereas the variables starting with all the other characters are real numbers or arrays. Original codes of the whole program are given in the Appendix a at the end of this pap Routine ugens This is the main subroutine, a fixed format of which has been provided in an ABAQUS user's manual. One can only meet one's own purpose by suitably specifying the relevant parameters and arrays in the format Analysis of laminate Structures by ABAQUS 1143 UGENS provides two arrays for the user to define a problem. One is State, a state variable array, and another is PROPS, which has been used to store original parameters necessary for running the UGENS. The size of the ProPs is specified by an integer variable, NPROPS, whereas the size of the state is defined by nstatv. both of these two parameters, dependent on the maximum possible layers of a laminate to be analyzed, should be defined in an input file(e.g,b inp,) for the ABAQUs according to the following format * SHELL GENERAL SECTION. ELSET=- PLATE. USER VARIABLES MSTATV, PROPERTIES NPROPS ome other controlling parameters are illustrated below NL= the number of layers in an element, which cannot be greater than 200(otherwise, proper changes must be made for the relevant working arrays MSEG= the number of linear segments to form a stress-strain curve of the matrix, which must be less than or equal to 20. For instance, if the matrix stress-strain curve is bilinear. MSEG=2. It is noted that all of the matrix stress-strain curves used in a specific analysis (i.e, at different temperatures, both at tension and compression) must consist of the same MSEG Segments NEM=MSEG, if the matrix material is temperature independent the number of different temperatures at which different matrix stress-strain curves are provided(25), if mechanical behavior of the matrix varies with temperature NEF=7, if the fiber material is temperature independent; the number of different temperatures at which different fiber properties are provided(=25), if mechanical behavior of the fiber is a function of temperature NPROPS=11*NL+13+NI+N2+N3+N4, where (refer to the next section for more detail) NI=7, if the fiber properties are temperature independent, or 6 NEF+5, if the fiber properties are dependent on temperature N2=4 NEM +3, if the matrix properties are temperature independent or (4 MSEG+5)*NEM, if the matrix properties are dependent on temperature. N3=0, if the maximum normal stress criterion is used to detect composite failure, or 5, if the Tsai-Wu or hashin failure criterion is used to detect composite failure N4=3, if the laminate is subjected to a 3-point or 4-point bending load, or 0, if the laminate is subjected to other kind of loads For instance, if the fiber properties are temperature independent and those of the matrix are dependent on temperature, Nl=20, MSeG=4, NEM=2, the Tsai-Wu criterion is used, and the laminate is subjected to an in-plane load, one has NPrOPS=287 NSTATV=7*NL +14(e.g, if NL=20, one specifies NSTATV= 154) Definitions for the other variables can refer to the aBaQus user's manual