Composite Structures 95 (2013) 53–62 Contents lists available at SciVerse ScienceDirect Composite Structures journal homepage: www.elsevier.com/locate/compstruct Progressive failure analysis of thin-walled composite structures Diego Cárdenas a, Hugo Elizalde b, Piergiovanni Marzocca c,⇑, Frank Abdi d, Levon Minnetyan e, Oliver Probst a a Physics Department, Instituto Tecnológico y de Estudios Superiores de Monterrey, Campus Monterrey, Eugenio Garza Sada 2501 Sur, Monterrey, N.L., Mexico CP 64849, Mexico School of Engineering, Instituto Tecnológico y de Estudios Superiores de Monterrey, Campus Ciudad de México, Eugenio Garza Sada 2501 Sur, Monterrey, N.L., Mexico CP 64849, Mexico c Mechanical and Aeronautical Engineering Department, Clarkson University, Potsdam, NY, USA d Alpha Star Corporation, 5150 East Pacific Coast Highway Ste. 650, Long Beach, CA, USA e Civil and Environmental Engineering Department, Clarkson University, Potsdam, NY, USA b a r t i c l e i n f o Article history: Available online 27 May 2012 Keywords: Progressive Failure Analysis (PFA) Damage model Failure criteria Thin-Wall Beam (TWB) Finite Element Model (FEM) Wind turbine blade a b s t r a c t A reduced-order finite-element model suitable for Progressive Failure Analysis (PFA) of composite structures under dynamic aeroelastic conditions based on a Thin-Walled Beam (TWB) formulation is presented. Validation of the PFA-TWB against an integrated PFA model based on a shell formulation and implemented in the commercial software tool GENOA is conducted for static load conditions. A helicopter blade made from composite material and previously used in literature for the discussion of damage propagation is used as the reference case. The failure criteria for the different layers of the composite material used in the PFA-TWB model have been formulated in analogy with the corresponding criteria implemented in the shell formulation. Comparisons between the predictions of both models for progressively increasing load have been conducted in terms of the cumulative overall damage volume in the thinwalled structure, the layer-resolved cumulative damage volume, as well as through spatially resolved damage maps for both models. A strikingly similar damage topology has been found from both models up to load values close to final failure, in spite of the restraining assumptions of the TWB formulation. In terms of damage volume the PFA-TWB models predicts slightly higher values which can be traced back to the inevitable differences in the failure criteria formulation in the one-dimensional and the shell model, respectively. It is shown that a good agreement with the predictions of the shell model in terms of the cumulative damage volume is obtained if the strength values of the composite material are adjusted upwards in a uniform manner by about 10%. Considering the common safety factors usually applied in the design process of composite material the agreement of the TWB and the shell model in terms of damage propagation is considered excellent, allowing for the PFA-TWB to be used in systematic design studies. Ó 2012 Elsevier Ltd. All rights reserved. 1. Introduction Thin-walled composite beams (TWBs hereafter) are widely applied in many fields of structural engineering. Helicopter and wind turbine blades are examples of flexible slender structures that can be modeled as TWB. The most critical parts of helicopters and wind turbines are the rotors, which provide thrust and lift (in the case of the helicopter) or allow extracting power from the wind. Rotor blade design is a complex process usually involving several disciplines such as aerodynamics and structural analysis, acoustics, Abbreviations: 1D, 2D, 3D, one-, two-, and three-dimensional, respectively; CFRP, composite fiber-reinforced polymer; DOF, degrees of freedom; FE, FEM, finite-element model; PFA, Progressive Failure Analysis; TWB, Thin-Walled Beam. ⇑ Corresponding author. E-mail address: pmarzocc@clarkson.edu (P. Marzocca). 0263-8223/$ - see front matter Ó 2012 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.compstruct.2012.05.022 and dynamics [1,2]. Recent investigations [3–11] describe the current state-of-the-art in aeroelastic analysis, the success of which highly depends on the structural model’s ability to reproduce static and dynamic behavior under complex aerodynamic flow conditions. To this end, Finite Element (FE) and modal models are the two main approaches. Due to the very extensive computational resources required in aeroelastic and rotor-dynamics coupled simulations, finite-element representations are limited in practice to 1D linear or non-linear beam models, with the potential of reducing the size of the analyzed system down to a few dozens of DOF, albeit at the expense of accuracy and loss of detail due to the simplification of geometry and material layup [12]. On the other hand, modal approaches rely on linear mode superposition in order to represent the overall structural behavior, and their accuracy depends on generating sufficiently refined mode shapes for different rotor speeds [14]. To accurately evaluate the structural integrity of a rotor blade 54 D. Cárdenas et al. / Composite Structures 95 (2013) 53–62 the practitioner should perform an assessment by aeroelastic simulations, in combination with failure models obtained experimentally or theoretically [3,15–19]. However, it is currently unfeasible to use this approach to account for the interaction between the damage progression on the blade and the aeroelastic response. It is well known that the presence of damage in a composite structure changes the stiffness of the structure [20–29], requiring a continuous update of the structural properties of the model as damage progresses. Current approaches for modeling composite blade dynamics in aeroelastic codes (based on classical 1D FE-beam and modal models) use condensed mechanical and geometrical properties; it is therefore of a paramount difficulty to track damage progression on a layerby-layer basis. Thin-Walled Beam (TWB) models, on the other hand, provide an effective one-dimensional representation of the dynamics of a composite structure, while allowing recovering the layerresolved strain and stress fields to evaluate structural damage progression. There is little work documented in literature on the combination of beam FE models, recovery of shell strains/stresses and damage tracking. Li [3] in her doctoral thesis uses the VABS methodology, developed in [13], in combination with a damage model to optimize the structural properties of a helicopter rotor blade. She first calculates the strain/stress distributions in the blade for a typical helicopter manoeuvre and then performs an offline fatigue damage analysis for the complete aeroelastic time series generated from a multi-body dynamic software package. In the case a non-acceptable damage occurs the design is revised and newly subjected to a failure analysis; otherwise the design is deemed to be appropriate. To the best of the authors’ knowledge, while this approach allows for blade optimization in complex aeroelastic situations it does not provide a platform for integrated damage progression analysis, where the structural properties of the blade are continuously updated while damage progresses. Pawar and Ganguli [32] described an integrated damage model obtained by combining a matrix-cracking damage model developed by Gudmundson et al. [34,35], ply de-bonding/de-lamination and fiber-breakage models developed by Shahid and Chang [36], and a TWB model developed by Chandra and Chopra [37]. In their work the authors obtained relationships between the blade response and damage densities. A static load is applied to the rotor blade and the structural response of the blade as a function of the damage level is calculated. In a follow-up paper [33] Pawar and Ganguli study the effect of the damage level on the natural frequencies of a rotating thin-walled beam. In both papers the stated objective of the authors lies with structural health monitoring, requiring a relationship between the global damage level and the static and dynamic response of a blade. The authors do not provide, however, a model capable of predicting damage progression as a function of applied load or spatial distributions of the damage in response to a given load. To fill this gap in knowledge and to provide the author’s view on how an integrated damage progression analysis can be performed, in the present work a Progressive Failure Analysis (PFA) Thin-Wall Beam (TWB) Finite Element (FE) model designed to predict the progression of damage in slender complex structures, like rotor blades, is described. As opposed to the work reviewed above, this model is capable of continuously updating the structural properties of the blade as damage progresses, thereby providing an integrated description of failure propagation. The selected TWB model was originally developed by Librescu and co-workers [38–41] and subsequently discretized via FE by Vo and Lee [42–44]; an application of the Librescu–Vo–Lee model to a realistic turbine blade was described in [45]. The TWB model allows including material anisotropy, arbitrary laminate layups and shear deformability and has shown to reproduce the fundamental structural behavior of 3D shell models with reasonable accuracy [45]. Stress/strain fields for individual layers can be recovered based solely on the knowledge of the nodal displacements obtained from the effective 1D finite-element model, thereby allowing for Progressive Failure Analysis (PFA) at any layer and position of the structure. A further advantage is the availability of analytical expressions for the stiffness of arbitrary cross-sections, making damage tracking a more straightforward task during a displacement-based FE analysis. For purposes of validation of the PFA-TWB model described here, the failure criteria proposed by Chamis and Minnetyan [46] were used, since these criteria are implemented in the GENOA package, a commercial software tool designed for composite failure analysis in complex structures [47]. As a case study the methodology presented was applied to the helicopter blade described by Pawar and Ganguli [32,33] for the case of static loads. Application of the model to damage progression under dynamic aeroelastic conditions will be reported in a follow-up paper. 2. Methodology The integrated TWB-PFA model described in the preceding section was developed and implemented for the purposes of real-time dynamic failure analysis under realistic load conditions; however, the objective of the present paper is to present a validation of the model against a static model implemented in the GENOA package. In this section the integrated approach is discussed in detail and an application to a reference case is provided. 2.1. Integrated Thin-Walled Beam Progressive Failure Analysis (TWB-PFA) model The Thin-Walled Beam (TWB) model used in the present work is based on the formulations given in Refs. [38–44]; an application of the TWB formulation to a realistic wind turbine blade was has been described in [45]. The main aspects of the TWB model will be briefly reviewed below. The TWB model is an effective one-dimensional representation of beam-type structures such as helicopter or wind turbine blades, capable of recovering full 3D strain and stress information based on the knowledge of the nodal displacements of the 1D model alone. Evidently, a certain set of conditions has to be obeyed in order to allow such a reconstruction of information. These conditions are the following [43]: (1) The structure is restricted to small strain values, (2) the beam cross sections remain undeformed in their own plane for all load conditions applied, (3) both the transverse shear strains c0xz ; c0yz and warping shear c0x are uniform over the cross section, (4) the Kirchhoff–Love assumption in classical plate theory remains valid for laminated composites. Based on these assumptions, a double integration of geometric and materials properties of the cross-section (through the wall and along the contour) can be conducted, yielding semi-analytical expressions for the stiffness tensor (Eij) of arbitrary cross-sections. Since the coefficients Eij contain detailed information about the stiffness of any material point in the shell (as specified by its angular and thickness location), the full deformation and stress fields can be recovered for post-processing purposes, such as for the application of failure criteria for composites. If a failure criterion is met at a given location in the structure, the corresponding elastic moduli at this location may be degraded according to the failure model chosen. The stiffness tensor can then be updated and a new set of nodal displacements may be calculated. The traditional approach is to conduct the elastic calculations and the failure analysis separately [3,15–19]; then a progressive failure analysis often becomes a tedious process, even in the case of static loads. PFA for dynamic loads, on the other hand, is practically impossible with this approach. The use of reduced-order 55 D. Cárdenas et al. / Composite Structures 95 (2013) 53–62 model such as TWB, on the other hand, allows implementing a coupled procedure potentially capable of analyzing the progression of damage under dynamic load conditions in real-time. Reduced-order models necessarily come at the expense of accuracy, so it is important to evaluate the validity of the assumption of the model (1)–(4) under the load conditions studied and identify a range of applicability of the model. One of the key assumptions which may be violated under severe load conditions is the hypothesis that the contours of the beam sections remain constant (assumption 2 of the TWB); the validity of this assumption will be further examinated in the results section. A Timoshenko formulation is used with seven independent variables given by the node translations U(z), V(z), W(z) where z is longitudinal or beam axis, the angular displacement /(z) of the cross sections around the z axis, and the angular displacements wx(z), wy(z), wx(z) around the x, y, and warping directions, respectively. Once the FE displacement field (axial direction and shear) has been calculated, the deformation field is recovered as follows: ez ¼ e0z þ ðx þ n sin hÞjy þ ðy  n cos hÞjx þ ðx  nqÞjx

csz ¼ c0xz cos h þ c0yz sin h þ c0x r þ 2tc þ n þ 2tc jsz ð1Þ where x, y are in-plane coordinates defining the cross-section’s midsurface, and n is the thickness coordinate measured from the midsurface; e0z is linear strain along the blade axis; cxz, cyz, cx are transverse shear strain and warping strain; x is the warping function; jx, jy, jx, jsz are curvatures in the x (flapwise bending), y (edgewise bending), x(warping) and shear (torsional) directions; t is the wall thickness while k is the St. Venant circuit shear flow. In addition, r and h are geometric variables relating the local coordinate system of the blade section with a global coordinate system fixed at the blade root. The stress field can be calculated for arbitrary materials layups, where each lamina k with the composite material obeys a constitutive law given by:  rz rsz  ¼ " Q 11 Q 16 Q 16 Q 66 # ez csz  ð2Þ where the coefficients Q ij are obtained from the lamina’s stiffness coefficients, as defined in a local reference frame given by the materials principal axes, in the following manner: first the 3D orthotropic law is reduced to plane stress/strain conditions [49], and subsequently a transformation from the local to the global coordinates is performed. Once the finite-element stress field has been calculated, it must be transformed back to the local (material) coordinates in order to evaluate the chosen failure criteria at each layer:  r11 r12 r21 r22 K ¼ ½T  rz K T ½T rsz  SS12(+), finally, are the in-plane shear strength values. The reason for choosing the Chamis–Minnetyan criteria over other, often more elaborate, failure criteria lies with the fact that these criteria have been implemented in the GENOA package [47] designed for progressive failure in a shell-based environment, thereby providing a useful reference tool against which the reduced-order model presented here can be compared (see Section 2.3 for a brief introduction into GENOA). 2.2. Implementation of the integrated TWB-PFA model As stated above, for the purposes of the present work only static load are considered. Fig. 1 shows the flowchart of this implementation, as described in the following. The general procedure consists in an iterative process where the materials properties are progressively degraded as the load is gradually increased. The initial step of each loop of a static TWB-PFA simulation consists in the definition of a nodal load vector; {f} = {f0} at the first iteration. Then, the TWB finite element model is solved for nodal displacements, where the current materials properties of the blade are used; the nodal displacements are stored in the vector {u}. Subsequently, the strain and stress field tensors expressed in the (n, s, z) coordinate system used by the TWB formulation are obtained for all layers and all azimuthal and spanwise locations. The stress tensor is then transformed at each layer into the corresponding materials coordinate system by means of a passive rotation (see Eq. (3)); this coordinate system is understood as the one defined by the principal axes of the corresponding material at this layer. Once the stresses in a given layer are known, the failure criteria given by Eq. (4) can be evaluated. If no increase in damage is found, the nodal load vector {f} is increased by a small step Df1, and the analysis is repeated. If, on the other hand, the increase in damaged volume is higher than a percentage pre-defined by the user, the load is reduced by Df2 (|Df2| < |Df1|), i.e. to an intermediate value between the last two. This maximum tolerable increase in damage volume for one simulation step can be adjusted by the user but it typically of the order of 1%. The objective of this procedure is to avoid too large increases in damage in one simulation step, which would lead to a loss of resolution. If a non-zero permissible damage level is observed in this step, then the load is maintained at a constant level for several simulation steps until the damage stabilizes. Only after stabilization the load is increased again. Progressive Failure Analysis (PFA) takes place whenever an acceptable (e.g. sufficiently small) ð3Þ where the subscripts 1 and 2 identify the directions along and transverse to the fiber, respectively, and [T] collects the director cosines between the global and local coordinate systems. The local stress values calculated in this manner can then be compared with the selected failure criteria. In the present case, the failure criteria proposed by Chamis and Minnetyan [46] have been chosen. No failure occurs while the local stresses at each lamina remain within the limits shown below: SL11C < r11 < SL11T SL22C < r22 < SL22T ð4Þ Ss12C < r12 < Ss12ðþÞ where SL11C and SL11T are the strength values for the compressive and tensile stress, respectively, for the along-fiber direction. Similarly, SL22C and SL22T are the in-plane compressive/tensile strength values for the direction perpendicular to the fibers; SS12() and Fig. 1. PFA simulation cycle. 56 D. Cárdenas et al. / Composite Structures 95 (2013) 53–62 increase in damage volume occurs. In such a case, materials properties are updated according to the type of damage yielded by the evaluation of the failure criteria, which in turns requires updating the TWB model, leading to updated nodal displacements, strain/ stress fields and new damage assessment, all for the same amount of load {f}. This PFA cycle is iterated until damage stabilizes for a given load level {f}. Simulation ends when the percentage of damaged volume, as pre-defined by the user, is reached in one or more TWB elements. Fig. 2. Cross sectional area of the blade. 2.3. Finite-element shell model (reference model) The GENOA [47] package was used a reference tool for comparisons with the TWB-PFA model described in this work. GENOA combines a commercial Finite Element Analysis (FEA) package with providing Progressive Failure Analysis (PFA) on various scales [47]. GENOÁs approach for PFA is based on the Composite Durability Structural Analysis computer program, CODSTRAN [46], originally developed at NASA Glenn Research Center for polymer– matrix composite structures, and significantly extended in capabilities by the Alpha Star Corporation. It consists of three modules: composite mechanics, finite-element (FE) analysis, and damage progression modeling. The ICAN composite mechanics module [48] is recalled before and after each FE analysis. The module computes the composite properties from the fiber and matrix constituent characteristics and the composite layup. Then, the FE module accepts the composite properties at each node and performs the analysis at each load increment, computing generalized nodal forces and deformations. This information is subsequently supplied to the ICAN module which evaluates the amount and nature of local damage, if any, in all plies of the composite laminate. Individual ply failure modes are assessed by ICAN using failure criteria associated with the negative and positive limits of the six ply-stress components in the material directions. No failure occurs as longs as the stress values remain within the limits exhibited by the following equation: SL11C < r11 < SL11T SL22C < r22 < SL22t SL33C < r33 < SL33T SS12ðÞ < r12 < SS12ðþÞ SS23ðÞ < r23 < SS23ðþÞ Along-fiber ðlongitudinalÞtensile=compressive strength In-plane tensile=compressive strength Through-the-thickness tensile=compressive strength In-plane shear strength 1st through-the-thickness shear strength SS13ðÞ < r13 < SS13ðþÞ 2nd through-the-thickness shear strength ð5Þ where r and S are ply stress and strength in directions defined by numerical indexes (1, 2, 3) according to usual tensorial notation, respectively. Directions 11, 22, 33 are given in relation to fiber orientation: along, in-plane transversal and through-the-thickness transversal, while 12, 23, 13 represent in-plane, first and second out-of-plane shear directions. Strengths are calculated by ICAN based on constituent fiber and matrix strengths and micromechanics equations [48]. In addition to failure criteria based on the stress limits described above, a modified distortion energy (MDE) failure criterion that takes into account combined stresses can also be considered [46]. For the current analysis, however, this failure criterion was not considered as to facilitate comparison with the TWB-PFA model. 2.4. Description of case study Both the reduced-order TWB-PFA model presented in this work and the GENOA package were used to conduct a comparative study for the damage prediction in a structural model of a helicopter blade as described by Pawar and Ganguli [32–33]. The bladés external geometry is based on the NACA 0012 airfoil with a constant 305 mm chord and a two-cell cross sectional area, as shown in Fig. 2; the total blade length is 5080 mm. The D-spar and skin sections are divided at 35% of the chord. The blade reported in [32–33] does not specify the use of a web; however, in this work both webless and webbed cases were simulated in order to facilitate extending results to wind turbine blades. For the latter case, the web consisted of a 2 mm width solid cross-section the centreline of which is located at 35% of the chord, running throughout the entire length of the blade. Table 1 lists the elastic and strength materials properties for Carbon Fiber Reinforced Polymer (CFRP, used for the D-spar and Skin sections) and Balsa (used for the web) used in the blade. Table 2 details the material layup for the D-spar and Skin sections, specifying the thickness of each layer from the outmost to the innermost. Based on the databases corresponding to Tables 1 and 2, two FE models were generated: Model 1 consists of 13,600 4-node shell elements, 13,467 nodes and 80,802 degrees of freedom (DOFs) and was subjected to progressive failure analysis using GENOA. Model 2 consists of 100 TWB elements, 101 nodes and 707 DOFs. Both models were subject to equivalent boundary conditions (fixed at the root end) and loading (tip force applied in the out-of-plane direction). To assess the extent of cross-section deformation in the Shell model and therefore the range of applicability of the TWB-PFA model, two versions of model 1 (webless and webbed) were constructed. 3. Results and discussion Fig. 3 illustrates the overall blade tip response during PFA. It can be observed that the three models (shell webless, webbed and TWB) have identical flexural stiffness for small loading, the region where no or only negligible damage would be expected. It should be noted that the reduced-order (TWB) model provides an excellent reproduction of the predictions of the far more detailed shell model. For higher loadings, the webless shell model exhibits a rapid decrease in flexural stiffness, explained by the collapse of the blade’s cross section, with a deviation from the other models of about 10% at a tip load of 1000 N and a rapidly growing discrepancy for higher loads. This reduction in flexural stiffness can be traced back to a collapse of the cross-sectional area of the webless Table 1 Material elastic and strength properties. E11 (MPa) E22 (MPa) G12 (MPa) V12 SL11 T (MPa) SL11 C (MPa) SL22 T (MPa) SL22C (MPa) SL12S (MPa) CFRP (D-spar and Skin) Balsa wood (web) 206,000 20,700 8300 0.30 1979 1000 59 223 103 2070 2070 848 0.22 13 13 13 13 3 D. Cárdenas et al. / Composite Structures 95 (2013) 53–62 57 Table 2 Material layup for the airfoil profile, listed from the outmost to the innermost layer. # Layer 1 2 3 4 5 6 7 D-Spar Skin Thickness (mm) Fiber angle (°) Thickness (mm) Fiber angle (°) 0.762 0.508 0.508 1.016 0.508 0.508 0.762 0 45 45 90 45 45 0 0.381 0.508 0.508 0.254 0.508 0.508 0.381 0 45 45 90 45 45 0 blade, as illustrated in Fig. 4. In this figure the cross section measured at 450 mm from the root end, the location where the most notorious change in profile area occurs, has been plotted for different values of the tip load. Interestingly, a reduction in cross-sectional area can be observed for load values well before a difference in flexural response becomes conspicuous, with the onset of the cross-sectional collapse occurring for load values as low as 200 N. The webbed shell model, on the other hand, maintains its cross-sectional area up to load values well into the region where structural damage occurs, as shown below. For a load value of 1000 N the reduction of the cross section is barely 5%, and still for 1500 N the decrease is limited to about 10%. It is evident from these findings that the load carrying capabilities of the blade are widely extended by the use of the web, in accordance with good engineering practices in other fields, such as wind turbine blade engineering [50]. Apart from providing higher load carrying capacity, the use of the web ensures the approximate constancy of the cross-sectional area and justifies the use of the TWB model used in this work. After these introductory findings, the core of the result section will concentrate on the comparison of the TWB and the webbed shell model. In a first level of analysis the total damage volume of both the TWB and the shell models was calculated and plotted against the load value (Fig. 5). Total damage volume is defined as the sum of damaged TWB or shell elements, respectively, resolved both by layer and span/contour-wise location. It is conspicuous that the prediction of the onset of damage is very similar in the TWB and webbed shell model, with an initial damage occurring at a load of about 1250 N in the case of the TWB model and 1300 N for the shell model. This corresponds to an error of about 3.8% in the prediction of the load where initial damage occurs, which is a fairly reasonable agreement. It should be noted that the TWB errors on the safe (conservative) side. Consistent with the slightly earlier onset of damage in the TWB model is the larger rate of increase with load which amounts to about 1.5% per 1000 N of load for the TWB in the linear region of Fig. 5 versus 1.2% per 1000 N of load for the Fig. 3. Out-of-plane blade tip displacement versus applied load. Fig. 4. Evolution of the blade cross-section at 450 mm from the root end as a function of the applied load. (webbed) shell model. While this difference in slope is far more substantial (25%) than the difference in the damage onset values, the subsequent layer-resolved damage analysis (below) will demonstrate that this discrepancy is actually much smaller if more appropriate error metrics are used. For the moment it can be stated that the general trend predicted by the TWB and the webbed shell model beyond initial failure is very similar, including the change of slope occurring at about 2000 N. Not unexpectedly, the webless model shows an onset of damage at about 1050 N, well before the TWB and the webbed shell model. At this load value the cross-section at 450 mm from the root end (Fig. 4) can be seen to be partially collapsed (at about 75% of its original area), leading to significant loss of the section’s structural properties. In the following, the comparative study will not further consider the webless model. As stated earlier, the main motivation for using the TWB-PFA model described in this work is its capability of providing detailed damage maps for different load values, thereby allowing the use of the model as a structural design tool, identifying and modifying critical areas in a sequence of rapid design cycles. In the following only the webbed shell and the TWB will be compared in detail. As described above, the damage maps are constructed from the nodal displacements of the low-DOF TWB model based on the analytical expressions (Eq. (1)) between the nodal displacements and local strains at arbitrary locations along the perimeter of a given contour and through the wall. In the following only the webbed shell and the TWB will be compared in detail. Only layers 4 (90°) and 7 (0°) of the thin-walled structure present damage at any of the load values studied. This finding is consistently predicted by both the TWB and the shell model. In order to compare the predictions of both models in a detailed manner, a damage map for layer 4 has been plotted for four different load values for both the TWB and the shell model; the results are shown in Fig. 6. Damaged cells are shown in red (dark) on top of the backdrop of a three-dimensional visualization of the Fig. 5. Percentage of damage volume vs. load for the three blade models. 58 D. Cárdenas et al. / Composite Structures 95 (2013) 53–62 (a) (b) Fig. 6. Damage progression of layer 4 (90°). (a) Plot showing location of damaged cells in layer 4 for different values of the external load as predicted by the GENOA model. (b) Corresponding plot built from the results obtained with the TWB-PFA model. Note that different scales have been used for visualization the deflection of the beam in the GENOA model (a) and the TWB-PFA model (b). blade, where undamaged cells are displayed in light blue. Please note that due to the different visualization tools used the scales of the blade displacements are different in the TWB and the shell part, respectively. Please also refer to Fig. 3 showing that the blade displacements of both models are essentially identical for identical loads. It should also be noted that due to the algorithmic nature of the load variation in both PFA tools the load value sequences corresponding the two models in Fig. 6 are not identical, although reasonably similar. It is conspicuous from Fig. 6 that the damage maps for corresponding loads at layer 4 are strikingly similar between the two models. Most importantly, both models predict the occurrence of damage under tensile stress only. Secondly, in both cases the damage zone is a tongue-like structure which propagates from the root zone towards the spanwise direction. While the damage region at the lowest load value shown in Fig. 6 is a little larger in the TWB model than the corresponding region in the shell model, the increase in size is practically identical in both cases. Taking the total length of the damage size as one possible metric, a linear increase as a function of load can be observed for load values higher than about 1600 N (the three highest load cases for both models). This increase is 1.84 mm per applied N of load for both models. If the first load point is included (leading to a somewhat less linear relationship), then the incremental length of the damage zone is 2.2 mm/N for the shell model and 2.3 mm/N for the TWB model, still only a 5.5% difference. The striking similarity of the topology of the damage zones for the two models is further illustrated by a plane map for the four load cases studied (Fig. 7). The similarities include the contour-wise location of the damage zone which is centered at about 60 mm for both models. As expected from its more detailed structure, the shell model shows more irregular damage patterns, but the overall topology is very similar. Damage in layer 4 is caused on the tension side of the blade and occurs because of the 90° orientation of the fibers with respect to the blade axis, leading to a situation where much of the load has to be carried by the matrix. Layer 7, on the other hand, fails on the compressed side of the blade, which can be explained by the fact that the stresses are highest in this layer because of the 0° fiber orientation and the strength for compression loads is relatively low compared to its tensile strength. Additional insight into the progression of damage at layer 4 can be obtained by studying the cumulative damage along the blade, as shown in Fig. 8. Three curves have been plotted in the sequence of four subplots, each corresponding to one load situation: the predic- (a) (b) Fig. 7. Top view of the damage progression of layer 4 (90° fiber orientation). (a) Damage map calculated from GENOA for different values of the tip load. (b) Corresponding damage map obtained with the TWB-PFA model. D. Cárdenas et al. / Composite Structures 95 (2013) 53–62 59 tions directly obtained from both the TWB and the webbed shell model, respectively, and an additional TWB model output curve with adjusted failure strength values. The adjustment factor was equal for all strength values at a given load, but was allowed to vary among load values. In all load cases a relatively small adjustment of the strength values used in conjunction with Eq. (4) is sufficient to reproduce the predictions of the more detailed (shell) model. The criterion use to fit the cumulative damage curves was the requirement that the total damage along the blade for a given load value should be identical. In the case of the lowest load value this was achieved by increasing the strength values used in the TWB-PFA model by 6.5%, whereas in the higher load cases the increase was 9.5%, 12%, and 14.5%, respectively. A still fairly acceptable fit is achieved if a uniform increase of 10% is used throughout. Even in the cases of the best fits for each load case the TWB-PFA can still be seen to slightly overpredict the damage at the root section, but the error is typically only of the order of a few percent. Moreover, the error is conservative in all cases, i.e. a designer relying on the TWB-PFA model will always err on the conservative side. Fig. 9. Damage progression of layer 7 (0°) as shown as damage maps for the two load cases where damage was observed. Left: Predictions of the shell model. Right: Predictions of the TWB-PFA model. Fig. 8. Cumulative damage volume as a function of the spanwise coordinate for layer 4 (90° fiber orientation) and four load cases. Continuous curves: Predictions of the shell model. Dashed curves: Predictions of the unadjusted TWB-PFA model. Fine-dashed curve: Predictions of the TWB-PFA model with adjusted strength values. While the agreement in the detailed predictions of the damage progression between the TWB-PFA and the GENOA (shell) model is not perfect, it is still surprisingly good, particularly if the huge reduction in degrees of freedom (DOF) is considered. Whereas the shell model needs 80,202 DOFs to accurately represent the geometry and materials layup of the blade, the reduced-order model uses only 707 DOFs, a mere 0.88% of the DOFs of the shell model. Consequently, the TWB model is generally more rigid and less capable of accommodating the externally imposed load by straining its internal DOFs. Evidently, the TWB-PFA model cannot account for the interaction between a web and the blade shell such as in the case of the webbed shell model used here in conjunction with GENOA. As shown in Figs. 9 and 10 for the case of the 1800 N load case a needle-like longitudinal damage feature can be seen to emerge in the shell model which is not predicted by the TWB model. This feature can be attributed to the presence of the web in the shell model, leading to a localized concentration of stress in the interface with the web. Apart from this feature, the agreement between the damage topologies predicted by both models is strikingly similar. Regarding the total damage volume at the two load values where damage arises (about 1800 and 2000 N), the agreement between the predictions between the TWB and GENOA is remarkably good (about 0.025% at 1800 N and 0.3% at 2000 N for both models, Fig. 11). It should be emphasized that layer 7 represents only a small fraction of the total damage in the blade (0% at the lowest two load values, about 0.025% of the total volume versus about 1% in layer 4 for the 1800 N load case, and about 0.3% versus 1.2% in layer 4 for the 2000 N load case), so discrepancies in this layer are of no practical concern at the design stage. Even so, the TWB correctly predicts the onset of damage, the total damage volume and (roughly) even the shape of the damage zone, with the exception of the needle-like feature caused by the presence of the web. This is quite an achievement for a reduced-order model. For forensic analyses, a more detailed model may be indispensible, but at the design stage the TWB-PFA model may be good enough. A few remarks on computational economy: setting up the shell model took about 3 person-h, while the corresponding task for the TWB model was estimated at 0.5 person-h. Execution time was some 4 h for the shell model, compared to 5 min for the reducedorder model. This difference in set-up and computation time is likely to increase for more complex structures than the helicopter blade studied in the present work. 60 D. Cárdenas et al. / Composite Structures 95 (2013) 53–62 Fig. 10. Top view of the damage progression of layer 7 (0° fiber orientation). Upper graph: Predictions of the shell model. Lower graph: Predictions of the TWB-PFA model. Fig. 11. Cumulative damage volume as a function of the spanwise coordinate for layer 7 (0°). Continuous curve: Predictions of the shell model. Dashed curve: Predictions of the TWB-PFA model. 4. Summary and conclusions In the present work, a reduced-order Thin-Wall Beam (TWB) Progressive Failure Analysis (PFA) model has been presented and a detailed validation of its predictions for a helicopter blade subject to static tip loads against a more elaborate shell model, built in the commercial PFA software suite GENOA, has been conducted. The TWB-PFA represents an equivalent one-dimensional model with a dramatically reduced number of degrees of freedom (DOF) compared to the detailed shell model. In spite of being a one-dimen- sional representation of the blade, the TWB-PFA model is capable of generating detailed maps of the progressively occurring damage in response to increasing external loads by taking advantage of the analytical relationships between the node displacements and the longitudinal and shear strains which, together with the constitutive laws for each lamina of the composite material, allowing applying failure criteria to the blade resolved both by layer and by contour-/spanwise position. Different levels of analysis were carried out in course of the study. At the first level the total damage volume as a function of the external load was calculated for both approaches; a similar value of the onset of damage was predicted by the shell (GENOA) and the TWB-PFA model if the shell model of the blade was equipped with an internal web, thereby avoiding a collapse of the cross sections near the root and ensuring the validity of one of the key assumptions of the TWB model used. The damage-onset value was predicted with an accuracy of about 5.5%, indicating a promising initial evidence for the suitability of the reduced-order model. While the damage after the onset progresses at a rate about 25% higher in the TWB-PFA model (1.5% of damage per 1000 N of load vs. 1.2%/1000 N for the shell model), this difference can be substantially reduced by a small adjustment of the strength values used with the failure criteria for the TWB-PFA model, as discussed below. In a second level of analysis the spatial distribution of the damage was studied for different load levels. Both in the TWB-PFA and the GENOA model damage is predicted in layers 4 and 7 only, and the topology of the damage maps is strikingly similar for both models. In the case of layer 4, where most of the damage occurs, the shape and the contour-wise location of the damage zone predicted by the reduced-order model is almost identical to the predictions of the more complex model, although the total length of the failure zone is somewhat larger for all load cases, which can be traced back to an almost constant offset between the two predictions, since the rate of growth of the length of the failure region D. Cárdenas et al. / Composite Structures 95 (2013) 53–62 is almost identical for both models. In the case of layer 7, where very little damage occurs and no damage is observed until fairly high load values (about 1800 N of tip load for both models) the shape of the damage zone is well predicted by the TWB-PFA model except for a needle-like damage feature which can be traced back to the presence of the web in the shell model. In a third level of analysis the layer-resolved cumulative damage volume as a function of the spanwise coordinate has been calculated for both models as a function of the applied tip load and the predictions have been compared. Consistently with the findings described above, the cumulative damage volume found from the TWB-PFA model is somewhat higher than in the GENOA model. A possible explanation for this behavior lies with the locally stiffer structure of the reduced-order model given its much smaller number of degrees of freedom (DOFs) and the correspondingly smaller ability of distributing the strain energy among the DOFs. Another issue refers to the appropriateness of the failure criteria, which may have to be modified for an effective reduced-order model. In order to explore this line of thought the strength values used in the TWBPFA model were increased by a constant factor and the simulations were re-run for layer 4 and all load cases, matching the total cumulative damage volume. It was possible to quite accurately reproduce the cumulative damage curves of the GENOA model by a relatively small adjustment in materials strength, with the adjustment factor ranging from 1.065 (lowest load case) to 1.14 (highest load case); a constant adjustment factor of 1.1 (10% higher) still yielded acceptable results. The increase in the value of the adjustment factor from 1.065 to 1.14 can be explained by the fact that the damage progressively weakens the contour of the affected sections, thereby allowing for a progressive reduction of the section stiffness due to the decrease in cross-sectional area. As discussed above, this behavior is correctly modeled by the shell but not the TWB model. In conclusion, the TWB-PFA model presented in this work has shown its capability of correctly predicting both the onset and the propagation of damage in the composite blade studied for the case of static loads, providing reliable information on the location of the damage zone (resolved both by layer and the contourand spanwise location of the damage). A small overprediction of the damage compared to the more detailed shell model can be corrected by a small upward adjustment of the (equivalent) strength values to be used in the model. In all cases, the reduced-order model proved to be more conservative that the shell model, so it should be possible to safely use it for design purposes, especially if a series of iteration with quick turnaround times is desired. While the final design of a rotor blade will still remain confined to the domain of higher resolution models such as shell and volume element models, the TWB-PFA approach is a promising tool for pre-design and aeroelastic design simulation, where the use of high-resolution models would be prohibitive. Acknowledgements Financial support from Tecnológico de Monterrey, including funds obtained through the Research Chair for Wind Energy (Contract No. CAT158), are gratefully acknowledged. One of the authors (PM) acknowledges support from the New York State Energy Research and Development Authority (NYSERDA) Project No. 18812, and the National Science Foundation, Grant No. CMMI-1031036, for partially funding this research project. Finally, the support of Alpha Star Corporation, provided in term of software and training, is greatly appreciated. References [1] Hau E. Wind turbines: fundamentals, economics. Berlin: Springer-Verlag; 2006. technologies, application, 61 [2] Sobieszczanski-Sobieski J, Haftka, R. Multidisciplinary aerospace design optimization: Survey of recent developments. In: American institute of aeronautics and astronautics, proceedings of the 34th aerospace sciences meeting and exhibit, Reno; 15–16 January, 1996. p. 1–33. [3] Li L. Structural design of composite rotor blades with consideration of manufacturability, durability and manufacturing uncertainties. Doctoral thesis in school of aerospace engineering. Georgia Institute of Technology; 2008. [4] Shahverdi H, Nobari AS, Behbahani-Nejad M, Haddadpour H. Aeroelastic analysis of helicopter rotor blade in hover using an efficient reduced-order aerodynamic model. J Fluids Struct 2009;25(8):1243–57. [5] Friedmann P, Glaz B, Palacios R. A moderate deflection composite helicopter rotor blade model with an improved cross-sectional analysis. Int J Solids Struct 2009;46(10):2186–200. [6] Murugan S, Chowdhury R, Adhikari S, Friswell MI. Helicopter aeroelastic analysis with spatially uncertain rotor blade properties. Aerosp Sci Technol 2011. http://dx.doi.org/10.1016/j.ast.2011.02.004. [7] Hansen MOL, Sorensen JN, Voutsinas S, Sorensen N, Madsen HA. State of the art in wind turbine aerodynamics and aeroelasticity. Prog Aerosp Sci 2006;42(4):285–330. [8] Jonkman JM, Sclavounos PD. Development of fully coupled aeroelastic and hydrodynamic models for offshore wind turbines. In: ASME, proceedings of the wind energy symposium, Reno; 10–12 January, 2006. [9] Ahlstrom A. Aeroelastic simulation of wind turbine dynamics. Doctoral thesis in structural mechanics. KTH, Sweden; 2005. [10] Buhl ML, Manjock A. A comparison of wind turbine aeroelastic codes used for certification. In: 44th AIAA aerospace sciences meeting and exhibit, Reno, Nevada; 9–12 January, 2006. [11] Passon P, KÜhn M, Butterfield D, Jonkman J, Camp T, Larsen TJ. Benchmark exercise of aero-elastic offshore wind turbine codes. In: Proceedings of EAWE conference, Lyngby, Denmark; 28–31 August, 2007. [12] Volovoi V, Hodges DH, Cesnik CES, Popescu B. Assessment of beam modelling methods for rotor blade applications. Math Comput Model 2001;33:1099–112. [13] Hodges D, Atilgan AR, Cesnik CES, Fulton MV. On a simplified strain energy function for geometrically nonlinear behaviour of anisotropic beams. Compos Eng 1992;2(5):513–26. [14] Jonkman JM, Marshall LB. FAST user’s guide. Technical, report NREL/EL-50038230; August 2005. [15] Noda M, Flay RGJ. A simulation model for wind turbine blade fatigue loads. J Wind Eng Ind Aerodynam 1999;83(1–3):527–40. [16] White D. New method for dual-axis fatigue testing of large wind turbine blades using resonance excitation and spectral loading. Technical, report NREL/TP-500-35268; April 2004. [17] Nijssen RPL. Fatigue life prediction and strength degradation of wind turbine rotor blade composites. Doctoral thesis in design and production of composite structures. Delft University of Technology; 2006. [18] Epaarachchia JA, Clausen PD. The development of a fatigue loading spectrum for small wind turbine blades. J Wind Eng Ind Aerodynam 2006;94(4):207–23. [19] Marin JC, Barroso A, Parıs F, Canas J. Study of fatigue damage in wind turbine blades. Eng Fail Anal 2009;16(2):656–68. [20] DimarogonaS AD. Vibration of cracked structures: a state of the art review. Eng Fract Mech 1996;55(5):831–57. [21] Gadelra RM. The effect of delamination on the natural frequencies of a laminated composite beam. J Sound Vib 1996;197(3):283–92. [22] Yan YJ, Chen L, Wu ZY, Yam LH. Development in vibration-based structural damage detection technique. Mech Syst Signal Process 2007;21(5):2198–211. [23] Zou Y, Tong L, Steven GP. Vibration-based model-dependent damage (delamination) identification and health monitoring for composite structures. J Sound Vib 2000;230(2):357–78. [24] Dilena M, Morassi A. Identification of crack location in vibrating beams from changes in node positions. J Sound Vib 2002;255(5):915–30. [25] Nahvi H, Jabbari M. Crack detection in beams using experimental modal data and finite element model. Int J Mech Sci 2005;2005(47):1477–97. [26] Alvandi A, Cremona C. Assessment of vibration-based damage identification techniques. J Sound Vib 2006;47(10):179–202. [27] Hu H, Wang BT, Lee CH, Su JS. Damage detection of surface cracks in composite laminates using modal analysis and strain energy method. Compos Struct 2006;74(4):399–405. [28] Karthikeyan M, Tiwari R, Talukdar S. Crack localization and sizing in a beam based on the free and forced response measurements. Mech Syst Signal Process 2007;21(3):1362–85. [29] Shih HW, Thambiratnam DP, Chan THT. Vibration based structural damage detection in flexural members using multi-criteria approach. J Sound Vib 2009;323(3–5):645–61. [32] Pawar PM, Ganguli R. Modeling progressive damage accumulation in thin walled composite beams for rotor blade applications. Compos Sci Technol 2006;66(13):2337–49. [33] Pawar PM, Ganguli R. On the effect of progressive damage on composite helicopter rotor system behavior. Compos Struct 2007;78(3):410–23. [34] Adolfsson E, Gudmundson P. Matrix crack induced stiffness reduction in [(hm/ 90n/+hp/hq)s]m composite laminates. Compos Eng 1995;5(1):107–23. [35] Gudmundson P, Zang W. An analytic model for thermoelastic properties of composite laminates containing transverse matrix cracks. Int J Solids Struct 1993;30(23):3211–31. [36] Shahid IS, Chang FK. An accumulative damage model for tensile and shear failure of laminated composite plates. J Compos Mater 1995;29:926–81. 62 D. Cárdenas et al. / Composite Structures 95 (2013) 53–62 [37] Chandra R, Chopra I. Structural response of composite beams and blades with elastic couplings. Compos Eng 1992;2(5–6):347–74. [38] Song O, Librescu L, Jeong NH. Static response of thin-walled composite I-beams loaded at their free-end cross-section: analytical solution. Compos Struct 2001;52(1):55–65. [39] Librescu L, Na S. Active vibration control of doubly tapered thin-walled beams using piezoelectric actuation. Thin-Walled Struct 2001;39(1):65–82. [40] Qin Z, Librescu L. On a shear-deformable theory of anisotropic thin-walled beams: further contribution and validations. Compos Struct 2002;56(1):345–58. [41] Librescu L, Song O. Thin-walled composite beams: theory and application. Springer; 2006. [42] Vo T, Lee J. Flexural torsional behavior of thin walled closed-section composite box beams. Eng Struct 2007;29(8):1774–82. [43] Lee J, Vo T. Flexural–torsional behavior of thin walled composite box beams using shear-deformable beam theory. Eng Struct 2008;30(7):1958–68. [44] Vo T, Lee J. Free vibration of thin-walled composite box beams. Compos Struct 2008;90(2):11–20. [45] Cardenas D, Elizalde H, Marzocca P, Probst O. Numerical validation of a finite element thin-walled beam model of a composite wind turbine blade. Wind Energy 2011. http://dx.doi.org/10.1002/we.462. [46] Chamis CC, Minnetyan L. Defect/damage tolerance of pressurized fiber composite shells. Compos Struct 2001;51(2):159–68. [47] http://www.ascgenoa.com/main/. [48] Murthy PLN, Chamis CC. ICAN: integrated composite analyzer user’s and programmer’s manual. Technical paper NASA TP-2515. National Aeronautics and Space Administration, Washington, DC; 1985. [49] Tuttle M. Structural analysis of polymeric composite materials. New York: Marcel Dekke; 2004. [50] Lindenburg C, Winkel GD. State of the art of rotor blade buckling tools. Technical report ECN-C-05-054. Energy Center of the Netherlands; 2005.