2. Different 3D finite element and analytical approaches incorporating effects of fiber misalignment have been ..... Mode II matrix fracture energy, GIIC (lbf/in).
Compressive failure modeling of composites with defects Dinh Chi Pham, Xiaodong Cui, Anand Karuppiah, and Jim Lua Global Engineering and Materials Inc. 1 Airport Place, Princeton, NJ 08540, USA
Fabrication and service induced defects have substantial influence on the performance of composite materials and structures. Characterization of their effects is essential for structural design, certification, and sustainment. In optimization of a fabrication process and risk informed design of composite structures with defects, high-fidelity modeling toolkits are essential to reduce testing matrix and design iterations for composite structures with defects. Among the manufacturing defects, ply waviness is one of the most important types of defects, which can significantly reduce the compressive strength of composite structures. To study the physical relation between the ply waviness and the compressive strength, a continuum damage mechanics (CDM) based progressive failure model is adopted. The initiation and evolution of matrix damage and fiber failure is predicted under compressive loading along with the formation of fiber kink bands. The effects of two and three dimensional wavy profiles are included in the progressive failure analysis by considering elevated degrees of fiber misalignments. The prediction from two dimensional (2D) ply waviness models using different misalignment angles are verified against the analytical solutions. Effects of hydrostatic pressure modeling on the compressive strength prediction are also examined. This approach is further extended to explore the effects of three dimensional (3D) ply waviness on the compressive strength of composite laminates.
I. Introduction Various forms of defects are found in composite structures including ply waviness, porosity, ply-thickness variation, etc. Among these defects, ply waviness [1-3] is the most common defect from a manufacturing process. Its distribution and geometric configuration are strongly reliant on the tooling and fabrication parameters such as vacuum pressure, curing temperature, resin viscosity, flow path, etc. [1]. The presence of these defects has a significant impact on the performance of composite structures. The effects of ply waviness on the mechanical behavior become profound when it is loaded in compression. Matrix failure and fiber kinking can occur under compressive loading. To analyze the effects of ply waviness, the inclusion of wavy geometries into finite element models can be made and characterization of their effects on the damage initiation and progression can be achieved through a progressive failure analysis [3-5]. Current researches have been focused mainly on the characterization of waviness effects along one direction (or 2D waviness) on the compressive response of unidirectional composites [1-3, 6-9]. Limited work has been performed by incorporating 3D wavy plies in a progressive failure analysis model. The effects of ply waviness along one direction or 2D ply waviness have been studied in Ref. [1]. It is shown that both the axial stiffness and compressive strength are degraded with increasing fiber misalignment angle. Based on the measurement of fiber or ply waviness angles, a conforming mesh can be generated at either micro- or meso- scale with fiber orientation dependent properties assigned for each element. As shown in Refs [6, 7], both in-plane and out-of-plane ply wavy angles are observed in composite wind turbine blades along with their uncertainty characterization. The influence of gross misalignment of plies on the compressive strength is investigated in Ref. [6] using a 2D finite element model with wavy plies characterized by a cubic spline interpolation function. As concluded in Ref. [6], moving to 3D and modeling each ply independently can improve the solution accuracy to capture the kinking and splitting failure.
1
Different 3D finite element and analytical approaches incorporating effects of fiber misalignment have been carried out to predict the compressive failure and fiber kinking mechanisms in composites. Argon [8] first developed a model for fiber kink-band formation in unidirectional composites. An interlaminar shear stress was reported in Ref. [8] due to initial fiber misalignment which cause a sliding and rotation of the fibers in a band width. Budiansky [9] considered the increased fiber rotation angle during the kinking and showed that the compressive strength is sensitive to the fiber misalignment. In both Ref. [8] and Ref. [9], the composite compressive strength was related to the matrix in-plane shear stress and the initial misalignment angle. Vogler et al. [10] studied the initiation of kink bands through 2D and 3D micromechanical models of unidirectional composites regarding local and global imperfections. They found that the kink band width increased with the fiber diameter and fiber volume fraction, whereas the band inclination angle is insensitive to those parameters. Yerramalli and Waas [11] presented a 3D micro model to demonstrate the increase in compressive strength with an increase in fiber diameter. Bishara et al. [12] revealed the complex aspects of composite compressive failure through fiber kinking at microscale. A 3D micro model was presented in Ref [12] to predict the compressive strengths considering fiber waviness, its statistical distribution and the angle of kink band. The kink band angle was shown to be dependent on the fiber tensile strength, and the compressive strength was sensitive to ratio of the maximum wave amplitude and the wave length of waviness. In an extended work by Bishara [13], the compressive failure of composite laminates were illustrated as an interaction of four failure modes: fiber kinking, fiber splitting, matrix cracking and delamination. A homogenized material model of unidirectional ply and associated micro mechanical model were then developed in Ref. [13] to simulate the kinking process. Tsuyoshi et at [14] illustrated the compressive failure mechanism and strength of unidirectional thermoplastic composites based on a modified kink band model by applying a quadratic yield function accounting for the effects of transverse tensile and shear stress in the kink band. The proposed model in Ref. [14] has been validated against series of compressive failure tests at different temperatures. Sun et al. [15] carried out the experiment and simulation of failure mechanisms in notched CFRP under compression by considering the effects of manufacturing defects. The results suggested that the voids in the resins lead to the initiation of micro matrix cracks. These matrix cracks reduced the fiber/matrix bonding shear strength, thus causing fiber splitting and inducing the kink-band formation. The failure of unidirectional laminates under the combined compression-transverse tension and compression-shear loads were investigated numerically in Ref. [15]. The transverse tension was shown not to influence the kink-band failure mode while an increase in the in-plane shear decreased the compressive strength. In summary, most of the studies on compressive modeling of composites with defects have been focused on the effects of fiber waviness along the fiber direction. General 3D waviness as a manufacturing defect can be seen in both the fiber and transverse-to-fiber direction. Such information for a 3D wavy surface can be obtained from X-ray scanning and optical microscopy of composite specimens [3]. A 3D finite element based progressive analysis with a mapping of 3D waviness from X-ray CT has not been fully explored to explicitly evaluate their effects on damage initiation and damage progression in laminated composite structures under compression. In this study, compressive responses of composite materials considering both 2D and 3D ply waviness are investigated by applying a continuum damage mechanics (CDM) approach. An automatic 3D finite element mesh generation method for composites with ply waviness along two directions is developed through Python scripting. Assignment of fiber orientation is made through a user-defined orient subroutine and cohesive interfaces between plies are automatically embedded using a Python script. The intra-ply damage is predicted using a developed CDM model and inter-ply delamination is predicted using Abaqus’ cohesive model. The predictive capability of our CDM model is demonstrated through its application to of 2D and 3D wavy specimens with different initial fiber misalignments. The numerical results of 2D wavy plies are verified with the corresponding analytical solutions from Arkon [8] and Buldansky [9]. The effects of hydrostatic pressures on the compressive strength are also investigated. Further studies using 3D wavy models based on the true configuration of 3D waviness from scanning images of an interlaminar tensile (ILT) specimen are performed to demonstrate the applicability of our CDM model in capturing ply waviness effects on the induced damage initiation and failure progression.
II. Overview of Continuum Damage Mechanics (CDM) Model In the CDM approach, mechanism driven failure criteria are used for damage initiation while energy driven es are used for damage propagation. A set of LaRC04 failure criteria following the method of Maimi et al. [16] is implemented to determine the constitutive failure induced composite damage initiation. The LaRC04 failure criteria distinguish four fracture mechanisms as expressed in Eq. (1).
2
t 11 12 22 Ff XT m m L 22 F c 12 f SL 2 2 22 22 12 t 1 g Y T g Y T S L Fm 1 L L 12 22 S 2 T L 2 eff eff t F f ST S L T
T
where X , Y , S
L
11 0
Fiber tension
11 0
Fiber compression
22 0
Matrix tension
22 0
Matrix compression
22 0
Matrix compression
(1)
T
and S are the longitudinal tensile strength, transverse tensile strength, longitudinal shear
strength and transverse shear strength, respectively. is the longitudinal friction coefficient. 11 , 22 and 12 are L
components of the undamaged effective stress tensor. 12 and 22 are the components of the effective stress tensor m
m
T L in the coordinate system associated with the rotation of the misaligned fiber. eff and eff are the effective stresses
calculated according to the sliding angle and compressive fracture angle 0 53 . To incorporate nonlinear shear behavior, the Ramberg-Osgood formulation is used to characterize the constitutive behavior under either an in-plane or out-of-plane shear. The nonlinear shear stress is given by:
/ G0 / K
1/ n
,
(2)
where G0 is the elastic shear modulus; K and n are material constants. Once an initial failure mode is detected, its subsequent damage accumulation is governed by an exponential damage evolution law given by, 1 (3) dM 1 exp AM 1 f N rN M 1, 1-, 2 , 2-, 6; N 1,1-, 2 , 2 f N rN where AM is a softening law parameter, rN and f N rN are the elastic domain thresholds and the functions that govern the softening behavior of a constitutive relation, respectively. In addition to the LARC04 criteria, the 2D Hashin failure criteria were also employed. The 2D Hashin criteria describe the matrix and fiber failure modes under tension and compression loading cases by the following equations:
Tensile fiber failure (FT) for σ11 ≥ 0: 2
FI FT
11 XC
2
(6)
Tensile matrix failure (MT) for σ22 ≥ 0: FI MT
(5)
Compressive fiber failure (FC) for σ11 < 0: FI FC
2 11 122 X T S12
22 2
YT2 Compressive matrix failure (MC) for σ22 < 0:
122 S122
(7)
Y 2 2 (8) FI MC C 1 22 122 4 S 232 2 S 23 YC S12 Post-peak failure in each failure mode of the 2D Hashin criteria is modeled by a crack band approach. In the crack band approach, it is assumed that distributed cracks are smeared out over a certain width within a finite element such
22
2
3
that the effect of progressive cracking is represented by linear strain softening in a continuum scheme. For composites, the crack band approach assumes that the crack evolution under each failure mode is governed by SERR in the respective mode [5, 17]. For an element of its characteristic length of Lc, the energy dissipation due to each failure mode is dictated by Gi as:
i eq
i c d ( eq L ) Gi
(i FT , FC , MT , MC )
(9)
i i i i where f , eq , eq and eq are the total crack length and equivalent stress, strain and crack displacement of each
failure mode.
III. Fiber Misalignment Dependent Compressive Strength 3D unidirectional composite models with uniform 2D wavy profile have been constructed through Python scripts considering various initial fiber misalignment angle 𝜃, ranging from 0 to 10 degrees (Fig. 1). The wavy profile is characterized by the initial wave angle, wavelength and amplitude. The 2D wavy models are subject to compressive loading at one end while the other end is constrained, as shown in Fig. 2. The selected material is AS4/8552 with its mechanical properties listed in Tables 1 and 2 and the layup is [0]5. The LARC04 and Hashin failure criteria were implemented in an user defined subroutine (UMAT) and the analyses are carried out using Abaqus’ implicit solver. Fig. 3 shows the predicted compressive strength using LaRC04 failure criteria in comparison with the analytical solutions proposed by Argon [8] and Budiansky [9]. Argon’s model relates the applied compressive stress to the induce shear stress y and initial fiber misalignment 𝜃 by Eq. (10).
XC y
(10)
Budiansky [9] extended Argon’s analysis and assumed that the transverse stress does not affect the compressive failure and the kink band forms at yield stress y . The compressive stress is then related to shear y , fiber misalignment 𝜃 and the shear strain at yielding y as in Eq. (11).
XC y
y
(11)
As can be seen in Fig. 3, the compressive strength predictions of unidirectional 2D wavy plies agree well with the analytical results from Argon [8] and Budiansky [9], when the misalignment angle is larger than 3º. The developed CDM model successfully predicts the reduction trend in compressive strengths with increasing fiber alignment angles. It is noted that the Argon model yields infinite compressive strength for the perfect case (𝜃 = 0°) while the Budianky model gives a very high compressive strength. The compressive strength prediction by LaRC04 for the case without fiber waviness matches the tested composite strength Xc [18]. Table 1 Mechanical properties of AS4/8552 composites [18]. Description Longitudinal tensile modulus, E11 (Msi) Transverse tensile modulus, E22 (Msi) In-plane shear modulus, G12 (Msi) Major Poisson’s ratio, v12 Longitudinal tensile strength, XT (psi) Longitudinal compressive strength, XC (psi) Transverse tensile strength, YT (psi) Transverse compressive strength, YC (psi) Longitudinal shear strength, SL (psi)
Value 19.09 1.34 0.7 0.3 299220 215290 9270 38850 13280
Table 2 AS4/8552 fracture energies. Description Tensile fiber fracture energy, GFT (lbf/in) Compressive fiber fracture energy, GFC (lbf/in) Mode I matrix fracture energy, GIC (lbf/in) Mode II matrix fracture energy, GIIC (lbf/in)
4
Value 485 485 0.94 4.71
00
10
20
30
60
40
80
50
100
Fig. 1 FE Models with different misalignment angles.
Fig. 2 Boundary conditions for 2D wavy laminates.
Compressive Strength Xc (psi)
8.E+05 7.E+05 Budiansky 6.E+05 Argon
5.E+05
FEMLaRC04
4.E+05 3.E+05 2.E+05 1.E+05 0.E+00
0
1
2
3
4
5
6
7
8
9
10
Misalignment angle (0)
Fig. 3 Comparison of the compressive strength predictions between the CDM model and analytical solutions from Argon [8] and Budiansky [9].
5
Effects of 2D and 3D waviness FE models considering 3D waviness have been built assuming the wavy profile in the y-direction is the same as the profile in the x-direction. The case of initial misalignment of 2ºwas selected. LaRC04 failure criteria were used to predict the compressive strength. Fig. 4 shows the 3D waviness FE model with the waviness along the x- and ydirections. The comparison on the predicted load-displacement curves between the 2D and 3D waviness models were plotted in Fig. 5. The predicted failure load by the 3D waviness model is much lower than the failure load predicted by the 2D waviness FEM with the initial ply misalignment of 2º. It can be seen that correctly characterizing and modeling 3D waviness effects in composites is important to accurately predict the compressive strength of composites.
Fig. 4 2D waviness model along x- and y- directions with initial misalignment of 2º.
Compressive strength Xc (ksi)
140
20
120 100
2D waviness FEM
80
3D waviness FEM
60 40
20 0 0
0.05
0.1 0.15 Displacement (in)
Fig. 5 Comparison on the predicted compressive strength by the 2D and 3D waviness model.
Effects of hydrostatic pressure on the composite compressive strength To examine the effects of hydrostatic pressure on the compressive strength of unidirectional composite, the composite has been further tested under both longitudinal compression and hydrostatic pressure. The FE model for the combined load case is shown in Fig. 6. Both the LaRC04 and Hashin failure criteria were used to analyze the effects of hydrostatic pressure on the compressive strength prediction. The comparison on the predicted results by the two failure criteria are provided in Fig. 7. The 2D Hashin model does not predict an obvious change in compressive strength with increasing hydrostatic pressure, while the LaRC04 model predicts a substantial increase in compressive strength. Therefore, the LaRC04 model is capable of capturing the compressive failure behavior of composite under hydrostatic pressure, which is consistent with the experimental observation [19].. Furthermore, it can be seen that effects of hydrostatic pressure on the compressive strength is not negligible and should be taken into consideration for applications with the presence of the hydrostatic pressure such as those resulted from the bolt preloading in composite bolted joints.
6
Fig. 6 FE compression model considering hydrostatic pressure effects. 1.8
Normalized (Xc/Xc0)
1.6 1.4 1.2 1 0.8
FEM LaRC04
0.6
FEM Hashin
0.4
0.2 0 0
5
10 15 Hydrostatic pressure (ksi)
20
Fig. 7 FE compression model considering hydrostatic pressure effects.
IV.
Model Predictions and Verifications
In addition to 2D ply waviness models, a representative 3D ply waviness model is extracted from an ILT beam specimen as shown in Fig. 8. The ply waviness along the length (x) direction is digitized from a cross-sectional view of photomicrographs while the ply waviness along the width (z) direction is pre-assumed. A layup of [45/90/-45/0]3 is used with a generic material of AS4/8552. Fig. 9 shows the FE model including 11 cohesive layers at each ply interface. Compressive failure analysis of the 3D waviness laminates is modeled using LaRC04. In addition to the intra-ply damage prediction by the LaRC04 failure criteria described above, inter-ply delamination is characterized using an existing Abaqus’ cohesive element. Delamination onset is assumed to follow a quadratic failure criterion:
n
2
N2
s2 S2
t2 T2
1
(12)
where n , s , t are the normal and two shear tractions and N, S, T are the corresponding normal and shear interfacial strength. For mixed mode delamination, equivalent traction Teq and displacement eq are introduced, and delamination propagation is indicated by a mixed-mode energy Benzeggagh and Kenane law:
GII GIII GC GIC (GIIC GIC ) GI GII GIII
T
G
eq d eq
(13)
(14)
The progressive damage results of the [45/90/-45/0]3 laminate with 3D waviness are presented in Fig. 10. As can be seen in Fig. 10, the predicted initial compressive matrix and fiber damage causes nonlinear behavior of the 3D laminate prior to the peak load (Point a). With the increase of the applied load, interface delamination takes places and compressive fiber failure is initiated in the 0º plies (Point b), leading to the main load drop. Further loading causes extensive delamination at various ply interfaces and additional fiber breakage (Point c, d). The final damage patterns of the [45/90/-45/0]3 wavy laminate is provided in Fig. 11. Numerous compressive matrix failure, fiber failure and delamination are found, and fiber kink banding is triggered after the compression fiber failure.
7
Note: Given ply waviness along z-direction is arbitrary for the demonstration of 3D waviness.
L-Beam Model
Wireframe of L-Beam
Fig. 8 Generation of a 3D waviness local model using the extracted ply waviness data from an ILT specimen.
Fig. 9 A local model with 3D wavy plies (left) and cohesive interfaces at each wavy ply interface (right).
8
(b) Compressive fiber failure (a) Compressive fiber damage 6,000
(a)
(b)
5,000
(c)
Axial Load (lbf)
4,000
(d) 3,000
(c) Compressive fiber failure and delamination
2,000 1,000
0 0
0.005
0.01
0.015
Axial Displacement (in)
Load displacement curve for [45/90/-45/0]3 laminate
(d) Delamination extension and fiber failure
Fig. 10 Compressive damage prediction of 3D waviness [45/90/-45/0]3 laminate. Compression
Compressive matrix failure
Compressive fiber failure
Fig. 11 Predicted compressive matrix and fiber failure and fiber-kink banding.
V. Conclusions The compressive responses and failure prediction of composite laminates with imperfection have been studied, focusing on the effects of 2D and 3D ply waviness, which is known to have detrimental effects on the compressive behavior of composite laminates. These effects have been examined on unidirectional composite laminates considering ply waviness with different fiber misalignment angles. It is found that the predicted composite compressive failure stress is very sensitive to the initial fiber misalignment. The CDM model implementation and numerical predictions have been verified using the analytical results reported in the literature and good agreement is achieved. In addition, the effects of hydrostatic pressure on the compressive strength have been addressed, showing an increase in composite strength with an increase in pressure application. Besides, the FE results for the 3D waviness local model extracted from photomicrograph of an ILT specimen show complicated failure modes including matrix cracks, fiber compressive failure, extensive delamination and kink band formation. Further studies will be performed for 3D waviness to capture the effects with different fiber alignment angles.
9
References [1] Hsiao, H. M., and Daniel, I. M. "Effect of Fiber Waviness on Stiffness and Strength Reduction of Unidirectional Composites Under Compressive Loading," Composites Science and Technology Vol. 56, No. 5, 1996, pp. 581-593. doi: https://doi.org/10.1016/0266-3538(96)00045-0 [2] Nikishkov, G., Nikishkov, Y., and Makeev, A. "Finite Element Mesh Generation for Composites with Ply Waviness Based on X-ray Computed Tomography," Advances in Engineering Software Vol. 58, 2013, pp. 35-44. doi: https://doi.org/10.1016/j.advengsoft.2013.01.002 [3] Lua, J., Pham, D. C., Sadeghirad, A., Karuppiah, A., Cui, X., Seneviratne, W., and Phan, N. "Characterization of Fabrication Induced Defects and Assessment of Their Effects on Integrity of Composite Structures," AHS International 74th Annual Forum & Technology Display. Phoenix, Arizona, 2018. [4] Pham, D. C., Cui, X., Lua, J., and Zhang, D. "A Continuum Damage Description for A Discrete Crack Modeling Approach for Delamination Migration in Composite Laminates," 2018 AIAA/ASCE/AHS/ASC Structures, Structural Dynamics, and Materials Conference. American Institute of Aeronautics and Astronautics, 2018. [5] Pham, D. C., Sun, X. S., Tan, V. B. C., Chen, B., and Tay, T. E. "Progressive Failure Analysis of Scaled Double-Notched Carbon/Epoxy Composite Laminates," International Journal of Damage Mechanics Vol. 21, No. 8, 2011, pp. 1154-1185. doi: 10.1177/1056789511430415 [6] Cairns, D. S., Nelson, J. W., Woo, K., and Miller, D. "Progressive Damage Analysis and Testing of Composite Laminates with Fiber Waves," Composites Part A: Applied Science and Manufacturing Vol. 90, 2016, pp. 51-61. doi: https://doi.org/10.1016/j.compositesa.2016.03.005 [7] Riddle, T., Cairns, D. S., and Nelson, J. "Effects of Defects Part A: Treatment of Manufacturing Defects as Uncertainty Variables in a Wind Blade Probabilistic Design Framework," 32nd ASME Wind Energy Symposium. American Institute of Aeronautics and Astronautics, 2014. [8] Argon, A. S. Treatise on Materials Science and Technology. New York: Academic Press, 1972. [9] Budiansky, B. "Micromechanics," Computers & Structures Vol. 16, No. 1, 1983, pp. 3-12. doi: https://doi.org/10.1016/0045-7949(83)90141-4 [10] Vogler, T. J., Hsu, S. Y., and Kyriakides, S. "On The Initiation and Growth of Kink Bands in Fiber Composites. Part II: Analysis," International Journal of Solids and Structures Vol. 38, No. 15, 2001, pp. 2653-2682. doi: https://doi.org/10.1016/S0020-7683(00)00175-X [11] Yerramalli, C. S., and Waas, A. M. "The Effect of Fiber Diameter on the Compressive Strength of Composites - A 3D Finite Element Based Study," Computer Modeling in Engineering & Sciences Vol. 6, No. 1, 2004, pp. 1-16. [12] Bishara, M., Rolfes, R., and Allix, O. "Revealing Complex Aspects of Compressive Failure of Polymer Composites – Part I: Fiber Kinking at Microscale," Composite Structures Vol. 169, 2017, pp. 105-115. doi: https://doi.org/10.1016/j.compstruct.2016.10.092 [13] Bishara, M., Vogler, M., and Rolfes, R. "Revealing Complex Aspects of Compressive Failure of Polymer Composites – Part II: Failure Interactions in Multidirectional Laminates and Validation," Composite Structures Vol. 169, 2017, pp. 116-128. doi: https://doi.org/10.1016/j.compstruct.2016.10.091 [14] Matsuo, T., and Kageyama, K. "Compressive Failure Mechanism and Strength of Unidirectional Thermoplastic Composites Based on Modified Kink Band Model," Composites Part A: Applied Science and Manufacturing Vol. 93, 2017, pp. 117-125. doi: https://doi.org/10.1016/j.compositesa.2016.11.018 [15] Sun, Q., Guo, H., Zhou, G., Meng, Z., Chen, Z., Kang, H., Keten, S., and Su, X. "Experimental and Computational Analysis of Failure Mechanisms in Unidirectional Carbon Fiber Reinforced Polymer Laminates under Longitudinal Compression Loading," Composite Structures Vol. 203, 2018, pp. 335-348. doi: https://doi.org/10.1016/j.compstruct.2018.06.028 [16] Maimí, P., Camanho, P. P., Mayugo, J. A., and Dávila, C. G. "A Continuum Damage Model for Composite Laminates: Part I – Constitutive Model," Mechanics of Materials Vol. 39, No. 10, 2007, pp. 897-908. doi: https://doi.org/10.1016/j.mechmat.2007.03.005 [17] Pham, D. C., and Sun, X. "Experimental And Computational Studies on Progressive Failure Analysis of Notched Cross-Ply CFRP Composite," International Journal of Computational Materials Science and Engineering Vol. 01, No. 03, 2012, p. 1250023. doi: 10.1142/s2047684112500236 [18] Marlett, K., Ng, Y., and Tomblin, J. "Hexcel 8552 AS4 Unidirectional Prepreg 190 gsm & 35% RC Qualification Material Property Data Report," National Center for Advanced Materials Performance, Wichita, Kansas. Test Report CAM-RP-2010002, Rev. A, 2011, pp. 1-278. [19] Pinho, S. T., Davila, C. G., Camanho, P. P., Iannucci, L., and Robinson, P. "Failure Models and Criteria for FRP under Inplane or Three-dimensional Stress States Including Shear Non-linearity." NASA/TM-2005-213530, 2005.
10
