Journal Search Engine
Search Advanced Search Adode Reader(link)
Download PDF Export Citaion korean bibliography PMC previewer
ISSN : 2093-5145(Print)
ISSN : 2288-0232(Online)
Journal of the Korean Society for Advanced Composite Structures Vol.5 No.4 pp.11-17

A Progressive Failure Analysis Procedure for Composite Laminates Ⅱ - Nonlinear Predictive Finite Element Analysis

Gyu-Sei Yi
Professor, Department of Civil Engineering, Sun Moon University, ChungNam, Korea
Corresponding Author : Yi, Gyu-Sei Department of Civil Engineering, Sun Moon University, ChungNam, Korea Tel: +82-41-530-2324, Fax: +82-41-530-2926,
December 22, 2014 December 29, 2014 December 30, 2014


A progressive failure analysis procedure for composite laminates is completed in here. An anisotropic plastic constitutive model for fiber-reinforced composite material is implemented into computer program for a predictive analysis procedure of composite laminates. Also, in order to describe material behavior beyond the initial yield, the anisotropic work-hardening model and subsequent yield surface are implemented into a computer code, which is Predictive Analysis for Composite Structures (PACS). The accuracy and efficiency of the anisotropic plastic constitutive model and the computer program PACS are verified by solving a number of various fiber-reinforced composite laminates with and without geometric discontinuity. The comparisons of the numerical results to the experimental and other numerical results available in the literature indicate the validity and efficiency of the developed model.

복합재료 거동특성의 파괴해석 Ⅱ-비선형 유한요소해석

이 규세
선문대학교 토목공학과 교수



    Composite materials partly behave in a nonlinear fashion, although composite materials generally have been modeled as linear elastic material. The nonlinearity of composite materials can be attributed to inherent material nonlinearity of individual constituents and to micromechanical failures such as fiber or matrix microcracking and interfacial debonding (Yener & Wolcott, 1988; Drucker, 1969). Furthermore, the development of advanced composite material (metal matrix composites) and the complex design of composite structures make it more complicated to model material behavior of composites (Shih & Lee, 1978).

    Generally, the literature has assumed linear stress-strain relationships for composite materials.

    However, it has also been shown that composite materials partly behave in a nonlinear fashion (Yener & Yi, 1989; Petit & Waddoups, 1969; Hahn & Tsai, 1973). The nonlinearity of composite materials can be attributed to inherent material nonlinearity of the individual constituents.

    As mentioned in other paper, the most constitutive models are often used in the practical design field for their simplicity and in spite of their limitations (Sandhus, 1974). Their unavoidable shortcoming for a predictive analysis procedure, however, is that they cannot predict the material behavior of permanent strain accumulation by large deformation. That is the reason why the anisotropic plastic constitutive model is developed in the other paper, which includes the anisotropic yield criterion, a non-proportional hardening rule for changing anisotropic parameters. Hence, highly anisotropic material such as unidirectional composite materials, and their model can be analyzed with account for differential between the tensile and compressive yield strengths (i.e., Bauschinger effect) (Hill, 1948, 1950; Valliappan, 1972).

    As mentioned earlier, general response prediction of composite structures becomes possible by developing a realistic and comprehensive analysis procedure for general loading conditions. Such an undertaking, among other considerations, requires very efficient constitutive model which can predict realistically nonlinear material behavior. Hence, an anisotropic plasticity constitutive model for fiber-reinforced composite laminates is implemented into a computer program for a predictive analysis procedure of composites (PACS) (Yener & Yi, 1994).


    2.1.Characteristics of Predictive Analysis

    The proposed predictive nonlinear analysis procedure is presented in here. In a predictive analysis procedure, load is applied to a structure incrementally. At each prescribed step of loading, the desired information describing the behavior of the structure are computed and kept to detect the progressive structural behavior. These behavioral characteristics generally necessitate the utilization of numerical methods. Clearly, because of its lengthy incremental nature, a predictive analysis requires the development of a general purpose computer program.

    2.2.A Nonlinear Analysis Procedure

    In Ref. 15, the finite element governing equations have been formulated, which contain the displacement, strain-displacement matrices and the constitutive matrix of the material. For a formulation to be applicable to general nonlinear problems, however, the constitutive description must be appropriate. The anisotropic elastic-plastic constitutive model, which has been presented in the previous section, is employed. By including the anisotropic elastic-plastic constitutive matrix into the governing equation, the current formulation is applicable to the geometric and material nonlinear problem. The details of finite element formulation and numerical implementation of the total Lagrangian approach is presented in the textbook (Chen & Han, 1988), so is omitted in here.

    The nature of the nonlinear equation requires the use of an incremental solution procedure. For each load increment, iterations are generally necessary to satisfy internal equilibrium. When appropriate, the behavior may be assumed to be linear for each iteration within the specified load increment. Successive iterations within each increment are performed to avoid divergence of the solution. The process, here referred to as the modified Newton-Raphson method, is schematically illustrated in Fig. 1. For each load increment, unbalanced load is calculated, which then is used to compute the displacement increment in the proceeding iteration step. This internal iteration process is repeated until the equilibrium is approximated to some acceptable tolerance.

    2.3.Incremental Elastic-Plastic Constitutive Relationship

    This section presents numerical calculations of stress increments. A proper elastic-plastic constitutive relationship integration algorithm is very important in maintaining numerical accuracy, stability, and efficiency for the whole solution procedure (owen & Figueiras, 1983). Several integration techniques for incremental elastic-plastic constitutive relationships are available in the literature (Yener & Yi, 1990, 1992; Tsai & Hahn, 1980). This study uses the stress integration scheme, which is a modified version of the elastic predictor-radial return method (Schreyer, Kulak & Kramer, 1979; Ortiz & Popov, 1985). The elastic predictor-radial return method is first due to Mendelson (1986). During application of an load increment, an element or part of an element can yield. The yield condition is checked to see whether or not plastic deformation has occurred at the Gaussian points. As a result, an element can behave elastically in some parts and plastically in other parts.

    If the material behaves plastically for each iteration, the stress should be calculated by an elastic-plastic constitutive relationship. The incremental constitutive relationship for elastic-plastic material is written in matrix form, and then the numerical implementation of the stress computation is followed.

    In matrix form, the stress increment {} can be expressed in terms of the elastic strain increment {e}, or the total strain increment {} as

    d σ = D e d e = D e d d p
    d σ = D ep d

    where [De] and [Dep] are the elastic and elastic-plastic stress-strain relationship matrix, respectively. The elastic-plastic constitutive relationship was already given as Eq. 36 in a tensor form, rewritten here in matrix form as

    D ep = D e D e f σ f σ T D e σ ¯ p + f σ T D e f σ

    Stress calculation is performed at all Gaussian integration points.

    Now, the first step of the numerical algorithm is to define a trial stress increment {Δσe}, assuming an elastic behavior,

    Δ σ e = D e Δ

    Accumulate the trial total stress by adding the assumed elastic stress increment.

    In a predictive analysis procedure, the stress state has to be monitored at every step. At this point, three possible conditions can be used to calculate the stress increment. Assuming that the stress state from the previous step was in elastic state, the initial yield condition must be checked. If the trial total stress vector does not violate the initial yield criterion, the elastic material behavior is presumed. In the second case, the trial stress increment has passed the initial yield surface, as shown in Fig. 2. Then, the strain increment is divided into two parts, such as the pure elastic response R{Δϵ}, and the elastic-plastic response (1 - R ){Δϵ}, where R is called a scaling factor. Hence, the stress increment is also divided into two parts that can be integrated as

    Δ σ = R Δ e + n 1 + R Δ n 1 + Δ D ep d

    If, as in the third case, the stress state at the end of n - 1th load step is already in a plastic state and in a plastic loading state, the stress increment is obtained by integrating the stress increment vector, Eq. 4 with R = 0 as

    Δ σ = R Δ e + n 1 n 1 + Δ D ep d

    Numerical implementations of integrations appearing in Eqs. 4 and 5 are detailed in Ref. 31.

    If the stress is an unloading state, an elastic constitutive relation should be used as

    Δ σ = Δ σ e

    where the stress increment vector is negative corresponding to negative strain increment.

    2.4.Computational Procedure for Predictive Analysis

    Based on the preceding discussion, a predictive finite element analysis procedure, is presented through the flow chart of Yener and Yi (1994). Also, it is implemented in the computer program, PACS which stands for Predictive Analysis of Composite Structures. As shown in the flow chart of Yener and Yi (1994), the essential steps required in a predictive solution algorithm for nonlinear problems are presented to provide a physical insight into the viability of the proposed algorithm.


    Anisotropic plasticity and a predictive analysis procedure have been presented for analysis of fiber-reinforced composite structures. Also, their models are successfully implemented into the computer program PACS (Yener & Yi, 1994). In order to verify the validity of the problem formulation and the efficiency of numerical implementation, several example problems are analyzed by using the computer program PACS. Numerical results are compared with various experimental and other numerical results.

    3.1.Elastic-Plastic Analysis of Composite Laminates

    In this section, the anisotropic elastic-plastic analysis capability of the computer program is verified with several composite laminates on different boundary and loading conditions. The elastic-plastic material behavior is observed and compared. Also, the behaviors of isotropic and anisotropic materials are compared.

    3.2.Clamped Square Plate

    A clamped square plate under uniformly distributed loading is considered. Only a quarter of the plate has been discretized because of symmetry considerations (Fig. 3). Eight-noded isoparametric elements with four-point (2x2) Gaussian integration in the plane are used, and eight equal thickness layers are taken through thickness direction. The plate has been solved with both isotropic and anisotropic materials; material properties are given in Table 1.

    The thickness and span of plate is h = 0.2 m and L = 6.0 m, respectively.

    In Fig. 4, vertical displacements at the plate's center are plotted for different load levels for both isotropic and anisotropic materials. The clear difference between isotropic and anisotropic material is observed. Fig. 4 also includes comparison with other numerical results, and agreement with present results is rather close.

    3.3.Simple Supported Square Plate

    This section analyzes the simple supported plate on four sides of the plate with geometry of plate and finite element discretization given in Fig. 5. The purpose of this problem is to see effects of different boundary conditions on the plastic behavior of the plate. In order to compare with the clamped plate, the same geometry as the clamped plate is used with a uniformly loaded case. The comparison is performed only on isotropic material. The development of plastic regions for both clamped and simple supported plates are presented in Fig. 6.

    Yielding starts at the corners in the simple supported plate, while starting in the middle of the sides of the clamped plate. Next, the center of both plates yields. Then the plastic regions propagate from the vicinity of the corners and the plate centers. The pattern of the plastic region growth agrees with that of Dodds (1987).

    3.4.Clamped Quadratic Square Shell

    An elastic-plastic analysis of a clamped spherical shell under self-weight condition is undertaken. The geometry of the spherical shell is given in Fig. 7. The coordinates of z-direction are obtained as

    Z = C L 2 2 x 2 + y 2

    where L is the width of shell and c is equal to L/10. Material characteristics are the same as employed in Table 1 of the plate problem. Finite element discretization is same as that of the plate problem in Fig. 5, where nine-nodes Heterosis element is used.

    Displacements at center of shell for both isotropic and anisotropic materials are plotted in Fig. 8, which show numerical results from Owen and Figueiras (1983) and good agreements.

    The spherical shell is again solved under a concentrated load at center. Vertical displacement at center is plotted against the applied load in Fig. 9.


    The development and use of an analytical procedure which has capability of predicting the progressive material behavior of structures, named as a predictive analysis procedure in here, is developed for more accurate assessment of structural safety and efficiency of composite structures.A quadratic anisotropic yield criterion in stress space is developed for general use with unidirectional and bidirectional composite lamina. The developed anisotropic work-hardening model allows for a nonproportional change of the yield values so that the subsequent yield surface can be a distorted shape. The predictive finite element analysis procedure is developed for accurate and efficient analysis in predicting progressive behavior of composite structures. As a result, a computer code, PACS (Predictive Analysis of Composite Structures) is developed, which adopts the abovementioned anisotropic plasticity model.

    The accuracy and efficiency of the computer code PACS are verified with various benchmark problems. Numerical predictions of the computer program PACS compare very well with available experimental, analytical and other numerical results. Comparisons illustrate the capability of the constitutive model and the computer program PACS. The developed constitutive and predictive analysis model predicts progressive nonlinear behavior from the beginning of loading, in plane and through thickness direction of composite laminates. The capability of PACS to predict nonlinear material behavior of composites can be used as a helpful device in parametric studies for design of fiber-reinforced laminates.



    Incremental-iteration Solution Procedure Shown on the Load-Deflection Curve [20].


    Schematic Illustration of Entering a Plastic State at the First Iteration Step in Proceeding from the nth to n + 1th Load Increment Step.


    The Geometry of Clamped Plate and Finite Element Discretization.


    Load-Displacement Curves at the Center of Clamped Plate for Isotropic and Anisotropic Material.


    (a) The geometry of simple supported plate and (b) finite element discretization.


    Development of plastic region for (a) clamped and (b) simple supported plate with uniform load.


    The Geometry of Clamped Quadratic Shell.


    The Curve for Displaceement at Center Versus Self-Weight of Spherical Shell


    The Load Versus Vertical Displaceement at Center of Shell.


    Material constants used in clamped plate (Unit: MN/m2


    1. Chen W.F , Han D.J (1988) Plasticity for Structural Engineering, Springer-Verlag,
    2. Herakovich CT Inelastic Behavior of Composite Materials, AMD, ASME, Vol.13; pp.1-15
    3. Dodds R.H Jr (1987) “Numerical Techniques for Plasticity Computations in Finite Element Analysis” , Computer and Structure, Vol.26; pp.767-779
    4. Hahn H , Tsai S.W (1973) “Nonlinear Elastic Behavior of Unidirectional Composite Laminae” , Journal of Composite Materials, Vol.7; pp.102-118
    5. Hill R (1948) “A Theory of Yielding and Plastic Flow of Anisotropic Metals” , Proc. Roy. Soc. Lond. Ser. A, pp.193-0
    6. Hill R (1950) The Mathematical Theory of Plasticity, Oxford University Press,
    7. Mendelson A (1986) Plasticity: Theory and Application, Krieger,
    8. Ortiz M , Popov E.P (1985) “Accuracy and Stability of Integration Algorithms for Elastoplastic Constitutive Relations” , International Journal for Numerical Methods in Engineering, Vol.21; pp.1561-1576
    9. Owen D.R.J , Figueiras J.A (1983) “Anisotropic Elasto-Plastic Finite Element Analysis of Thick and Thin Plates and Shells” , International Journal for Numerical Methods in Engineering, Vol.19; pp.541-566
    10. Petit P.H , Waddoups M.E (1969) “A Method of Predicting the Nonlinear Behavior of Laminated Composites” , Journal of Composite Materials, Vol.3; pp.2-19
    11. Sandhu R.S (1974) “Nonlinear Response of Unidirectional and Angle Ply Laminates”,
    12. Shreyer H.L , Kulak R.F , Kramer J.M (1979) “Accurate Numerical Solutions for Elastic-Plastic Models” , ASME Journal of Pressure Vessel Technology, Vol.101; pp.226-234
    13. Shih C.F , Lee D (1978) “Further Developments in Anisotropic Plasticity” , Journal of Engineering Material and Technology, Vol.100; pp.294-302
    14. Tsai S.W , Hahn H.T (1980) “Introduction to Composite Materials”, Technomic Publishing Company Inc, pp.0-0
    15. Valliappan S (1972) “Elasto-Plastic Analysis of Anisotropic Work-Hardening Materials” , Archives of Mechanics, Vol.24 (3) ; pp.465-481
    16. Yener M , Wolcott E (1988) “Damage Assessment Analysis of Composite Pressure Vessels Subjected to Random Impact Loading” , ASME, ASME,
    17. Yener M , Yi G.S (1990) “On the General Characteristics of Composite Materials” , Utah State University, Structural Engineering and Mechanics Division Report, CEE - SEMD - 90-07, Logan, UT,
    18. Yener M , Yi G.S (1992) “Constitutive Relationships and Failure models for Composite Materials” , Utah State University, Structural Engineering and Mechanics Division Report, CEE-SEMD-92-11, Logan, UT,
    19. Yener M , Yi G.S (1994) “User's Manual for Predictive Analysis of Composite Structures, Utah State Univ” , Structural Engineering and Mechanics Division Report, CEE - SEMD - 94-10, Logan. UT,
    20. Yi G.S (2014) “A New Orthotropic Yield Criteria For Fiber-Reinforced Composite Materials” , The 4th International Conference on Concergence Technology 2014, Vol.1; pp.313-316
    21. Yi G.S , Yoo H.J , Choi J.Y , Lee J.O , Lim N.H (2012) “Elastic Stability of Circular arches with the Open Thin-Walled Monosymmetric Section Considering the Prebuckling Deformation” , Open Civil Engineering Journal, Vol.6; pp.87-97