PROGRESSIVE FAILURE ANALYSIS OF HELICOPTER ROTOR BLADE UNDER AEROELASTIC LOADING

Unlike metal structure, composite structures don’t give any clue till the fatal final collapse. The problem is more complicated when applied load on the structure is aeroelastic in nature. Under such loading, composite laminate experiences stresses. The first layer failure happens when stresses in the weakest ply exceed the allowable strength of the laminate. This initial layer-based failure changes overall material characteristics. It is important now to degrade the composite laminate characteristics for the subsequent failure prediction. The constitutive relations are required to be updated by the reduction in stiffness. The rest of the undamaged laminates continue to take the load till the updated strength is reached. In the present work, layer wise progressive failure analysis under aeroelastic loading has been performed by the inclusion of different failure criteria which allow for the identification of the location of the failure. ANSYS APDL environment has been used to model geometry of helicopter rotor. Under the loading conditions, stresses are calculated in the blade. Using stress tensor and failure criteria, failure location and modes have been predicted. It has been found that failure starts at higher speeds and failure starts from the root chord and tend towards the tip chord.


Introduction
Composite materials have shown distinct characteristics in term of strength and weight in such a way that all industries in general and aerospace industry in particular are drifting towards modifying their manufacturing processes and using composite materials in their product development. Helicopter blades are also made of composites (Alkahe & Rand, 2000;Azzam & Andrew, 1992;Morozov, Sylantiev, & Evseev, 2003;Shahid & Chang, 1995). Operating conditions of helicopter rotor are unsteady aerodynamic and acceleration/deceleration loads leading to harsh loading conditions. This may lead to induce damage in the composite material of the rotor. Therefore, progressive failure of helicopter rotor under dynamic loading conditions has attained significant attention of researchers. Because of the interdisciplinary nature of the work, few authors have investigated failure analysis. Modal data has been used by some authors to check damage in rotor blades (Pawar & Ganguli, 2003;Alkahe, Rand, & Oshman, 2003). Experimental studies, even more difficult aspect, have been investigated by few authors (Cattarius & Inman, 2000;Gho-shal et al., 2001;Kiddy & Pines, 2001). One way of composite damage modeling is by some percentage reduction in the bending and torsional stiffnesses while some have gone through crack modeling based on fracture mechanics approach (Ganguli, 2002). With the availability of light weight, high strength composite materials, it was felt to perform detailed modeling of damage modes in composite rotor blades. Notable efforts have been done by Pawar and Ganguli (2005) considering ply laminates and matrix cracking was the initial damage observed. Gudmundson and Weilin (1993), Adolfsson and Gudmundson (1997) have investigated generalized composite laminates under combined loading consisting of bending and axial loading. When the matrix cracking effect starts saturating, a local de-lamination can initiate from the matrix crack tips. The de-lamination induced due to matrix cracking was modeled in the literature using various approaches such as classical plate theory (Obrien, 1991), first-order shear deformation laminate plate theory (Armanios, Sriram, & Badir, 1991), variation techniques (Nairn & Hu, 1992), shear lag techniques (Dharani & Tang, 1990) and 3-D finite element analysis of 90° lamina (Salpekar & O'Brien, 1991). Most recently, the shear lag method was applied by Berthelot and Corre (2000) to analyze the stress fields in cross-ply laminates containing transverse cracks and crack tip de-lamination with shear friction between the delaminated plies. While several researchers have modeled separate damage modes as discussed above, few researchers have modeled the progressive damage accumulation in composite structures (Allen, Harris, & Groves, 1987). Most of these models are focused on characterization of cracking in a single sub-laminate with 90° plies. However, under realistic loading condition matrix cracks can appear in all plies and in a multi-directional laminate and may lead to laminate failure by de-lamination or fiber breakage. Such a model was proposed by Shahid and Chang (1995) for predicting the accumulated damage and studied the effect of such damage on the in-plane response of laminated composite plates subjected to tensile and shear loads.
Progressive damage of moderately thick composite laminated plates has been investigated by  considering in-plane compressive load and uniform lateral pressure. Under same loading conditions, Ghannadpour  presented analysis of laminate with different shapes and size cut-outs. The concept of minimum potential energy was used for the formulations. Hashin failure criteria have been used to predict the location of failure. Post buckling analysis of cracked imperfect composite plates has been explored by Ghannadpour and Karimi (2018). Progressive damage and nonlinear geometric analysis of relatively thick laminates has been presented by , . Two different versions of finite strip method have been used by  to predict geometrically nonlinear behavior of composite laminate under progressive end-shortening and to pressure loading. Three geometric degradation models have been employed to estimate the degradation zone around the failure location. Using classical plate theory and first order shear deformation theory, Dawe, Lam, and Azizian (1993) has described the analysis geometrically nonlinear composite laminate behavior to predict local buckling.  and (Ovesy, Ghannadpour, & Zia-Dehkordi), have used exact finite strip method to predict buckling analysis of moderately thick composite laminate. In their further work, Ghannadpour (Ovesy, Totounferoush, & Ghannadpour, 2015;Ghannadpour & Ovesy, 2009) and Ghannadpour et al. (2014) have presented de-lamination and dynamic buckling of composite laminate using semi analytical finite strip and exact finite strip methods. Other authors (Oktay & Sal, 2015, 2016Oktay & Sultan, 2013) have also explained active and passive control and helicopter rotor blade morphing for energy efficiency.
From the presented literature review, it is revealed that noticeable work has been done by several authors on the geometry and control of the helicopter rotor blade, how-ever, aeroelastic loading, centrifugal loads and their effect on the composite laminate structure along with failure mode prediction is the inter-disciplinary area in which not enough work has been done. The same has been explored in the current work.
In the present work, a composite helicopter rotor blade (Pawar & Ganguli, 2006) using shell and solid elements has been used. Progressive failure analysis (PFA) has been performed under static aeroelastic loading. Aerodynamic loads are calculated using strip theory. As the velocity is increased so does the load and hence damage is modeled with progressive increase in the load. (Helicopter blade is rotating and hence its vibration modes will be different with that at wind off conditions this effect has been accounted by axial loading of the blade and calculating pre-stressed modes). The key damage modes in composite structure such as matrix cracking, de-lamination and fiber breakage are modeled. The effect of progressive damage accumulation on the composite helicopter rotor blade in forward flight is investigated using a two-cell airfoil section-beam with stiffness properties representing a stiff inplane helicopter rotor blade.

Variation formulation
Helicopter blade undergoes flap (out-of-plane) bending, lag (in-plane) bending, elastic twist and axial displacement Using a generalized Hamilton's principal applicable to non-conservative systems (Pawar & Ganguli, 2007) is given as The δU and δT include energy contributions from components that are attached to the blade, e.g., pitch link, lag damper, etc. After finite element discretization, Hamilton's principle is written as Assembling the blade finite element equations and applying boundary conditions in Eq results in the final form of equation (3) (4)

Aerodynamic, inertia and centrifugal forces
Helicopter rotor blade is subjected to lift, gravity and centrifugal forces making its structural and aeroelastic analysis difficult. Gravity force acts downward, lift upward while centrifugal force acts radially outward as shown in Figure 1. The lifting force on a rotating blade is related to its local angle of attack and local dynamic pressure. There is velocity distribution on the blade, zero flow velocity at the rotational axis and the maximum at the tip. Moreover, geometrical twist angle increases from tip to root. Although, the velocity is maximum at the tip, lift is not maximum at the tip because of tip vortices. Lift produced by rotating blade may be described by Blade Element Theory (BET) as Eq.
here, the velocity V, and angle of attack α, are the function of blade radius r. BET assumes no influence of adjacent blade elements section. Moreover, non-uniform induced flow across the blade is obtained through a distribution of geometrical twist along the blade span. Linear geometrical twist has been assumed as given in Figure 2 in the current study. Assumed velocity and lift distributions at 15 rad/ sec are shown in Figure 3. Composite materials are used for the manufacturing of helicopter rotor blade. As the density is assigned to a particular laminate, its mass is calculated by the finite element software through its volume. Gravity force can then be applied globally and the same has been considered in current work. As the helicopter rotor blade rotates, centrifugal force act on the blade radially whose magnitude, mω 2 r, is much more than the gravity and lift force even at low rotation speed. Centrifugal force has been applied using the parameter of rotation speed, ω.

Progressive damage model
The progressive damage model involves the steps of stress analysis, failure analysis and material properties degradation. These steps, which operate in an iterative procedure following the mentioned order, were incorporated into ANSYS code using APDL. Stress analysis was performed numerically using the commercial Finite Element code ANSYS which can provide the ply-by-ply calculated stresses to be used as input to progressive model. Failure analysis was performed by implementing the commonly used Hashin polynomial failure criteria. These criteria were selected because they can distinguish between different failure modes and may be applied with ease in a finite element formulation. Once failure is detected in a ply, the material properties of the ply are degraded by implementing a sudden degradation rule. The scope of using degradation rules is to disable the ply from carrying a certain load. Thus, for each mode of failure there exists a corresponding degradation rule. The available ply-b-ply material property degradation rules are empirical and include assumptions resulting from engineering constraints in the properties of composite materials.

FEA calculation process
The previously described steps of progressive model were programmed into an ANSYS parametric macro-routine, which works iteratively. The macro-routine involves the following steps.
1. Generating 3-D model of rotor blade by giving as input the initial material properties, specimen geometry, boundary conditions and initial loads. 2. Performing non-linear stress analysis and calculation of the on-axis stresses for each ply. 3. Performing failure analysis by applying the failure criteria. 4. Checking for ply failures. 5. Increasing the applied load by an increment If no failure is detected. Then the program returns back to step 2 for stress analysis. 6. Continuing the program to the next step If any mode of failure is detected. Non-linear stress analysis with the same load is performed again in order to calculate the redistributed stresses. Convergence at a load step is assumed when no additional failures are detected. 11. Repeating the procedure until final failure occurs. The progressive damage model was implemented for rotor blade in order to evaluate its residual strength under incremental tensile loading.

Failure models
Failure criteria compare the loading state at a point with a set of values reflecting the strength of the material at that point (often referred to as the material allowable). For unidirectional materials, this is typically in the direction of the fibers.
In general, the load is represented by a full stress or strain tensor having six independent components. By convention, for lamina materials the material X axis lies in the direction of the warp fibers while the Z axis lies in the through-thickness direction of the sheet.

Damage criteria and degradation rules
Hashin damage criteria in the destruction of the fiber pull-off formula is given in (Tserpes et al., 2002). Damage criteria for fiber pull-out failure mode: Fiber tensile failure is given as: Fiber compressive failure is given as: Matrix tensile failure is given as: Matrix compressive failure is given as: Fiber matrix shear failure is given as: Delamination failure mode under tensile loading is given as: Delamination failure mode under compression load is given as: here, σ xx , σ yy σ zz are the normal stresses in uni-directional composites in 1, 2 and 3 directions respectively. σ yz and σ xz represent the shear stresses; X T and X C is the tensile and compressive strength in longitudinal direction Y T and Y C is the tensile and compressive strength in transverse direction respectively; S xy , S xz , S yz are the shear strengths in three planes. Z T and Z C are the tensile and compression strengths in thickness direction. It can be generally approximated that Z T and Z C are equal to the matrix tensile and compressive strength.
When the elements are damaged, there is a need for guideline rules related to the degradation of material properties. The stiffness degradation rules used in this research work are given in Table 1.

Computational model
In this study, structural model of a hypothetical helicopter blade as described by Pawar and Ganguli (2006), Pawar and Ganguli (2007) has been used. The blade outer geometry is based on the NACA 0012 airfoil with a constant 305 mm chord and a two-cell cross sectional area, as shown in Figure 4. The blade span is 5080 mm. The D-spar and skin sections are divided at 35% of the chord. 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). Table 2 and table 3 details the material layup for the D-spar and Skin sections, specifying the thickness of each layer from the outmost to the inner most. The same blade has been analyzed by D Cardenas et al. (Cárdenas et al., 2013) along with 2 mm balsa web placed at 35% of the chord as shown in Figure 4. Fixed boundary conditions were applied at the root end. Keeping in view the magnitude of centrifugal force which is much higher than the rest of the two forces, failure modes and the first ten modes of vibration (six flap, two lag, two torsion and no axial) have been evaluated with and without centrifugal forces. Figure 5 shows D-spar, shear web and the skin in the finite element model.

Finite element model validation
In order to validate the blade structural model, nonlinear results have been compared with those presented by D Cardenas et al. (2013) in Figure 6. For the validation purpose, only tip load was applied (no gravity force considered) and corresponding tip deflection data were plotted and the result shows good agreement between the two.

Results and discussion
An APDL code has been developed in ANSYS environment putting finite element structural model, lift, gravity, centrifugal forces and the composite failure criteria at one place making the code attractive for the solution of any practical problem. To demonstrate the capability of the code, benchmark hypothetical helicopter blade geometry was chosen and its structural model was validated in the last section. Nonlinear tip deflection matches closely to that of the literature value. There are three types of loadings a helicopter rotor blade experiences during its normal operation which are gravity, lift and centrifugal forces. Failure modes have been predicted first considering gravity and lift forces, later the centrifugal force is also added to see the combined effect.
In the first case, gravity and lift forces were applied as both act in the opposite direction. Lift force was varied through rotational speed. Five speeds were considered i.e., 11-15 rad/sec. As the speed is changed, lift force computed through BET is changed; however, gravity force remains constant. Failure in the composite materials is predicted through the failure criteria. Corresponding mode changes have also been computed. No failure is observed as the rotation speed of 11 rad/sec and 12 rad/sec is kept, however, failure is observed near the root section at 13 rad/sec as shown in Figure 7. Nodes plotted in black color experience no failure while matrix cracking was computed at the red plotted nodes. The composite damage progressively increases near the root section as the speed is increased to 14 rad/sec and 15 rad/sec as shown in Figures 8 and Figure 9. As the tip deflection increases with the increase of rotational speed, there is stress stiffening effect in the blade. This behavior is obvious from the plot of first 10 modes in Figure 10. As the rotational speed is set to 11 rad/sec, frequencies of all modes increase; however, there is no further stress stiffening as the speed is further increased from11 rad/sec to 15 rad/sec. In the next computational scheme, centrifugal force is also added along with lift and gravity forces. Out of plane deflections and twist about span axis with and without centrifugal force are shown in Figures 11 and Figure 12. It is obvious from the figures that centrifugal force acts as the stabilizing force and deflection/twist values have reduced in presence of centrifugal force. Moreover, no failure is observed as the rotational speeds are changed from 11 rad/sec to 15 rad/sec. Modal data is plotted in Figure 13. The trend of all modes remains the same as in absence of centrifugal force, however, there is more stress stiffening and hence frequencies at 11 rad/sec is more than those in absence of centrifugal force. Moreover, as expected there is axial and chord wise deformations, although smaller in magnitude, in the presence of centrifugal force as shown in Figures 14 and Figure 15. It can be observed from these figures that magnitude of centrifugal force generated at such low rotational speed, i.e., 15 rad/sec, is quite high as compared to lift, gravity and chord wise force. In actual operation, rotational speed is quite high than this value and hence centrifugal force consideration plays central role in the design of rotating machinery blade.

Conclusions
Helicopter rotor blade is subjected to complex aerodynamic loading and therefore, its structural design is of utmost importance. To get high performance of helicopter, rotors are being manufactured from composite materials. Unlike metal structures, composite materials failure modes are complex and no clue is given by these structures before their failure. Hence, investigation of composite structures failure modes has gained popularity in the recent research. In the present work, an APDL code in ANSYS environment has been generated encompassing structural model of composite laminate and failure criteria of composites along with gravity, lift and centrifugal forces at one place making the code very attractive for the nonlinear static solution of practical problems. To demonstrate the code capability, a hypothetical helicopter rotor blade from the literature was selected and structural model was validated through the comparison of deflections. It is subjected to gravity and lift forces and failure modes were observed because of out of plane deflections. It was realized that damage starts near the root section and it was matrix cracking. Modal data predicted a stress stiffening effect and hence frequency values of first 10 modes increased. Centrifugal force was further added by rotating the structures at different speeds. It was realized that centrifugal force acted as a stabilizing agent as it restricted the out of plane deflections during the axial rotations. Hence no failure was observed in the presence of centrifugal force on the given speeds. Modal results indicated further stress stiffening and frequency values for first 10 modes further increased. Additionally, the trend of increasing frequency values for bending and torsion with increasing speed is same. In the future work, dynamic analysis will be incorporated to make the code more realistic for practical problems of rotating machinery. It is pertinent to mention that although all possible failure models have been incorporated in APDL code, it is not necessary that all failure modes have been experienced. It depends on the loading conditions and their magnitude. However, it is clarified that as the matrix failure is observed in specific region, stiffness of that region goes to zero and it is not further able to withstand the load.

Funding
This work was supported by the Pakistan Science Foundation under Grant No: PSF/P&D/TG-ND (463)/19.