- Email: [email protected]
Impact response of laminated composite plates: Prediction and verification
Composite Structures 28 (1994) 61-72
' ~ ; ' "': L
Impact response of laminated composite plates: Prediction and verification H. V. Lakshminarayana, R. Boukhili & R. Gauvin Department of Mechanical Engineering, Ecole Polytechnique de Montreal, Canada Methods and procedures for predicting the impact response of laminated composite plates using a commercial finite element system are described. Results of element evaluation, procedure verification and a correlation study are presented and discussed. The need for a hybrid experimental-numerical approach and combined geometric and material nonlinear finite element analysis is identified. A methodology for the prediction of delamination (onset and growth) is outlined.
1 INTRODUCTION
while few studies have included nonlinear effects. Almost all of them make use of special-purpose programs. Such programs, valuable in their own right, are difficult to access, are not well documented and hence are difficult to use by the practitioners. (2) Assessment of convergence and accuracy of the deflection history, stress distributions and damage zones calculated by the Finite Element Method (FEM) have not been given explicit attention. (3) Experimental investigations that document the response recorded during the impact test are rather limited. The majority of them provide a qualitative/quantitative description of the accumulated damage but not its growth. This information is vital to perform correlation studies. (4) Few investigators have considered full-scale components. The majority of the reported studies use generic structures such as beams, circular plates, square plates and cylindrical panels. There is a real danger in extrapolating conclusions drawn from studies on test specimens to real life components. It is more appropriate to use commercial FEM systems. They are widely distributed, well documented and user friendly. However, there is a need to verify the accuracy of material models, finite elements and analysis procedures in such systems before using them for the intended application, namely, prediction of the impact response of composite plates. The predictability issue itself demands a correlation study. This, in fact, is the aim and scope of the present study. A brief description of the specific methods and procedures used is given in the next section. In
Laminated composite plates are easily damaged by impacts, especially those normal to the plane of the laminate. Systematic study of such problems can be divided into three distinct areas: structural mechanics, damage mechanics and residual strength prediction. A study in that order will provide a mechanistic basis for both design and assessment of damage tolerance. Reliable and accurate prediction of the impact response of multilayered anisotropic plates covering a wide range of parameters is the focus of this study. Impact response means contact force history, deformation history, surface strain history, stress distribution across the laminate thickness (including interlaminar stresses), nature and extent of damage and stiffness and strength loss associated with that damage. Parameters significantly influencing the impact response include impact velocity/energy, material system, ply orientation and stacking sequence, plate geometry (shape, size and thickness), wall construction (solid laminate, integrally stiffened, and sandwich), support conditions, pre-stress state and initial curvature (curved panels). A review of previous work is available in a monograph L2 and two review papers, 3,4 which provide a background to the present study. A critical assessment of this vast literature, with particular reference to the focus of the present study, can be summarized as follows. (1) Finite element analysis of linear and transient structural behaviour has been the focus of intense research, 61
Composite Structures 0263-8223/94/S07.00 © 1994 Elsevier Science Limited, England. Printed in Great Britain
62
H. V. Lakshminarayana, R. Boukhili, R. Gauvin
succeeding sections, results of element evaluation, procedure verification and correlation study are presented and discussed. The presentation concludes by identifying directions for further work.
2 M E T H O D S AND P R O C E D U R E S
Numerical results for this study were generated using ABAQUS -- a general-purpose finite element code with emphasis on nonlinear applications. ~ This program is capable of modelling multilayered anisotropic materials. It provides elements suitable for dynamic analysis of composite plates and shells taking into account bendingmembrane coupling and transverse shear deformation effects. Among these, the S8R, an isoparametric quadrilateral plate/shell element, is employed in the present study. This element has eight nodes and six engineering degrees of freedom at each node. The user can specify within each element an arbitrary number of layers, each with its own thickness, ply orientation and orthotropic elastic properties. The formulation uses a value of 5/6 for the shear correction factors as default. However, the user has the option to use any other value through independent input of transverse shear stiffness. The element output includes membrane stress resultants, bending stress resultants and transverse shear-stress resultants, either at the nodes or at the four integration points. The user can also request ply-byply stresses at the integration points and at a maximum of 3 section points within each ply. It should be noted that the formulation of S8R does not ensure the continuity of interlaminar normal and shear stresses at the interfaces between plies. ABAQUS provides two procedures for calculating the response of structures subjected to impulsive loads. They are modal analysis and dynamic analysis by direct implicit integration of the equations of motion. The modal analysis is limited to linear transient behaviour. Its convergence and accuracy is dependent on the number of natural modes considered in the analysis. The pulse shape and the duration (t¢~) for which the contact force acts (in comparison with the period T~ of the fundamental mode of free vibration of the structure) strongly influence the response calculated by the modal method. For complex pulse shapes and for (t,,/Ti)'~ 1, a very large number of modes is required for convergence. Unfortunately, accurate determination of the eigenvalues and eigenfunctions associated with higher modes
places very heavy demands on computational resources by the FEM, Incidentally, eigenvalue extraction in ABAQUS is done using the subspace iteration method. The transient response calculation by the direct integration of the equations of motion is applicable to linear as well as nonlinear structural behaviour. When applying this step-by-step method, the time step At should be selected with caution, because a system of nonlinear algebraic equations must be solved at each time increment. This is done in ABAQUS iteratively by using Newton's method or, if preferred, the quasi-Newton's method. This time stepping, nonlinear equation-solving procedure is computationally expensive. The principal advantage of this procedure is that it is unconditionally stable, which means that there is no mathematical limit on the size of the time increment that can be used. In practice, At should be small enough to adequately define the history of excitation, its value being chosen on the basis of the shortest period which corresponds to the highest natural mode likely to contribute to the response. A recommended methodology(' for determining the contact force history is to measure the local contact stiffness of the composite plate in static tests and use this in conjunction with a finite element model. This approach, however, could not be considered truly predictive since it requires fabrication of the plate and static indentation tests for every impactor under consideration. Alternatively, it can be measured during impact tests using suitable instrumentation and data acquisition system. This approach is applicable to the nonlinear response also. The nature and extent of impact-induced damage is estimated by first calculating the dyamic stresses and their spatial distribution in the laminate and incorporating these in appropriate failure criteria. Multiple matrix cracks, delaminations and fiber breaks are the failure modes observed after impact tests. These failure modes and complex interactions between them complicate the prediction of damage. The tensor polynomial failure criterion proposed by Tsai and Wu 7 is employed in the present study to calculate the failure index (FI) given by F I = k ' l o , + F2o~ + F~o~ + F11o~ + F2:o! + F(~,o~ + _F~2o 102
where k~j are the strength tensors; 7 o~, o 2 are the lamina stresses in the fiber direction and transverse direction respectively; and o~, is the in-plane
Impact responseof laminated compositeplates shear stress. Loci of points at which FI = 1 define the damage zone. The failure modes are identified using the maximum stress criteria. 7 The stiffness and strength loss due to impactinduced damage are not at present truly predictable. A suggested approach involves flexure tests. ~ Evaluation of the accuracy of the S8R element for the analysis of composite plates in general is presented in the next section. We identify three distinct procedures: (1) modal analysis (procedure # 1); (2) linear and transient response analysis (TRA) (procedure #2); and (3) nonlinear and transient response analysis (NLTRA) (procedure # 3). The convergence and accuracy of each one of these is verified in a section entitled procedure verification. Finally, a critical assessment of their predictability with particular reference to the impact response of laminated composite plates is presented in Section 5.
3 ELEMENT EVALUATION Application of the full set of test problems proposed in Ref. 9 to the S8R element is not pre-
I I
6 _
W= 5111111
Ill Ill
W(O,O) E,h~ 10s Pa" I
t
4 PROCEDURE VERIFICATION
CPT Solution
1
-1
10 Fig. 1.
sented here. The results showed that a particular problem, a homogeneous, anisotropic, clamped, square plate under uniform pressure, provided an intensive measure of element performance. This problem, shown as the insert in Fig. 1, was chosen to evaluate the combined effect of material anisotropy and shear deformation on the accuracy of the S8R element. Numerical results were obtained for a unidirectional laminate, made of a highmodulus graphite/epoxy composite (EI/E 2= 40, E2 = 5.17 GPa (0.75 x 106 psi), Gi2--- 3.10 GPa (0"45 x 10 ~' psi), G23/GI2=0"8, vl2=0.25 ), for a ply orientation of 45 and for various values of the thickness ratio (h/a). Computed results are compared with the converged solutions given in Ref. 9. The agreement is very good for displacements as well as stress resultants provided (h/a)> 0.01. The effect of h/a on the predicted central deflection is shown in Fig. 1, which also has the classical plate theory (CPT) solution for comparison. For a given thickness ratio, the inaccuracy associated with the omission of transverse shear deformation effects is obtained from this figure. Obviously, the accuracy of the S8R element deteriorates for (hi a)<0"01, indicating that the formulation is not free from the so-called shear-locking phenomenon. There is scope to improve the performance of the SSR element for applications to thin and highly anisotropic situations. An approach to accomplish this is demonstrated in Ref. 10.
II
IIIIII
0
63
-2
10 h/a
-3
10
Effect of (h/a) on the central deflection predicted using the S8R element (8 x 8 mesh, whole plate).
The chosen elements and procedures should be capable of providing converged and accurate solutions to standard test problems for which alternative solutions are available for comparison. The specific problem used in this study is a crossply laminated, graphite/epoxy composite, square plate with all edges simply supported and subjected to impact at its center, which was investigated by Sun and Chen. ~ The geometry, ply orientation/stacking sequence, finite element mesh and contact force history (corresponding to impact by a 12.7 mm diameter steel ball at velocity V0 -- 3 m/s) used in the computations are shown in Fig. 2. Representative displacement time histories calculated using procedure # 2 (linear and transient response by direct implicit integration) are given in Fig. 3. These results are in excellent agreement with those reported in Ref. 11, in spite of the larger time step (6t=0.05 ms) used. No significant
64
H. V. Lakshminarayana, R. Boukhili, R. Gauvin
P(t)
h
w(t) P(t)
a
Pm a x
0°
h = 5.1 m m a = 203.2 mm
\
0° 90 °
Pr.~x= 27.4 Kg to = 0.2 ms
I I I
90'
/
I
0° 0 °
i
',\
t 1 = 0.25 t o
90 °
i
0 ° 90 °
t t
t
1
0°
0
Contact Force History
Ply Orientation/Stacking S e q u e n c e
/;
j/-
/ i
J
/
J
Finite Element Model (Quarter Plate)
Fig. 2.
Test problem for procedure verification.
improvement is noticed using a finer mesh and a smaller time step, indicating that a converged and accurate solution has been obtained. Figure 3 also indicates that two pulse shapes (triangular and half sine-wave) produce almost identical results. Since the half sine-wave pulse can be represented by an analytical expression, it is ideally suited to create benchmarks for the impact response of composites. For the same problem, Fig. 4 shows displacement-time histories calculated using procedure # 1 (modal analysis). There is no sign of convergence as the number of modes used is increased from 10 to 15. These results do not show any comparison with the reference solution in Ref. 11. Not only is the displacement history different, the maximum amplitudes differ by one order of
magnitude. This observation is substantiated by the findings of a recent round robin study.~2 It is noted that whenever the duration of impact is very small (in comparison with the period of the fundamental mode of free vibration of the target), the number of modes to be considered for convergence may be so large as to be computationally inefficient. It appears that reliable prediction of the impact response is outside the domain of the modal analysis. For the problem specified in Fig. 2, a comparison of the displacement histories calculated using procedure # 3 (nonlinear and transient response by direct implicit integration) and procedure # 2 is made in Fig. 5. The two solutions are in fact identical because at this low load level the resulting deflections are very small (in comparison
Impact response of laminated composite plates 0.24
0.22
D
0.22
0.2 0.18
0.2
0.16
0.18
0.14
0.16
~
0.12
¢"- 0.14
d
o.1
"*"= 0.12 0 0.1 v0
r
0.08 I~1 0.06 5 ¢G)
II
0.04
I~:
0.02
c-
E
0GI -0,02
0
Pmax= 27.4 Kg to= 0.2 ms t, = 0.25 t o
(3. .~_ -o.04
-o.os
O
006
0.08
0.04
E 0.02
o,
0
65
P(t
oo -0.02
=~
a
-0.04
\,
-0.08 -0.1 -0.12 0
11
I
0.5
1
0
; f/
i
-o ~ ~_o
' i (;t~
0
0.5
-0.06
t
-0,08
0
I
I, t¢ l
!1 Prnax= 27.4 Kg ' 1 1, to= 0"2,ms
-0.1
1.5
2
1
1.5
Time t(miliseconds)
Time t(miliseconds) Fig. 3.
\ fltl=
Central displacement history calculated using procedure # 2.
2
2
1.8
1.8
N u m b e r of m o d e s = 1 0
1.6
1.6
1.4
1.4
1.2
1.2 t-
1
1
d o.8
"~. o.e o <:5 o.6
0
v
v
0.4
II
0.2
15
I~
o,
•~ e-
!P(t)
c
E
0.6 0.4
o.2 0
Pr..x P(t)
E
8
o
e~
h
e~
w,t,
._~ a
P(t) = P S i n
i
°u
a
(nt/to)
Pmax= 27.4 Kg
0
to
0
t = 0.2 ms o
I
2
0
I
I
I
1
2
3
1
Time t(miliseconds) Fig. 4.
Time t(miliseconds)
Central displacement history calculated using procedure # 1.
66
H. V. Lakshminarayana, R. Boukhili, R. Gauvin 0.24
0.24
0.22
Procedure # 2 6t = 0.05 ms
0.22 0.2 0.18
0.18
7
0.16 e.-
0.16
0.14
~
0.14
C~ 0.12 0 0.1
d
o.12
II
0.1
0.08
0.08
15 ¢.. 006 E
0.04
0
0.02
Procedure # 3 6t = 0.05 ms
0.2
P(t)
/
to= 0.2 ms
_ g 0.04
Pmax
0.02
..........
~.
~. 0 .oo E3 -o.o2
o?,,J,,2t
-0.04 -0.06 -0.08
0
-0.1 0
f
I
I
0.5
1
Time Fig. 5.
"o
I
1.5
t(miliseconds)
o
-0.02 P(t) = Pmeeqin (rrt/t o)
-0.04 --
Pm~ = 27.4 Kg VO 3 m/s O= 0.2 ms
-0.06 I -0.08
-0.1 J 0
I
0.5
l
1
4
1.5
T i m e t(miliseconds)
Comparison of displacement-time histories calculated by procedures # 2 and # 3 (impact velocity = 3 m/s).
with the thickness h) and hence nonlinear effects do not show up. Impact velocity is a very important parameter controlling the response. Figures 6 and 7 provide a comparison of displacement time histories calculated using procedure # 2 and procedure # 3 at V~=10 m/s and V0=30 m/s, respectively. Obviously, the responses predicted by the two procedures do not agree with one another. Both the displacement-time history as well as the maximum amplitudes predicted differ rather significantly. Basically, the load levels are such that resulting deflections are of the order of the plate thickness and therefore procedure # 3, which includes nonlinear effects, is more appropriate. Unfortunately, the accuracy of the numerical solutions presented in Figs 6 and 7 could not be verified due to the nonavailability of reference solutions. Procedure # 3 is therefore more appropriate for numerical solutions of impact tests. In fact, for complex structures, this procedure may provide more cost-effective information than experimentation. However, there is a need to further validate the predictability of this procedure for the impact
response. This aspect is addressed in the next section. 5 CORRELATION STUDY The predictability aspect of the procedures verified in the previous section, with particular reference to the impact response of composite plates is assessed using a bench mark. Bench marks are fully specified and standard problems, which resemble instances found in practical applications, and for which reference solutions have been obtained using both analytical/numerical and experimental methods. The bench mark chosen for this study was created by Aggour and Sun. j3 It consists of a cross-ply laminated, E-glass/epoxy composite (Ej=38"6 GPa (5"6x106 psi), E2=I0.34 GPa (1.5×106 psi), G2=4"14 GPa (0.6x106 psi), v,2 =0.25) square plate with all edges clamped. The geometric parameters, ply orientation/stacking sequence, and the finite element discretization used in the computations are shown in Fig. 8. The contact force history corresponds to an impact by
Impact response of laminated composite plates
67
1 0.9
0.9 f
Procedure #2 6t = 0.05 ms
0.8
0.7
0.7
0.6
0.6
0.5
~, ~d
t-
0.5
0.4
0.4
0.3
C) 0.3 O 0.2
0.2
II 0.1
0.1 0
t
e(9
~
-0.1
~
-0.2
E -0.1 (9 -0.2
~.
-0.3
~ . -0.3
i'~
-(I.4
u) ~,~ -0.4
-0.5
-0.5
-0.6
-0.6
-0.7
-0.8
Procedure # 3 5t = 0.05 ms
0.8
0
0.5
1
~W(t) ~ Ih P(t)= PmaxSin(~t/t 0) V 0 = 10 m / s P m ~ = 9 1 . 2 Kg
-0.7
to
o
~-~
1.5
-0.8 0
2
t 0= 0.2
I
I
1
0.5
Time t(miliseconds) Fig. 6.
ms
I
1.5
2
Time t(miliseconds)
Comparison of displacement-time histories calculated by procedures # 2 and # 3 (impact velocity = 10 m/s).
4
2
Procedure #2 : LTRA
Procedure # 3 : NLTRA
1.8
6t = 0.05 ms
1.6 1.4 1.2
J::
¢-
1 0.8
c:) 0.6
o 0
0
II
O.4
II 0.2
,4..I
t5
0[
(9
E
P(t)
1
o,
c (9 -0.2
E
(9 -0.4 o ¢0 -0,6
o
P(t) = Pm=Sin (nl/t o ) Prr,= = 364.8 Kg
-2
-1.2
Vo= 30 m/s t -- 0.2 ms
-1.4
Vo= 30 m/s Pro= = 364.8 Kg to= 0.2 ms
0
-3 0
0.5
1
Time t(miliseconds) Fig. 7.
a•-~ -0.8 -1
1.5
-1.6
0.5
1
1.5
2
Time t(miliseconds)
Comparison of displacement-time histories calculated using procedures # 2 and # 3 (impact velocity = 30 m/s).
68
H.V. Lakshminarayana, R. Boukhili, R. Gauvin z [ P(t)
r
--.
~
X
a P(t)
P ,x ....
/
/
\ /
'\
'\ ',,
I
ho o
/iiiiiiiiiiiiiiii
9o!
I!iii!i!~iiii!iiii!iiiiiiliiiiiiiiiii!iiiiii!i ilili ~i!i!i!i!~!i!iiiii~ii~iililiiililili9o i!iiiiiJo
\\\
~_t
t °°
0 °
t Contact Force Histor~ P(t) = Pm~xSin(nUt o ) P~x = 310.1 Kg to= 0.25 ms
Ply Orientation/Stacking Sequence a = 139.7 mm h=4.1
J
mm I
A
B
E
H
Finite Element Model (Quater Plate).
Fig. 8. Testproblemfor correlation study.
a steel cylinder (diameter 9.5 mm and length 25.4 mm) at a velocity E, = 22.6 m/s. The calculated central displacement history using procedure # 3 is presented in Fig. 9. Results using procedure # 2 and test data taken from Ref. 13 are also included in the same figure to enable a three-way correlation. Results obtained using linear finite element analysis by Aggour and Sun ~3 are in close agreement with those obtained in the present study using procedure # 2. Furthermore, the responses calculated using procedure # 3 and procedure # 2 are identical, indicating that nonlinear effects are negligible. It is gratifying to note that computed results closely follow the central deflection measured during an impact test for t-<0"5t0, where to is the duration of impact. How-
ever, for higher values of t there is considerable difference between the computed and measured results. Therefore, to achieve a satisfactory correlation between the two sets of results over a larger duration, further refinement of the procedures used appears to be essential. Basically, the effects of impact-induced damage has to be included in the analysis. More specifically, a combined geometric and material nonlinear dynamic finite element analysis is indispensable for this purpose. However, this is more easily said than actually performed, because a priori knowledge of the nature and extent of impact-induced damage and accurate constitutive modelling of failed plies and damaged laminates will be required after each time step. At present,
Impact response of laminated compositeplates
69
1 Procedure # 2
0.9
6t = 0.025 ms
0.8
0.6
Testdata
0.4
~,
0.4
0.3
~
0.3
oc
o.1
g
~ -0.1
0{
~ -0.1
~ -0.2
~ -0.2
~
-0.3
~
~
-o.g
i~ -0.4
-0.5
-0.5
-0.6
-0.6
-0.7
-0.7
-0.8
-0.3
-0.8 0.5
l 0
Time t(miliseconds) Fig. 9.
Comparison of measured and predicted central displacement histories.
1.8
Procedure # 2
l'k 1.6
1.4
1.61.4-
VO= 40.0 m/s []
Vo-- 2226 m/s
¢...
0.6
Vo= 40.0 m/s
1
Vo= 30.0 m/s ©
0.8
Vo= 22•6 m/s
0,6 oo-
\
ji
0.4
/'//
0.2
II -0.2
i
.,.., e.- -0.2 1
(1) -0.4
E -o.4 ~
-0.6
\'~ "t~r~ \\
(:~ -0.6
-o8
~.
-0.8
-1
i~
-1
-1.2
-1.2
-1.4
~N
-1.6
-1.4 -1.6
-1.8 0.5
Time t(miliseconds) Fig. 10.
Procedure # 3
1.2
Vo= 30.0 m/s ©
0.8
0 ~
0.5
Time t(miliseconds)
1.8
!
]
0.2
I
0.1
c5 0.4
Studyl
Testdata
0.5
0.2
[
Prese~
0.7 0.6
0.5
c~
#3
6t = 0.05 ms
0.8
Present Study
0.7
Procedure
0.9
-1.8 0.5
Time t(miliseconds)
Comparison of predicted displacement-time histories by linear and nonlinear analysis procedures.
70
H. V. Lakshminarayana, R. Boukhili, R. Gauvin 1.6
Procedure #2 V = 40 m/s
1.4
0
t
= 0~25 ms
t=
Procedure #3 V = 40 m/s
14 ,2
o2oms
0
[
i
I t = 0z15 ms
x::
[t = o
, ~f i
ms/
/t = 0.10 msl
"~'I~ 0.8! ..
0.6 I
.6
0.4
0,4
0.2 O"I
-02
0.1
0.2
0.3
0.4
05
x/a Fig. 11.
-02
0
0.1
0.2
0.3
0.4
0.5
x/a
Comparisonof predicted displacementdistributionsby linear and nonlinearanalysisprocedures.
the former cannot be truly calculated while information on the latter is not available. The progressive failure finite element analysis outlined is in fact a topic for further research.~4 Reliable prediction of displacement history as a function of impact velocity (Fig. 10) constitutes only a fraction of the impact response story. Spatial distribution of displacement (Fig. 11 ), stress distribution across the laminate thickness (Fig. 12) and damage growth (Fig. 13) are in fact even more important. The F E M in general and A B A Q U S in particular is able to provide such results. Indeed a major problem is to find experimental techniques that can provide data that are equally detailed to validate the results presented here.
6 CONCLUDING REMARKS A hybrid experimental-numerical approach is necessary to predict the impact response of laminated composite plates covering a wide range of the parameters involved. Experimental determination of contact force history is essential if the structural behaviour is nonlinear. Nondestructive test methods are indispensable for characterisation of damage. Post-impact tests are needed to determine residual stiffness and strength. The combined geometric and material nonlinear finite element analysis capability required for this purpose is available in commercial finite element systems. However, accurate constitutive models
for composite laminates with multiple matrix cracks, fiber breaks and delaminations are not yet available. Delaminations are often the primary, life-limiting failure modes. The task of developing methods and data for predicting the onset and growth of delamination due to impact was not considered in the present study. Accurate evaluation of interlaminar normal and shear stresses a verified interface criteria and test methods to measure material properties associated with such a criterion are the prerequisites for prediction of delamination onset. It is to be noted that an appropriate finite element model for this purpose should be based on a lamination theory that enforces the continuity of interlaminar normal and shear stresses at ply interfaces,~5 and these are not yet implemented in commercial FEM systems such as ABAQUS. For the prediction of delamination growth, the fracture mechanics concept ~6is indispensable. This procedure involves calculation of energy release rates associated with delamination growth, the use of a mixed-mode fracture criterion and test methods to measure interlaminar fracture toughness. A drastically different finite element modelling approach is necessary for the numerical determination of energy release rates. These are topics deserving further research. This study sets the stage for confident application of a commercial FEM system ABAQUS for numerical simulation of impact tests conducted on composite plate and shell structures.
Impact response of laminated compositeplates Procedure #3
71
]
z/h
at t = 0 . 2 ms, V o = 4 0 m / s __
X = 1022.3 MPa
0 °
0'
•i•i•i:!•i•i:•i•i:•i•!i:•!i:•i:•i:•i•i:•i•i:!•i:•i:•i
90'
?::::..:.v....:.v..::.v.~..`..v:.......v.....v.v...v.`....:..:......`.:.v..::...
I1
~ : ~ : ~ : ~ : ~ : ~ : ~
90'
I
~......................................................... :~:~:~:~:~::~:~:~:~
90'
~::..`.:......`.:.v`.::~.:~.....:.v:.....v.....v`.`.`.`.:~.~.`.:::......::.v..`.`.`.. i..:.v..`.:....::.`....:::......`.:............:.v...v....:•..`.........v.....v:...v::.v :::~::~::~:~:~:~:~::~:~:~:~::~:~:~::~::~:~:~::~:~:~:~:~ 90' O'
--
O*
i
I
I
I
I
-1
1
o
01/X O° 0°
Y = 27.3 MPa
•!i•i•i•i:•i•i•i•i•i•ii:•i•i:i:•i•i•i•i•i:•i•ii:i:•i•i!•i!:•i•ii(.:.~...~.~.~.~.~::..~..~..~...~..~...~......`..~.~.::..`.9o' ~.`.`.~`~...~.`..~ •i•!i•!:!•i:!•:!•i•:i•!:i•!i•!h•i!•i:!•i•!i•!i:!•i•!:!•i 9O'
[
~.~.):.~.~.~.)):~)~.)))~.~.~`~.~.~.33)&:.:.:.~.))3~`~.~.)~.)~.)90' )~.:.:.:.~.~.333~:.~ "
0' 0° q
-10
lO
~/y z/h S = 40.9 MPa
•i•i:•i•i•i•i•i:•i•ii:•!•i!•ih:i:i:•i!•i•i•i:•i!:i:!•i!
9o'
:......:.....v:.v...v.......v~..`:..:.v...~`v.......v.v......:..::..:.v....`....`....: ~0' 0°
.1
o
otis Fig. 12.
Stress distribution across laminate thickness.
ACKNOWLEDGEMENTS t =0.05 ms ~]
t = 0.15 ms
~
t=O.2Oms
900
i 0°
This research was supported by the Natural Sciences and Engineering Research Council of Canada and la Direction de la Recherche de l'Ecole Polytechnique de Montr6al.
REFERENCES
V0= 40 m/s Procedure #3 Fig. 13.
Growth of impact damage.
1. Greszczuk, L. B., Damage in composite materials due to low-velocity impact. In h n p a c t D y n a m i c s , ed. J. A. Zukas et al. John Wiley & Sons, New York, USA, 1982. 2. Zukas, J. A., Numerical simulation of impact phenomena. In I m p a c t D y n a m i c s , ed. J. A. Zukas et al. John Wiley & Sons, New York, USA, 1982.
72
H. V. Lakshminarayana, R. Boukhili, R. Gauvin
3. Abrate, S., Impact on laminated composite material. Applied Mechanics Reviews, 44 ( 1991 ) 155-90. 4. Cantwell, W. J. & Morton, J., The impact resistance of composite materials -- Review. Int. J. Composites, 22 ( 1991 ) 347-62. 5. ABAQUS -A general purpose finite element code, users manual, version 4-7. Hibbit, Carlson and Sorensen Inc.. Providence, RI, USA. 6. Yan, T. M. & Sun, C. T., Use of statistical indentation laws in the impact analysis of laminated composite plates. J. Appl. Mech., 52 (1985)6-12. 7. S. W. Tsai, & Hahn, H. T., brtroduction to Composite Materials. Technomic, Lancaster, PA, USA, 1982. 8. Nemes, J. A. & Simmonds, K. E., Low-velocity impact response of foam-core sandwich composites. J. C~mp. Mat., 26 (1992) 500-19. 9. kakshminarayana. H. V. & Sridhara Murthy, S., A shear - - flexible triangular finite element model for laminated composite plates. Int. J. Nttmer. Meth. In Eng., 20 11984) 591-623. 10. Lakshminarayana, H. V. & Kailash, K., A shear -deformable curved shell element of quadrilateral shape.
Computers and Structures, 33 (1989) 987-1001. 11. Sun, C. T. & C h e n . J. K., On the impact of initially stressed composite laminates. J. (~mp. Materials. 19 (1985) 490-504. 12. Ewins, D. J. & lmregun, M., A survey to assess structural dynamic response prediction capabilities. In Quulit3, Assurance in FEM Technology, ed. J. Robinson. Robinson Associates, UK, 1987. 13. Aggour, H. & Sun, C. T., Finite element analysis of laminated composite plate subjected to circularly distributed central impact loading. Computers & Structures, 28 (1988) 729-36. 14. Yener, M. & Wolcott, E., Damage assessment analysis of composite pressure vessels subjected to random impact loading. J. Pressure Vessel Technolog3', 111 (1989) 124-9. 15. kec, C. Y. & kiu, D., lnterlaminar stress continuity theory for laminated composite analysis. AIAA ,/.. 29 1 9 9 ( ) 2010-12. 16. Friedrich. K. (ed.). Application of kracture Mechanics to ('omposite Materials. Elsevier Applied Science Publishers, Amsterdam, The Netherlands, 1989.
' ~ ; ' "': L
Impact response of laminated composite plates: Prediction and verification H. V. Lakshminarayana, R. Boukhili & R. Gauvin Department of Mechanical Engineering, Ecole Polytechnique de Montreal, Canada Methods and procedures for predicting the impact response of laminated composite plates using a commercial finite element system are described. Results of element evaluation, procedure verification and a correlation study are presented and discussed. The need for a hybrid experimental-numerical approach and combined geometric and material nonlinear finite element analysis is identified. A methodology for the prediction of delamination (onset and growth) is outlined.
1 INTRODUCTION
while few studies have included nonlinear effects. Almost all of them make use of special-purpose programs. Such programs, valuable in their own right, are difficult to access, are not well documented and hence are difficult to use by the practitioners. (2) Assessment of convergence and accuracy of the deflection history, stress distributions and damage zones calculated by the Finite Element Method (FEM) have not been given explicit attention. (3) Experimental investigations that document the response recorded during the impact test are rather limited. The majority of them provide a qualitative/quantitative description of the accumulated damage but not its growth. This information is vital to perform correlation studies. (4) Few investigators have considered full-scale components. The majority of the reported studies use generic structures such as beams, circular plates, square plates and cylindrical panels. There is a real danger in extrapolating conclusions drawn from studies on test specimens to real life components. It is more appropriate to use commercial FEM systems. They are widely distributed, well documented and user friendly. However, there is a need to verify the accuracy of material models, finite elements and analysis procedures in such systems before using them for the intended application, namely, prediction of the impact response of composite plates. The predictability issue itself demands a correlation study. This, in fact, is the aim and scope of the present study. A brief description of the specific methods and procedures used is given in the next section. In
Laminated composite plates are easily damaged by impacts, especially those normal to the plane of the laminate. Systematic study of such problems can be divided into three distinct areas: structural mechanics, damage mechanics and residual strength prediction. A study in that order will provide a mechanistic basis for both design and assessment of damage tolerance. Reliable and accurate prediction of the impact response of multilayered anisotropic plates covering a wide range of parameters is the focus of this study. Impact response means contact force history, deformation history, surface strain history, stress distribution across the laminate thickness (including interlaminar stresses), nature and extent of damage and stiffness and strength loss associated with that damage. Parameters significantly influencing the impact response include impact velocity/energy, material system, ply orientation and stacking sequence, plate geometry (shape, size and thickness), wall construction (solid laminate, integrally stiffened, and sandwich), support conditions, pre-stress state and initial curvature (curved panels). A review of previous work is available in a monograph L2 and two review papers, 3,4 which provide a background to the present study. A critical assessment of this vast literature, with particular reference to the focus of the present study, can be summarized as follows. (1) Finite element analysis of linear and transient structural behaviour has been the focus of intense research, 61
Composite Structures 0263-8223/94/S07.00 © 1994 Elsevier Science Limited, England. Printed in Great Britain
62
H. V. Lakshminarayana, R. Boukhili, R. Gauvin
succeeding sections, results of element evaluation, procedure verification and correlation study are presented and discussed. The presentation concludes by identifying directions for further work.
2 M E T H O D S AND P R O C E D U R E S
Numerical results for this study were generated using ABAQUS -- a general-purpose finite element code with emphasis on nonlinear applications. ~ This program is capable of modelling multilayered anisotropic materials. It provides elements suitable for dynamic analysis of composite plates and shells taking into account bendingmembrane coupling and transverse shear deformation effects. Among these, the S8R, an isoparametric quadrilateral plate/shell element, is employed in the present study. This element has eight nodes and six engineering degrees of freedom at each node. The user can specify within each element an arbitrary number of layers, each with its own thickness, ply orientation and orthotropic elastic properties. The formulation uses a value of 5/6 for the shear correction factors as default. However, the user has the option to use any other value through independent input of transverse shear stiffness. The element output includes membrane stress resultants, bending stress resultants and transverse shear-stress resultants, either at the nodes or at the four integration points. The user can also request ply-byply stresses at the integration points and at a maximum of 3 section points within each ply. It should be noted that the formulation of S8R does not ensure the continuity of interlaminar normal and shear stresses at the interfaces between plies. ABAQUS provides two procedures for calculating the response of structures subjected to impulsive loads. They are modal analysis and dynamic analysis by direct implicit integration of the equations of motion. The modal analysis is limited to linear transient behaviour. Its convergence and accuracy is dependent on the number of natural modes considered in the analysis. The pulse shape and the duration (t¢~) for which the contact force acts (in comparison with the period T~ of the fundamental mode of free vibration of the structure) strongly influence the response calculated by the modal method. For complex pulse shapes and for (t,,/Ti)'~ 1, a very large number of modes is required for convergence. Unfortunately, accurate determination of the eigenvalues and eigenfunctions associated with higher modes
places very heavy demands on computational resources by the FEM, Incidentally, eigenvalue extraction in ABAQUS is done using the subspace iteration method. The transient response calculation by the direct integration of the equations of motion is applicable to linear as well as nonlinear structural behaviour. When applying this step-by-step method, the time step At should be selected with caution, because a system of nonlinear algebraic equations must be solved at each time increment. This is done in ABAQUS iteratively by using Newton's method or, if preferred, the quasi-Newton's method. This time stepping, nonlinear equation-solving procedure is computationally expensive. The principal advantage of this procedure is that it is unconditionally stable, which means that there is no mathematical limit on the size of the time increment that can be used. In practice, At should be small enough to adequately define the history of excitation, its value being chosen on the basis of the shortest period which corresponds to the highest natural mode likely to contribute to the response. A recommended methodology(' for determining the contact force history is to measure the local contact stiffness of the composite plate in static tests and use this in conjunction with a finite element model. This approach, however, could not be considered truly predictive since it requires fabrication of the plate and static indentation tests for every impactor under consideration. Alternatively, it can be measured during impact tests using suitable instrumentation and data acquisition system. This approach is applicable to the nonlinear response also. The nature and extent of impact-induced damage is estimated by first calculating the dyamic stresses and their spatial distribution in the laminate and incorporating these in appropriate failure criteria. Multiple matrix cracks, delaminations and fiber breaks are the failure modes observed after impact tests. These failure modes and complex interactions between them complicate the prediction of damage. The tensor polynomial failure criterion proposed by Tsai and Wu 7 is employed in the present study to calculate the failure index (FI) given by F I = k ' l o , + F2o~ + F~o~ + F11o~ + F2:o! + F(~,o~ + _F~2o 102
where k~j are the strength tensors; 7 o~, o 2 are the lamina stresses in the fiber direction and transverse direction respectively; and o~, is the in-plane
Impact responseof laminated compositeplates shear stress. Loci of points at which FI = 1 define the damage zone. The failure modes are identified using the maximum stress criteria. 7 The stiffness and strength loss due to impactinduced damage are not at present truly predictable. A suggested approach involves flexure tests. ~ Evaluation of the accuracy of the S8R element for the analysis of composite plates in general is presented in the next section. We identify three distinct procedures: (1) modal analysis (procedure # 1); (2) linear and transient response analysis (TRA) (procedure #2); and (3) nonlinear and transient response analysis (NLTRA) (procedure # 3). The convergence and accuracy of each one of these is verified in a section entitled procedure verification. Finally, a critical assessment of their predictability with particular reference to the impact response of laminated composite plates is presented in Section 5.
3 ELEMENT EVALUATION Application of the full set of test problems proposed in Ref. 9 to the S8R element is not pre-
I I
6 _
W= 5111111
Ill Ill
W(O,O) E,h~ 10s Pa" I
t
4 PROCEDURE VERIFICATION
CPT Solution
1
-1
10 Fig. 1.
sented here. The results showed that a particular problem, a homogeneous, anisotropic, clamped, square plate under uniform pressure, provided an intensive measure of element performance. This problem, shown as the insert in Fig. 1, was chosen to evaluate the combined effect of material anisotropy and shear deformation on the accuracy of the S8R element. Numerical results were obtained for a unidirectional laminate, made of a highmodulus graphite/epoxy composite (EI/E 2= 40, E2 = 5.17 GPa (0.75 x 106 psi), Gi2--- 3.10 GPa (0"45 x 10 ~' psi), G23/GI2=0"8, vl2=0.25 ), for a ply orientation of 45 and for various values of the thickness ratio (h/a). Computed results are compared with the converged solutions given in Ref. 9. The agreement is very good for displacements as well as stress resultants provided (h/a)> 0.01. The effect of h/a on the predicted central deflection is shown in Fig. 1, which also has the classical plate theory (CPT) solution for comparison. For a given thickness ratio, the inaccuracy associated with the omission of transverse shear deformation effects is obtained from this figure. Obviously, the accuracy of the S8R element deteriorates for (hi a)<0"01, indicating that the formulation is not free from the so-called shear-locking phenomenon. There is scope to improve the performance of the SSR element for applications to thin and highly anisotropic situations. An approach to accomplish this is demonstrated in Ref. 10.
II
IIIIII
0
63
-2
10 h/a
-3
10
Effect of (h/a) on the central deflection predicted using the S8R element (8 x 8 mesh, whole plate).
The chosen elements and procedures should be capable of providing converged and accurate solutions to standard test problems for which alternative solutions are available for comparison. The specific problem used in this study is a crossply laminated, graphite/epoxy composite, square plate with all edges simply supported and subjected to impact at its center, which was investigated by Sun and Chen. ~ The geometry, ply orientation/stacking sequence, finite element mesh and contact force history (corresponding to impact by a 12.7 mm diameter steel ball at velocity V0 -- 3 m/s) used in the computations are shown in Fig. 2. Representative displacement time histories calculated using procedure # 2 (linear and transient response by direct implicit integration) are given in Fig. 3. These results are in excellent agreement with those reported in Ref. 11, in spite of the larger time step (6t=0.05 ms) used. No significant
64
H. V. Lakshminarayana, R. Boukhili, R. Gauvin
P(t)
h
w(t) P(t)
a
Pm a x
0°
h = 5.1 m m a = 203.2 mm
\
0° 90 °
Pr.~x= 27.4 Kg to = 0.2 ms
I I I
90'
/
I
0° 0 °
i
',\
t 1 = 0.25 t o
90 °
i
0 ° 90 °
t t
t
1
0°
0
Contact Force History
Ply Orientation/Stacking S e q u e n c e
/;
j/-
/ i
J
/
J
Finite Element Model (Quarter Plate)
Fig. 2.
Test problem for procedure verification.
improvement is noticed using a finer mesh and a smaller time step, indicating that a converged and accurate solution has been obtained. Figure 3 also indicates that two pulse shapes (triangular and half sine-wave) produce almost identical results. Since the half sine-wave pulse can be represented by an analytical expression, it is ideally suited to create benchmarks for the impact response of composites. For the same problem, Fig. 4 shows displacement-time histories calculated using procedure # 1 (modal analysis). There is no sign of convergence as the number of modes used is increased from 10 to 15. These results do not show any comparison with the reference solution in Ref. 11. Not only is the displacement history different, the maximum amplitudes differ by one order of
magnitude. This observation is substantiated by the findings of a recent round robin study.~2 It is noted that whenever the duration of impact is very small (in comparison with the period of the fundamental mode of free vibration of the target), the number of modes to be considered for convergence may be so large as to be computationally inefficient. It appears that reliable prediction of the impact response is outside the domain of the modal analysis. For the problem specified in Fig. 2, a comparison of the displacement histories calculated using procedure # 3 (nonlinear and transient response by direct implicit integration) and procedure # 2 is made in Fig. 5. The two solutions are in fact identical because at this low load level the resulting deflections are very small (in comparison
Impact response of laminated composite plates 0.24
0.22
D
0.22
0.2 0.18
0.2
0.16
0.18
0.14
0.16
~
0.12
¢"- 0.14
d
o.1
"*"= 0.12 0 0.1 v0
r
0.08 I~1 0.06 5 ¢G)
II
0.04
I~:
0.02
c-
E
0GI -0,02
0
Pmax= 27.4 Kg to= 0.2 ms t, = 0.25 t o
(3. .~_ -o.04
-o.os
O
006
0.08
0.04
E 0.02
o,
0
65
P(t
oo -0.02
=~
a
-0.04
\,
-0.08 -0.1 -0.12 0
11
I
0.5
1
0
; f/
i
-o ~ ~_o
' i (;t~
0
0.5
-0.06
t
-0,08
0
I
I, t¢ l
!1 Prnax= 27.4 Kg ' 1 1, to= 0"2,ms
-0.1
1.5
2
1
1.5
Time t(miliseconds)
Time t(miliseconds) Fig. 3.
\ fltl=
Central displacement history calculated using procedure # 2.
2
2
1.8
1.8
N u m b e r of m o d e s = 1 0
1.6
1.6
1.4
1.4
1.2
1.2 t-
1
1
d o.8
"~. o.e o <:5 o.6
0
v
v
0.4
II
0.2
15
I~
o,
•~ e-
!P(t)
c
E
0.6 0.4
o.2 0
Pr..x P(t)
E
8
o
e~
h
e~
w,t,
._~ a
P(t) = P S i n
i
°u
a
(nt/to)
Pmax= 27.4 Kg
0
to
0
t = 0.2 ms o
I
2
0
I
I
I
1
2
3
1
Time t(miliseconds) Fig. 4.
Time t(miliseconds)
Central displacement history calculated using procedure # 1.
66
H. V. Lakshminarayana, R. Boukhili, R. Gauvin 0.24
0.24
0.22
Procedure # 2 6t = 0.05 ms
0.22 0.2 0.18
0.18
7
0.16 e.-
0.16
0.14
~
0.14
C~ 0.12 0 0.1
d
o.12
II
0.1
0.08
0.08
15 ¢.. 006 E
0.04
0
0.02
Procedure # 3 6t = 0.05 ms
0.2
P(t)
/
to= 0.2 ms
_ g 0.04
Pmax
0.02
..........
~.
~. 0 .oo E3 -o.o2
o?,,J,,2t
-0.04 -0.06 -0.08
0
-0.1 0
f
I
I
0.5
1
Time Fig. 5.
"o
I
1.5
t(miliseconds)
o
-0.02 P(t) = Pmeeqin (rrt/t o)
-0.04 --
Pm~ = 27.4 Kg VO 3 m/s O= 0.2 ms
-0.06 I -0.08
-0.1 J 0
I
0.5
l
1
4
1.5
T i m e t(miliseconds)
Comparison of displacement-time histories calculated by procedures # 2 and # 3 (impact velocity = 3 m/s).
with the thickness h) and hence nonlinear effects do not show up. Impact velocity is a very important parameter controlling the response. Figures 6 and 7 provide a comparison of displacement time histories calculated using procedure # 2 and procedure # 3 at V~=10 m/s and V0=30 m/s, respectively. Obviously, the responses predicted by the two procedures do not agree with one another. Both the displacement-time history as well as the maximum amplitudes predicted differ rather significantly. Basically, the load levels are such that resulting deflections are of the order of the plate thickness and therefore procedure # 3, which includes nonlinear effects, is more appropriate. Unfortunately, the accuracy of the numerical solutions presented in Figs 6 and 7 could not be verified due to the nonavailability of reference solutions. Procedure # 3 is therefore more appropriate for numerical solutions of impact tests. In fact, for complex structures, this procedure may provide more cost-effective information than experimentation. However, there is a need to further validate the predictability of this procedure for the impact
response. This aspect is addressed in the next section. 5 CORRELATION STUDY The predictability aspect of the procedures verified in the previous section, with particular reference to the impact response of composite plates is assessed using a bench mark. Bench marks are fully specified and standard problems, which resemble instances found in practical applications, and for which reference solutions have been obtained using both analytical/numerical and experimental methods. The bench mark chosen for this study was created by Aggour and Sun. j3 It consists of a cross-ply laminated, E-glass/epoxy composite (Ej=38"6 GPa (5"6x106 psi), E2=I0.34 GPa (1.5×106 psi), G2=4"14 GPa (0.6x106 psi), v,2 =0.25) square plate with all edges clamped. The geometric parameters, ply orientation/stacking sequence, and the finite element discretization used in the computations are shown in Fig. 8. The contact force history corresponds to an impact by
Impact response of laminated composite plates
67
1 0.9
0.9 f
Procedure #2 6t = 0.05 ms
0.8
0.7
0.7
0.6
0.6
0.5
~, ~d
t-
0.5
0.4
0.4
0.3
C) 0.3 O 0.2
0.2
II 0.1
0.1 0
t
e(9
~
-0.1
~
-0.2
E -0.1 (9 -0.2
~.
-0.3
~ . -0.3
i'~
-(I.4
u) ~,~ -0.4
-0.5
-0.5
-0.6
-0.6
-0.7
-0.8
Procedure # 3 5t = 0.05 ms
0.8
0
0.5
1
~W(t) ~ Ih P(t)= PmaxSin(~t/t 0) V 0 = 10 m / s P m ~ = 9 1 . 2 Kg
-0.7
to
o
~-~
1.5
-0.8 0
2
t 0= 0.2
I
I
1
0.5
Time t(miliseconds) Fig. 6.
ms
I
1.5
2
Time t(miliseconds)
Comparison of displacement-time histories calculated by procedures # 2 and # 3 (impact velocity = 10 m/s).
4
2
Procedure #2 : LTRA
Procedure # 3 : NLTRA
1.8
6t = 0.05 ms
1.6 1.4 1.2
J::
¢-
1 0.8
c:) 0.6
o 0
0
II
O.4
II 0.2
,4..I
t5
0[
(9
E
P(t)
1
o,
c (9 -0.2
E
(9 -0.4 o ¢0 -0,6
o
P(t) = Pm=Sin (nl/t o ) Prr,= = 364.8 Kg
-2
-1.2
Vo= 30 m/s t -- 0.2 ms
-1.4
Vo= 30 m/s Pro= = 364.8 Kg to= 0.2 ms
0
-3 0
0.5
1
Time t(miliseconds) Fig. 7.
a•-~ -0.8 -1
1.5
-1.6
0.5
1
1.5
2
Time t(miliseconds)
Comparison of displacement-time histories calculated using procedures # 2 and # 3 (impact velocity = 30 m/s).
68
H.V. Lakshminarayana, R. Boukhili, R. Gauvin z [ P(t)
r
--.
~
X
a P(t)
P ,x ....
/
/
\ /
'\
'\ ',,
I
ho o
/iiiiiiiiiiiiiiii
9o!
I!iii!i!~iiii!iiii!iiiiiiliiiiiiiiiii!iiiiii!i ilili ~i!i!i!i!~!i!iiiii~ii~iililiiililili9o i!iiiiiJo
\\\
~_t
t °°
0 °
t Contact Force Histor~ P(t) = Pm~xSin(nUt o ) P~x = 310.1 Kg to= 0.25 ms
Ply Orientation/Stacking Sequence a = 139.7 mm h=4.1
J
mm I
A
B
E
H
Finite Element Model (Quater Plate).
Fig. 8. Testproblemfor correlation study.
a steel cylinder (diameter 9.5 mm and length 25.4 mm) at a velocity E, = 22.6 m/s. The calculated central displacement history using procedure # 3 is presented in Fig. 9. Results using procedure # 2 and test data taken from Ref. 13 are also included in the same figure to enable a three-way correlation. Results obtained using linear finite element analysis by Aggour and Sun ~3 are in close agreement with those obtained in the present study using procedure # 2. Furthermore, the responses calculated using procedure # 3 and procedure # 2 are identical, indicating that nonlinear effects are negligible. It is gratifying to note that computed results closely follow the central deflection measured during an impact test for t-<0"5t0, where to is the duration of impact. How-
ever, for higher values of t there is considerable difference between the computed and measured results. Therefore, to achieve a satisfactory correlation between the two sets of results over a larger duration, further refinement of the procedures used appears to be essential. Basically, the effects of impact-induced damage has to be included in the analysis. More specifically, a combined geometric and material nonlinear dynamic finite element analysis is indispensable for this purpose. However, this is more easily said than actually performed, because a priori knowledge of the nature and extent of impact-induced damage and accurate constitutive modelling of failed plies and damaged laminates will be required after each time step. At present,
Impact response of laminated compositeplates
69
1 Procedure # 2
0.9
6t = 0.025 ms
0.8
0.6
Testdata
0.4
~,
0.4
0.3
~
0.3
oc
o.1
g
~ -0.1
0{
~ -0.1
~ -0.2
~ -0.2
~
-0.3
~
~
-o.g
i~ -0.4
-0.5
-0.5
-0.6
-0.6
-0.7
-0.7
-0.8
-0.3
-0.8 0.5
l 0
Time t(miliseconds) Fig. 9.
Comparison of measured and predicted central displacement histories.
1.8
Procedure # 2
l'k 1.6
1.4
1.61.4-
VO= 40.0 m/s []
Vo-- 2226 m/s
¢...
0.6
Vo= 40.0 m/s
1
Vo= 30.0 m/s ©
0.8
Vo= 22•6 m/s
0,6 oo-
\
ji
0.4
/'//
0.2
II -0.2
i
.,.., e.- -0.2 1
(1) -0.4
E -o.4 ~
-0.6
\'~ "t~r~ \\
(:~ -0.6
-o8
~.
-0.8
-1
i~
-1
-1.2
-1.2
-1.4
~N
-1.6
-1.4 -1.6
-1.8 0.5
Time t(miliseconds) Fig. 10.
Procedure # 3
1.2
Vo= 30.0 m/s ©
0.8
0 ~
0.5
Time t(miliseconds)
1.8
!
]
0.2
I
0.1
c5 0.4
Studyl
Testdata
0.5
0.2
[
Prese~
0.7 0.6
0.5
c~
#3
6t = 0.05 ms
0.8
Present Study
0.7
Procedure
0.9
-1.8 0.5
Time t(miliseconds)
Comparison of predicted displacement-time histories by linear and nonlinear analysis procedures.
70
H. V. Lakshminarayana, R. Boukhili, R. Gauvin 1.6
Procedure #2 V = 40 m/s
1.4
0
t
= 0~25 ms
t=
Procedure #3 V = 40 m/s
14 ,2
o2oms
0
[
i
I t = 0z15 ms
x::
[t = o
, ~f i
ms/
/t = 0.10 msl
"~'I~ 0.8! ..
0.6 I
.6
0.4
0,4
0.2 O"I
-02
0.1
0.2
0.3
0.4
05
x/a Fig. 11.
-02
0
0.1
0.2
0.3
0.4
0.5
x/a
Comparisonof predicted displacementdistributionsby linear and nonlinearanalysisprocedures.
the former cannot be truly calculated while information on the latter is not available. The progressive failure finite element analysis outlined is in fact a topic for further research.~4 Reliable prediction of displacement history as a function of impact velocity (Fig. 10) constitutes only a fraction of the impact response story. Spatial distribution of displacement (Fig. 11 ), stress distribution across the laminate thickness (Fig. 12) and damage growth (Fig. 13) are in fact even more important. The F E M in general and A B A Q U S in particular is able to provide such results. Indeed a major problem is to find experimental techniques that can provide data that are equally detailed to validate the results presented here.
6 CONCLUDING REMARKS A hybrid experimental-numerical approach is necessary to predict the impact response of laminated composite plates covering a wide range of the parameters involved. Experimental determination of contact force history is essential if the structural behaviour is nonlinear. Nondestructive test methods are indispensable for characterisation of damage. Post-impact tests are needed to determine residual stiffness and strength. The combined geometric and material nonlinear finite element analysis capability required for this purpose is available in commercial finite element systems. However, accurate constitutive models
for composite laminates with multiple matrix cracks, fiber breaks and delaminations are not yet available. Delaminations are often the primary, life-limiting failure modes. The task of developing methods and data for predicting the onset and growth of delamination due to impact was not considered in the present study. Accurate evaluation of interlaminar normal and shear stresses a verified interface criteria and test methods to measure material properties associated with such a criterion are the prerequisites for prediction of delamination onset. It is to be noted that an appropriate finite element model for this purpose should be based on a lamination theory that enforces the continuity of interlaminar normal and shear stresses at ply interfaces,~5 and these are not yet implemented in commercial FEM systems such as ABAQUS. For the prediction of delamination growth, the fracture mechanics concept ~6is indispensable. This procedure involves calculation of energy release rates associated with delamination growth, the use of a mixed-mode fracture criterion and test methods to measure interlaminar fracture toughness. A drastically different finite element modelling approach is necessary for the numerical determination of energy release rates. These are topics deserving further research. This study sets the stage for confident application of a commercial FEM system ABAQUS for numerical simulation of impact tests conducted on composite plate and shell structures.
Impact response of laminated compositeplates Procedure #3
71
]
z/h
at t = 0 . 2 ms, V o = 4 0 m / s __
X = 1022.3 MPa
0 °
0'
•i•i•i:!•i•i:•i•i:•i•!i:•!i:•i:•i:•i•i:•i•i:!•i:•i:•i
90'
?::::..:.v....:.v..::.v.~..`..v:.......v.....v.v...v.`....:..:......`.:.v..::...
I1
~ : ~ : ~ : ~ : ~ : ~ : ~
90'
I
~......................................................... :~:~:~:~:~::~:~:~:~
90'
~::..`.:......`.:.v`.::~.:~.....:.v:.....v.....v`.`.`.`.:~.~.`.:::......::.v..`.`.`.. i..:.v..`.:....::.`....:::......`.:............:.v...v....:•..`.........v.....v:...v::.v :::~::~::~:~:~:~:~::~:~:~:~::~:~:~::~::~:~:~::~:~:~:~:~ 90' O'
--
O*
i
I
I
I
I
-1
1
o
01/X O° 0°
Y = 27.3 MPa
•!i•i•i•i:•i•i•i•i•i•ii:•i•i:i:•i•i•i•i•i:•i•ii:i:•i•i!•i!:•i•ii(.:.~...~.~.~.~.~::..~..~..~...~..~...~......`..~.~.::..`.9o' ~.`.`.~`~...~.`..~ •i•!i•!:!•i:!•:!•i•:i•!:i•!i•!h•i!•i:!•i•!i•!i:!•i•!:!•i 9O'
[
~.~.):.~.~.~.)):~)~.)))~.~.~`~.~.~.33)&:.:.:.~.))3~`~.~.)~.)~.)90' )~.:.:.:.~.~.333~:.~ "
0' 0° q
-10
lO
~/y z/h S = 40.9 MPa
•i•i:•i•i•i•i•i:•i•ii:•!•i!•ih:i:i:•i!•i•i•i:•i!:i:!•i!
9o'
:......:.....v:.v...v.......v~..`:..:.v...~`v.......v.v......:..::..:.v....`....`....: ~0' 0°
.1
o
otis Fig. 12.
Stress distribution across laminate thickness.
ACKNOWLEDGEMENTS t =0.05 ms ~]
t = 0.15 ms
~
t=O.2Oms
900
i 0°
This research was supported by the Natural Sciences and Engineering Research Council of Canada and la Direction de la Recherche de l'Ecole Polytechnique de Montr6al.
REFERENCES
V0= 40 m/s Procedure #3 Fig. 13.
Growth of impact damage.
1. Greszczuk, L. B., Damage in composite materials due to low-velocity impact. In h n p a c t D y n a m i c s , ed. J. A. Zukas et al. John Wiley & Sons, New York, USA, 1982. 2. Zukas, J. A., Numerical simulation of impact phenomena. In I m p a c t D y n a m i c s , ed. J. A. Zukas et al. John Wiley & Sons, New York, USA, 1982.
72
H. V. Lakshminarayana, R. Boukhili, R. Gauvin
3. Abrate, S., Impact on laminated composite material. Applied Mechanics Reviews, 44 ( 1991 ) 155-90. 4. Cantwell, W. J. & Morton, J., The impact resistance of composite materials -- Review. Int. J. Composites, 22 ( 1991 ) 347-62. 5. ABAQUS -A general purpose finite element code, users manual, version 4-7. Hibbit, Carlson and Sorensen Inc.. Providence, RI, USA. 6. Yan, T. M. & Sun, C. T., Use of statistical indentation laws in the impact analysis of laminated composite plates. J. Appl. Mech., 52 (1985)6-12. 7. S. W. Tsai, & Hahn, H. T., brtroduction to Composite Materials. Technomic, Lancaster, PA, USA, 1982. 8. Nemes, J. A. & Simmonds, K. E., Low-velocity impact response of foam-core sandwich composites. J. C~mp. Mat., 26 (1992) 500-19. 9. kakshminarayana. H. V. & Sridhara Murthy, S., A shear - - flexible triangular finite element model for laminated composite plates. Int. J. Nttmer. Meth. In Eng., 20 11984) 591-623. 10. Lakshminarayana, H. V. & Kailash, K., A shear -deformable curved shell element of quadrilateral shape.
Computers and Structures, 33 (1989) 987-1001. 11. Sun, C. T. & C h e n . J. K., On the impact of initially stressed composite laminates. J. (~mp. Materials. 19 (1985) 490-504. 12. Ewins, D. J. & lmregun, M., A survey to assess structural dynamic response prediction capabilities. In Quulit3, Assurance in FEM Technology, ed. J. Robinson. Robinson Associates, UK, 1987. 13. Aggour, H. & Sun, C. T., Finite element analysis of laminated composite plate subjected to circularly distributed central impact loading. Computers & Structures, 28 (1988) 729-36. 14. Yener, M. & Wolcott, E., Damage assessment analysis of composite pressure vessels subjected to random impact loading. J. Pressure Vessel Technolog3', 111 (1989) 124-9. 15. kec, C. Y. & kiu, D., lnterlaminar stress continuity theory for laminated composite analysis. AIAA ,/.. 29 1 9 9 ( ) 2010-12. 16. Friedrich. K. (ed.). Application of kracture Mechanics to ('omposite Materials. Elsevier Applied Science Publishers, Amsterdam, The Netherlands, 1989.