- Email: [email protected]
Numerical simulation of the crushing process of a corrugated composite plate
Composites: Part A 42 (2011) 1119–1126
Contents lists available at ScienceDirect
Composites: Part A journal homepage: www.elsevier.com/locate/compositesa
Numerical simulation of the crushing process of a corrugated composite plate Vladimir S. Sokolinsky ⇑, Kyle C. Indermuehle, Juan A. Hurtado Dassault Systèmes SIMULIA Corp., Providence, RI 02909, USA
a r t i c l e
i n f o
Article history: Received 2 December 2010 Received in revised form 6 April 2011 Accepted 27 April 2011 Available online 1 May 2011 Keywords: A. Fabrics/textiles C. Computational modelling C. Damage mechanics C. Finite element analysis (FEA) Crushing response
a b s t r a c t The significant potential of composite materials as highly effective energy absorbers has contributed to their increasing usage in the aerospace, automotive, and railway industries. Traditionally, the development of crashworthiness composites design relies upon laborious and costly experimental testing. Therefore, the ability to simulate accurately the crushing response of composites and their energy absorption mechanisms can significantly reduce the product development cycle and cost. The present work presents a physics-based finite element model of a corrugated carbon–epoxy fabric composite plate subject to quasi-static crushing. The model accounts for both intralaminar (in-plane) and interlaminar (delamination) failure mechanisms. The in-plane response of the fabric plies is modeled using a proposed constitutive model that has been implemented in the commercial finite element code Abaqus/Explicit; and the delamination response of the composite plate is simulated using the cohesive surface capability in Abaqus. The results of the Abaqus/Explicit simulations show very good quantitative and qualitative agreement with the experimental data, thus demonstrating that the proposed methodology and tools can be used for realistic crush simulations of composite structures. Ó 2011 Elsevier Ltd. All rights reserved.
1. Introduction Composites exhibit a significantly higher energy absorption capacity per unit weight than metallic structures. Therefore, the use of composite structures as energy absorbers is rapidly growing in many industries. Energy-absorbing composite components, which are usually designed to dissipate kinetic energy under controlled collapse or disintegration, substantially increase structural crashworthiness in exchange for a reasonable increase in overall vehicle cost. Traditionally, the development of crashworthiness design relies upon experimental testing. When several loading directions are involved, the experimental testing of the crushing process in composites is rather costly. In addition, it is difficult to carry out a meaningful experimental program during the initial design stages. Thus, the ability to simulate accurately the crushing response of composites and their energy absorption mechanisms can significantly reduce the product development cycle and cost. However, the physics of composites crushing is very complicated because composites absorb crash energy through a combination of interlaminar (delamination) and several intralaminar (in-plane) fracture mechanisms. The experimental characterization of specific energy absorption (SEA) of composite materials is frequently performed using thinwalled tubular specimens of various geometries. Less frequently ⇑ Corresponding author. Tel.: +1 401 276 8480; fax: +1 401 276 4408. E-mail address: [email protected] (V.S. Sokolinsky). 1359-835X/$ - see front matter Ó 2011 Elsevier Ltd. All rights reserved. doi:10.1016/j.compositesa.2011.04.017
used methods employ: (a) simple flat plate specimens in conjunction with very involved anti-buckling supporting fixtures or, alternatively, (b) self-supporting specimens of more complicated geometries that do not require the use of dedicated test fixtures. Self-supporting specimens in the form of corrugated web structures show robustness to both geometric shape parameters and the type of crush initiator while requiring a minimal amount of supporting fixtures and offering lower manufacturing costs than other test methods used for SEA characterization [1]. The work presented here examines the numerical simulation of the crushing response of a carbon–epoxy fabric corrugated plate specimen using the Abaqus/Explicit finite element code. The geometry and lay-up of the corrugated plate are presented first. Then a proposed constitutive model for fabric-reinforced composites is described, followed by a description of a physics-based finite element model that accounts for both delamination and in-plane failure of fabric-reinforced composite plies. Finally, the numerical results are compared with the experimental data and conclusions are presented. 2. Description of the test specimen The corrugated plate discussed here was developed and tested by the Crashworthiness Working Group of the CMH-17 organization within the framework of a recent round-robin effort during which a series of crush tests on a corrugated plate specimen were performed. The test results were made available to the authors who participated in the round-robin effort [2].
1120
V.S. Sokolinsky et al. / Composites: Part A 42 (2011) 1119–1126
2.1. Geometry and composite lay-up The cross-section of the corrugated plate specimen can be obtained by repeating a semicircular segment of radius 6.35 mm (0.25 in.) three times at alternating sides with respect to the midplane (see Fig. 1). The specimen has straight edges of material on each side of the corrugation (end-lips) for additional stability. The test specimens were manufactured of a carbon–epoxy TORAYCA fabric (see Table 1) by press-molding through a set of aluminum matching tools [1]. The [0/90]2S composite lay-up was 2 mm (0.079 in.) thick. To facilitate stable crushing, a single-sided 45-degree chamfer was machined at the bottom of the specimen end after cure. Each specimen was 76.2 mm long (3 in.) and 50.80 mm wide (2 in.).
2.2. Experimental set-up The specimens were tested in the vertical configuration standing free on an as-milled dry steel surface. The crushing plate of the experimental rig was free to slide along four posts. A self-aligning sphere was used to introduce the crushing load from the test frame onto the specimen. Because the corrugated specimens are self-supporting, no supports were used on the specimens. The testing was performed at a quasi-static speed of 0.2 mm/s (0.5 in./ min). Fig. 2 shows typical views of the crushed corrugated fabric plate.
3. Constitutive model for fabric-reinforced composites This section describes a proposed in-plane constitutive model for a fabric composite ply that has been implemented in Abaqus/ Explicit. The fabric-reinforced composite ply is modeled as a homogeneous orthotropic material that is capable of sustaining progressive stiffness degradation due to fiber/matrix cracking and plastic deformation under shear loading. Fig. 2. Crushed corrugated plate specimen: (a) side view; (b) top view (image courtesy of Prof. Paolo Feraboli, the University of Washington, Seattle). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
3.1. Elastic continuum damage mechanics model
Fig. 1. Schematic of a corrugated fabric panel (all dimensions are in millimeters). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
The fiber directions of woven fabric reinforcement are assumed to be orthogonal. The constitutive stress–strain relations are formulated in a local Cartesian coordinate system with the base vectors aligned with the fiber directions. The in-plane elastic stress–strain relations are given according to orthotropic damaged elasticity [3]:
8 9 > < e11 > =
Variable
1 ð1d1 ÞE1
6 m21 e22 ¼ 6 4 E2 > > : el ; e12 0
Table 1 Material properties of the TORAYCA T700/2510 fabric [2]. Description
2
mE121
0
1 ð1d2 ÞE2
0
0
1 ð1d12 Þ2G12
Value T
Longitudinal modulus (GPa) Transverse moduli (GPa) Principal Poisson’s ratio Shear moduli (GPa) Longitudinal tensile strength (MPa) Longitudinal compressive strength (MPa) Transverse tensile strength (MPa) Transverse compressive strength (MPa) In-plane shear strength (MPa) Mode-I fracture toughness (J/m2) Mode-II fracture toughness (J/m2)
E11 E22 = E33
m12 G12 = G23 = G31 X1+ X1 Y2+ Y2 S GIC GIIC
55.9 54.4 0.042 4.2 911.3 704.0 770.1 698.2 131.6 504 1566
38 9 > r11 > = 7< 7 r22 5> > : ;
ð1Þ
r12
where e ¼ f e11 ; e22 ; eel12 g is the elastic strain vector; {r = r11, r22, r12}T is the stress vector; E1 and E2 are the Young’s moduli in the principal orthotropic directions; G12 is the in-plane shear modulus; m12 and m21 are the principal Poisson’s ratios; d1 and d2 are the damage parameters associated with fiber fracture along the principal orthotropic directions; and d12 is the damage parameter associated with the matrix micro-cracking due to the in-plane shear deformation. The damage parameters vary between 0 and 1 and represent stiffness degradation caused by microdamage in the material.
1121
V.S. Sokolinsky et al. / Composites: Part A 42 (2011) 1119–1126
Note that there is no damage parameter associated with the degradation in Poisson’s ratios. As a result, m12 and m21 are reduced by factors (1 d1) and (1 d2), respectively, with progression of damage. To differentiate between the tensile and compressive fiber failure modes, the corresponding damage variables are calculated based on the stress state in the fiber directions as follows (i = 1, 2):
di ¼ diþ
hrii i hrii i þ di jrii j jrii j
ð2Þ
where di+ and di represent the tensile and compressive damage variables of the fibers in the ith direction, respectively. The damage evolution equations are assumed in the following form:
~ 1þ Þ; d1 ¼ d1 ðr ~ 1 Þ d1þ ¼ d1þ ðr ~ 2þ Þ d2 ¼ d2 ðr ~ 2 Þ d2þ ¼ d2þ ðr ~ 12 Þ d12 ¼ d12 ðr
ð3Þ
/12 ¼
~ iþ and r ~ i (i = 1, 2) are defined as: where the effective stresses r
r~ iþ ¼ r~ 12 ¼
hrii i 1 diþ
r~ i ¼
r12
hrii i 1 di
ð4Þ
pffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffi r~ iþ ¼ 2Ei Y iþ r~ i ¼ 2Ei Y i pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi r~ 12 ¼ 2G12 Y 12
ð5Þ
where the thermodynamic forces Yi+, Yi, and Y12 are responsible for damage progression. Thus, it can be seen that fiber and shear damage modes are uncoupled. At any given time, the elastic domain is defined in terms of the damage activation functions Fi+, Fi, and F12 as:
F iþ ¼ /iþ r iþ 6 0 F i ¼ /i r i 6 0 F 12 ¼ /12 r 12 6 0
ð6Þ
here ri+, ri, and r12 are the tensile, compressive, and shear damage thresholds, respectively; and the functions ui+ and ui provide a criterion for fiber failure and are given by (i = 1, 2):
r~ iþ X iþ
r~ 12
ð10Þ
S
where S is the shear strength; and, after shear damage initiation (u12 = 1), the damage threshold, r12, is given by
r12 ðtÞ ¼ max /12 ðsÞ
ð11Þ
s6t
1 d12
The effective stresses are directly related to the thermodynamic forces, which are conjugate to the damage variables, through the following expression (i = 1, 2):
/iþ ¼
where Lc is the characteristic length of the element; Gaf is the fracture energy per unit area under uniaxial tensile or compressive loading; and g a0 is the elastic energy per unit volume (elastic energy density) at the point of damage initiation. The formulation of the damage evolution law ensures that: (a) the damage variables are monotonically increasing quantities (i.e., d_ a P 0) and (b) the correct amount of energy is dissipated when the lamina is subjected to uniaxial loading along the fiber directions irrespective of the element size. The latter is achieved by formulating the damage evolution law in Eq. (9) using the characteristic element length, Lc. Thus, the use of Lc is critical to ensure objectivity of the numerical results with respect to mesh refinement. The formulation imposes an upper bound, Lmax ¼ Gaf =g a0 , on the characteristic element length, Lc, of the finite element mesh. Similarly, the function u12, which provides a criterion for the initiation of shear damage in the matrix, is assumed to be
/i ¼
r~ i X i
ð7Þ
where Xi+ and Xi are the tensile and compressive strengths, respectively, for uniaxial loading along the fiber directions. The damage thresholds ri+ and ri are initially set to unity. After damage initiation, i.e., when ui+ = 1 or/and ui = 1, they increase with damage progression according to:
riþ ðtÞ ¼ max /iþ ðsÞ s6t
ri ðtÞ ¼ max /i ðsÞ
ð8Þ
s6t
The definition ensures that the damage thresholds are nondecreasing quantities, i.e., that r_ a ðtÞ P 0 (a = ±1, ±2). The damage thresholds are assumed to obey the Kuhn–Tucker complementary conditions and the consistency condition ðr_ a F_ a ¼ 0Þ. Following the work in Ref. [4], the evolution of the damage variables is given by the following expression (a = ±1, ±2):
" # 1 2g a0 Lc da ¼ 1 exp a ðra 1Þ ra Gf g a0 Lc
ð9Þ
Furthermore, in accordance with Ref. [3], it is assumed that the shear damage variable increases with the logarithm of r12 until a maximum value of d12 is reached, i.e., max
d12 ¼ minða12 lnðr 12 Þ; d12 Þ
ð12Þ
max d12
where a12 > 0 and 6 1 are material properties. Please note that max the values of S, a12, and, d12 can be determined using a calibration procedure presented in the Appendix. 3.2. Elastic–plastic shear model The in-plane shear response is dominated by the non-linear behavior of the matrix that exhibits both stiffness degradation due to matrix microcracking and plasticity. The elastic part of the shear response was discussed previously. Here, the plastic shear response of the material is considered. The matrix response may be inelastic due to extensive cracking or plasticity. This leads to permanent deformations in the ply upon unloading. To account for these effects, a classical plasticity model with an elastic domain function and a hardening law, which is applied to the effective stresses in the damaged material, is used. The elastic domain function, F, is given by:
~ 12 j r ~ 0 ðepl Þ 6 0 F ¼ jr
ð13Þ
and the hardening law is taken to be of the form:
r~ 0 ðepl Þ ¼ r~ y0 þ Cðepl Þp
ð14Þ
~ y0 is the initial effective shear yield stress; C and p are coefwhere r ficients; and epl is the equivalent plastic strain due to shear defor~ y0 , C, and p can be determined using the mation. The values of r calibration procedure described in the Appendix. 4. Description of the finite element model In this section, a finite element model of the corrugated specimen described in Section 2 is detailed. Both intralaminar (in-plane) and interlaminar (delamination) failure mechanisms are considered in order to simulate the real physics of composites crushing. The finite-element model of the corrugated composite fabric plate was created in Abaqus/CAE [5]. The crush initiator in the form of a single-sided 45-degree chamfer was accounted for in the model (see Fig. 3). Each of the eight fabric plies was individually meshed
1122
V.S. Sokolinsky et al. / Composites: Part A 42 (2011) 1119–1126
with continuum shell elements (SC8R). A uniform mesh with an element size of 1 mm by 0.5 mm (0.04 in. by 0.02 in.) was used for each ply. Continuum shell elements in Abaqus have the geometry of a three-dimensional solid element but their kinematic and constitutive behaviors are similar to those of conventional shell elements. This type of element allows for the accurate modeling of stacked composite plies. The finite element model used for this simulation consisted of approximately one-half million degrees of freedom.
4.1. Intralaminar response of the plate The in-plane response of the fabric plies was modeled using the Abaqus constitutive model described previously. This model requires the input of the: (a) elasticity constants; (b) damage initiation coefficients; (c) damage evolution coefficients; and (d) shear plasticity coefficients. The elasticity constants include the Young’s moduli in the fiber directions, the principal Poisson’s ratio, and the in-plane shear modulus. The damage initiation coefficients include the tensile and compressive strengths along the fiber directions and the shear strength at the onset of shear damage. The values of the user material constants and damage initiation coefficients are given in Table 1. Please note that the strength and toughness properties of the TORAYCA fabric were adopted directly from Table 1 without any further adjustments.
The damage evolution coefficients are represented by the tensile and compressive fracture energies per unit area along the fiber max directions and by the parameters a12 and d12 in Eq. (12). The fracture energies can be determined through the testing procedure described in Ref. [6]. However, for the present work this testing was not performed, and average values of the fracture energies G1þ and f G1 for similar material systems from Ref. [6] were utilized instead f (see Table 2). The critical values of the fracture energies along the 1þ fiber direction 2 were calculated so that the ratios G2þ and f =Gf 1 G2 =G were approximately equal to the ratios Y /X and Y2/ 2+ 1+ f f X1, respectively. max The parameters a12 and d12 should be determined using the calibration procedure as described in the Appendix. For the present work, however, the cyclic tensile tests were not conducted so the max values of a12 = 0.15 and d12 ¼ 1 were used in the analysis based on the calibration results for similar material systems. The shear plastic coefficients include the initial effective shear ~ y0 , and the coefficients C and p (see Eq. (14)). These yield stress, r values are also determined as shown in the Appendix. In the absence of the cyclic tensile tests, these coefficients were also adopted based on the calibration results for similar material systems (see Table 3). Please note that shearing of the fabric ply is not a dominant failure mode for this analysis. Therefore, the use of typical values for the damage and plasticity shear material parameters do not have a strong influence on the results presented here. This statement was numerically verified by perturbing the discussed shear parameters over a significantly wide range of values without changes in the results. 4.2. Interlaminar response of the plate In order to account for the adhesive bonds between adjacent plies and model the delamination phenomenon, a generalized traction-separation law was used. Damage initiation at the interfaces was defined using the quadratic nominal stress criterion according to which damage is assumed to initiate when the quadratic interaction function involving the stress ratios reaches a value of one [5]:
( )2 ( )2 ( )2 ht n i ts tt þ 0 þ 0 ¼1 t 0n ts tt
ð15Þ
where t = {tn, ts, tt}T is the nominal traction vector and the interface damage initiation properties used in the present analysis are listed in Table 4. These properties are usually difficult to determine experimentally with good accuracy; therefore, they can be used as calibration parameters, if required. In the present work, the variation of the damage initiation parameters within 15% from the typical values given in Table 4 did not influence the analysis results. A 15% variation was used because it represents the typical maximum deviation of the material properties in manufactured composite parts. Progression of damage at the interfaces was modeled using the critical mixed mode energy behavior based on the Benzeggagh–Kenane (BK) law [5]:
Table 2 Fracture energies per unit area along the fiber directions.
Fig. 3. Abaqus/CAE model of the corrugated composite fabric plate: (a) top view; (b) side view (the finite element mesh is not shown). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
Description
Variable
Value (kJ/m2)
Tensile fracture along fiber direction 1
G1þ f G1 f G2þ f G2 f
125
Compressive fracture along fiber direction 1 Tensile fracture along fiber direction 2 Compressive fracture along fiber direction 2
250 95 245
1123
V.S. Sokolinsky et al. / Composites: Part A 42 (2011) 1119–1126 Table 3 Shear plasticity coefficients. Description
Variable
Value
Initial effective shear yield stress (MPa) Coefficient in the hardening equation (Eq. (14)) Power term in the hardening equation (Eq. (14))
r~ y0
195 1053 0.41
G g S GCn þ GCs GCn ¼ GC GT
C p
ð16Þ
where GS = Gs + Gt; GT = Gn + GS; Gn, Gs, and Gt refer to the work done by the normal and the first and the second shear forces acting in the interface, respectively; GCn and GCs are critical fracture energies required to cause failure in the normal and shear directions, respectively ðGCt ¼ GCs Þ; GC is the total critical mixed-mode fracture energy; and g is the power coefficient taken to be 2.284. The values of the required normal ðGCn Þ and shear ðGCs Þ mode critical fracture energies were taken from the last two rows of Table 1. After damage initiation, the evolution of the damage variable until final failure was modeled using an exponential softening law [5].
(S4R). The crushing load was applied to the composite plate through another rigid steel plate moving with a constant velocity of 200 mm/s (470 in./min). This loading rate is approximately eight times higher than the maximum quasi-static test speed reported in Ref. [1] and was used to achieve a reasonable run time for the explicit dynamics analysis. This artificial increase of the loading velocity has no effect on the analysis results because according to previous work (see Ref. [1]), dynamic effects become significant beyond 500–1000 mm/s (1200–2400 in./min). As pointed out in Ref. [7], the friction between the surfaces of the fronds and the crushing surface and between the fronds and the debris wedge plays an important role in the analysis. Here, the friction coefficient between the surfaces of the debonded composite plies was assumed to be 0.3 [7]. The friction coefficient between the composite plate and the rig was taken to be 0.12 based on the data available in the literature. The interactions between the debonded composite plies and between the composite plies and the rig were modeled using the general contact capability of Abaqus/Explicit. After delamination between any two adjacent plies occurs, the general contact definition is automatically updated to account for possible post-delamination contact.
4.3. Element deletion
4.5. Mass scaling
Failed continuum shell elements were removed from the model to prevent their excessive distortion and premature termination of the analysis. In the present work, damage-based element deletion is combined with a deformation-based deletion criterion. The damage-based element deletion is activated when any damage variable along the fiber directions or the equivalent plastic strain due to shear deformation reaches a maximum specified value. The deformation-based deletion criterion is used when either the tensile or compressive principal logarithmic strain reaches its maximum or minimum specified value, respectively (see Table 5). Additionally, the detached composite fragments are deleted when they move far enough away from the crush zone. Conceptually, the zone is defined by a three-dimensional box that encompasses the crush area. With the aid of Abaqus/Explicit user subroutine VUSDFLD, the detached parts are detected and deleted when they move out of the range of this box. This is done to prevent non-physical numerical distortions of the detached fragments that may cause premature termination of the analysis.
The mass of the composite plate was increased by three orders of magnitude at the beginning of the analysis using the fixed mass scaling capability of Abaqus/Explicit. Similar to the artificial increase of the loading velocity discussed previously, this numerical device was employed to achieve a reasonable run time. A comparison between the kinetic and strain energies confirmed that the quasi-static character of the problem was not changed as a result.
4.4. Loads and boundary conditions The composite plate was positioned vertically on a rigid steel plate which was modeled using conventional shell elements
5. Simulation results and discussion The simulation required 9 h of wall time to run on a 16-core Linux cluster. The deformed shape of the corrugated fabric plate at the end of the crushing simulation is shown in Fig. 4. Note that the deletion of damaged elements in the plies creates detached parts that are removed from the model as explained previously. A comparison between Figs. 2 and 4 demonstrates very good qualitative correlation between the experimental and numerical results. Specifically, it can be seen from this comparison that the Abaqus/Explicit simulation reproduced important characteristic features of the crushing response of composites such as the frond formation in the composite plate, the propagation of the main delamination through the plate, and the accumulation of debris between the loading rig and the composite plate.
Table 4 Interface damage initiation properties. Description
Variable
Value (MPa)
Maximum nominal stress in the normal-only mode
t 0n
54
Maximum nominal stress in the first shear direction (for a mode that involves separation only in this direction)
t 0S
70
Maximum nominal stress in the second shear direction (for a mode that involves separation only in this direction)
t 0t
70
Table 5 Element deletion controls. Description
Variable
Maximum value of the damage variable in the fiber directions used in the damage-based deletion criterion Maximum value of the equivalent plastic strain used in the damage-based deletion criterion
dmax
Maximum (positive) principal logarithmic strain used in the deformation-based deletion criterion Minimum (negative) principal logarithmic strain used in the deformation-based deletion criterion
epl max ^emax ^emin
Value 0.9999 0.5 0.1 0.1
1124
V.S. Sokolinsky et al. / Composites: Part A 42 (2011) 1119–1126
Fig. 4. Simulated deformed shapes of the corrugated composite plate: (a) front view (XOZ-plane); (b and c) side views (YOZ-plane); and (d) isometric view. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
Moreover, a comparison of Fig. 2b with Fig. 4b–d reveals that the simulation was capable of correctly predicting an important qualitative difference in the delamination propagation through the end-lips of the composite plate. Namely, one of the end-lips was split in half by the propagation of the main delamination through its height all the way from top to bottom. In the other end-lip, the main delamination traveled through approximately one-third of its height. It should be emphasized that the use of a uniform mesh as described in Section 4 plays an important role in the accurate simulation of delamination propagation through the end-lips of the composite plate. A close match between the experimental (Ref. [2]) and simulated load displacement curves is demonstrated in Fig. 5. Both the peak and average crush forces were predicted accurately. It can be seen from the figure that the height of the crushed plate was reduced by half with respect to its undeformed configuration. Please note that both the experimental results and the simulation
results were filtered using a SAE filter. However, the frequency cut-off used for the experimental results was 600 Hz, whereas for the simulation results, 100 Hz frequency cut-off was used. The lower cut-off frequency was used to remove the numerical noise characteristic of an explicit analysis crushing simulation. The numerical results for the energy absorbed (EA), which is calculated as the total area under the load–displacement curve in Fig. 5, are shown versus stroke in Fig. 6. The term stroke here is used in reference to the length of the plate that has been destroyed during crushing [1]. As can be seen from the figure, the experimental and simulated curves nearly coincide with each other. Finally, a comparison between the experimental and simulated specific energy absorption (SEA) of the corrugated plate is shown in Fig. 7. SEA shows the ability of a material to dissipate energy; it is measured in Joules per gram and usually varies for composites between 15 and 80 J/g [1]. SEA is computed by dividing EA by the crushed mass of the plate, which can be expressed as the product
Fig. 5. A comparison between the experimental and numerically-predicted load–displacement curves. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
V.S. Sokolinsky et al. / Composites: Part A 42 (2011) 1119–1126
1125
Fig. 6. A comparison between the experimental and numerical results for the absorbed energy versus stroke. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
of the stroke, cross-sectional area of the plate, and its density. As can be seen from Fig. 7, the numerical and experimental results agree well, especially in the second half of the crushing process. The largest discrepancy between the results is observed at the early loading stages when the chamfer part of the plate is crushed. Taking into consideration the extremely complicated nature of the crushing process, the overall correspondence between the measured and simulated SEA curves is judged to be acceptable. Some comments on the robustness of the presented methodology are in order. As has been mentioned previously, the strength and interface toughness properties of the TORAYCA fabric were adopted directly from Ref. [2] without any further adjustments, and average values for similar material systems from Ref. [6] were used for the critical fracture energies per unit area along the fiber directions. Furthermore, the interface damage initiation properties were based on the data available in the literature, and the model response did not change as a result of a realistic deviation from these properties. Moreover, the perturbations of the shear damage and plasticity parameters over a wide range of values did not have any appreciable impact on the results. Therefore, it can be
concluded that the presented methodology exhibits robustness with respect to the input material properties. 6. Conclusions A finite element model of a corrugated fabric composite plate subject to quasi-static crushing has been presented. A physicsbased approach that accounts for both intralaminar (in-plane) and interlaminar (delamination) failure mechanisms has been used to model the crushing behavior of the composite plate. The in-plane response of the fabric plies has been modeled with the aid of the proposed constitutive model implemented in Abaqus/Explicit, whereas delamination has been simulated using the Abaqus cohesive surface capability. The results of an Abaqus/Explicit simulation show very good quantitative and qualitative agreement with experimental data and exhibit robustness with respect to the input material properties. Thus, the proposed methodology shows considerable promise for realistic crush simulations of composite structures. This may permit substitution of costly experimental testing with numerical simulation for the
Fig. 7. A comparison between the experimental and numerical results for the specific energy absorption. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
1126
V.S. Sokolinsky et al. / Composites: Part A 42 (2011) 1119–1126
Fig. A2. Schematic of the calibration of the: (a) shear damage parameters; (b) shear hardening curve. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) Fig. A1. Schematic of typical response of a fabric-reinforced composite. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
crashworthiness design of composites. Further research is needed to investigate the limits of the proposed numerical methodology by applying it to a wide variety of composites crush simulations. This will be the focus of future work. Acknowledgements The first author would like to express his sincere gratitude to Qingping Huang and Fernando Carranza of SIMULIA whose technical advice enhanced the quality of the present work. Appendix. Calibration procedure for the shear response parameters used in the constitutive model for fabric-reinforced composites The calibration of the shear response parameters uses the results of cyclic tensile tests for ±45° laminates. An idealized representation of the outcome of such tests (the actual response exhibits hysteretic behavior upon unloading and reloading) appears in Fig. A1. First, the level of damage is calculated from the ratio of the unloading stiffness to the initial (undamaged) elastic stiffness (see Fig. A1). This generates pairs of stress–damage values, (r12, d12), for each unloading branch of the cyclic shear stress– ~ 12 ; d12 Þ as strain curve. These are represented as coordinates ðln r
shown in Fig. A2a. Then the slope of a linear fit curve for these pairs provides the value of a12, and its intersection with the abscissa, the value of S. If the damage data indicate saturation (see Fig. A2a), the max corresponding value of d12 (d12 6 1) should be used as d12 ; othermax wise, d12 ¼ 1 is adopted. Furthermore, for each unloading branch of the cyclic shear stress–strain curve (Fig. A1), the plastic strain epl 12 at the onset of unloading is determined from the value of residual deformation in the unloaded state. Then the coefficients of a nonlinear fitting ~ 12 ; epl curve for the points ðr 12 Þ at the onset of unloading are used as the plastic coefficients of the fabric material as illustrated in Fig. A2b. References [1] Feraboli P. Development of a corrugated test specimen for composite materials energy absorption. J Compos Mater 2008;42(3). [2] CMH-17 meeting results. 6th CMH-17 Crashworthiness Working Group Meeting, Ottawa, CA; 5 August 2008. [3] Johnson AF. Modelling fabric reinforced composites under impact loads. Composites Part A 2001;32:1–2. [4] Maimí P, Camanho PP, Mayugo J-A, Dávila CG. A thermodynamically consistent damage model for advanced composites. Technical Report NASA/TM-2006214282, NASA; 2006. Ò [5] Abaqus 6.10 Analysis user’s manual. Dassault Systèmes; 2010. [6] Pinho ST, Robinson P, Iannucci L. Fracture toughness of tensile and compressive fibre failure modes in laminated composites. Compos Sci Technol 2006;66(13). [7] Pinho ST, Camanho PP, De Moura MF. Numerical simulation of the crushing process of composite materials. Int J Crashworthiness 2004;9(3).
Contents lists available at ScienceDirect
Composites: Part A journal homepage: www.elsevier.com/locate/compositesa
Numerical simulation of the crushing process of a corrugated composite plate Vladimir S. Sokolinsky ⇑, Kyle C. Indermuehle, Juan A. Hurtado Dassault Systèmes SIMULIA Corp., Providence, RI 02909, USA
a r t i c l e
i n f o
Article history: Received 2 December 2010 Received in revised form 6 April 2011 Accepted 27 April 2011 Available online 1 May 2011 Keywords: A. Fabrics/textiles C. Computational modelling C. Damage mechanics C. Finite element analysis (FEA) Crushing response
a b s t r a c t The significant potential of composite materials as highly effective energy absorbers has contributed to their increasing usage in the aerospace, automotive, and railway industries. Traditionally, the development of crashworthiness composites design relies upon laborious and costly experimental testing. Therefore, the ability to simulate accurately the crushing response of composites and their energy absorption mechanisms can significantly reduce the product development cycle and cost. The present work presents a physics-based finite element model of a corrugated carbon–epoxy fabric composite plate subject to quasi-static crushing. The model accounts for both intralaminar (in-plane) and interlaminar (delamination) failure mechanisms. The in-plane response of the fabric plies is modeled using a proposed constitutive model that has been implemented in the commercial finite element code Abaqus/Explicit; and the delamination response of the composite plate is simulated using the cohesive surface capability in Abaqus. The results of the Abaqus/Explicit simulations show very good quantitative and qualitative agreement with the experimental data, thus demonstrating that the proposed methodology and tools can be used for realistic crush simulations of composite structures. Ó 2011 Elsevier Ltd. All rights reserved.
1. Introduction Composites exhibit a significantly higher energy absorption capacity per unit weight than metallic structures. Therefore, the use of composite structures as energy absorbers is rapidly growing in many industries. Energy-absorbing composite components, which are usually designed to dissipate kinetic energy under controlled collapse or disintegration, substantially increase structural crashworthiness in exchange for a reasonable increase in overall vehicle cost. Traditionally, the development of crashworthiness design relies upon experimental testing. When several loading directions are involved, the experimental testing of the crushing process in composites is rather costly. In addition, it is difficult to carry out a meaningful experimental program during the initial design stages. Thus, the ability to simulate accurately the crushing response of composites and their energy absorption mechanisms can significantly reduce the product development cycle and cost. However, the physics of composites crushing is very complicated because composites absorb crash energy through a combination of interlaminar (delamination) and several intralaminar (in-plane) fracture mechanisms. The experimental characterization of specific energy absorption (SEA) of composite materials is frequently performed using thinwalled tubular specimens of various geometries. Less frequently ⇑ Corresponding author. Tel.: +1 401 276 8480; fax: +1 401 276 4408. E-mail address: [email protected] (V.S. Sokolinsky). 1359-835X/$ - see front matter Ó 2011 Elsevier Ltd. All rights reserved. doi:10.1016/j.compositesa.2011.04.017
used methods employ: (a) simple flat plate specimens in conjunction with very involved anti-buckling supporting fixtures or, alternatively, (b) self-supporting specimens of more complicated geometries that do not require the use of dedicated test fixtures. Self-supporting specimens in the form of corrugated web structures show robustness to both geometric shape parameters and the type of crush initiator while requiring a minimal amount of supporting fixtures and offering lower manufacturing costs than other test methods used for SEA characterization [1]. The work presented here examines the numerical simulation of the crushing response of a carbon–epoxy fabric corrugated plate specimen using the Abaqus/Explicit finite element code. The geometry and lay-up of the corrugated plate are presented first. Then a proposed constitutive model for fabric-reinforced composites is described, followed by a description of a physics-based finite element model that accounts for both delamination and in-plane failure of fabric-reinforced composite plies. Finally, the numerical results are compared with the experimental data and conclusions are presented. 2. Description of the test specimen The corrugated plate discussed here was developed and tested by the Crashworthiness Working Group of the CMH-17 organization within the framework of a recent round-robin effort during which a series of crush tests on a corrugated plate specimen were performed. The test results were made available to the authors who participated in the round-robin effort [2].
1120
V.S. Sokolinsky et al. / Composites: Part A 42 (2011) 1119–1126
2.1. Geometry and composite lay-up The cross-section of the corrugated plate specimen can be obtained by repeating a semicircular segment of radius 6.35 mm (0.25 in.) three times at alternating sides with respect to the midplane (see Fig. 1). The specimen has straight edges of material on each side of the corrugation (end-lips) for additional stability. The test specimens were manufactured of a carbon–epoxy TORAYCA fabric (see Table 1) by press-molding through a set of aluminum matching tools [1]. The [0/90]2S composite lay-up was 2 mm (0.079 in.) thick. To facilitate stable crushing, a single-sided 45-degree chamfer was machined at the bottom of the specimen end after cure. Each specimen was 76.2 mm long (3 in.) and 50.80 mm wide (2 in.).
2.2. Experimental set-up The specimens were tested in the vertical configuration standing free on an as-milled dry steel surface. The crushing plate of the experimental rig was free to slide along four posts. A self-aligning sphere was used to introduce the crushing load from the test frame onto the specimen. Because the corrugated specimens are self-supporting, no supports were used on the specimens. The testing was performed at a quasi-static speed of 0.2 mm/s (0.5 in./ min). Fig. 2 shows typical views of the crushed corrugated fabric plate.
3. Constitutive model for fabric-reinforced composites This section describes a proposed in-plane constitutive model for a fabric composite ply that has been implemented in Abaqus/ Explicit. The fabric-reinforced composite ply is modeled as a homogeneous orthotropic material that is capable of sustaining progressive stiffness degradation due to fiber/matrix cracking and plastic deformation under shear loading. Fig. 2. Crushed corrugated plate specimen: (a) side view; (b) top view (image courtesy of Prof. Paolo Feraboli, the University of Washington, Seattle). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
3.1. Elastic continuum damage mechanics model
Fig. 1. Schematic of a corrugated fabric panel (all dimensions are in millimeters). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
The fiber directions of woven fabric reinforcement are assumed to be orthogonal. The constitutive stress–strain relations are formulated in a local Cartesian coordinate system with the base vectors aligned with the fiber directions. The in-plane elastic stress–strain relations are given according to orthotropic damaged elasticity [3]:
8 9 > < e11 > =
Variable
1 ð1d1 ÞE1
6 m21 e22 ¼ 6 4 E2 > > : el ; e12 0
Table 1 Material properties of the TORAYCA T700/2510 fabric [2]. Description
2
mE121
0
1 ð1d2 ÞE2
0
0
1 ð1d12 Þ2G12
Value T
Longitudinal modulus (GPa) Transverse moduli (GPa) Principal Poisson’s ratio Shear moduli (GPa) Longitudinal tensile strength (MPa) Longitudinal compressive strength (MPa) Transverse tensile strength (MPa) Transverse compressive strength (MPa) In-plane shear strength (MPa) Mode-I fracture toughness (J/m2) Mode-II fracture toughness (J/m2)
E11 E22 = E33
m12 G12 = G23 = G31 X1+ X1 Y2+ Y2 S GIC GIIC
55.9 54.4 0.042 4.2 911.3 704.0 770.1 698.2 131.6 504 1566
38 9 > r11 > = 7< 7 r22 5> > : ;
ð1Þ
r12
where e ¼ f e11 ; e22 ; eel12 g is the elastic strain vector; {r = r11, r22, r12}T is the stress vector; E1 and E2 are the Young’s moduli in the principal orthotropic directions; G12 is the in-plane shear modulus; m12 and m21 are the principal Poisson’s ratios; d1 and d2 are the damage parameters associated with fiber fracture along the principal orthotropic directions; and d12 is the damage parameter associated with the matrix micro-cracking due to the in-plane shear deformation. The damage parameters vary between 0 and 1 and represent stiffness degradation caused by microdamage in the material.
1121
V.S. Sokolinsky et al. / Composites: Part A 42 (2011) 1119–1126
Note that there is no damage parameter associated with the degradation in Poisson’s ratios. As a result, m12 and m21 are reduced by factors (1 d1) and (1 d2), respectively, with progression of damage. To differentiate between the tensile and compressive fiber failure modes, the corresponding damage variables are calculated based on the stress state in the fiber directions as follows (i = 1, 2):
di ¼ diþ
hrii i hrii i þ di jrii j jrii j
ð2Þ
where di+ and di represent the tensile and compressive damage variables of the fibers in the ith direction, respectively. The damage evolution equations are assumed in the following form:
~ 1þ Þ; d1 ¼ d1 ðr ~ 1 Þ d1þ ¼ d1þ ðr ~ 2þ Þ d2 ¼ d2 ðr ~ 2 Þ d2þ ¼ d2þ ðr ~ 12 Þ d12 ¼ d12 ðr
ð3Þ
/12 ¼
~ iþ and r ~ i (i = 1, 2) are defined as: where the effective stresses r
r~ iþ ¼ r~ 12 ¼
hrii i 1 diþ
r~ i ¼
r12
hrii i 1 di
ð4Þ
pffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffi r~ iþ ¼ 2Ei Y iþ r~ i ¼ 2Ei Y i pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi r~ 12 ¼ 2G12 Y 12
ð5Þ
where the thermodynamic forces Yi+, Yi, and Y12 are responsible for damage progression. Thus, it can be seen that fiber and shear damage modes are uncoupled. At any given time, the elastic domain is defined in terms of the damage activation functions Fi+, Fi, and F12 as:
F iþ ¼ /iþ r iþ 6 0 F i ¼ /i r i 6 0 F 12 ¼ /12 r 12 6 0
ð6Þ
here ri+, ri, and r12 are the tensile, compressive, and shear damage thresholds, respectively; and the functions ui+ and ui provide a criterion for fiber failure and are given by (i = 1, 2):
r~ iþ X iþ
r~ 12
ð10Þ
S
where S is the shear strength; and, after shear damage initiation (u12 = 1), the damage threshold, r12, is given by
r12 ðtÞ ¼ max /12 ðsÞ
ð11Þ
s6t
1 d12
The effective stresses are directly related to the thermodynamic forces, which are conjugate to the damage variables, through the following expression (i = 1, 2):
/iþ ¼
where Lc is the characteristic length of the element; Gaf is the fracture energy per unit area under uniaxial tensile or compressive loading; and g a0 is the elastic energy per unit volume (elastic energy density) at the point of damage initiation. The formulation of the damage evolution law ensures that: (a) the damage variables are monotonically increasing quantities (i.e., d_ a P 0) and (b) the correct amount of energy is dissipated when the lamina is subjected to uniaxial loading along the fiber directions irrespective of the element size. The latter is achieved by formulating the damage evolution law in Eq. (9) using the characteristic element length, Lc. Thus, the use of Lc is critical to ensure objectivity of the numerical results with respect to mesh refinement. The formulation imposes an upper bound, Lmax ¼ Gaf =g a0 , on the characteristic element length, Lc, of the finite element mesh. Similarly, the function u12, which provides a criterion for the initiation of shear damage in the matrix, is assumed to be
/i ¼
r~ i X i
ð7Þ
where Xi+ and Xi are the tensile and compressive strengths, respectively, for uniaxial loading along the fiber directions. The damage thresholds ri+ and ri are initially set to unity. After damage initiation, i.e., when ui+ = 1 or/and ui = 1, they increase with damage progression according to:
riþ ðtÞ ¼ max /iþ ðsÞ s6t
ri ðtÞ ¼ max /i ðsÞ
ð8Þ
s6t
The definition ensures that the damage thresholds are nondecreasing quantities, i.e., that r_ a ðtÞ P 0 (a = ±1, ±2). The damage thresholds are assumed to obey the Kuhn–Tucker complementary conditions and the consistency condition ðr_ a F_ a ¼ 0Þ. Following the work in Ref. [4], the evolution of the damage variables is given by the following expression (a = ±1, ±2):
" # 1 2g a0 Lc da ¼ 1 exp a ðra 1Þ ra Gf g a0 Lc
ð9Þ
Furthermore, in accordance with Ref. [3], it is assumed that the shear damage variable increases with the logarithm of r12 until a maximum value of d12 is reached, i.e., max
d12 ¼ minða12 lnðr 12 Þ; d12 Þ
ð12Þ
max d12
where a12 > 0 and 6 1 are material properties. Please note that max the values of S, a12, and, d12 can be determined using a calibration procedure presented in the Appendix. 3.2. Elastic–plastic shear model The in-plane shear response is dominated by the non-linear behavior of the matrix that exhibits both stiffness degradation due to matrix microcracking and plasticity. The elastic part of the shear response was discussed previously. Here, the plastic shear response of the material is considered. The matrix response may be inelastic due to extensive cracking or plasticity. This leads to permanent deformations in the ply upon unloading. To account for these effects, a classical plasticity model with an elastic domain function and a hardening law, which is applied to the effective stresses in the damaged material, is used. The elastic domain function, F, is given by:
~ 12 j r ~ 0 ðepl Þ 6 0 F ¼ jr
ð13Þ
and the hardening law is taken to be of the form:
r~ 0 ðepl Þ ¼ r~ y0 þ Cðepl Þp
ð14Þ
~ y0 is the initial effective shear yield stress; C and p are coefwhere r ficients; and epl is the equivalent plastic strain due to shear defor~ y0 , C, and p can be determined using the mation. The values of r calibration procedure described in the Appendix. 4. Description of the finite element model In this section, a finite element model of the corrugated specimen described in Section 2 is detailed. Both intralaminar (in-plane) and interlaminar (delamination) failure mechanisms are considered in order to simulate the real physics of composites crushing. The finite-element model of the corrugated composite fabric plate was created in Abaqus/CAE [5]. The crush initiator in the form of a single-sided 45-degree chamfer was accounted for in the model (see Fig. 3). Each of the eight fabric plies was individually meshed
1122
V.S. Sokolinsky et al. / Composites: Part A 42 (2011) 1119–1126
with continuum shell elements (SC8R). A uniform mesh with an element size of 1 mm by 0.5 mm (0.04 in. by 0.02 in.) was used for each ply. Continuum shell elements in Abaqus have the geometry of a three-dimensional solid element but their kinematic and constitutive behaviors are similar to those of conventional shell elements. This type of element allows for the accurate modeling of stacked composite plies. The finite element model used for this simulation consisted of approximately one-half million degrees of freedom.
4.1. Intralaminar response of the plate The in-plane response of the fabric plies was modeled using the Abaqus constitutive model described previously. This model requires the input of the: (a) elasticity constants; (b) damage initiation coefficients; (c) damage evolution coefficients; and (d) shear plasticity coefficients. The elasticity constants include the Young’s moduli in the fiber directions, the principal Poisson’s ratio, and the in-plane shear modulus. The damage initiation coefficients include the tensile and compressive strengths along the fiber directions and the shear strength at the onset of shear damage. The values of the user material constants and damage initiation coefficients are given in Table 1. Please note that the strength and toughness properties of the TORAYCA fabric were adopted directly from Table 1 without any further adjustments.
The damage evolution coefficients are represented by the tensile and compressive fracture energies per unit area along the fiber max directions and by the parameters a12 and d12 in Eq. (12). The fracture energies can be determined through the testing procedure described in Ref. [6]. However, for the present work this testing was not performed, and average values of the fracture energies G1þ and f G1 for similar material systems from Ref. [6] were utilized instead f (see Table 2). The critical values of the fracture energies along the 1þ fiber direction 2 were calculated so that the ratios G2þ and f =Gf 1 G2 =G were approximately equal to the ratios Y /X and Y2/ 2+ 1+ f f X1, respectively. max The parameters a12 and d12 should be determined using the calibration procedure as described in the Appendix. For the present work, however, the cyclic tensile tests were not conducted so the max values of a12 = 0.15 and d12 ¼ 1 were used in the analysis based on the calibration results for similar material systems. The shear plastic coefficients include the initial effective shear ~ y0 , and the coefficients C and p (see Eq. (14)). These yield stress, r values are also determined as shown in the Appendix. In the absence of the cyclic tensile tests, these coefficients were also adopted based on the calibration results for similar material systems (see Table 3). Please note that shearing of the fabric ply is not a dominant failure mode for this analysis. Therefore, the use of typical values for the damage and plasticity shear material parameters do not have a strong influence on the results presented here. This statement was numerically verified by perturbing the discussed shear parameters over a significantly wide range of values without changes in the results. 4.2. Interlaminar response of the plate In order to account for the adhesive bonds between adjacent plies and model the delamination phenomenon, a generalized traction-separation law was used. Damage initiation at the interfaces was defined using the quadratic nominal stress criterion according to which damage is assumed to initiate when the quadratic interaction function involving the stress ratios reaches a value of one [5]:
( )2 ( )2 ( )2 ht n i ts tt þ 0 þ 0 ¼1 t 0n ts tt
ð15Þ
where t = {tn, ts, tt}T is the nominal traction vector and the interface damage initiation properties used in the present analysis are listed in Table 4. These properties are usually difficult to determine experimentally with good accuracy; therefore, they can be used as calibration parameters, if required. In the present work, the variation of the damage initiation parameters within 15% from the typical values given in Table 4 did not influence the analysis results. A 15% variation was used because it represents the typical maximum deviation of the material properties in manufactured composite parts. Progression of damage at the interfaces was modeled using the critical mixed mode energy behavior based on the Benzeggagh–Kenane (BK) law [5]:
Table 2 Fracture energies per unit area along the fiber directions.
Fig. 3. Abaqus/CAE model of the corrugated composite fabric plate: (a) top view; (b) side view (the finite element mesh is not shown). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
Description
Variable
Value (kJ/m2)
Tensile fracture along fiber direction 1
G1þ f G1 f G2þ f G2 f
125
Compressive fracture along fiber direction 1 Tensile fracture along fiber direction 2 Compressive fracture along fiber direction 2
250 95 245
1123
V.S. Sokolinsky et al. / Composites: Part A 42 (2011) 1119–1126 Table 3 Shear plasticity coefficients. Description
Variable
Value
Initial effective shear yield stress (MPa) Coefficient in the hardening equation (Eq. (14)) Power term in the hardening equation (Eq. (14))
r~ y0
195 1053 0.41
G g S GCn þ GCs GCn ¼ GC GT
C p
ð16Þ
where GS = Gs + Gt; GT = Gn + GS; Gn, Gs, and Gt refer to the work done by the normal and the first and the second shear forces acting in the interface, respectively; GCn and GCs are critical fracture energies required to cause failure in the normal and shear directions, respectively ðGCt ¼ GCs Þ; GC is the total critical mixed-mode fracture energy; and g is the power coefficient taken to be 2.284. The values of the required normal ðGCn Þ and shear ðGCs Þ mode critical fracture energies were taken from the last two rows of Table 1. After damage initiation, the evolution of the damage variable until final failure was modeled using an exponential softening law [5].
(S4R). The crushing load was applied to the composite plate through another rigid steel plate moving with a constant velocity of 200 mm/s (470 in./min). This loading rate is approximately eight times higher than the maximum quasi-static test speed reported in Ref. [1] and was used to achieve a reasonable run time for the explicit dynamics analysis. This artificial increase of the loading velocity has no effect on the analysis results because according to previous work (see Ref. [1]), dynamic effects become significant beyond 500–1000 mm/s (1200–2400 in./min). As pointed out in Ref. [7], the friction between the surfaces of the fronds and the crushing surface and between the fronds and the debris wedge plays an important role in the analysis. Here, the friction coefficient between the surfaces of the debonded composite plies was assumed to be 0.3 [7]. The friction coefficient between the composite plate and the rig was taken to be 0.12 based on the data available in the literature. The interactions between the debonded composite plies and between the composite plies and the rig were modeled using the general contact capability of Abaqus/Explicit. After delamination between any two adjacent plies occurs, the general contact definition is automatically updated to account for possible post-delamination contact.
4.3. Element deletion
4.5. Mass scaling
Failed continuum shell elements were removed from the model to prevent their excessive distortion and premature termination of the analysis. In the present work, damage-based element deletion is combined with a deformation-based deletion criterion. The damage-based element deletion is activated when any damage variable along the fiber directions or the equivalent plastic strain due to shear deformation reaches a maximum specified value. The deformation-based deletion criterion is used when either the tensile or compressive principal logarithmic strain reaches its maximum or minimum specified value, respectively (see Table 5). Additionally, the detached composite fragments are deleted when they move far enough away from the crush zone. Conceptually, the zone is defined by a three-dimensional box that encompasses the crush area. With the aid of Abaqus/Explicit user subroutine VUSDFLD, the detached parts are detected and deleted when they move out of the range of this box. This is done to prevent non-physical numerical distortions of the detached fragments that may cause premature termination of the analysis.
The mass of the composite plate was increased by three orders of magnitude at the beginning of the analysis using the fixed mass scaling capability of Abaqus/Explicit. Similar to the artificial increase of the loading velocity discussed previously, this numerical device was employed to achieve a reasonable run time. A comparison between the kinetic and strain energies confirmed that the quasi-static character of the problem was not changed as a result.
4.4. Loads and boundary conditions The composite plate was positioned vertically on a rigid steel plate which was modeled using conventional shell elements
5. Simulation results and discussion The simulation required 9 h of wall time to run on a 16-core Linux cluster. The deformed shape of the corrugated fabric plate at the end of the crushing simulation is shown in Fig. 4. Note that the deletion of damaged elements in the plies creates detached parts that are removed from the model as explained previously. A comparison between Figs. 2 and 4 demonstrates very good qualitative correlation between the experimental and numerical results. Specifically, it can be seen from this comparison that the Abaqus/Explicit simulation reproduced important characteristic features of the crushing response of composites such as the frond formation in the composite plate, the propagation of the main delamination through the plate, and the accumulation of debris between the loading rig and the composite plate.
Table 4 Interface damage initiation properties. Description
Variable
Value (MPa)
Maximum nominal stress in the normal-only mode
t 0n
54
Maximum nominal stress in the first shear direction (for a mode that involves separation only in this direction)
t 0S
70
Maximum nominal stress in the second shear direction (for a mode that involves separation only in this direction)
t 0t
70
Table 5 Element deletion controls. Description
Variable
Maximum value of the damage variable in the fiber directions used in the damage-based deletion criterion Maximum value of the equivalent plastic strain used in the damage-based deletion criterion
dmax
Maximum (positive) principal logarithmic strain used in the deformation-based deletion criterion Minimum (negative) principal logarithmic strain used in the deformation-based deletion criterion
epl max ^emax ^emin
Value 0.9999 0.5 0.1 0.1
1124
V.S. Sokolinsky et al. / Composites: Part A 42 (2011) 1119–1126
Fig. 4. Simulated deformed shapes of the corrugated composite plate: (a) front view (XOZ-plane); (b and c) side views (YOZ-plane); and (d) isometric view. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
Moreover, a comparison of Fig. 2b with Fig. 4b–d reveals that the simulation was capable of correctly predicting an important qualitative difference in the delamination propagation through the end-lips of the composite plate. Namely, one of the end-lips was split in half by the propagation of the main delamination through its height all the way from top to bottom. In the other end-lip, the main delamination traveled through approximately one-third of its height. It should be emphasized that the use of a uniform mesh as described in Section 4 plays an important role in the accurate simulation of delamination propagation through the end-lips of the composite plate. A close match between the experimental (Ref. [2]) and simulated load displacement curves is demonstrated in Fig. 5. Both the peak and average crush forces were predicted accurately. It can be seen from the figure that the height of the crushed plate was reduced by half with respect to its undeformed configuration. Please note that both the experimental results and the simulation
results were filtered using a SAE filter. However, the frequency cut-off used for the experimental results was 600 Hz, whereas for the simulation results, 100 Hz frequency cut-off was used. The lower cut-off frequency was used to remove the numerical noise characteristic of an explicit analysis crushing simulation. The numerical results for the energy absorbed (EA), which is calculated as the total area under the load–displacement curve in Fig. 5, are shown versus stroke in Fig. 6. The term stroke here is used in reference to the length of the plate that has been destroyed during crushing [1]. As can be seen from the figure, the experimental and simulated curves nearly coincide with each other. Finally, a comparison between the experimental and simulated specific energy absorption (SEA) of the corrugated plate is shown in Fig. 7. SEA shows the ability of a material to dissipate energy; it is measured in Joules per gram and usually varies for composites between 15 and 80 J/g [1]. SEA is computed by dividing EA by the crushed mass of the plate, which can be expressed as the product
Fig. 5. A comparison between the experimental and numerically-predicted load–displacement curves. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
V.S. Sokolinsky et al. / Composites: Part A 42 (2011) 1119–1126
1125
Fig. 6. A comparison between the experimental and numerical results for the absorbed energy versus stroke. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
of the stroke, cross-sectional area of the plate, and its density. As can be seen from Fig. 7, the numerical and experimental results agree well, especially in the second half of the crushing process. The largest discrepancy between the results is observed at the early loading stages when the chamfer part of the plate is crushed. Taking into consideration the extremely complicated nature of the crushing process, the overall correspondence between the measured and simulated SEA curves is judged to be acceptable. Some comments on the robustness of the presented methodology are in order. As has been mentioned previously, the strength and interface toughness properties of the TORAYCA fabric were adopted directly from Ref. [2] without any further adjustments, and average values for similar material systems from Ref. [6] were used for the critical fracture energies per unit area along the fiber directions. Furthermore, the interface damage initiation properties were based on the data available in the literature, and the model response did not change as a result of a realistic deviation from these properties. Moreover, the perturbations of the shear damage and plasticity parameters over a wide range of values did not have any appreciable impact on the results. Therefore, it can be
concluded that the presented methodology exhibits robustness with respect to the input material properties. 6. Conclusions A finite element model of a corrugated fabric composite plate subject to quasi-static crushing has been presented. A physicsbased approach that accounts for both intralaminar (in-plane) and interlaminar (delamination) failure mechanisms has been used to model the crushing behavior of the composite plate. The in-plane response of the fabric plies has been modeled with the aid of the proposed constitutive model implemented in Abaqus/Explicit, whereas delamination has been simulated using the Abaqus cohesive surface capability. The results of an Abaqus/Explicit simulation show very good quantitative and qualitative agreement with experimental data and exhibit robustness with respect to the input material properties. Thus, the proposed methodology shows considerable promise for realistic crush simulations of composite structures. This may permit substitution of costly experimental testing with numerical simulation for the
Fig. 7. A comparison between the experimental and numerical results for the specific energy absorption. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
1126
V.S. Sokolinsky et al. / Composites: Part A 42 (2011) 1119–1126
Fig. A2. Schematic of the calibration of the: (a) shear damage parameters; (b) shear hardening curve. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) Fig. A1. Schematic of typical response of a fabric-reinforced composite. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
crashworthiness design of composites. Further research is needed to investigate the limits of the proposed numerical methodology by applying it to a wide variety of composites crush simulations. This will be the focus of future work. Acknowledgements The first author would like to express his sincere gratitude to Qingping Huang and Fernando Carranza of SIMULIA whose technical advice enhanced the quality of the present work. Appendix. Calibration procedure for the shear response parameters used in the constitutive model for fabric-reinforced composites The calibration of the shear response parameters uses the results of cyclic tensile tests for ±45° laminates. An idealized representation of the outcome of such tests (the actual response exhibits hysteretic behavior upon unloading and reloading) appears in Fig. A1. First, the level of damage is calculated from the ratio of the unloading stiffness to the initial (undamaged) elastic stiffness (see Fig. A1). This generates pairs of stress–damage values, (r12, d12), for each unloading branch of the cyclic shear stress– ~ 12 ; d12 Þ as strain curve. These are represented as coordinates ðln r
shown in Fig. A2a. Then the slope of a linear fit curve for these pairs provides the value of a12, and its intersection with the abscissa, the value of S. If the damage data indicate saturation (see Fig. A2a), the max corresponding value of d12 (d12 6 1) should be used as d12 ; othermax wise, d12 ¼ 1 is adopted. Furthermore, for each unloading branch of the cyclic shear stress–strain curve (Fig. A1), the plastic strain epl 12 at the onset of unloading is determined from the value of residual deformation in the unloaded state. Then the coefficients of a nonlinear fitting ~ 12 ; epl curve for the points ðr 12 Þ at the onset of unloading are used as the plastic coefficients of the fabric material as illustrated in Fig. A2b. References [1] Feraboli P. Development of a corrugated test specimen for composite materials energy absorption. J Compos Mater 2008;42(3). [2] CMH-17 meeting results. 6th CMH-17 Crashworthiness Working Group Meeting, Ottawa, CA; 5 August 2008. [3] Johnson AF. Modelling fabric reinforced composites under impact loads. Composites Part A 2001;32:1–2. [4] Maimí P, Camanho PP, Mayugo J-A, Dávila CG. A thermodynamically consistent damage model for advanced composites. Technical Report NASA/TM-2006214282, NASA; 2006. Ò [5] Abaqus 6.10 Analysis user’s manual. Dassault Systèmes; 2010. [6] Pinho ST, Robinson P, Iannucci L. Fracture toughness of tensile and compressive fibre failure modes in laminated composites. Compos Sci Technol 2006;66(13). [7] Pinho ST, Camanho PP, De Moura MF. Numerical simulation of the crushing process of composite materials. Int J Crashworthiness 2004;9(3).