Non-linear strain invariant failure approach for fibre

0 downloads 0 Views 500KB Size Report
Abstract: The strain invariant failure theory (Gosse and Christensen, 2001) is a ..... engineering stiffness constants and three non-linear elastic material ...
284

Int. J. Materials and Structural Integrity, Vol. 6, Nos. 2/3/4, 2012

Non-linear strain invariant failure approach for fibre reinforced composite materials Kumar D. Mishra and Rani F. El-Hajjar* Engineering Mechanics and Composites Laboratory, Department of Civil Engineering and Mechanics, College of Engineering and Applied Science, University of Wisconsin-Milwaukee, Milwaukee, WI 53211, USA E-mail: [email protected] E-mail: [email protected] *Corresponding author Abstract: The strain invariant failure theory (Gosse and Christensen, 2001) is a micromechanics-based failure theory that depends on computation of critical strain invariants for the fibre and matrix phases. In this paper, a novel approach is developed for including the non-linear effects of the matrix through the implementation of the J2 flow theory for the matrix response. The non-linear isotropic hardening is introduced for the analysis of a non-linear matrix in square and hexagonal fibre packaging geometries. The results are presented for a common type of structural composite made from carbon fibres and an epoxy matrix. The proposed approach is used to predict the stress and strain concentration factors in a multi-cell array and shows the necessity of including the non-linear response in the analytical approach in this theory. Keywords: micromechanical approaches; unit cell; finite element method; strain invariant failure theory; SIFT; fibre reinforced composite; p-FEM. Reference to this paper should be made as follows: Mishra, K.D. and El-Hajjar, R.F. (2012) ‘‘Non-linear strain invariant failure approach for fibre reinforced composite materials’’, Int. J. Materials and Structural Integrity, Vol. 6, Nos. 2/3/4, pp.284––296. Biographical notes: Kumar D. Mishra completed his MS degree from the University of Texas at Austin after which he became a Research Engineer at the Engineering Mechanics and Composites Laboratory at the University of Wisconsin at Milwaukee. Currently, he is pursuing doctoral studies at Colorado State University while researching at Articulated Motion Laboratory. Rani F. El-Hajjar earned his PhD at the Georgia Institute of Technology and his Master’’s degree in Civil Engineering from the University of Texas at Austin. He is a licensed Professional Engineer in the State of Washington and was previously a Structural Analyst at The Boeing Company and a Senior Process Engineer at Intel Corporation.

Copyright © 2012 Inderscience Enterprises Ltd.

Non-linear strain invariant failure approach

1

285

Introduction

Several researchers have investigated the use of the strain invariant failure theory (SIFT) for predicting failure in composite structures. The SIFT (Gosse and Christensen, 2001) is a micromechanics-based failure theory that depends on computation of critical strain invariants for the fibre and matrix phases. The strain invariants achieved from analysis are amplified using results from micromechanics-based models that take into account the interaction of the fibre and matrix phases and account for the thermal stresses associated with the manufacturing process. The amplification factors are used to modify the strain values obtained from homogenised solutions obtained by using classical methods. Physically, the first strain invariant is related to the dilational failure associated with crazing in polymers or interlaminar failure of the matrix in laminated composites whereas the equivalent strain is related to the distortional failure where no volume change occurs. The critical values for these quantities are obtained from transverse tension and off-axis tests, respectively. The world wide failure exercise (Soden et al., 2004) investigated 19 failure theories for predicting the deformation and failure response of polymer composite laminates. The effort, while comprehensive in many aspects, did not explain the failure predictions using interactive theories that are based on semi-empirical methods for determining the fitting factors. The SIFT has been used by several researchers to predict failure of composite structures (Gosse and Christensen, 2001; Gosse, 2001). Tay et al. (2005) obtained good correlation between experimental results when using the SIFT for composite damage progression analysis on a three-point bend specimen, using an element failure method that partially fails elements when the failure criteria are met. Li et al. (2003) used the first strain invariant to successfully predict the matrix-dominated failure in I-beams and T-cleat structures. The studies considered have typically assumed linear material properties for the matrix materials. Material non-linearity in the matrix is known to have an effect on the failure response of composite materials. Odom and Adams (1992) have studied the applicability of Hooke’’s law to polymer materials. Using a cylinder-type specimen and tests on different epoxy polymers, they conclude that beyond 0.5% strain, the stress and strain relationship is non-linear and that the linear assumptions of Hooke’’s law are no longer valid. Several micromechanics approaches that account for material non-linearity have been proposed previously, such as those based on using the Aboudi (1989) unit cell method. Recently, Haj-Ali (2009) introduced a new approach using a nine-cell model with cohesive layers to predict failure in composite laminates with a notch. In this approach, cohesive elements are embedded between the fibre-fibre, fibre-matrix, and matrix-matrix subcells. Incorporating cohesive models and micromechanics was also performed in Haj-Ali and El-Hajjar (2003) where failure was propagated along predicted planes concurrently with a micro-mechanics-based approach for the non-linear behaviour of the pultruded composite material. Mechanical and thermal amplification factors were studied by Yudhanto (2005) and Yudhanto et al. (2006) using a 3D finite element approach for various materials and fibre packing geometries using linear properties. Despite increased use of SIFT, a limited number of studies have reported the amplification factors obtained from strain analysis of a micromechanical unit cell with non-linear materials. Closed form solutions have been performed before, but their predictions did not match finite element predictions, especially when non-linear matrix properties are used. An example of one of these

286

K.D. Mishra and R.F. El-Hajjar

approaches used a stress function approach with linear properties to study the effect of fibre volume fractions (FVF) on the stress concentrations in the matrix around the fibres (Upadhyay and Lyons, 1996). This study uses the p-version approach to verify the proposed method. In the p-version implementation, increasing the polynomial degree of elements controls the error of approximation. This makes it convenient to obtain the error estimates in terms of the variable of interest. In this study, mechanical amplification factors are determined using the p-version finite element approach for multi-fibre hexagonal and square packing arrays. The study examined the effects of the FVF of 0.3, 0.5 and 0.7 respectively. The unit cell models are loaded under transverse tension and shear loadings to study the effect of the fibre content on the amplification factors. A Ramberg Osgood approximation for the matrix is proposed to predict the non-linear material effects on the amplification factors.

2

Modelling approach

A multi-fibre cell approach to extract the amplification factors has been performed to study the effects of fibre interactions with the non-linear matrix. The amplification factors are determined by comparing the ratio of the local invariant strain within the micromechanical model to the average strain applied on the model boundary. The SIFT is formulated based on the first strain invariant J1 and the von Mises or equivalent strain, İvm. Failure in the matrix is therefore considered to occur when İvm • İvm––Crit or J1 • J1––Crit where İvm––Crit is the critical von Mises strain and J1––Crit is the critical first strain invariant. The critical values for the matrix are determined from transverse tension and off-axis tests, respectively. The first strain invariant, J1 is determined from the cubic characteristic equation determined from the strain tensor (Ford and Alexander, 1977): İ 3  J1İ 2  J 2 İ  J 3

0

(1)

where J1, is given by: J1

İx  İ y  İz

(2)

where İx, İy and İz are the normal components of the strain vector in the general Cartesian coordinates. The J1 criterion is believed to be applicable to interlaminar failure under the tension load cases where the failure is dominated by volume change in the matrix phase. The von Mises strain criterion is proposed for cases where the matrix failure is dominated by distortion. The second deviatoric strain invariant, J 2c is defined by:

J 2c

1ª 1 İx  İ y 2  İ y  İz 2  İx  İz 2 º¼  İxy2  İ 2yz  İxz2 6¬ 4

(3)

The von Mises (or equivalent) strain is related to the second deviatoric strain invariant, J 2c and the principal strains, İ1, İ2 and İ3 by: İvm

3J 2c

1ª 2 2 2 ¬ İ1  İ2  İ1  İ3  İ2  İ3 º¼ 2

(4)

Non-linear strain invariant failure approach

3

287

Models

A p-version finite element approach is used to model the fibre and matrix multi-cell models. This is accomplished by having the order of shape functions increased by increasing the nodes in the element. The refinement in these models is based on increasing the polynomial degree of the elements, which is directly related to degrees of freedom in the model. The primary advantage of this approach is the ability to check the convergence of the solutions with increasing element order. Increasing the polynomial degree of the elements is related to the error of approximation. These errors are reduced exponentially as the number of degrees of freedom is increased. The global energy norm is calculated for each polynomial order as an indicator of the convergence of the solution. In addition to the global energy norm, the quantity being investigated (J1 or İvm) is also checked for convergence with each polynomial order solution. The boundary and loading conditions used for the transverse tension and shear loadings are shown in Table 1. The selected loading cases were based on the examination of the amplification factors and the controlling cases obtained from Yudhanto et al. (2006). The material modelled for this study was selected because it is widely used and has been extensively examined. The composite used consists of AS4 carbon fibres and the Hexcel produced 3501-6 epoxy. The material properties used for the carbon fibre and epoxy matrix are shown in Tables 2 and 3 respectively. StressCheck (2009) FE software was used to perform the analysis of the micromechanical models. The typical unit cell model used for the hexagonal packing contained 198 elements consisting of 72 hexahedral and 106 pentahedral elements, whereas the square packing contained 256 elements consisting of 128 hexahedral and 128 pentahedral elements. Table 1

Boundary conditions for p-version micromechanical models

Direction of applied displacement

Boundary condition

Transverse tension (x-direction)

İxx = 1.0, İyy = İzz = Ȗxy = Ȗxz = Ȗyz = 0

In-plane shear loading (Txy)

Ȗxy = 1.0, İxx = İyy = İzz = Ȗyz = Ȗzx = 0

Shear loading (Tyz)

Ȗyz = 1.0, İxx = İyy = İzz = Ȗxy = Ȗzx = 0

Shear loading (Tzx)

Ȗzx = 1.0, İxx = İyy = İzz = Ȗxy = Ȗyz = 0

Table 2

Material properties for AS4 fibre

Fibre diameter microns (in)

E1 GPa (psi)

7.11 (0.28 × 10––3)

E2 GPa (psi)

E3 GPa (psi)

v12

v13

v23

G12 GPa (psi)

G23 GPa (psi)

G13 GPa (psi)

234.4 15.2 15.2 0.20 0.013 0.013 27.6 6.89 6.89 (34.0 × 106) (2.2 × 106) (2.2 × 106) (4.0 × 106) (1.0 × 106) (1.0 × 106)

Source: Daniel and Ishai (2006) Table 3

3501-6 epoxy material properties and Ramberg-Osgood parameters

E GPa (psi) 5

6.2 × 10

v

n

ıo MPa (ksi)

D

Max tensile strain

0.35

5

172.4

3/7

2% to 5% (Daniel and Ishai, 2006)

(25)

3.875% (Odom and Adams, 1992)

288

K.D. Mishra and R.F. El-Hajjar

The non-linear material properties are accounted for by using a Ramberg-Osgood approach to capture the non-linear stress strain curve for the 3501-6 epoxy resin. The Ramberg-Osgood plastic deformation model represents the total strain as the sum of the elastic and plastic strains (Dowling, 1999). The polymeric matrix is assumed as a non-linear (J2 –– deformation plasticity) isotropic material, which requires two engineering stiffness constants and three non-linear elastic material parameters for the 1D Ramberg-Osgood stress-strain of the matrix. The matrix non-linear behaviour was characterised by the modulus of elasticity, E, the Poisson’’s ratio, v, and the parameters ıo, D, and n:

İ

ı Įı o § ı ·  E E ¨© ı o ¸¹

n

(5)

The values chosen for the Ramberg-Osgood fitting parameters are shown in Table 3. Figure 1 shows the non-linear material response of the 3501-6 epoxy obtained from experimental data in the literature. The figure also shows the fitting of the non-linear behaviour using the Ramberg-Osgood relations. The stress strain response was selected for the compressive properties because the behaviour captures the greater range of non-linearity not seen in tension tests. It is believed that the tension testing results in premature failure largely as a result of the experimental difficulties in performing tension tests. The models, however, are then compared with tensile failure strains available in the literature. Figure 1

Stress strain response of the Hexcel 3501-6 epoxy and Ramberg-Osgood approximation 30

200

25 150

15

100

10 Odom & Adams, 1992 Proposed Ramberg Osgood Model Linear Behavior

5

0

0

0.02

0.04

Strain

0.06

0.08

50

0

Stress (MPa)

Stress, (ksi)

20

Non-linear strain invariant failure approach

4

289

Results and discussion

Error analysis is performed on the results by examining the strain invariant results and analysing the global energy norm versus the polynomial order. Figure 2 shows the error in the global energy norm versus the element polynomial order for the square and hexagonal models for the transverse tension models. It can be seen that the convergence is achieved rapidly in all cases with the error less than 1%. The strain invariant quantities were also converged in between the different models executed with the increasing polynomial order. Figure 2

Error in global energy norm versus element polynomial order for square and hexagonal models

Error in Global Energy Norm, %

20

15 FVF 0.7 FVF 0.5 FVF 0.3 FVF 0.7 FVF 0.5 FVF 0.3

10

Square Square Square Hex Hex Hex

5

0

0

1

2

3

4

5

Element Polynomial Order

6

7

8

The equivalent strain distribution is shown for the modelled unit cells in Figures 3 and 4 for the hexagonal and square arrays, respectively. The volume fraction for the models shown is for a fibre volume fraction of 0.7. The models were designed such that 16 fibres are used in the square array model and 9 for the hexagonal model. The reasoning behind this approach is to isolate any boundary effects that may occur if a single fibre unit cell is used. The locations examined in this study are the strains in the matrix around the fibre and the strains in the region between the fibres or the inter-fibre spacing. The strain values in Figure 5 are the normalised J1 matrix strains versus the angle around a single fibre for the transverse tension loading case. The results show the sensitivity of the J1 strain invariant to the fibre volume fraction. The invariants for the hexagonal and square arrays converge as the volume fraction is reduced. This result illustrates the dominance of the matrix on the amplification factors for the lower volume fractions. The amplification factor for the equivalent strain is lower than the J1 value and occurs away from the fibre/matrix interface for the case of linear material properties (Figure 6). The strains in the inter-fibre region are shown in Figures 7 and 8 for the transverse tension case. For the equivalent strain, the maximum strains do not occur on the interface, with the results between interface and interfibre becoming even more

290

K.D. Mishra and R.F. El-Hajjar

exaggerated with the higher volume fractions. The results for the square and hexagonal arrays converge for lower values of fibre volume fraction (Vf values of 0.3) illustrating that the amplification factors at that range are largely dominated by the matrix phase and there is less interaction or influence of the fibres. An investigation of the shear effects on the micromodels was also conducted. For the cases considered, the results indicate that the equivalent strains were much more affected by the loading case than the J1 invariant. Figure 9 shows the amplification factors for the models under the different shear loading scenarios. The Tzx loading was found to result in the highest amplification factors representing the case where the fibres are sheared against loading. Figure 3

Equivalent strain distribution of hexagonal array of an AS4/3501-6 composite under transverse tension loading (Vf = 0.7) (see online version for colours)

Figure 4

Equivalent strain distribution of square array of an AS4/3501-6 composite under transverse tension loading (Vf = 0.7) (see online version for colours)

Non-linear strain invariant failure approach Figure 5

291

Normalised J1 in the matrix versus angle around a single fibre for square and hexagonal array models of an AS4/3501-6 composite under transverse tension loading (linear)

Normalized Strain Invariant, J1

3

2.5

FVF 0.5 Sq FVF 0.7 Sq FVF 0.3 Hex FVF 0.5 Hex

1.5

FVF 0.7 Hex

1

0.5

0

Figure 6

FVF 0.3 Sq

2

0

50

100

Angle,T (relative to x-axis, degrees)

150

Normalised equivalent strain in the matrix versus angle around a single fibre for square and hexagonal array models of an AS4/3501-6 composite under transverse tension loading (linear) 2 FVF 0.3 Hex

Normalized Equivalent Strain, Heq

1.8

FVF 0.5 Hex FVF 0.7 Hex

1.6

FVF 0.3 Sq FVF 0.5 Sq FVF 0.7 Sq

1.4 1.2 1 0.8 0.6 0.4 0.2 0

0

50

100

Angle, T, (relative to x-axis, degrees)

150

292 Figure 7

K.D. Mishra and R.F. El-Hajjar Normalised J1 strain between fibres for square and hexagonal array models of an AS4/3501-6 composite under transverse tension loading (linear)

Normalized Strain Invariant, J1

3

2.5

2

1.5

1

0.5

0

0

0.2

0.4

Hex Hex Hex Sq Sq Sq

0.6

0.8

Normalized Inter Fiber Spacing (x-direction)

1

Normalised equivalent strain between fibres for square and hexagonal array models of an AS4/3501-6 composite under transverse tension loading (linear) 2 1.8

Normalized Equivalent Strain, Heq

Figure 8

FVF 0.3 FVF 0.5 FVF 0.7 FVF 0.3 FVF 0.5 FVF 0.7

1.6 1.4 1.2 1 0.8 FVF FVF FVF FVF FVF FVF

0.6 0.4

0.3 0.5 0.7 0.3 0.5 0.7

Hex Hex Hex Sq Sq Sq

0.2 0

0

0.2

0.4

0.6

0.8

Normalized Inter Fiber Spacing (x-direction)

1

Non-linear strain invariant failure approach Figure 9

293

Normalised equivalent strain around a fibre for square and hexagonal array models of an AS4/3501-6 composite under various shear loadings with linear properties (FVF = 0.7)

Normalized Equivalent Strain, Heq

3.5 Hex Txy Hex Tyz Hex Tzx Sq Txy Sq Tyz Sq Tzx

3

2.5

2

1.5

1

0.5

0

0

50

100

150

Angle, T (relative to x-axis, degrees)

A series of non-linear analyses were also performed to determine the effect of the non-linear properties on the use of micromechanical models for analysis with the SIFT approach. The results in Figure 10 indicates that for the range of strain associated with a 5% strain value, the hexagonal models with the highest fibre volume fraction showed the largest sensitivity to non-linear material behaviour of the epoxy. The von Mises strain values did not dominate the behaviour in the linear response under transverse tension, but the values exceed the behaviour when nonlinear behaviour is concerned. The linear predictions are shown to highlight the effect of the nonlinear material behaviour for the higher strain values. The prediction for the hexagonal and square models yields similar results at lower volume fractions, indicating the effect of the packing arrangement to be diminished for those cases. For the transverse tension case, Figure 11 shows the predictions from non-linear and non-linear analysis cases. The linear results together with non-linear analysis for 2% and 3.5% maximum average strain are shown related to values expected for failure of the epoxy matrix. The equivalent strain is seen to peak at a different location around the fibre to that where the maximum J1 value occurred in the linear case. This result is attributed to the large shear type deformation occurring near the fibres during the transverse tension loading when non-linear properties are included.

294

K.D. Mishra and R.F. El-Hajjar

Figure 10

Non-linear stress-strain predictions in hexagonal and square multi-fibre models of an AS4/3501-6 composite under transverse tension loading 80 70 FVF 0.3 Hex NL FVF 0.5 Hex NL FVF 0.7 Hex NL FVF 0.3 Sq NL FVF 0.5 Sq NL FVF 0.7 Sq NL Linear Predictions

50

400

40 30

200

Stress (MPa)

Stress (ksi)

60

20 10 0

0

0.01

0.02

Strain

0.03

0 0.05

0.04

Effect of non-linear analyses on the J1 strain invariant and equivalent strain for hexagonal packing arrangements under transverse tension loading (Vf = 0.7)

Figure 11

5

Normalized Strain Invariant

4.5 4

J1 FVF 0.7 Hex (Linear) J1 FVF 0.7 Hex (Nonlinear, Max Strain 2%) J1 FVF 0.7 Hex (Nonlinear, Max Strain 3.5%) EEQ FVF 0.7 Hex (Linear)

3.5

EEQ FVF 0.7 Hex (Nonlinear, Max Strain 2%) EEQ FVF 0.7 Hex (Nonlinear, Max Strain 3.5%)

3 2.5 2 1.5 1 0.5 0

5

0

50

100

Angle,T (relative to x-axis, degrees)

150

Conclusions

A series of linear and non-linear modelling approaches were used to study the strain amplification factors due to mechanical loading in square and hexagonal packing arrays. The fully parameterised models were used to study the effects of the non-linear material properties and the fibre volume fraction on the overall amplification factors. The

Non-linear strain invariant failure approach

295

approach used a p-version finite element approach to model the multi-cell array and to study the errors of approximation associated with the analysis. The models were capable of capturing the non-linear matrix behaviour by representing the matrix with a non-linear Ramberg Osgood relationship. For the transverse tension loading, the non-linear material behaviour of the epoxy matrix was found to significantly affect the von Mises strain around the fibre for the range of strain and loading conditions considered. However, it maybe that peaks occurring away from the fibres, may make the material non-linearity less of an issue. The results from the non-linear analysis show the effects of material non-linearity on the equivalent strain exceeding those seen on the J1 in determining the amplification factor. The results show that non-linear analysis may be necessary when considering failure analysis using the amplification factor approach in the SIFT. In addition, the study highlights the need to carefully consider the fibre volume fraction in the various areas of the manufactured components. The results of this study point to the need to understand how manufacturing defects such as resin migration or fibre waviness may affect the mechanical amplification factors. This conclusion is reached because these defects may lead to localised areas where the fibre volume fraction approaches the theoretical limits. The effect of strains due to curing will also need to be considered for complete consideration of thermal related amplification factors. The p-version finite element method offers an accurate way to implement fibres having a large aspect ratio than those available in the usual isoparametric approach within a non-linear matrix. A study of the errors of approximation associated with discretisation shows prompt convergence and the models were capable of capturing the non-linear matrix behaviour. Validation with experimental data at the fibre level will likely be challenging due to the scale involved; however, the results included can be used as a guide to understand how manufacturability of a composite can affect the interactions of stresses and strain at this level. Future studies will also need to investigate implementation of the unit cells and analysis under mixed mode loadings.

References Aboudi, J. (1989) ‘‘Micromechanical analysis of composites by the method of cells’’, Appl. Mech. Rev. Vol. 42, No. 7, pp.193––221. Daniel, I. and Ishai, O. (2006) Engineering Mechanics of Composite Materials, 2nd ed., Oxford University Press, New York, NY. Dowling, N.E. (1999) Mechanical Behavior of Materials, 2nd ed., Prentice Hall, New Jersey. Ford, H. and Alexander, J. (1977) Advanced Mechanics of Materials, 2nd ed., Ellis Horwood, Chichester, UK. Gosse, J. (2001) ‘‘Strain invariant failure criteria for polymers in composite materials’’, 42nd AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics, and Materials Conference, April, pp.16––19, Seattle, WA. Gosse, J.H. and Christensen, S. (2001) ‘‘Strain invariant failure criteria for polymers in composite materials’’, AIAA-2001-1184. Haj-Ali, R. (2009) ‘‘Cohesive micromechanics: a new approach for progressive damage modeling in laminated composites’’, Int. J. Dam. Mech., Vol. 18, No. 8, pp.691––719. Haj-Ali, R. and El-Hajjar, R. (2003) ‘‘Crack propagation analysis of mode-I fracture in pultruded composites using micromechanical constitutive models’’, Mech. Mater. Vol. 35, No. 9, pp.885––902.

296

K.D. Mishra and R.F. El-Hajjar

Li, R., Kelly, D. and Ness, R. (2003) ‘‘Application of a first invariant strain criterion for matrix failure in composite materials’’, J. Compos. Mater., Vol. 37, No. 22, pp.1977––2000. Odom, E.M. and Adams, D.F. (1992) ‘‘An investigation of the isotropy of epoxy polymers’’, J. Mat. Res., Vol. 7, No. 12, pp.3352––3358. Soden, P.D., Kaddour, A.S. and Hinton, M.J. (2004) ‘‘Recommendations for designers and researchers resulting from the world-wide failure exercise’’, Compos. Sci. Technol., Vol. 64, Nos. 3––4, pp.589––604. StressCheck (2009) Release 9, Engineering Software Research & Development, Inc., St. Louis, MO, USA. Tay, T.E., Tan, S.H., Tan, V.B. and Gosse, J.H. (2005) ‘‘Damage progression by the element-failure method (EFM) and strain invariant failure theory (SIFT)’’, Comp. Sci. Tech., Vol. 65, No. 6, pp.935––944. Upadhyay, P.C. and Lyons, D.W. (1996) ‘‘On the stress concentrations in the matrix around the fibers in laminae’’, J. Reinf. Plas. Comp., Vol. 15, No. 8, pp.818––826. Yudhanto, A. (2005) ‘‘Effects of micromechanical factors in strain invariant failure theory for composites’’, thesis, Department of Mechanical Engineering, National University of Singapore, Singapore. Yudhanto, A., Tay, T.E. and Tan, V.B.C. (2006) ‘‘Micromechanical characterization parameters for a new failure criterion for composite structures’’, Key Eng. Mater., Vols. 306––308, pp.781––786.