- Email: [email protected]
A multi-scale finite element approach for modelling damage progression in woven composite structures
Composite Structures 94 (2012) 977–986
Contents lists available at SciVerse ScienceDirect
Composite Structures journal homepage: www.elsevier.com/locate/compstruct
A multi-scale finite element approach for modelling damage progression in woven composite structures P. Römelt, P.R. Cunningham ⇑ Department of Aeronautical and Automotive Engineering, Loughborough University, Leicestershire LE11 3TU, United Kingdom
a r t i c l e
i n f o
Article history: Available online 6 November 2011 Keywords: Woven composites Multi-scale modelling Damage
a b s t r a c t A significant challenge in the numerical modelling of composite structures with a multi-axis fibre architecture is the reproducibility of the textile mechanics [1]. A numerical analysis procedure for woven composite structures using a multi-scale finite element approach has been developed, and is presented in this paper. The approach is demonstrated for a flat two-dimensional woven glass/epoxy laminate. Digital microscopy is used to estimate tow cross-section and path, and quantify the amount of variation of these parameters. This data is used to generate both a meso-scale model of a single unit cell as well as a macroscale model of the complete structure. Numerical results from the proposed approach are compared to experimental stress–strain data, which show good agreement in the lower strain range. Ó 2011 Elsevier Ltd. All rights reserved.
1. Introduction In order to be used efficiently, a design procedure for twodimensional woven composites needs to include a validated damage model. To correctly model damage initiation and propagation, models which take into account the different modes of failure, such as tow matrix de-bonding and crack propagation within the tows and matrix are required. Since these failure events can be best observed on the meso-scale, a multi-scale approach which uses meso-scale models to generate material input parameters for macro-cells would be the most suitable to model the behaviour of a woven macrostructure. Three different approaches to model multi-axial structures can be found in the literature, the full finite element approach [2], the voxel technique [3] and the binary model [4]. The full finite element approach [2] uses precise geometrical data to generate a high resolution model of a unit cell. This allows for an estimation of the stress distribution within the unit cell, which in turn allows for realistic modelling of damage initiation and progression within the cell’s limit. Lomov et al. have shown that this method has been used successfully to predict damage initiation [5]. However, models of a single unit cell require a high number of elements, resulting in high computational costs, which makes this method unsuitable for modelling structures larger than a few unit cells. Bogdanovich defines voxels as sub-meso-volumes [3] and Crookston et al. [6] have developed a modelling approach based on such voxels. A two-dimensional grid is laid over the top of a geometrical ⇑ Corresponding author. E-mail address: [email protected] (P.R. Cunningham). 0263-8223/$ - see front matter Ó 2011 Elsevier Ltd. All rights reserved. doi:10.1016/j.compstruct.2011.10.024
unit cell model and the architecture is analysed in the throughthickness direction. Average material properties are calculated depending on the through-thickness architecture, which are then assigned to three-dimensional elements generated at the twodimensional grid point positions. The method reduces both meshing effort and the size of the finite element model considerably, resulting in reduced computational effort [6]. However, by averaging material properties the stress distribution is also averaged over the voxel. This means that stress concentrations, especially on the sites of damage initiation like the tow/matrix interface, are averaged out. Crookston et al. successfully applied this method to reproduce the low strain linear behaviour of a plain weave reinforced glass/ polyester composite but could not reproduce the nonlinear behaviour for higher strain values [6]. The binary model [4] established by Cox and Yang utilises onedimensional elements, which run along the yarn centre line. These elements represent the axial properties of the yarns. The other yarn properties are incorporated in a three-dimensional representative matrix element, referred to as ‘‘effective medium elements’’, using micromechanics. As with the voxel technique, the binary model reduces meshing effort and computational cost while averaging out stress concentrations. Cox and Yang successfully used the binary model to predict ultimate stress of tri-axially braided carbon/ epoxy composites and reproduce the linear part of the stress– strain curve but could not reproduce the nonlinear part [4]. The aim of this paper is to demonstrate a multi-scale procedure using digital microscopy and finite elements to model damage in woven composite macro-structures. Making use of the advantages of the averaging techniques and the full finite element model, this procedure attempts to predict the stiffness degradation of macrostructures due to meso-scale damage. The approach taken is
978
P. Römelt, P.R. Cunningham / Composite Structures 94 (2012) 977–986
Nomenclature A a, b C d D, E Et KIc K L t t0 t x, y, z
r
cross-sectional area (mm2) major, minor elliptical half axes (mm) sine function amplitude (mm) damage parameter (unity) sine function parameters (mm) tangent modulus (GPa) fracture toughness (MPa (m)0.5) stiffness (N/mm) sinusoidal wavelength (mm) traction (MPa) allowable traction (MPa) traction in undamaged material (MPa) coordinate values tensile stress (MPa)
illustrated in Fig. 1. Digital Microscopy was used to find the required geometric input parameters for a meso-scale model of a unit cell of a woven macrostructure using a full finite element approach. The results from the full finite element model were then used to define the material properties for an equivalent binary model of a unit cell, which shows the same behaviour as the full finite element model. The macrostructure was then modelled using these equivalent binary unit cells and stress–strain predictions were made for a simple unidirectional tensile load case. The results from a series of tensile tests of glass–epoxy specimens were used to validate the finite element predictions.
s
shear stress (MPa)
Vectors t e
tangent vector directional vector
Subscripts 1D one dimensional 3D three dimensional n normal direction s, t transverse shear directions spring nonlinear spring element x, y, z coordinate directions
2. Tensile testing and microscopic analysis 2.1. Test specimen manufacture A 2.5 mm thick laminate plate measuring 330 mm 330 mm was manufactured from thirty plies of plain weave glass fabric (Gurit RE86P) with a layup of [0/90]30, which was infused with Gurit Prime 20-LV epoxy resin using a double bag vacuum infusion procedure to afford better consolidation and even thickness of the cured laminate [7]. Nine test specimens were cut from this laminate with the geometry and dimensions as shown in Fig. 2. A rectangular rosette strain gauge with a 5 mm gauge length and a resistance of 120 X was bonded to one side of each of the test specimens in the position shown in Fig. 2. 2.2. Tensile testing Tensile tests were conducted according to ASTM D3039. Although the use of strain gauges for mechanical testing of woven composites has been criticised in the literature [3], Lang and Chou [8] have demonstrated that good results can be achieved if the gauge length is longer than a single unit cell length. The resulting stress–strain curves from the tensile tests are shown in Fig. 3, with the mean stress–strain curve also marked. There is an average variation of about 7% about the mean in the resulting stress–strain data, with failure strengths ranging from 375 MPa to 450 MPa with a +9.8% to 8.5% variation at failure. This is consistent with findings from Daggumati et al. [9], who showed that damage initiation and progressive stiffness degradation is highly dependent on the amount of variation in tow path in both warp and weft direction and tow cross sectional area. 2.3. Microscopic analysis A total of 123 microscopic specimens were cut using the same laminate from which the tensile test specimens were obtained. A
Fig. 1. Flow chart of numerical procedure.
Fig. 2. Tensile test specimen, geometry and strain gauge position.
P. Römelt, P.R. Cunningham / Composite Structures 94 (2012) 977–986
979
Fig. 3. Tensile test results.
Fig. 4. Laminate cutting scheme.
Fig. 5. Magnified images from microscopic analysis: (a) sine fit and (b) ellipse fit.
typical cutting pattern for a portion of the laminate can be seen in Fig. 4. After cutting the specimens were cast in a dyed epoxy resin to provide contrast under the microscope, and the specimens were polished using a five-step process. Tow path data was generated using 20 times magnified microscopic images as shown in Fig. 5a. This was done to ensure that the averaged tow path and cross-section parameters were the same for the numerical models as for the experimental specimen. Lines were placed along the centre of the tow in the 20 times magnified
microscope image and the coordinates of the start and end points of these lines were recorded. A sine function was fitted through those points using a least-square fit algorithm as given in the following equation:
zðxÞ ¼ D sin where
2px 2px þ E cos L L
ð1Þ
980
P. Römelt, P.R. Cunningham / Composite Structures 94 (2012) 977–986
Fig. 6. Comparison of least-square fit results with microscopic data.
Table 1 Geometrical properties of two-dimensional woven laminate meso-structure. Parameter (unit)
Average value ± standard deviation
C (mm) L (mm) a (mm) b (mm) A (mm2)
0.025 ± 0.009 1.604 ± 0.164 0.187 ± 0.066 0.037 ± 0.007 0.022
Table 3 Published E-glass fibre properties [10]. E (GPa) G (GPa) m (–)
72.4 30.0 0.2
Table 4 Gurit Prime 20LV/epoxy properties [11].
P P 2 ni¼1 zi sin 2pLxi E ni¼1 sin 4pLxi D¼ P n ni¼1 cos 4pLxi
E (GPa) G (GPa) m (–) Ultimate tensile stress, rult (MPa) Ultimate shear stress, sult (MPa) KIc (MPa (m)0.5)
ð2Þ
3.5 1.3 0.35 73 137 3.69
and
P P 2 ni¼1 zi cos 2pLxi D ni¼1 sin 4pLxi E¼ P n ni¼1 cos 4pLxi
ð3Þ
Tow cross-section data was generated using 80 times magnified images of the same specimen as shown in Fig. 5b. Lines were placed along the tow cross-section in the 80 times magnified microscopic image and a standard ellipse was fitted through the start and end points using a least-square fit algorithm as given in Eq. (4). A total of 385 tow paths were analysed as well as 449 tow cross-sections.
rffiffiffiffiffiffiffiffiffiffiffiffiffiffi y2 a2
zðyÞ ¼ b 1
ð4Þ
where
a2 ¼
Pn
4 i¼1 yi
Pn
2 i¼1 yi
Pn i¼1
Pn
4 i¼1 zi
2
b ¼
ð5Þ
ðy2i z2i Þ 2
b
and
Pn
2 i¼1 zi
Pn i¼1
ð6Þ
ðy2i z2i Þ
a2
Fig. 6 shows a comparison of the measured points from an example specimen of the microscopic analysis and the curves fitted through them using the least square algorithms given by Eqs. (1) and (4). Fig. 6a shows a good correlation with exceptions at the extremes, especially at the beginning of the measurements. The correlation for the ellipse is reasonable as shown in Fig. 6b but gives a larger result for both half-axes of the ellipse compared with the measured data. The microscopic analysis showed variations between 10% and 100% for geometric parameters. Daggumati et al. [9] showed that these variations lead to variations in the point of damage initiation in a unit cell. For simplicity reasons only average values are used for geometric parameters in this paper. The average values as well as standard deviations for tow path amplitude C (mm) and wavelength L (mm) and the cross-section
Table 2 Weave properties. Warp
Weft
Filament type
Weave type
Nominal (g/m2)
Fibre (tex)
Count (ends/cm)
Theoretical (g/m2)
Count (ends/cm)
Theoretical (g/m2)
E-glass
Plain
85
34
12.0
41.5
12.5
43.3
P. Römelt, P.R. Cunningham / Composite Structures 94 (2012) 977–986
981
ellipse parameters a (mm) and b (mm) are listed in Table 1. The cross-sectional area A (mm2) for the tow ellipse was calculated using the averaged ellipse parameters a and b from the fit. As can be seen from the results, the standard deviation of the fit for the samples measured is relatively low, with a maximum standard deviation of 10% occurring for the tow path wavelength. 3. Material properties of constituents A number of parameters of the constituent materials, such as Young’s modulus and Poisson’s ratio, were required for the modelling approach. However, data on the constituents is very limited. Therefore a number of parameters had to be assumed with reference values found in the literature. Isotropic linear elastic material behaviour was assumed for both the matrix and the glass fibre tows. The two-dimensional glass fibre weave used in the laminate was Gurit RE86P E-glass plain weave, the properties for which are listed in Table 2. The material properties of the E-glass fibre, as shown in Table 3, were taken from the MATWEB database [10]. The resin used was Gurit Prime 20LV, and the material properties are listed in Table 4. Not all of the required resin properties were available from the manufacturer, therefore certain properties, such as fracture toughness, were also taken from the MATWEB database for epoxy materials [11]. Finally, assumptions had to be made for the interface strength between tows and matrix, on which numerous studies were
Fig. 9. Traction separation formulation.
available in the literature. The value of 62.5 MPa from Ref. [12] was chosen for this analysis. 4. Numerical modelling 4.1. Full finite element unit cell model A unit cell of the two-dimensional weave pattern, with geometrical parameters obtained from the microscopic analysis, was modelled using a full finite element approach [2]. PATRAN was used as
Fig. 7. Full finite element unit cell model and equivalent binary unit cell model.
Fig. 8. Unit cell boundary conditions.
982
P. Römelt, P.R. Cunningham / Composite Structures 94 (2012) 977–986
Fig. 10. Graphical representation of analytical approach.
Fig. 11. Eight-block equivalent binary unit cell model.
a pre-processor to model geometry and for meshing. An input file was generated and amended to include material and section properties, boundary conditions, solver parameter and requested output data. ABAQUS/Standard ver. 6–10 was used as the solver, with ABAQUS/CAE for post-processing. The major challenge when generating full finite element unit cell models of multi-axial composites is to achieve an accurate description of the unit cell while keeping the computational resources required to solve the model
to a minimum [13]. Due to the highly elliptical cross-section of the tow and the sinusoidal tow path, coarse meshes tend to have high aspect ratios and skew angles. This can lead to significant errors in predicting the local stress–strain field [14]. Mesh refinement is usually used to mitigate this problem, which results in higher computational costs. However, because of the complex geometry of the model the mesh would have to be refined to a level that computational costs become unreasonably high. Therefore, some distorted elements in the mesh are deemed to be acceptable [15]. Eight-node brick elements were used to generate a mesh with a total of 117,096 degrees of freedom. For this mesh, 3% of the elements had an aspect ratio greater than the recommended value of 5 or skew angles greater than the recommended value of 30°, which was deemed an acceptable level for a model of this type [15]. Glaessgen et al. [15] achieved convergence for a plain weave glass fibre epoxy unit cell model using 50,000 degrees of freedom. The mesh used is shown in Fig. 7a, matrix elements have been omitted for an easier portrayal of the mesh. It has been demonstrated by Ivanov et al. [16] that boundary conditions on the meso-scale model depend on the position of the unit cell in the structure. However periodic boundary conditions, derived by Whitcomb et al. [17], were assumed in this paper, again for simplicity reasons. A constant displacement was applied to all nodes on the y–z faces of the unit cell in the axial direction. The faces in the transverse directions of the unit cell, x–y and x–z faces, were kept straight to account for the effect of neighbouring cells in the weave while allowing for Poisson’s effect. This also constrains rigid body rotations around all axes. The centre node of the model was constrained in both y and z direction to prevent rigid body motion in these directions. The boundary conditions are shown diagrammatically in Fig. 8. Global strains of the unit cell were calculated by dividing the prescribed displacement by the length of the unit cell. Global stresses were calculated by the summation of axial forces of every node on the left y–z face, which was then divided by the cross sectional area of the same face. Cohesive elements with a traction–separation formulation (Fig. 9) similar to Ridha et al. [18] were used to model crack
Fig. 12. Critical cross section of tensile test specimen.
Fig. 13. Linear analysis stress distribution at 2.5% strain.
P. Römelt, P.R. Cunningham / Composite Structures 94 (2012) 977–986
Fig. 14. Damage parameter at various global tensile strain values: (a) 0.37% strain, (b) 0.75% strain, (c) 0.75% strain, (d) 1.66% strain and (e) 2.5% strain.
983
984
P. Römelt, P.R. Cunningham / Composite Structures 94 (2012) 977–986
initiation, crack propagation and tow matrix de-bonding. A maximum traction criterion, as shown in Eq. (7), was used in conjunction with a damage variable d, which reduced the resulting traction t for the undamaged material at a given strain level as shown in Eq. (8). In other words, ABAQUS/Standard first calculates the tractions for the undamaged material at a given strain level and then multiplies these tractions with the damage parameter d. The damage parameter d is calculated based on the combined effect of normal and shear tractions within the element. It reflects overall damage within the cell and is therefore used to reduce all local element stiffnesses, independent of the cohesive element orientation.
( ) jt n j ts t t ¼1 max ; ; t 0n t0s t 0t tn ¼
ð1 dÞt n ; tn > 0 tn
ð7Þ ð8Þ
where
t s ¼ ð1 dÞt s
ð9Þ
and
t t ¼ ð1 dÞt t
ð10Þ
4.2. Binary unit cell model In the next step a similar model of the weave was made using the binary model approach (Fig. 7b) [4]. The binary model used three-dimensional hex elements to model transverse tow and surrounding matrix properties. Axial stiffness of the tows was modelled using one-dimensional spring elements, which were embedded in the three-dimensional element mesh. The properties of elements, i.e. spring stiffness for one-dimensional elements and linear isotropic material properties for three-dimensional elements, were calculated using the tow geometrical properties measured during the microscopic analysis and the rule-of mixtures for an equivalent unidirectional laminate with a fibre volume fraction comparable to the multi-axial composite. Unlike the initially proposed model by Cox and Yang [4], nonlinear one-dimensional spring elements were used to model damage in this work. This means all damage modes observed in the full finite element unit cell model were accounted for in the load– displacement relationship of the one-dimensional elements. The required stiffness at different strain levels of the unit cell can be estimated using the tangent modulus Et derived from the resultant true stress–true strain curve of the full finite element model. The stiffness of the one-dimensional element at different strain levels was calculated using the results from the full finite element model. For this approach one-dimensional elements and effective medium elements were modelled as a series of springs, as shown in Fig. 10. Four elements were assumed to be sufficient to reproduce the tow path and therefore model the axial stiffness of the tow. These four elements (K1D) were assumed to be in series with each other (Kspring) and in parallel with the three-dimensional solid element (K3D). The cross-sectional area and length of the equivalent binary unit cell is denoted by A and L respectively. The local orientation of the one-dimensional spring elements, which is different to the axial direction of the unit cell due to tow crimp, had to be considered when calculating the required stiffness of the individual one-dimensional element since the resultant stress–strain curve of the full finite element unit cell model was derived for the axial direction of the unit cell. This was done using the scalar product of the normalised local tow path tangent vector t with the global axial unit vector ex as shown in the following equations:
K 1D ¼ K spring
4 X ti ex jt ij i¼1
K spring ¼ ðEt E3D Þ
A L
ð11Þ
ð12Þ
In order to check the force–displacement relationship generated by the analytical approach a simple binary model was generated. However, assigning the boundary conditions that have been used in the full finite element model to a single equivalent unit cell was not possible because only a single three-dimensional solid element was used. Therefore a block of eight equivalent binary unit cells was modelled and the same boundary conditions applied. This block of equivalent binary unit cells has a total of 594 degrees of freedom and is shown in Fig. 11. 4.3. Binary macro-scale model In the final step a binary model of the critical cross-section of the test specimen was modelled, as shown in Fig. 12. Nine equivalent binary unit cells are placed next to each other in the widthwise and axial direction respectively and in order to achieve the correct thickness, thirty of these plies are stacked giving the same number of plies used in the tensile test specimen. This results in a total of 540,186 degrees of freedom in the model. Boundary conditions were chosen as to best represent the conditions at the critical cross-section. A uniform displacement was applied to the y–z faces of the model. This condition not only constrains axial movement but also rotation around the z and y axes. The centre node of the model was fixed in both the y and z direction to constrain rigid body motion in both transverse directions. Finally, to constrain rotation around the x axis the two nodes in the centre of the x–z faces were constrained in the z direction. 5. Results and discussion 5.1. Full finite element model A linear analysis without damage was conducted first to find the stress concentration points in the unit cell. To account for damage, cohesive elements were placed at the point of highest stress in the full finite element unit cell model using the results from the linear analysis. The elements were placed along the line of highest stress and around the interface between tow and matrix. Fig. 13 shows stress in the axial direction in the linear analysis, with stress concentrations near the edge of the weft tows and at the tow/matrix interface. Daggumati et al. [9] and Gao et al. [19]
Fig. 15. One-dimensional spring load–displacement relationship.
P. Römelt, P.R. Cunningham / Composite Structures 94 (2012) 977–986
Fig. 16. Tow and matrix loading at different strain levels.
985
in the unit cell. Because of the complexity and large variation of local strain states within the unit cell, tensile strain over the entire unit cell is used as a reference to describe the state of the unit cell model. Damage initiates between the weft tow and the matrix at the point of highest stress concentration (Fig. 14a). From there damage instantaneously propagates into the weft tow and the matrix (Fig. 14c). On one side damage spans the whole width of the unit cell while on the other side only half the width is spanned by the damage (Fig. 14c–e). This has been observed by Gao et al. [19] for a laminate of 8-harness satin weave plies under uni-axial tension. Consequently, the matrix between warp and weft is unloaded in the numerical model due to damage propagating and tow matrix interface failure (Fig. 14d). The second mode of failure in the numerical model is tow/matrix de-bonding (Fig. 14b). Again this is consistent with experimental findings in [9,19]. 5.2. Unit cell binary model results
Fig. 17. Stress–strain results of full finite element unit cell model and equivalent binary unit cell model.
showed that damage initiates and propagates at similar points in woven laminates under uni-axial tension using digital microscopy. Fig. 14a–e shows the scalar value of damage parameter d at different strain levels. This gives an idea of the propagation of damage
The stress–strain curve obtained from the full finite element unit cell model was used, together with the analytical approach presented in Section 4.2, to generate the nonlinear force–displacement curve for the one-dimensional tow elements as shown in Fig. 15. At the start of this curve the spring stiffness, given by the slope of the force–displacement curve, is degrading with increasing displacement as expected. However, after a displacement of 0.004 mm, which equates to a strain of 0.5%, the spring stiffness increases and then remains constant. In the numerical model damage initiates at the weft yarn/matrix interface and progresses towards the top and bottom of the unit cell. As damage grows to the boundary of the unit cell, the increasing load is carried by the tow only (Fig. 16). Since the tow was assumed to be linear-elastic, the resulting stress–strain curve of the full finite element model becomes linear, requiring a linear increase in spring stiffness. In addition, Poisson’s effect was ignored in the analytical approach used to calculate the nonlinear spring stiffness, which would also account for a deviation of predicted and experimental results for higher strains when Poisson’s effect becomes significant. Fig. 17 shows a comparison of the response of the full finite element model to the response of the binary model in the axial direction. Good correlation is achieved for the majority of the strain
Fig. 18. Stress–strain results of macro-scale binary model and experimental results.
986
P. Römelt, P.R. Cunningham / Composite Structures 94 (2012) 977–986
range with a deviation of less than 5% for strains up to 1.5%. For larger strains the deviation is up to 7% at a maximum strain of 2.5%. This can be explained by the unusual shape of the load–displacement curve that was used as an input for the one-dimensional tow elements, as previously discussed. 6. Macro-scale modelling Fig. 18 shows the comparison of the final macro-scale model with the results from the tensile test. Initially the stress–strain curve predicted by the macro-scale follows the averaged experimental stress–strain curve up until to an axial strain of about 0.015. For higher strain levels the predicted stress–strain curve deviates from the average by about 6% at the extreme, but stays within the experimental bounds. This difference between the numerical macro-scale and experimental results is due to two reasons. In the numerical models damage could not progress across the unit cell boundary whereas in the experimental investigation damage progresses throughout the laminate on a scale exceeding the dimensions of a unit cell [9,19]. Also, the mechanical behaviour of both constituents in the numerical model was assumed to be linear with linear degradation of selected matrix elements while E-glass fibre tows show strain dependent material properties [20] and epoxies behave visco-elastically [21,22]. Despite this, good agreement with experimental results is obtained using the method outlined in this paper. 7. Conclusion A numerical procedure to simulate the macro-behaviour of twodimensional woven laminates has been suggested in this paper. A binary model of the critical area of the macro-structure was generated and nonlinear springs were used to model the stiffness degradation due to damage. The force–displacement relationship for the nonlinear springs was derived from a meso-scale analysis of a single unit cell of the composite. Cohesive elements were used in the full finite element analysis to model damage initiation and progression. The predicted stress–strain curve shows good agreement with the averaged experimental results to a strain level of about 1.5%. For higher strains the predicted curve deviates from the averaged experimental results but stays within the bounds of the experimental results. Possible reasons for the differences in the result at higher strain values has been attributed to damage not being able to grow beyond the unit cell boundary in the numerical model and the nonlinear behaviour of the constituent materials not taken into account by the full finite element unit cell model. Errors introduced due to the latter can easily be mitigated by use of nonlinear material models, if precise data on constituent material behaviour is available. Mitigating errors due to the limitation of damage growth can only be achieved by increasing the size of the unit cell in the model.
References [1] Rueckert C, Kolax M. New design approaches for composite fuselage structures – requirements for advanced materials, processes and manufacturing technologies. In: Proceedings of 14th international conference on composites materials (ICCM-14), San Diego, California, USA, 14–18 July; 2003. [2] Tserpes KI, Labeas G, Pantelakis S. Multi-scale modelling of the mechanical response of plain weave composites and cellular solids. Theor Appl Fract Mech 2010. doi:10.1016/j.tafmech.2010.10.01. [3] Bogdanovich AE. Multi-scale modeling, stress and failure analyses of 3-D woven composites. J Mater Sci 2006;41:6547–90. [4] Cox BN, Yang Q. Predicting failure in textile composites using the binary model with gauge averaging. Eng Fract Mech 2010. doi:10.1016/ j.engfracmech.2010.08.00. [5] Lomov S, Ivanov DS, Verpoest I, Zako M, Kurashiki T, Nakai H, et al. Meso-FE modelling of textile composites: road map, data flow and algorithms. Compos Sci Technol 2007;67:1870–91. [6] Crookston JJ, Ruijter W, Long AC, Jones IA. Modelling mechanical performance including damage development for textile composites using a grid-based finite element method with adaptive mesh refinement. In: Proceedings of 8th international conference on textile composites (TEXCOMP-8), Nottingham, UK, 16–18 October; 2006. [7] Williams C, Summerscales J, Grove S. Resin infusion under flexible tooling (RIFT): a review. Composites Part A 1996;27:517–24. [8] Lang EJ, Chou T-W. The effect of strain gage size on measurement errors in textile composite materials. Compos Sci Technol 1998;58:539–48. [9] Daggumati S, De Baere I, Van Paepegem W, Degrieck J, Xu J, Lomov SV, et al. Local damage in a 5 – harness satin weave composite under static tension: part I – experimental analysis. Compos Sci Technol 2010. doi:10.1016/j.compscitech. 2010.07.00. [10] Matweb Material Property Data Generic Glass Fibre; January, 2011.. [11] Matweb Material Property Data Epoxy Resin; January, 2011. . [12] Keusch S, Queck H, Gliesche K. Influence of glass fibre/epoxy resin interface on static mechanical properties of unidirectional composites and on fatigue performance of cross ply composites. Composites Part A 1998;29:701–5. [13] Crookston JJ, Long AC, Jones IA. A summery review of mechanical properties prediction methods for textile reinforced polymer composites. Proc IMechE Part L: J Mater: Des Appl 2005;219:91–109. [14] Kim HJ, Swan CC. Voxel-based meshing and unit-cell analysis of textile composites. Int J Numer Methods Eng 2003;56:977–1006. [15] Glaessgen EH, Pastore CM, Hayden Griffin O, Birger A. Geometrical and finite element modelling of textile composites. Composites Part B 1996;27:43–50. [16] Ivanov DS, Lomov SV, Ivanov SG, Verpoest I. Stress distribution in outer and inner plies of textile laminates and novel boundary conditions for unit cell analysis. Composites Part A 2010;41:571–80. [17] Whitcomb JD, Chapman CD, Tang X. Derivation of boundary conditions for micromechanics analyses of plain and satin weave composites. J Compos Mater 2000;34:724–47. [18] Ridha M, Tan VBC, Tay TE. Traction–separation laws for progressive failure of bonded scarf repair of composite panel. Compos Struct 2011;93:1239–45. [19] Gao F, Boniface L, Ogin SL, Smith PA, Greaves RP. Damage accumulation in woven-fabric CFRP laminates under tensile loading: Part 1. Observations of damage accumulation. Compos Sci Technol 1999;59:123–36. [20] Bazhenov SL, Rogozinskii AK, Evstiforov SS, Berlin AA. Abnormal effect of tensile strain on Young’s modulus of inorganic E-glass fibres. Doklady Phys Chem 2010;433:118–20. [21] Hu Y, Xia Z, Ellyin F. Deformation behavior of en epoxy resin subject to multiaxial loadings. Part I: experimental investigations. Polym Eng Sci 2003;43:721–33. [22] Seidel A, Lepenies I, Engler T, Cherif C, Zastrau BW. Aspect of creep behavior of textile reinforcements for composites materials. Open Mater Sci J 2009;3:67–79.
Contents lists available at SciVerse ScienceDirect
Composite Structures journal homepage: www.elsevier.com/locate/compstruct
A multi-scale finite element approach for modelling damage progression in woven composite structures P. Römelt, P.R. Cunningham ⇑ Department of Aeronautical and Automotive Engineering, Loughborough University, Leicestershire LE11 3TU, United Kingdom
a r t i c l e
i n f o
Article history: Available online 6 November 2011 Keywords: Woven composites Multi-scale modelling Damage
a b s t r a c t A significant challenge in the numerical modelling of composite structures with a multi-axis fibre architecture is the reproducibility of the textile mechanics [1]. A numerical analysis procedure for woven composite structures using a multi-scale finite element approach has been developed, and is presented in this paper. The approach is demonstrated for a flat two-dimensional woven glass/epoxy laminate. Digital microscopy is used to estimate tow cross-section and path, and quantify the amount of variation of these parameters. This data is used to generate both a meso-scale model of a single unit cell as well as a macroscale model of the complete structure. Numerical results from the proposed approach are compared to experimental stress–strain data, which show good agreement in the lower strain range. Ó 2011 Elsevier Ltd. All rights reserved.
1. Introduction In order to be used efficiently, a design procedure for twodimensional woven composites needs to include a validated damage model. To correctly model damage initiation and propagation, models which take into account the different modes of failure, such as tow matrix de-bonding and crack propagation within the tows and matrix are required. Since these failure events can be best observed on the meso-scale, a multi-scale approach which uses meso-scale models to generate material input parameters for macro-cells would be the most suitable to model the behaviour of a woven macrostructure. Three different approaches to model multi-axial structures can be found in the literature, the full finite element approach [2], the voxel technique [3] and the binary model [4]. The full finite element approach [2] uses precise geometrical data to generate a high resolution model of a unit cell. This allows for an estimation of the stress distribution within the unit cell, which in turn allows for realistic modelling of damage initiation and progression within the cell’s limit. Lomov et al. have shown that this method has been used successfully to predict damage initiation [5]. However, models of a single unit cell require a high number of elements, resulting in high computational costs, which makes this method unsuitable for modelling structures larger than a few unit cells. Bogdanovich defines voxels as sub-meso-volumes [3] and Crookston et al. [6] have developed a modelling approach based on such voxels. A two-dimensional grid is laid over the top of a geometrical ⇑ Corresponding author. E-mail address: [email protected] (P.R. Cunningham). 0263-8223/$ - see front matter Ó 2011 Elsevier Ltd. All rights reserved. doi:10.1016/j.compstruct.2011.10.024
unit cell model and the architecture is analysed in the throughthickness direction. Average material properties are calculated depending on the through-thickness architecture, which are then assigned to three-dimensional elements generated at the twodimensional grid point positions. The method reduces both meshing effort and the size of the finite element model considerably, resulting in reduced computational effort [6]. However, by averaging material properties the stress distribution is also averaged over the voxel. This means that stress concentrations, especially on the sites of damage initiation like the tow/matrix interface, are averaged out. Crookston et al. successfully applied this method to reproduce the low strain linear behaviour of a plain weave reinforced glass/ polyester composite but could not reproduce the nonlinear behaviour for higher strain values [6]. The binary model [4] established by Cox and Yang utilises onedimensional elements, which run along the yarn centre line. These elements represent the axial properties of the yarns. The other yarn properties are incorporated in a three-dimensional representative matrix element, referred to as ‘‘effective medium elements’’, using micromechanics. As with the voxel technique, the binary model reduces meshing effort and computational cost while averaging out stress concentrations. Cox and Yang successfully used the binary model to predict ultimate stress of tri-axially braided carbon/ epoxy composites and reproduce the linear part of the stress– strain curve but could not reproduce the nonlinear part [4]. The aim of this paper is to demonstrate a multi-scale procedure using digital microscopy and finite elements to model damage in woven composite macro-structures. Making use of the advantages of the averaging techniques and the full finite element model, this procedure attempts to predict the stiffness degradation of macrostructures due to meso-scale damage. The approach taken is
978
P. Römelt, P.R. Cunningham / Composite Structures 94 (2012) 977–986
Nomenclature A a, b C d D, E Et KIc K L t t0 t x, y, z
r
cross-sectional area (mm2) major, minor elliptical half axes (mm) sine function amplitude (mm) damage parameter (unity) sine function parameters (mm) tangent modulus (GPa) fracture toughness (MPa (m)0.5) stiffness (N/mm) sinusoidal wavelength (mm) traction (MPa) allowable traction (MPa) traction in undamaged material (MPa) coordinate values tensile stress (MPa)
illustrated in Fig. 1. Digital Microscopy was used to find the required geometric input parameters for a meso-scale model of a unit cell of a woven macrostructure using a full finite element approach. The results from the full finite element model were then used to define the material properties for an equivalent binary model of a unit cell, which shows the same behaviour as the full finite element model. The macrostructure was then modelled using these equivalent binary unit cells and stress–strain predictions were made for a simple unidirectional tensile load case. The results from a series of tensile tests of glass–epoxy specimens were used to validate the finite element predictions.
s
shear stress (MPa)
Vectors t e
tangent vector directional vector
Subscripts 1D one dimensional 3D three dimensional n normal direction s, t transverse shear directions spring nonlinear spring element x, y, z coordinate directions
2. Tensile testing and microscopic analysis 2.1. Test specimen manufacture A 2.5 mm thick laminate plate measuring 330 mm 330 mm was manufactured from thirty plies of plain weave glass fabric (Gurit RE86P) with a layup of [0/90]30, which was infused with Gurit Prime 20-LV epoxy resin using a double bag vacuum infusion procedure to afford better consolidation and even thickness of the cured laminate [7]. Nine test specimens were cut from this laminate with the geometry and dimensions as shown in Fig. 2. A rectangular rosette strain gauge with a 5 mm gauge length and a resistance of 120 X was bonded to one side of each of the test specimens in the position shown in Fig. 2. 2.2. Tensile testing Tensile tests were conducted according to ASTM D3039. Although the use of strain gauges for mechanical testing of woven composites has been criticised in the literature [3], Lang and Chou [8] have demonstrated that good results can be achieved if the gauge length is longer than a single unit cell length. The resulting stress–strain curves from the tensile tests are shown in Fig. 3, with the mean stress–strain curve also marked. There is an average variation of about 7% about the mean in the resulting stress–strain data, with failure strengths ranging from 375 MPa to 450 MPa with a +9.8% to 8.5% variation at failure. This is consistent with findings from Daggumati et al. [9], who showed that damage initiation and progressive stiffness degradation is highly dependent on the amount of variation in tow path in both warp and weft direction and tow cross sectional area. 2.3. Microscopic analysis A total of 123 microscopic specimens were cut using the same laminate from which the tensile test specimens were obtained. A
Fig. 1. Flow chart of numerical procedure.
Fig. 2. Tensile test specimen, geometry and strain gauge position.
P. Römelt, P.R. Cunningham / Composite Structures 94 (2012) 977–986
979
Fig. 3. Tensile test results.
Fig. 4. Laminate cutting scheme.
Fig. 5. Magnified images from microscopic analysis: (a) sine fit and (b) ellipse fit.
typical cutting pattern for a portion of the laminate can be seen in Fig. 4. After cutting the specimens were cast in a dyed epoxy resin to provide contrast under the microscope, and the specimens were polished using a five-step process. Tow path data was generated using 20 times magnified microscopic images as shown in Fig. 5a. This was done to ensure that the averaged tow path and cross-section parameters were the same for the numerical models as for the experimental specimen. Lines were placed along the centre of the tow in the 20 times magnified
microscope image and the coordinates of the start and end points of these lines were recorded. A sine function was fitted through those points using a least-square fit algorithm as given in the following equation:
zðxÞ ¼ D sin where
2px 2px þ E cos L L
ð1Þ
980
P. Römelt, P.R. Cunningham / Composite Structures 94 (2012) 977–986
Fig. 6. Comparison of least-square fit results with microscopic data.
Table 1 Geometrical properties of two-dimensional woven laminate meso-structure. Parameter (unit)
Average value ± standard deviation
C (mm) L (mm) a (mm) b (mm) A (mm2)
0.025 ± 0.009 1.604 ± 0.164 0.187 ± 0.066 0.037 ± 0.007 0.022
Table 3 Published E-glass fibre properties [10]. E (GPa) G (GPa) m (–)
72.4 30.0 0.2
Table 4 Gurit Prime 20LV/epoxy properties [11].
P P 2 ni¼1 zi sin 2pLxi E ni¼1 sin 4pLxi D¼ P n ni¼1 cos 4pLxi
E (GPa) G (GPa) m (–) Ultimate tensile stress, rult (MPa) Ultimate shear stress, sult (MPa) KIc (MPa (m)0.5)
ð2Þ
3.5 1.3 0.35 73 137 3.69
and
P P 2 ni¼1 zi cos 2pLxi D ni¼1 sin 4pLxi E¼ P n ni¼1 cos 4pLxi
ð3Þ
Tow cross-section data was generated using 80 times magnified images of the same specimen as shown in Fig. 5b. Lines were placed along the tow cross-section in the 80 times magnified microscopic image and a standard ellipse was fitted through the start and end points using a least-square fit algorithm as given in Eq. (4). A total of 385 tow paths were analysed as well as 449 tow cross-sections.
rffiffiffiffiffiffiffiffiffiffiffiffiffiffi y2 a2
zðyÞ ¼ b 1
ð4Þ
where
a2 ¼
Pn
4 i¼1 yi
Pn
2 i¼1 yi
Pn i¼1
Pn
4 i¼1 zi
2
b ¼
ð5Þ
ðy2i z2i Þ 2
b
and
Pn
2 i¼1 zi
Pn i¼1
ð6Þ
ðy2i z2i Þ
a2
Fig. 6 shows a comparison of the measured points from an example specimen of the microscopic analysis and the curves fitted through them using the least square algorithms given by Eqs. (1) and (4). Fig. 6a shows a good correlation with exceptions at the extremes, especially at the beginning of the measurements. The correlation for the ellipse is reasonable as shown in Fig. 6b but gives a larger result for both half-axes of the ellipse compared with the measured data. The microscopic analysis showed variations between 10% and 100% for geometric parameters. Daggumati et al. [9] showed that these variations lead to variations in the point of damage initiation in a unit cell. For simplicity reasons only average values are used for geometric parameters in this paper. The average values as well as standard deviations for tow path amplitude C (mm) and wavelength L (mm) and the cross-section
Table 2 Weave properties. Warp
Weft
Filament type
Weave type
Nominal (g/m2)
Fibre (tex)
Count (ends/cm)
Theoretical (g/m2)
Count (ends/cm)
Theoretical (g/m2)
E-glass
Plain
85
34
12.0
41.5
12.5
43.3
P. Römelt, P.R. Cunningham / Composite Structures 94 (2012) 977–986
981
ellipse parameters a (mm) and b (mm) are listed in Table 1. The cross-sectional area A (mm2) for the tow ellipse was calculated using the averaged ellipse parameters a and b from the fit. As can be seen from the results, the standard deviation of the fit for the samples measured is relatively low, with a maximum standard deviation of 10% occurring for the tow path wavelength. 3. Material properties of constituents A number of parameters of the constituent materials, such as Young’s modulus and Poisson’s ratio, were required for the modelling approach. However, data on the constituents is very limited. Therefore a number of parameters had to be assumed with reference values found in the literature. Isotropic linear elastic material behaviour was assumed for both the matrix and the glass fibre tows. The two-dimensional glass fibre weave used in the laminate was Gurit RE86P E-glass plain weave, the properties for which are listed in Table 2. The material properties of the E-glass fibre, as shown in Table 3, were taken from the MATWEB database [10]. The resin used was Gurit Prime 20LV, and the material properties are listed in Table 4. Not all of the required resin properties were available from the manufacturer, therefore certain properties, such as fracture toughness, were also taken from the MATWEB database for epoxy materials [11]. Finally, assumptions had to be made for the interface strength between tows and matrix, on which numerous studies were
Fig. 9. Traction separation formulation.
available in the literature. The value of 62.5 MPa from Ref. [12] was chosen for this analysis. 4. Numerical modelling 4.1. Full finite element unit cell model A unit cell of the two-dimensional weave pattern, with geometrical parameters obtained from the microscopic analysis, was modelled using a full finite element approach [2]. PATRAN was used as
Fig. 7. Full finite element unit cell model and equivalent binary unit cell model.
Fig. 8. Unit cell boundary conditions.
982
P. Römelt, P.R. Cunningham / Composite Structures 94 (2012) 977–986
Fig. 10. Graphical representation of analytical approach.
Fig. 11. Eight-block equivalent binary unit cell model.
a pre-processor to model geometry and for meshing. An input file was generated and amended to include material and section properties, boundary conditions, solver parameter and requested output data. ABAQUS/Standard ver. 6–10 was used as the solver, with ABAQUS/CAE for post-processing. The major challenge when generating full finite element unit cell models of multi-axial composites is to achieve an accurate description of the unit cell while keeping the computational resources required to solve the model
to a minimum [13]. Due to the highly elliptical cross-section of the tow and the sinusoidal tow path, coarse meshes tend to have high aspect ratios and skew angles. This can lead to significant errors in predicting the local stress–strain field [14]. Mesh refinement is usually used to mitigate this problem, which results in higher computational costs. However, because of the complex geometry of the model the mesh would have to be refined to a level that computational costs become unreasonably high. Therefore, some distorted elements in the mesh are deemed to be acceptable [15]. Eight-node brick elements were used to generate a mesh with a total of 117,096 degrees of freedom. For this mesh, 3% of the elements had an aspect ratio greater than the recommended value of 5 or skew angles greater than the recommended value of 30°, which was deemed an acceptable level for a model of this type [15]. Glaessgen et al. [15] achieved convergence for a plain weave glass fibre epoxy unit cell model using 50,000 degrees of freedom. The mesh used is shown in Fig. 7a, matrix elements have been omitted for an easier portrayal of the mesh. It has been demonstrated by Ivanov et al. [16] that boundary conditions on the meso-scale model depend on the position of the unit cell in the structure. However periodic boundary conditions, derived by Whitcomb et al. [17], were assumed in this paper, again for simplicity reasons. A constant displacement was applied to all nodes on the y–z faces of the unit cell in the axial direction. The faces in the transverse directions of the unit cell, x–y and x–z faces, were kept straight to account for the effect of neighbouring cells in the weave while allowing for Poisson’s effect. This also constrains rigid body rotations around all axes. The centre node of the model was constrained in both y and z direction to prevent rigid body motion in these directions. The boundary conditions are shown diagrammatically in Fig. 8. Global strains of the unit cell were calculated by dividing the prescribed displacement by the length of the unit cell. Global stresses were calculated by the summation of axial forces of every node on the left y–z face, which was then divided by the cross sectional area of the same face. Cohesive elements with a traction–separation formulation (Fig. 9) similar to Ridha et al. [18] were used to model crack
Fig. 12. Critical cross section of tensile test specimen.
Fig. 13. Linear analysis stress distribution at 2.5% strain.
P. Römelt, P.R. Cunningham / Composite Structures 94 (2012) 977–986
Fig. 14. Damage parameter at various global tensile strain values: (a) 0.37% strain, (b) 0.75% strain, (c) 0.75% strain, (d) 1.66% strain and (e) 2.5% strain.
983
984
P. Römelt, P.R. Cunningham / Composite Structures 94 (2012) 977–986
initiation, crack propagation and tow matrix de-bonding. A maximum traction criterion, as shown in Eq. (7), was used in conjunction with a damage variable d, which reduced the resulting traction t for the undamaged material at a given strain level as shown in Eq. (8). In other words, ABAQUS/Standard first calculates the tractions for the undamaged material at a given strain level and then multiplies these tractions with the damage parameter d. The damage parameter d is calculated based on the combined effect of normal and shear tractions within the element. It reflects overall damage within the cell and is therefore used to reduce all local element stiffnesses, independent of the cohesive element orientation.
( ) jt n j ts t t ¼1 max ; ; t 0n t0s t 0t tn ¼
ð1 dÞt n ; tn > 0 tn
ð7Þ ð8Þ
where
t s ¼ ð1 dÞt s
ð9Þ
and
t t ¼ ð1 dÞt t
ð10Þ
4.2. Binary unit cell model In the next step a similar model of the weave was made using the binary model approach (Fig. 7b) [4]. The binary model used three-dimensional hex elements to model transverse tow and surrounding matrix properties. Axial stiffness of the tows was modelled using one-dimensional spring elements, which were embedded in the three-dimensional element mesh. The properties of elements, i.e. spring stiffness for one-dimensional elements and linear isotropic material properties for three-dimensional elements, were calculated using the tow geometrical properties measured during the microscopic analysis and the rule-of mixtures for an equivalent unidirectional laminate with a fibre volume fraction comparable to the multi-axial composite. Unlike the initially proposed model by Cox and Yang [4], nonlinear one-dimensional spring elements were used to model damage in this work. This means all damage modes observed in the full finite element unit cell model were accounted for in the load– displacement relationship of the one-dimensional elements. The required stiffness at different strain levels of the unit cell can be estimated using the tangent modulus Et derived from the resultant true stress–true strain curve of the full finite element model. The stiffness of the one-dimensional element at different strain levels was calculated using the results from the full finite element model. For this approach one-dimensional elements and effective medium elements were modelled as a series of springs, as shown in Fig. 10. Four elements were assumed to be sufficient to reproduce the tow path and therefore model the axial stiffness of the tow. These four elements (K1D) were assumed to be in series with each other (Kspring) and in parallel with the three-dimensional solid element (K3D). The cross-sectional area and length of the equivalent binary unit cell is denoted by A and L respectively. The local orientation of the one-dimensional spring elements, which is different to the axial direction of the unit cell due to tow crimp, had to be considered when calculating the required stiffness of the individual one-dimensional element since the resultant stress–strain curve of the full finite element unit cell model was derived for the axial direction of the unit cell. This was done using the scalar product of the normalised local tow path tangent vector t with the global axial unit vector ex as shown in the following equations:
K 1D ¼ K spring
4 X ti ex jt ij i¼1
K spring ¼ ðEt E3D Þ
A L
ð11Þ
ð12Þ
In order to check the force–displacement relationship generated by the analytical approach a simple binary model was generated. However, assigning the boundary conditions that have been used in the full finite element model to a single equivalent unit cell was not possible because only a single three-dimensional solid element was used. Therefore a block of eight equivalent binary unit cells was modelled and the same boundary conditions applied. This block of equivalent binary unit cells has a total of 594 degrees of freedom and is shown in Fig. 11. 4.3. Binary macro-scale model In the final step a binary model of the critical cross-section of the test specimen was modelled, as shown in Fig. 12. Nine equivalent binary unit cells are placed next to each other in the widthwise and axial direction respectively and in order to achieve the correct thickness, thirty of these plies are stacked giving the same number of plies used in the tensile test specimen. This results in a total of 540,186 degrees of freedom in the model. Boundary conditions were chosen as to best represent the conditions at the critical cross-section. A uniform displacement was applied to the y–z faces of the model. This condition not only constrains axial movement but also rotation around the z and y axes. The centre node of the model was fixed in both the y and z direction to constrain rigid body motion in both transverse directions. Finally, to constrain rotation around the x axis the two nodes in the centre of the x–z faces were constrained in the z direction. 5. Results and discussion 5.1. Full finite element model A linear analysis without damage was conducted first to find the stress concentration points in the unit cell. To account for damage, cohesive elements were placed at the point of highest stress in the full finite element unit cell model using the results from the linear analysis. The elements were placed along the line of highest stress and around the interface between tow and matrix. Fig. 13 shows stress in the axial direction in the linear analysis, with stress concentrations near the edge of the weft tows and at the tow/matrix interface. Daggumati et al. [9] and Gao et al. [19]
Fig. 15. One-dimensional spring load–displacement relationship.
P. Römelt, P.R. Cunningham / Composite Structures 94 (2012) 977–986
Fig. 16. Tow and matrix loading at different strain levels.
985
in the unit cell. Because of the complexity and large variation of local strain states within the unit cell, tensile strain over the entire unit cell is used as a reference to describe the state of the unit cell model. Damage initiates between the weft tow and the matrix at the point of highest stress concentration (Fig. 14a). From there damage instantaneously propagates into the weft tow and the matrix (Fig. 14c). On one side damage spans the whole width of the unit cell while on the other side only half the width is spanned by the damage (Fig. 14c–e). This has been observed by Gao et al. [19] for a laminate of 8-harness satin weave plies under uni-axial tension. Consequently, the matrix between warp and weft is unloaded in the numerical model due to damage propagating and tow matrix interface failure (Fig. 14d). The second mode of failure in the numerical model is tow/matrix de-bonding (Fig. 14b). Again this is consistent with experimental findings in [9,19]. 5.2. Unit cell binary model results
Fig. 17. Stress–strain results of full finite element unit cell model and equivalent binary unit cell model.
showed that damage initiates and propagates at similar points in woven laminates under uni-axial tension using digital microscopy. Fig. 14a–e shows the scalar value of damage parameter d at different strain levels. This gives an idea of the propagation of damage
The stress–strain curve obtained from the full finite element unit cell model was used, together with the analytical approach presented in Section 4.2, to generate the nonlinear force–displacement curve for the one-dimensional tow elements as shown in Fig. 15. At the start of this curve the spring stiffness, given by the slope of the force–displacement curve, is degrading with increasing displacement as expected. However, after a displacement of 0.004 mm, which equates to a strain of 0.5%, the spring stiffness increases and then remains constant. In the numerical model damage initiates at the weft yarn/matrix interface and progresses towards the top and bottom of the unit cell. As damage grows to the boundary of the unit cell, the increasing load is carried by the tow only (Fig. 16). Since the tow was assumed to be linear-elastic, the resulting stress–strain curve of the full finite element model becomes linear, requiring a linear increase in spring stiffness. In addition, Poisson’s effect was ignored in the analytical approach used to calculate the nonlinear spring stiffness, which would also account for a deviation of predicted and experimental results for higher strains when Poisson’s effect becomes significant. Fig. 17 shows a comparison of the response of the full finite element model to the response of the binary model in the axial direction. Good correlation is achieved for the majority of the strain
Fig. 18. Stress–strain results of macro-scale binary model and experimental results.
986
P. Römelt, P.R. Cunningham / Composite Structures 94 (2012) 977–986
range with a deviation of less than 5% for strains up to 1.5%. For larger strains the deviation is up to 7% at a maximum strain of 2.5%. This can be explained by the unusual shape of the load–displacement curve that was used as an input for the one-dimensional tow elements, as previously discussed. 6. Macro-scale modelling Fig. 18 shows the comparison of the final macro-scale model with the results from the tensile test. Initially the stress–strain curve predicted by the macro-scale follows the averaged experimental stress–strain curve up until to an axial strain of about 0.015. For higher strain levels the predicted stress–strain curve deviates from the average by about 6% at the extreme, but stays within the experimental bounds. This difference between the numerical macro-scale and experimental results is due to two reasons. In the numerical models damage could not progress across the unit cell boundary whereas in the experimental investigation damage progresses throughout the laminate on a scale exceeding the dimensions of a unit cell [9,19]. Also, the mechanical behaviour of both constituents in the numerical model was assumed to be linear with linear degradation of selected matrix elements while E-glass fibre tows show strain dependent material properties [20] and epoxies behave visco-elastically [21,22]. Despite this, good agreement with experimental results is obtained using the method outlined in this paper. 7. Conclusion A numerical procedure to simulate the macro-behaviour of twodimensional woven laminates has been suggested in this paper. A binary model of the critical area of the macro-structure was generated and nonlinear springs were used to model the stiffness degradation due to damage. The force–displacement relationship for the nonlinear springs was derived from a meso-scale analysis of a single unit cell of the composite. Cohesive elements were used in the full finite element analysis to model damage initiation and progression. The predicted stress–strain curve shows good agreement with the averaged experimental results to a strain level of about 1.5%. For higher strains the predicted curve deviates from the averaged experimental results but stays within the bounds of the experimental results. Possible reasons for the differences in the result at higher strain values has been attributed to damage not being able to grow beyond the unit cell boundary in the numerical model and the nonlinear behaviour of the constituent materials not taken into account by the full finite element unit cell model. Errors introduced due to the latter can easily be mitigated by use of nonlinear material models, if precise data on constituent material behaviour is available. Mitigating errors due to the limitation of damage growth can only be achieved by increasing the size of the unit cell in the model.
References [1] Rueckert C, Kolax M. New design approaches for composite fuselage structures – requirements for advanced materials, processes and manufacturing technologies. In: Proceedings of 14th international conference on composites materials (ICCM-14), San Diego, California, USA, 14–18 July; 2003. [2] Tserpes KI, Labeas G, Pantelakis S. Multi-scale modelling of the mechanical response of plain weave composites and cellular solids. Theor Appl Fract Mech 2010. doi:10.1016/j.tafmech.2010.10.01. [3] Bogdanovich AE. Multi-scale modeling, stress and failure analyses of 3-D woven composites. J Mater Sci 2006;41:6547–90. [4] Cox BN, Yang Q. Predicting failure in textile composites using the binary model with gauge averaging. Eng Fract Mech 2010. doi:10.1016/ j.engfracmech.2010.08.00. [5] Lomov S, Ivanov DS, Verpoest I, Zako M, Kurashiki T, Nakai H, et al. Meso-FE modelling of textile composites: road map, data flow and algorithms. Compos Sci Technol 2007;67:1870–91. [6] Crookston JJ, Ruijter W, Long AC, Jones IA. Modelling mechanical performance including damage development for textile composites using a grid-based finite element method with adaptive mesh refinement. In: Proceedings of 8th international conference on textile composites (TEXCOMP-8), Nottingham, UK, 16–18 October; 2006. [7] Williams C, Summerscales J, Grove S. Resin infusion under flexible tooling (RIFT): a review. Composites Part A 1996;27:517–24. [8] Lang EJ, Chou T-W. The effect of strain gage size on measurement errors in textile composite materials. Compos Sci Technol 1998;58:539–48. [9] Daggumati S, De Baere I, Van Paepegem W, Degrieck J, Xu J, Lomov SV, et al. Local damage in a 5 – harness satin weave composite under static tension: part I – experimental analysis. Compos Sci Technol 2010. doi:10.1016/j.compscitech. 2010.07.00. [10] Matweb Material Property Data Generic Glass Fibre; January, 2011.