- Email: [email protected]
Modeling composite wing aeroelastic behavior with uncertain damage severity and material properties
Mechanical Systems and Signal Processing 32 (2012) 32–43
Contents lists available at SciVerse ScienceDirect
Mechanical Systems and Signal Processing journal homepage: www.elsevier.com/locate/ymssp
Modeling composite wing aeroelastic behavior with uncertain damage severity and material properties G. Georgiou a, A. Manan a, J.E. Cooper b,n a b
School of Engineering, University of Liverpool, Liverpool, L69 3GH, UK Department of Aerospace Engineering, University of Bristol, Bristol, BS8 1TH, UK
a r t i c l e in f o
abstract
Article history: Received 20 June 2011 Received in revised form 1 May 2012 Accepted 7 May 2012 Available online 16 June 2012
The effect of uncertain material properties and severity of damage on the aeroelastic behavior of a finite element composite wing model are predicted by applying the Polynomial Chaos Expansion method (PCE). Different damage modes, including the transverse matrix cracking and broken fibers, are induced into pre-defined locations in the laminates and the aeroelastic stability and dynamic response of the wing due to ‘‘1-cosine’’ vertical gusts are evaluated. For this purpose, PCE models that predict the variation due to uncertainty of the flutter speed and an ‘‘Interesting Quantity’’ (root shear force) of the wing box are developed based upon a small sample of observations, exploiting the efficient Latin Hypercube sampling technique. The uncertainty propagation on the output responses, in the form of probability density functions, is evaluated at low computational cost, implementing the PCE models and verified successfully against the actual results. & 2012 Elsevier Ltd. All rights reserved.
Keywords: Composite wing Damage severity Uncertainty propagation Aeroelastic behavior
1. Introduction Composite materials are being used increasingly in aerospace structures, primarily due to their attractive strength– weight ratios. A further advantage is the inherent anisotropic behavior which can be exploited to induce beneficial bending–torsion couplings in aircraft wings during flight, thus providing a greater design flexibility compared with metallic structures. A number of aeroelastic tailoring applications have sought to use composite materials for a range of applications, including the increase of divergence speed, the minimization of drag, the reduction of gust loads, etc.[1–8]. However, the complex manufacturing process of fibrous laminates induces an increased number of possible uncertain geometrical (thickness) and material (in-plane moduli, fiber orientation, etc.) parameters, which affect in a considerable degree the strength and the overall performance of composite structures. Additionally, under monotonic or cyclic loading the strength difference between the fibers and the matrix is usually related to damage mechanisms. The damage process is initiated from defects randomly distributed in a ply, leading to the development of transverse and longitudinal matrix cracks, inter-laminar delamination and fiber fracture [9]. The transverse matrix cracking, which starts at low loads and propagates in the fiber direction, is the predominant damage mode related to a significant reduction of the effective stiffness and strength for composite laminates [10]. Taking into account that the initiation and propagation of transverse matrix cracks is a stochastic procedure, the severity of damage in composite materials represents a considerable source of uncertainty. Moreover, at high loads, the matrix cracking triggers also the growth of resin-dominated damage modes, such as delamination and fiber breakage, which is
n
Corresponding author. Tel.: þ 44 117 954 5388; fax: þ 44 117 927 2771. E-mail addresses: [email protected] (G. Georgiou), [email protected] (J.E. Cooper).
0888-3270/$ - see front matter & 2012 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.ymssp.2012.05.003
G. Georgiou et al. / Mechanical Systems and Signal Processing 32 (2012) 32–43
33
important to be considered during the design cycle of composite structures. Various theories have been developed in an effort to predict numerically the effect of damage on the material properties of laminates, based upon continuum, micromechanics and macro-mechanics damage principles [10,11]. Mathematically, the uncertainty propagation is handled using several methods classified mainly into two categories, the probabilistic and non-probabilistic theories [12–17]. In probabilistic approaches, the probability density functions of the random input variables are properly defined, while in non-probabilistic methods, the interval of upper and lower bounds of the parameters are used to represent the uncertainty range. Systematic work regarding the influence of structural uncertainties on the behavior of aeroelastic structures has already been undertaken. More specifically, the quantification and tailoring of various aerodynamic problems under uncertainties were investigated [18–23]. Additionally, the robust design of aeroelastic systems subjected to parametric uncertainties has also been proposed [24]. Moreover, the Polynomial Chaos Expansion (PCE) method has been employed to explore the variability of response in control [25], computational fluid dynamics [26,27] and buckling problems [28]. There is a wide range of different parametric uncertainties that affects the unsteady aerodynamic loads of aerospace structures, for instance the mass/inertia distribution of the fuel, the atmospheric properties (temperature and air density) as well as the variations in the composite fiber lay-up (e.g. [18]). The aeroelastic stability and loads design process already includes the analysis of a wide range of different parameters such fuel mass/inertia, center of gravity position, and consequently these were not considered as uncertain variables in the present study. The atmospheric properties are a different matter, and for instance it has been argued that the variations in the nonlinear LCO behavior of F-16 aircraft is due to variations in the atmosphere. However, the authors’ belief is (shown later) that for the examined aeroelastic model the effect of atmospheric variation is much smaller than the damage, and consequently it was not the main thrust of this work. The motivation behind this effort was that research studies have been published on the variability of the aeroelastic stability and response of composite structures due to uncertain fiber lay-up [18,19,24], although without considering the inclusion of damage severity as uncertain variable in combination with stochastic material parameters. In this paper, the Polynomial Chaos Expansion (PCE) method [29–31] is applied to predict the uncertainty propagation on the aeroelastic stability and gust response of a finite element composite wing model. Two sources of uncertainty with Gaussian distribution were incorporated during the modeling process; the longitudinal Young’s modulus and the structural damping of all the components of the wing box along with the crack density of pre-defined layers of the wing structure, are considered as independent input variables. Additionally, the inclusion of fiber breakages at specific locations is investigated separately. Typical output responses, including the flutter speed and the maximum/minimum root shear force due to gust loads are used as metrics for the aeroelastic behavior of the composite wing. Exploiting an efficient Latin Hypercube sampling technique, PCE models of the output responses are built using a small sample of observations and the corresponding probability density functions (PDFs) are evaluated in a systematic way. The derived PDFs are correlated and verified against relevant results obtained by performing adequate number of simulations of the actual models. The paper is structured as follows. The finite element model of the composite wing is presented in Section 2, and the method for modeling damage is outlined in Section 3. In Section 4, the stochastic expansion approach is described, followed by numerical results for the deterministic and probabilistic aeroelastic models without and with damage (5.2 and 5.3) in Section 5. Finally, in Section 6, the conclusions from this work are drawn. 2. Composite wing modeling The geometrical configuration of the composite wing used in this study is shown in Fig. 1(a). The semi-span of the wing was 5 m, while the root and tip chord were 2 m and 1 m, respectively. Additionally, the leading and trailing edge of the wing structure was swept backwards 301 and 20.71, and the thickness-to-chord ratio was four percent. A structural wing box, consisting of top and bottom skins, front and rear spars and ribs (see Fig. 1(b)), was inserted into the wing and positioned such that the front spar lies at the quarter chord and the rear spar lies at the three quarter chord. Ribs were distributed along the span-wise direction of the wing in order to increase the chord-wise stiffness and suppress local panel deflections.
Fig. 1. (a) Geometrical and (b) structural definition of the composite wing.
34
G. Georgiou et al. / Mechanical Systems and Signal Processing 32 (2012) 32–43
Table 1 Material properties for the composite layers. Material properties E1 (GPa) E2 (GPa)
n12 G12 (GPa) Ply thickness (mm) Density (kg/m3)
136.52 10.0 0.3 4.83 0.125 1570.0
Table 2 Lay-up scheme for the composite wing.
Front spar Rear spar Top skin Bottom skin Ribs
Lay-up scheme
Total thickness (mm)
(( 451)4)s (( 451)4)s (( 451)2, (451)3,90)s (( 451)2, (451)3,90)s (( 451)2)s
1.00 1.00 1.50 1.50 0.50
The finite element model of the wing box was developed using shell elements. A graphite/epoxy laminate was assigned at the wing box using a two-dimensional orthotropic material model and the corresponding material properties are summarized in Table 1 [2]. The modeling of the orthotropic material was performed by defining a material coordinate system for each component (skins, spars and ribs) of the wing box (e.g. system xyz for the top skin in Fig. 1(b)). The ply angle of each layer was then defined with respect to the x-axis of the local coordinate system. A symmetric configuration of the laminate was selected for all the components of the wing box and the fiber orientations in the stacking direction with respect to the relevant material coordinate system are summarized in Table 2. 3. Damage modeling Unidirectional fiber composites are characterized by high material stiffness and strength in the fiber direction (longitudinal) and low mechanical properties in the matrix direction (transverse). Under monotonic or cyclic loading, the strength difference between the fibers and the matrix is usually related to damage mechanisms. The damage process is initiated from defects present in a ply, such as small fiber–resin debonds, resin-rich regions, resin voids, high fiber volume fraction areas etc., leading to the development of transverse and longitudinal matrix cracks, inter-laminar delamination and fiber fracture (Fig. 2(a)) [9]. The transverse matrix cracking, which starts at low loads and propagates in the fiber direction, is the predominant damage mode related to a significant reduction of the effective stiffness and strength for composite laminates [10]. At high loads, the matrix cracking also triggers the growth of resin-dominated damage modes, such as delamination and fiber breakage. A number of theories have been developed in an effort to predict the effect of damage on the material properties of laminates, based upon continuum, micro-mechanics and macro-mechanics damage principles [10,11]. Briefly, in the continuum mechanics approach, the damage parameters such as matrix cracks and inter-laminar debonding are characterized by a set of internal state variables, which physically represents damage magnitude and modes. According to the micro-mechanics theories, a detailed micro-field of stress–strain in the vicinity of the defect or crack is evaluated and the effective stiffness of the material is predicted. Whereas, in macro-mechanics methods, such as two-phase or threephase models, a composite layer is considered as an orthotropic homogeneous solid with cracks and the degraded material properties are calculated using analytical expressions. In the present study, the self-consistent method and particularly the two-phase model was incorporated in order to evaluate the elastic constants of unidirectional composite lamina with transverse matrix cracks. According to the twophase model, when the fiber diameter is much smaller than the crack length or the ply thickness, the fibers and the matrix in a fiber-reinforced material can be replaced by an effective homogeneous medium which contains aligned cracks of length 2a (Fig. 2(b)). Consequently, the stiffness and compliance changes of composite laminates due to the existence of open cracks can be predicted using analytical expressions described in the literature [10]. More specifically, the degraded material properties are calculated with respect to a metric which characterize the severity of damage in a cracked layer. The damage metric b, called crack density, is defined as the average distance between cracks and under the assumption that the cracks are equally spaced in a ply, the distance between successive cracks is 2a/b (Fig. 2(b)). Therefore, if b ¼0 the composite layer is uncracked, while if b ¼1 there is one crack in a square of side 2a. As a result, the in-plane elastic orthotropic properties of a unidirectional lamina, such as longitudinal Young’s modulus (fiber direction), transverse Young’s modulus (matrix direction) and the shear modulus can be efficiently predicted as a
G. Georgiou et al. / Mechanical Systems and Signal Processing 32 (2012) 32–43
35
Fig. 2. (a) Damage modes and (b) aligned slit cracks in an infinite fibrous media and a fiber lamina.
Fig. 3. Reduced material properties with respect to the crack density: (a) transverse Young’s modulus and (b) shear modulus.
function of crack density (Fig. 3). However, it is worth mentioning that according to the self-consistent method for the two-phase model, the longitudinal stiffness does not degrade due to the transverse matrix cracks. Additionally, it becomes apparent from Fig. 3 that the reduction rate for transverse Young’s modulus is higher than the corresponding rate for the shear modulus. On the other hand, in order to resemble the fiber fracture numerically, the elastic constants of the material were assumed to be close to zero in the region of damage. Typical locations on a structural wing box where high stress concentrations are usually observed, were selected to apply the degraded material properties due to the initiation of transverse matrix cracks or the existence of fiber breakages. Thus, the investigated damage modes were induced at the middle layers of the stacking sequence for the leading spar and at the outer layers for the bottom skin (Fig. 4), in an effort to examine the damage effect on the aeroelastic stability and the gust response of the composite wing model. The described damage modeling does not include any increase in damping as it was considered that the reduction in stiffness has the most significant effect on the aeroelastic response of the composite wing. Also, the damping behavior cannot be modeled accurately without first performing vibration testing and there is not going to be a simple relationship between the damage size and the damping level. Finally, any resulting increase in damping (it cannot decrease) is likely to be small and will be beneficial in terms of both the gust response and flutter speed. 4. Stochastic expansion Stochastic expansion techniques, such as the Polynomial Chaos Expansion (PCE) and Stochastic Collocation (SC), expand the output response in a series of random variables for uncertainty quantification and propagation at reduced computational effort. The development of the stochastic expansions started when Wiener [29] introduced a mathematical
36
G. Georgiou et al. / Mechanical Systems and Signal Processing 32 (2012) 32–43
Fig. 4. Damage locations for transverse matrix cracks or fiber breakages.
model, in an effort to describe the irregularities appeared in Brownian motion using a multiple stochastic integral with homogeneous chaos. Thereafter, Ito [30] modified Weiner’s work and showed that any stochastic process can be described as a Weiner process. Ghanem and Spanos [31] stated a simple definition of the Polynomial Chaos Expansion as a convergent series of the form u ¼ a0 H 0 þ
1 X
ai1 H1 ðxi1 Þ þ
i1 1 X X i1 ¼ 1 i2 ¼ 1
i¼1
i1 i2 1 X X X
ai1 i2 H2 ðxi1 , xi2 Þ þ
ai1 i2 i3 H3 ðxi1 , xi2 , xi3 Þ þ ,
ð1Þ
i1 ¼ 1 i2 ¼ 1 i3 ¼ 1
1
where fxi gi ¼ 1 represents a set of independent stochastic Gaussian variables, Hp ðxi1 ,. . ., xip Þ is a set of multidimensional Hermite polynomials of order p and ai1 ,. . .,aip are properly defined coefficients. The approximation of the output response u can be rewritten using different orthogonal polynomials which belong to the Askey scheme, such as the Laguerre, Jacobi and Legendre if the input random variables follow the gamma, beta and uniform distribution, respectively. 1 In the present study, the uncertain parameters fxi gi ¼ 1 were defined as normalized variables using the following expression:
xi ¼
xi mi
si
,
ð2Þ
where mi and s2i represent the mean value and the variance of the i-th random variable. The Hermite polynomials with respect to the independent variable x are defined by H0 ðxÞ ¼ 1, H1 ðxÞ ¼ x, 2
H2 ðxÞ ¼ x 1, ð3Þ
3
H3 ðxÞ ¼ x 3x, 4
2
H4 ðxÞ ¼ x 6x þ 3, ^ ^ As a result, in the case of a single input parameter, the Polynomial Chaos Expansion analytical formula truncated at the fourth order is written as 2
3
4
2
u ¼ a0 þa1 x þ a2 ðx 1Þ þ a3 ðx 3xÞ þa4 ðx 6x þ 3Þ, whereas, in the case of two uncertain variables, x1 and x2 the second order stochastic expansion is obtained as 2
2
u ¼ a0 þa1 x1 þa2 x2 þa3 ðx1 1Þ þ a4 x1 x2 þ a5 ðx2 1Þ: The unknown coefficients ai of the Polynomial Chaos Expansions can be calculated using a regression analysis or the statistically averaging method. Here, a linear regression process was used and the polynomial coefficients were predicted by solving the following system of equations for the case of a single uncertain variable: 2 32 3 H01 H11 Hp1 a0 2 3 u1 6H 76 a 7 H H 6 7 6 02 12 p2 17 6^ 7 6 76 7, ð4Þ 4 5¼6 76 ^ 7 ^ ^ ^ ^ 4 54 5 un H0n H1n Hpn ap where the term ui for i¼1,y,n is a discrete observation of the output response, while the element Hpi represents the value of the p-th order Hermite polynomial calculated for the i-th value of the input parameter. It becomes apparent that the
G. Georgiou et al. / Mechanical Systems and Signal Processing 32 (2012) 32–43
37
accuracy of the regression coefficients depends upon the quantity and the quality of the obtained observations and for this reason, an efficient sampling technique, such as the Latin Hypercube method, which ensures that the points of the input variables are selected with equal probability, was implemented prior to the regression analysis. After the evaluation of the unknown coefficients, the Polynomial Chaos Expansion models can be emulated for given values of the uncertain parameters and predict statistical properties of the output response, such as the mean value, the standard deviation and the probability density function at low computational cost. The Polynomial Chaos Expansion method was selected in the current study as a very efficient non-intrusive method used to produce probabilistic models for the output responses of complicated structures with many variables. Additionally, an important advantage of the Polynomial Chaos Expansion method is the simplistic formulation of the output responses in a polynomial form, making the technique suitable in a range of different applications (e.g. [18, 25–28]).
5. Numerical results In this section, representative numerical results on the aeroelastic stability and gust response of the composite wing model incorporating deterministic or probabilistic material and/or damage parameters are presented. Initially, the baseline aeroelastic model of the composite wing with no variation in material properties and without the inclusion of damage was investigated in a sequence of unsteady aerodynamic analyses, including flutter analysis and dynamic aeroelastic response analysis. The purpose of this baseline study was to identify the dynamic instabilities (flutter mechanisms) and the gust response history of typical ‘‘Interesting Quantities’’, such as the wing root shear force. Following this study, the uncertain nature of the mechanical properties was taken into account in a probabilistic manner; longitudinal Young’s modulus and the structural damping were varied in a stochastic manner. The probability density functions for the flutter speed, and also the maximum and minimum root shear force were calculated based upon initial Latin Hypercube samples and applying the Polynomial Chaos Expansion method. Finally, the effect of damage due to transverse matrix cracks and fiber breakage was introduced into the probabilistic aeroelastic model of the wing, and the severity of the damage on the aeroelastic stability and gust response was studied.
5.1. Deterministic aeroelastic models 5.1.1. Flutter analysis The ‘‘p-k’’ method embedded in MSC.NASTRAN software and in particular the PKNL variant which enables consistent ‘‘matched’’ air density, Mach number and air speeds, was applied to determine the frequency and damping ratio trends of the composite wing at subcritical speeds (Fig. 5) [32]. Inspection of the frequency and damping ratio curves in Fig. 5 shows that the first bending mode (Mode 1) enters a stable non-oscillatory motion after 65 m/s. Additionally, the frequency values of the first torsion mode (Mode 5) exhibit a significant reduction with respect to the air speed, and approximately at air speed 114 m/s the frequency curve intersects the corresponding curve of the third bending mode of the structure (Mode 4), without although affecting the flutter stability. However, Mode 4 starts to interact with the first sway mode (Mode 3) at around 120 m/s, leading to a flutter instability of Mode 3 at an air speed of 167.53 m/s.
Fig. 5. Frequency and damping ratio trends versus air speed for the composite wing model.
38
G. Georgiou et al. / Mechanical Systems and Signal Processing 32 (2012) 32–43
5.1.2. Gust analysis The unsteady aerodynamic behavior of the composite wing was also investigated in a dynamic aeroelastic analysis which included a discrete ‘‘1-cosine’’ vertical gust. The maximum gust amplitude wg0 was calculated using the following expression [33]: 1=6 H ð5Þ wg0 ¼ wref 106:17 where H represents the gust gradient (half the gust wavelength Lg), while the reference gust velocity wref reduces linearly from 17.07 m/s EAS at sea level to 13.41 m/s EAS at 4752 m, and then again to 6.36 m/s EAS at 18288 m. Correspondingly, in the time domain the gust velocity is defined by 8 < wg0 1cos 2pV t if 0 r t r t g 2 Lg ð6Þ wg ðtÞ ¼ :0 if t Z t g where V is the air speed, while tg ¼Lg/V. Three different wavelengths Lg ¼[18 m, 116 m, 214 m](Fig. 6), at sea level with air speed V ¼174.85 m/s were selected as flight conditions in order to evaluate the time history of the root shear force developed at the bottom and at the top skin (Fig. 7). 5.2. Probabilistic aeroelastic models without damage In an attempt to predict the behavior of the composite wing under uncertainty, probabilistic aeroelastic models of the wing box were built and simulated. Longitudinal Young’s modulus and the structural damping of all the components of the wing box were considered as Gaussian parameters and the aeroelastic stability and gust response of the wing were evaluated. The uncertain longitudinal Young’s modulus was assumed to have a mean value 136.52 GPa and coefficient of variance 0.03, whereas the structural damping was chosen to have mean value 0.02 and coefficient of variance 0.2.
Fig. 6. Gust velocity for different wavelengths.
Fig. 7. Time histories of the root shear force for different wavelengths obtained: (a) at the bottom skin and (b) at the top skin.
G. Georgiou et al. / Mechanical Systems and Signal Processing 32 (2012) 32–43
39
The Polynomial Chaos Expansion method [29–31] was implemented in a non-intrusive way in order to predict the uncertainty propagation for the flutter speed, the maximum and minimum root shear force of the composite wing. The PCE models of the aforementioned quantities were constructed by exploiting the statistical properties of a small sample of observations (size 50), which was obtained using the Latin Hypercube sampling technique. A large number of emulations (100,000) for the PCE models were conducted at low computational cost to produce the probability density functions for the flutter speed, the maximum and the minimum root shear force. The PDFs obtained from the PCE models were compared to those found from 5000 simulations of the actual models (a Latin Hypercube sample) and an excellent agreement was achieved. 5.2.1. Flutter analysis Several case studies were conducted in an attempt to determine the variability on the aeroelastic stability of the composite wing with respect to uncertain material properties applied to all the components of the wing box. The probability density function of the flutter speed was predicted using the Polynomial Chaos Expansion method described above, for variations in longitudinal Young’s modulus and the structural damping. The obtained PDFs are presented in Fig. 8(a) and (b), and it becomes apparent that the inclusion of the structural damping as independent variable has little effect upon the shape and the range of the derived distribution. Direct comparison of the PDFs obtained by simulation of the actual models and the corresponding approximated PDFs revealed the accuracy of developed PCE models. 5.2.2. Gust analysis In an effort to predict the uncertainty propagation of material properties on the gust response of the composite wing, several representative scenarios were performed and a typical ‘‘Interesting Quantity’’ of the wing structure was evaluated. The probability density function of the maximum and minimum root shear force was predicted using the Polynomial Chaos Expansion method, selecting as uncertain variables the longitudinal Young’s modulus and the structural damping. The derived PDFs are included in Fig. 9(a) and (b), and it is worth to mention that the shape as well as the dispersion of the values for the maximum and minimum root shear force presented significant differences. However, the PDFs approximations of the output responses were verified by the comparison with simulations of the actual models performed using the Latin Hypercube sampling technique and a very good agreement was achieved. 5.3. Probabilistic aeroelastic models with damage The effect of damage on the aeroelastic stability and response of the composite wing box was examined by inducing two damage modes, the transverse matrix cracking as well as the fiber breakage in a separate sequence at the outer layers of the bottom skin and at the middle layers of the leading spar. The severity of damage due to the transverse matrix cracking was investigated, considering the crack density as an uncertain Gaussian variable with mean value 0.3 and coefficient of variance 0.25. According to the self-consistent method, only the stiffness and compliance of composite laminates are changed due to the existence of open cracks, while the structural damping remains unaffected. However, previous studies (e.g. [34]) show that the damping in composite materials does tend to increase with the severity of damage. Thus, a combined case study was conducted in which the crack density and material properties, including longitudinal Young’s modulus and the structural damping were varied in a stochastic manner with mean value 136.52 GPa and coefficient of variance 0.03, and mean value 0.02 and coefficient of variance 0.2, respectively. Accurate predictions of the flutter speed and the maximum/minimum root shear force were achieved using the Polynomial Chaos Expansion method. A regression method was implemented in order to calculate the coefficients of the Hermite polynomials, using a small sample of observations (size 50) obtained by applying the Latin Hypercube sampling
Fig. 8. Probability density functions of the flutter speed: (a) for uncertain longitudinal Young’s modulus and (b) for uncertain longitudinal Young’s modulus and structural damping.
40
G. Georgiou et al. / Mechanical Systems and Signal Processing 32 (2012) 32–43
Fig. 9. Probability density functions of the wing root shear force for uncertain longitudinal Young’s modulus and structural damping: (a) maximum case and (b) minimum case.
Fig. 10. Probability density function of the flutter speed: (a) for uncertain crack density, (b) for uncertain longitudinal Young’s modulus, structural damping and crack density and (c) for uncertain longitudinal Young’s modulus, structural damping and inclusion of fiber breakages.
technique. The uncertainty quantification of the aeroelastic behavior for the composite wing was accomplished by performing 100,000 emulations of the developed PCE models. The statistical properties of the PCE models were verified against the relevant characteristics of the actual models by conducting 5000 simulations and through the evaluation of the probability density functions for the flutter speed, the maximum and minimum root shear force. 5.3.1. Flutter analysis Three discrete case scenarios were conducted in an effort to predict the uncertainty propagation of the damage metric (crack density) on the aeroelastic stability of the composite wing taking into account the variability of the crack density (Fig. 10(a)) as well as the Gaussian distribution of the longitudinal Young’s modulus, the structural damping and the crack density (Fig. 10(b)). It is worth to mention that the PDF of the flutter speed due to the uncertainty of the crack density presented a skewed distribution with small dispersion of the values. Consequently, the variability of the output response
G. Georgiou et al. / Mechanical Systems and Signal Processing 32 (2012) 32–43
41
Fig. 11. Probability density functions of the wing root shear force for uncertain crack density: (a) maximum case and (b) minimum case.
in the case where longitudinal Young’s modulus and the structural damping were considered as additional random variables (Fig. 10(b)) did not present significant differences in comparison with the corresponding results included in Fig. 8(b). A separate scenario was performed at which the longitudinal Young’s modulus and the structural damping were selected as uncertain parameters along with the inclusion of fiber breakages in pre-defined composite layers (Fig. 10(c)). The relevant shape of the derived PDF followed the normal distribution, although the mean value of the flutter speed was significantly reduced. In all the examined cases, the developed PCE models provided accurate PDFs predictions in comparison with the corresponding results obtained through the simulation of the actual models. Up to the edge of the stratosphere (11 km), the temperature decreases linearly (6.5 1C per km) and this has a major effect on the aerodynamic forces at a given altitude. The International Standard Atmosphere sets the sea level temperature at 15 1C, so it is of interest to note the effect of varying the effective sea level temperature between the possible extremes that might be encountered. Fig. A1 in Appendix 1 shows how the flutter speed varies over a range of about 4 m/s between the extreme ranges of 30 1C and þ40 1C. This range is much larger that might be expected in reality, so this variation in flutter speed can be considered to be of less consequence than the structural damage. This is an area that requires further detailed examination. 5.3.2. Gust analysis The severity of damage on the gust response of the composite wing was predicted through various case studies in which the probability density function of the maximum and minimum root shear force was evaluated using the Polynomial Chaos Expansion method. Initially, the crack density (Fig. 11(a) and (b)) was considered as the uncertain variable. It is worth to mention that the obtained distributions displayed significant dispersion of the values as well as negative and positive skewness for the maximum and minimum root shear force, respectively. This investigation was extended by incorporating as additional stochastic parameters the longitudinal Young’s modulus and the structural damping. The derived PDFs are presented in Fig. 12(a) and (b) in superposition with the corresponding results gathered for the uncracked material case (Fig. 9(a) and (b)). Direct comparison of the PDFs for the two discrete cases (uncracked and cracked material) revealed that the distribution for the maximum root shear force of the composite wing with cracked material was shifted towards higher values and the shape was changed considerably. However, the shape and dispersion of the values for the distribution of the minimum root shear force for the uncracked and cracked material cases did not present significant differences and the relevant PDF was only shifted towards lower values. Finally, similar behavior was observed in the scenario in which fiber breakages were included in the damage modeling of the composite layers (Fig. 12(c)). Particularly, the obtained PDF of the maximum root shear force was shifted towards higher values in comparison with the undamaged material case and the distribution presented a wider dispersion of the values. 6. Conclusions The Polynomial Chaos Expansion (PCE) method has been applied to predict the uncertainty propagation on the aeroelastic stability and gust response of a composite wing, due to uncertainties related to the material properties and the severity of damage, in an accurate and computationally efficient way. Initially, the material properties of the wing structure were used as stochastic variables with Gaussian distribution and characteristic output responses, such as the flutter speed, the maximum and minimum root shear force were considered as metrics for the aeroelastic behavior of the composite wing. The statistical properties of the output responses, in the form of probability density functions (PDF), were accurately determined by implementing PCE models, built from a small sample of observations. The PDF of the flutter speed was predicted using the PCE method, for variations in longitudinal Young’s modulus and the structural damping, where it became apparent that the inclusion of the structural damping as independent variable had little effect upon the
42
G. Georgiou et al. / Mechanical Systems and Signal Processing 32 (2012) 32–43
Fig. 12. Probability density functions of the wing root shear force for uncertain longitudinal Young’s modulus, structural damping and crack density: (a) maximum case, (b) minimum case and (c) maximum case with inclusion of fiber breakages.
shape and the range of the derived distribution. Concurrently, the PDFs of the maximum and minimum root shear force were predicted using the PCE method and it was found that the shape as well as the dispersion of the values for the interesting quantities presented significant variability due to the manufacturing uncertainties. The severity of damage in pre-defined locations of the composite layers was also considered as stochastic variable with Gaussian distribution and the effect of damage due to transverse matrix cracking was incorporated into the laminates using the self-consistent method. The extreme case of fiber fracture was also investigated. The PDF of the flutter speed was predicted using the PCE method, for variations in the material properties and taking into account the damage severity. The results did not present significant differences in comparison with the corresponding findings for the uncracked material, due mainly to the assumption for the damage locations. On the other hand, the mean value of the flutter speed was reduced when the fiber breakages were included in pre-defined composite layers. Different trends were observed for the maximum and minimum root shear force after the superposition of the results gathered for the uncracked and cracked material cases, where the shape and dispersion of the values for the distribution of the interesting quantities presented significant differences. All the derived PDFs were correlated against corresponding results obtained by performing an adequate number of simulations of the actual models, providing an excellent agreement and the future possibility to apply the developed methodology in order to efficiently investigate the effect of damage severity in robust design of composites wings.
Acknowledgments This work is supported through the Virtual Engineering Centre (VEC), which is a University of Liverpool initiative in partnership with the Northwest Aerospace Alliance, the Science and Technology Facilities Council (Daresbury Laboratory), BAE Systems, Morson Projects and Airbus (UK). The VEC is funded by the Northwest Regional Development Agency (NWDA) and European Regional Development Fund (ERDF) to provide a focal point for virtual engineering research, education and skills development, best practice demonstration, and knowledge transfer to the aerospace sector. Appendix 1 The effect of sea level temperature on the flutter speed for the composite wing model is depicted in Fig. A1.
G. Georgiou et al. / Mechanical Systems and Signal Processing 32 (2012) 32–43
43
Fig. A1. Flutter speed trend versus air temperature for the composite wing model.
References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28] [29] [30] [31] [32] [33] [34]
T.A. Weisshaar, Aeroelastic tailoring of forward swept composite wings, J. Aircr. 18 (8) (1981) 669–676. F.E. Eastep, V.A. Tischler, V.B. Venkayya, N.S. Khot, Aeroelastic tailoring of composite structures, J. Aircr. 36 (6) (1999) 1041–1047. T.A. Weisshaar, D.K. Duke, Induced drag reduction using aeroelastic tailoring with adaptive control surfaces, J. Aircr. 43 (1) (2006) 157–164. H. Arizono, K. Isogai, Application of genetic algorithm for aeroelastic tailoring of a cranked-arrow wing, J. Aircr. 42 (2) (2005) 493–499. M. Kameyama, H. Fukunaga, Optimum design of composite plate wings for aeroelastic characteristics using lamination parameters, Comput. Struct. 85 (2007) 213–224. S. Guo, Aeroelastic optimization of an aerobatic aircraft wing structure, Aerosp. Sci. Technol. 11 (2007) 396–404. C.L. Pettit, R.V. Grandhi, Optimization of a wing structure for gust response and aileron effectiveness, J. Aircr. 40 (6) (2003) 1185–1191. T.U. Kim, I.H. Hwang., Optimal design of composite wing subjected to gust loads, Comput. Struct. 83 (2005) 1546–1554. J.-M. Berthelot, Transverse cracking and delamination in cross-ply glass-fiber and carbon-fiber reinforced plastic laminates: static and fatigue loading, Appl. Mech. Rev. 56 (2003) 111–147. G.J. Dvora, N. Laws, M. Hejazi, Analysis of progressive matrix cracking in composite laminates. I. Thermoelastic properties of a ply with cracks, J. Compos. Mater. 19 (1985) 216–234. R. Talreja, Stiffness properties of composite laminates with matrix cracking and interior delamination, Eng. Fract. Mech. 25 (1986) 751–762. P. Soundappan, E. Nikolaidis, R.T. Haftka, R.V. Grandhi, R.A. Canfield, Comparison of evidence theory and Bayesian theory for uncertainty modeling, Reliab. Eng. Syst. Saf. 85 (2004) 295–311. K.L. Wood, K.N. Otto, E.K. Antonsson, Engineering design calculations with fuzzy parameters, Fuzzy Sets Syst. 52 (1992) 1–20. H. Bae, R.V. Grandhi, R.A. Canfield, Uncertainty quantification of structural response using evidence theory, AIAA J. 41 (10) (2003) 2026–2068. A. Hunter, S. Parsons, Applications of Uncertainty Formalisms, Springer, New York, 1998, pp. 8–32. G. Shafer, A Mathematical Theory of Evidence, Princeton University Press, Princeton, NJ, 1976. Y. Ben-Haim, Information-Gap Decision Theory: Decisions under Severe Uncertainty, Academic Press, 2001. A. Manan, J.E. Cooper, Design of composite wings including uncertainties: a probabilistic approach, J. Aircr. 46 (2) (2009) 601–607. A. Manan, J.E. Cooper, Robust Design of Composite Wings for Gust Response, Proceedings of the 50th AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics, and Materials Conference, Palm Springs, California, May 4–7, 2009. C.L. Pettit, Uncertainty quantification in aeroelasticity: recent results and research challenges, J. Aircr. 41 (5) (2004) 1217–1229. P.S. Beran, C.L. Pettit, D.R. Millman, Uncertainty quantification of limit-cycle oscillations, J. Comput. Phys. 217 (2006) 217–247. J. Kuttenkeuler, U. Ringertz, Aeroelastic tailoring considering uncertainties in material properties, Struct. Optim. 15 (1998) 157–162. S.C. Catravete, R.A. Ibrahim, Effect of stiffness uncertainties on the flutter of cantilever wing, AIAA J. 46 (4) (2008) 925–935. L.G. Crespo, D.P. Giesy, S.P. Kenny, Robustness analysis and robust design of uncertain systems, AIAA J. 46 (2) (2008) 388–396. F.S. Hover, M.S. Triantafyllou, Application of polynomial chaos in stability and control, Automatica 42 (2006) 789–795. L. Mathellin, M.Y. Hussaini, T.A. Zang, Stochastic approaches to uncertainty quantification in CFD simulations, Numer. Algorithms 38 (2005) 209–236. O.M. Knioa, O.P. Le Maˆıtre, Uncertainty propagation in CFD using polynomial chaos decomposition, J. Fluid Dyn. Res. 38 (2006) 616–640. S.K. Choi, R.A. Grandhi, R.A. Canfield, C.L. Pettit, Polynomial chaos expansion with Latin hypercube sampling for estimating response variability, AIAA J. 42 (6) (2004) 1191–1198. N. Wiener, The homogeneous chaos, Am. J. Math. 60 (4) (1938) 897–936. K. Ito, Multiple Wiener integrals, J. Math. Soc., Jpn 3 (1951) 157–169. R.G. Ghanem, P.D. Spanos, Stochastic Finite Elements: A Spectral Approach, Springer-Verlag Inc, New York, 1991. W.P. Rodden, E.H. Johnson, MSC.Nastran Aeroelastic Analysis User’s Guide v68, The MacNeal-Schwendler Corporation, 1994. J.R. Wright, J.E. Cooper, Introduction to Aircraft Aeroelasticity and Loads, John Wiley and Sons Ltd, 2007. Z. Zhang, G. Hartwig, Relation of damping and fatigue damage of unidirectional fibre composites, Int. J. Fatigue 24 (2002) 713–718.
Contents lists available at SciVerse ScienceDirect
Mechanical Systems and Signal Processing journal homepage: www.elsevier.com/locate/ymssp
Modeling composite wing aeroelastic behavior with uncertain damage severity and material properties G. Georgiou a, A. Manan a, J.E. Cooper b,n a b
School of Engineering, University of Liverpool, Liverpool, L69 3GH, UK Department of Aerospace Engineering, University of Bristol, Bristol, BS8 1TH, UK
a r t i c l e in f o
abstract
Article history: Received 20 June 2011 Received in revised form 1 May 2012 Accepted 7 May 2012 Available online 16 June 2012
The effect of uncertain material properties and severity of damage on the aeroelastic behavior of a finite element composite wing model are predicted by applying the Polynomial Chaos Expansion method (PCE). Different damage modes, including the transverse matrix cracking and broken fibers, are induced into pre-defined locations in the laminates and the aeroelastic stability and dynamic response of the wing due to ‘‘1-cosine’’ vertical gusts are evaluated. For this purpose, PCE models that predict the variation due to uncertainty of the flutter speed and an ‘‘Interesting Quantity’’ (root shear force) of the wing box are developed based upon a small sample of observations, exploiting the efficient Latin Hypercube sampling technique. The uncertainty propagation on the output responses, in the form of probability density functions, is evaluated at low computational cost, implementing the PCE models and verified successfully against the actual results. & 2012 Elsevier Ltd. All rights reserved.
Keywords: Composite wing Damage severity Uncertainty propagation Aeroelastic behavior
1. Introduction Composite materials are being used increasingly in aerospace structures, primarily due to their attractive strength– weight ratios. A further advantage is the inherent anisotropic behavior which can be exploited to induce beneficial bending–torsion couplings in aircraft wings during flight, thus providing a greater design flexibility compared with metallic structures. A number of aeroelastic tailoring applications have sought to use composite materials for a range of applications, including the increase of divergence speed, the minimization of drag, the reduction of gust loads, etc.[1–8]. However, the complex manufacturing process of fibrous laminates induces an increased number of possible uncertain geometrical (thickness) and material (in-plane moduli, fiber orientation, etc.) parameters, which affect in a considerable degree the strength and the overall performance of composite structures. Additionally, under monotonic or cyclic loading the strength difference between the fibers and the matrix is usually related to damage mechanisms. The damage process is initiated from defects randomly distributed in a ply, leading to the development of transverse and longitudinal matrix cracks, inter-laminar delamination and fiber fracture [9]. The transverse matrix cracking, which starts at low loads and propagates in the fiber direction, is the predominant damage mode related to a significant reduction of the effective stiffness and strength for composite laminates [10]. Taking into account that the initiation and propagation of transverse matrix cracks is a stochastic procedure, the severity of damage in composite materials represents a considerable source of uncertainty. Moreover, at high loads, the matrix cracking triggers also the growth of resin-dominated damage modes, such as delamination and fiber breakage, which is
n
Corresponding author. Tel.: þ 44 117 954 5388; fax: þ 44 117 927 2771. E-mail addresses: [email protected] (G. Georgiou), [email protected] (J.E. Cooper).
0888-3270/$ - see front matter & 2012 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.ymssp.2012.05.003
G. Georgiou et al. / Mechanical Systems and Signal Processing 32 (2012) 32–43
33
important to be considered during the design cycle of composite structures. Various theories have been developed in an effort to predict numerically the effect of damage on the material properties of laminates, based upon continuum, micromechanics and macro-mechanics damage principles [10,11]. Mathematically, the uncertainty propagation is handled using several methods classified mainly into two categories, the probabilistic and non-probabilistic theories [12–17]. In probabilistic approaches, the probability density functions of the random input variables are properly defined, while in non-probabilistic methods, the interval of upper and lower bounds of the parameters are used to represent the uncertainty range. Systematic work regarding the influence of structural uncertainties on the behavior of aeroelastic structures has already been undertaken. More specifically, the quantification and tailoring of various aerodynamic problems under uncertainties were investigated [18–23]. Additionally, the robust design of aeroelastic systems subjected to parametric uncertainties has also been proposed [24]. Moreover, the Polynomial Chaos Expansion (PCE) method has been employed to explore the variability of response in control [25], computational fluid dynamics [26,27] and buckling problems [28]. There is a wide range of different parametric uncertainties that affects the unsteady aerodynamic loads of aerospace structures, for instance the mass/inertia distribution of the fuel, the atmospheric properties (temperature and air density) as well as the variations in the composite fiber lay-up (e.g. [18]). The aeroelastic stability and loads design process already includes the analysis of a wide range of different parameters such fuel mass/inertia, center of gravity position, and consequently these were not considered as uncertain variables in the present study. The atmospheric properties are a different matter, and for instance it has been argued that the variations in the nonlinear LCO behavior of F-16 aircraft is due to variations in the atmosphere. However, the authors’ belief is (shown later) that for the examined aeroelastic model the effect of atmospheric variation is much smaller than the damage, and consequently it was not the main thrust of this work. The motivation behind this effort was that research studies have been published on the variability of the aeroelastic stability and response of composite structures due to uncertain fiber lay-up [18,19,24], although without considering the inclusion of damage severity as uncertain variable in combination with stochastic material parameters. In this paper, the Polynomial Chaos Expansion (PCE) method [29–31] is applied to predict the uncertainty propagation on the aeroelastic stability and gust response of a finite element composite wing model. Two sources of uncertainty with Gaussian distribution were incorporated during the modeling process; the longitudinal Young’s modulus and the structural damping of all the components of the wing box along with the crack density of pre-defined layers of the wing structure, are considered as independent input variables. Additionally, the inclusion of fiber breakages at specific locations is investigated separately. Typical output responses, including the flutter speed and the maximum/minimum root shear force due to gust loads are used as metrics for the aeroelastic behavior of the composite wing. Exploiting an efficient Latin Hypercube sampling technique, PCE models of the output responses are built using a small sample of observations and the corresponding probability density functions (PDFs) are evaluated in a systematic way. The derived PDFs are correlated and verified against relevant results obtained by performing adequate number of simulations of the actual models. The paper is structured as follows. The finite element model of the composite wing is presented in Section 2, and the method for modeling damage is outlined in Section 3. In Section 4, the stochastic expansion approach is described, followed by numerical results for the deterministic and probabilistic aeroelastic models without and with damage (5.2 and 5.3) in Section 5. Finally, in Section 6, the conclusions from this work are drawn. 2. Composite wing modeling The geometrical configuration of the composite wing used in this study is shown in Fig. 1(a). The semi-span of the wing was 5 m, while the root and tip chord were 2 m and 1 m, respectively. Additionally, the leading and trailing edge of the wing structure was swept backwards 301 and 20.71, and the thickness-to-chord ratio was four percent. A structural wing box, consisting of top and bottom skins, front and rear spars and ribs (see Fig. 1(b)), was inserted into the wing and positioned such that the front spar lies at the quarter chord and the rear spar lies at the three quarter chord. Ribs were distributed along the span-wise direction of the wing in order to increase the chord-wise stiffness and suppress local panel deflections.
Fig. 1. (a) Geometrical and (b) structural definition of the composite wing.
34
G. Georgiou et al. / Mechanical Systems and Signal Processing 32 (2012) 32–43
Table 1 Material properties for the composite layers. Material properties E1 (GPa) E2 (GPa)
n12 G12 (GPa) Ply thickness (mm) Density (kg/m3)
136.52 10.0 0.3 4.83 0.125 1570.0
Table 2 Lay-up scheme for the composite wing.
Front spar Rear spar Top skin Bottom skin Ribs
Lay-up scheme
Total thickness (mm)
(( 451)4)s (( 451)4)s (( 451)2, (451)3,90)s (( 451)2, (451)3,90)s (( 451)2)s
1.00 1.00 1.50 1.50 0.50
The finite element model of the wing box was developed using shell elements. A graphite/epoxy laminate was assigned at the wing box using a two-dimensional orthotropic material model and the corresponding material properties are summarized in Table 1 [2]. The modeling of the orthotropic material was performed by defining a material coordinate system for each component (skins, spars and ribs) of the wing box (e.g. system xyz for the top skin in Fig. 1(b)). The ply angle of each layer was then defined with respect to the x-axis of the local coordinate system. A symmetric configuration of the laminate was selected for all the components of the wing box and the fiber orientations in the stacking direction with respect to the relevant material coordinate system are summarized in Table 2. 3. Damage modeling Unidirectional fiber composites are characterized by high material stiffness and strength in the fiber direction (longitudinal) and low mechanical properties in the matrix direction (transverse). Under monotonic or cyclic loading, the strength difference between the fibers and the matrix is usually related to damage mechanisms. The damage process is initiated from defects present in a ply, such as small fiber–resin debonds, resin-rich regions, resin voids, high fiber volume fraction areas etc., leading to the development of transverse and longitudinal matrix cracks, inter-laminar delamination and fiber fracture (Fig. 2(a)) [9]. The transverse matrix cracking, which starts at low loads and propagates in the fiber direction, is the predominant damage mode related to a significant reduction of the effective stiffness and strength for composite laminates [10]. At high loads, the matrix cracking also triggers the growth of resin-dominated damage modes, such as delamination and fiber breakage. A number of theories have been developed in an effort to predict the effect of damage on the material properties of laminates, based upon continuum, micro-mechanics and macro-mechanics damage principles [10,11]. Briefly, in the continuum mechanics approach, the damage parameters such as matrix cracks and inter-laminar debonding are characterized by a set of internal state variables, which physically represents damage magnitude and modes. According to the micro-mechanics theories, a detailed micro-field of stress–strain in the vicinity of the defect or crack is evaluated and the effective stiffness of the material is predicted. Whereas, in macro-mechanics methods, such as two-phase or threephase models, a composite layer is considered as an orthotropic homogeneous solid with cracks and the degraded material properties are calculated using analytical expressions. In the present study, the self-consistent method and particularly the two-phase model was incorporated in order to evaluate the elastic constants of unidirectional composite lamina with transverse matrix cracks. According to the twophase model, when the fiber diameter is much smaller than the crack length or the ply thickness, the fibers and the matrix in a fiber-reinforced material can be replaced by an effective homogeneous medium which contains aligned cracks of length 2a (Fig. 2(b)). Consequently, the stiffness and compliance changes of composite laminates due to the existence of open cracks can be predicted using analytical expressions described in the literature [10]. More specifically, the degraded material properties are calculated with respect to a metric which characterize the severity of damage in a cracked layer. The damage metric b, called crack density, is defined as the average distance between cracks and under the assumption that the cracks are equally spaced in a ply, the distance between successive cracks is 2a/b (Fig. 2(b)). Therefore, if b ¼0 the composite layer is uncracked, while if b ¼1 there is one crack in a square of side 2a. As a result, the in-plane elastic orthotropic properties of a unidirectional lamina, such as longitudinal Young’s modulus (fiber direction), transverse Young’s modulus (matrix direction) and the shear modulus can be efficiently predicted as a
G. Georgiou et al. / Mechanical Systems and Signal Processing 32 (2012) 32–43
35
Fig. 2. (a) Damage modes and (b) aligned slit cracks in an infinite fibrous media and a fiber lamina.
Fig. 3. Reduced material properties with respect to the crack density: (a) transverse Young’s modulus and (b) shear modulus.
function of crack density (Fig. 3). However, it is worth mentioning that according to the self-consistent method for the two-phase model, the longitudinal stiffness does not degrade due to the transverse matrix cracks. Additionally, it becomes apparent from Fig. 3 that the reduction rate for transverse Young’s modulus is higher than the corresponding rate for the shear modulus. On the other hand, in order to resemble the fiber fracture numerically, the elastic constants of the material were assumed to be close to zero in the region of damage. Typical locations on a structural wing box where high stress concentrations are usually observed, were selected to apply the degraded material properties due to the initiation of transverse matrix cracks or the existence of fiber breakages. Thus, the investigated damage modes were induced at the middle layers of the stacking sequence for the leading spar and at the outer layers for the bottom skin (Fig. 4), in an effort to examine the damage effect on the aeroelastic stability and the gust response of the composite wing model. The described damage modeling does not include any increase in damping as it was considered that the reduction in stiffness has the most significant effect on the aeroelastic response of the composite wing. Also, the damping behavior cannot be modeled accurately without first performing vibration testing and there is not going to be a simple relationship between the damage size and the damping level. Finally, any resulting increase in damping (it cannot decrease) is likely to be small and will be beneficial in terms of both the gust response and flutter speed. 4. Stochastic expansion Stochastic expansion techniques, such as the Polynomial Chaos Expansion (PCE) and Stochastic Collocation (SC), expand the output response in a series of random variables for uncertainty quantification and propagation at reduced computational effort. The development of the stochastic expansions started when Wiener [29] introduced a mathematical
36
G. Georgiou et al. / Mechanical Systems and Signal Processing 32 (2012) 32–43
Fig. 4. Damage locations for transverse matrix cracks or fiber breakages.
model, in an effort to describe the irregularities appeared in Brownian motion using a multiple stochastic integral with homogeneous chaos. Thereafter, Ito [30] modified Weiner’s work and showed that any stochastic process can be described as a Weiner process. Ghanem and Spanos [31] stated a simple definition of the Polynomial Chaos Expansion as a convergent series of the form u ¼ a0 H 0 þ
1 X
ai1 H1 ðxi1 Þ þ
i1 1 X X i1 ¼ 1 i2 ¼ 1
i¼1
i1 i2 1 X X X
ai1 i2 H2 ðxi1 , xi2 Þ þ
ai1 i2 i3 H3 ðxi1 , xi2 , xi3 Þ þ ,
ð1Þ
i1 ¼ 1 i2 ¼ 1 i3 ¼ 1
1
where fxi gi ¼ 1 represents a set of independent stochastic Gaussian variables, Hp ðxi1 ,. . ., xip Þ is a set of multidimensional Hermite polynomials of order p and ai1 ,. . .,aip are properly defined coefficients. The approximation of the output response u can be rewritten using different orthogonal polynomials which belong to the Askey scheme, such as the Laguerre, Jacobi and Legendre if the input random variables follow the gamma, beta and uniform distribution, respectively. 1 In the present study, the uncertain parameters fxi gi ¼ 1 were defined as normalized variables using the following expression:
xi ¼
xi mi
si
,
ð2Þ
where mi and s2i represent the mean value and the variance of the i-th random variable. The Hermite polynomials with respect to the independent variable x are defined by H0 ðxÞ ¼ 1, H1 ðxÞ ¼ x, 2
H2 ðxÞ ¼ x 1, ð3Þ
3
H3 ðxÞ ¼ x 3x, 4
2
H4 ðxÞ ¼ x 6x þ 3, ^ ^ As a result, in the case of a single input parameter, the Polynomial Chaos Expansion analytical formula truncated at the fourth order is written as 2
3
4
2
u ¼ a0 þa1 x þ a2 ðx 1Þ þ a3 ðx 3xÞ þa4 ðx 6x þ 3Þ, whereas, in the case of two uncertain variables, x1 and x2 the second order stochastic expansion is obtained as 2
2
u ¼ a0 þa1 x1 þa2 x2 þa3 ðx1 1Þ þ a4 x1 x2 þ a5 ðx2 1Þ: The unknown coefficients ai of the Polynomial Chaos Expansions can be calculated using a regression analysis or the statistically averaging method. Here, a linear regression process was used and the polynomial coefficients were predicted by solving the following system of equations for the case of a single uncertain variable: 2 32 3 H01 H11 Hp1 a0 2 3 u1 6H 76 a 7 H H 6 7 6 02 12 p2 17 6^ 7 6 76 7, ð4Þ 4 5¼6 76 ^ 7 ^ ^ ^ ^ 4 54 5 un H0n H1n Hpn ap where the term ui for i¼1,y,n is a discrete observation of the output response, while the element Hpi represents the value of the p-th order Hermite polynomial calculated for the i-th value of the input parameter. It becomes apparent that the
G. Georgiou et al. / Mechanical Systems and Signal Processing 32 (2012) 32–43
37
accuracy of the regression coefficients depends upon the quantity and the quality of the obtained observations and for this reason, an efficient sampling technique, such as the Latin Hypercube method, which ensures that the points of the input variables are selected with equal probability, was implemented prior to the regression analysis. After the evaluation of the unknown coefficients, the Polynomial Chaos Expansion models can be emulated for given values of the uncertain parameters and predict statistical properties of the output response, such as the mean value, the standard deviation and the probability density function at low computational cost. The Polynomial Chaos Expansion method was selected in the current study as a very efficient non-intrusive method used to produce probabilistic models for the output responses of complicated structures with many variables. Additionally, an important advantage of the Polynomial Chaos Expansion method is the simplistic formulation of the output responses in a polynomial form, making the technique suitable in a range of different applications (e.g. [18, 25–28]).
5. Numerical results In this section, representative numerical results on the aeroelastic stability and gust response of the composite wing model incorporating deterministic or probabilistic material and/or damage parameters are presented. Initially, the baseline aeroelastic model of the composite wing with no variation in material properties and without the inclusion of damage was investigated in a sequence of unsteady aerodynamic analyses, including flutter analysis and dynamic aeroelastic response analysis. The purpose of this baseline study was to identify the dynamic instabilities (flutter mechanisms) and the gust response history of typical ‘‘Interesting Quantities’’, such as the wing root shear force. Following this study, the uncertain nature of the mechanical properties was taken into account in a probabilistic manner; longitudinal Young’s modulus and the structural damping were varied in a stochastic manner. The probability density functions for the flutter speed, and also the maximum and minimum root shear force were calculated based upon initial Latin Hypercube samples and applying the Polynomial Chaos Expansion method. Finally, the effect of damage due to transverse matrix cracks and fiber breakage was introduced into the probabilistic aeroelastic model of the wing, and the severity of the damage on the aeroelastic stability and gust response was studied.
5.1. Deterministic aeroelastic models 5.1.1. Flutter analysis The ‘‘p-k’’ method embedded in MSC.NASTRAN software and in particular the PKNL variant which enables consistent ‘‘matched’’ air density, Mach number and air speeds, was applied to determine the frequency and damping ratio trends of the composite wing at subcritical speeds (Fig. 5) [32]. Inspection of the frequency and damping ratio curves in Fig. 5 shows that the first bending mode (Mode 1) enters a stable non-oscillatory motion after 65 m/s. Additionally, the frequency values of the first torsion mode (Mode 5) exhibit a significant reduction with respect to the air speed, and approximately at air speed 114 m/s the frequency curve intersects the corresponding curve of the third bending mode of the structure (Mode 4), without although affecting the flutter stability. However, Mode 4 starts to interact with the first sway mode (Mode 3) at around 120 m/s, leading to a flutter instability of Mode 3 at an air speed of 167.53 m/s.
Fig. 5. Frequency and damping ratio trends versus air speed for the composite wing model.
38
G. Georgiou et al. / Mechanical Systems and Signal Processing 32 (2012) 32–43
5.1.2. Gust analysis The unsteady aerodynamic behavior of the composite wing was also investigated in a dynamic aeroelastic analysis which included a discrete ‘‘1-cosine’’ vertical gust. The maximum gust amplitude wg0 was calculated using the following expression [33]: 1=6 H ð5Þ wg0 ¼ wref 106:17 where H represents the gust gradient (half the gust wavelength Lg), while the reference gust velocity wref reduces linearly from 17.07 m/s EAS at sea level to 13.41 m/s EAS at 4752 m, and then again to 6.36 m/s EAS at 18288 m. Correspondingly, in the time domain the gust velocity is defined by 8 < wg0 1cos 2pV t if 0 r t r t g 2 Lg ð6Þ wg ðtÞ ¼ :0 if t Z t g where V is the air speed, while tg ¼Lg/V. Three different wavelengths Lg ¼[18 m, 116 m, 214 m](Fig. 6), at sea level with air speed V ¼174.85 m/s were selected as flight conditions in order to evaluate the time history of the root shear force developed at the bottom and at the top skin (Fig. 7). 5.2. Probabilistic aeroelastic models without damage In an attempt to predict the behavior of the composite wing under uncertainty, probabilistic aeroelastic models of the wing box were built and simulated. Longitudinal Young’s modulus and the structural damping of all the components of the wing box were considered as Gaussian parameters and the aeroelastic stability and gust response of the wing were evaluated. The uncertain longitudinal Young’s modulus was assumed to have a mean value 136.52 GPa and coefficient of variance 0.03, whereas the structural damping was chosen to have mean value 0.02 and coefficient of variance 0.2.
Fig. 6. Gust velocity for different wavelengths.
Fig. 7. Time histories of the root shear force for different wavelengths obtained: (a) at the bottom skin and (b) at the top skin.
G. Georgiou et al. / Mechanical Systems and Signal Processing 32 (2012) 32–43
39
The Polynomial Chaos Expansion method [29–31] was implemented in a non-intrusive way in order to predict the uncertainty propagation for the flutter speed, the maximum and minimum root shear force of the composite wing. The PCE models of the aforementioned quantities were constructed by exploiting the statistical properties of a small sample of observations (size 50), which was obtained using the Latin Hypercube sampling technique. A large number of emulations (100,000) for the PCE models were conducted at low computational cost to produce the probability density functions for the flutter speed, the maximum and the minimum root shear force. The PDFs obtained from the PCE models were compared to those found from 5000 simulations of the actual models (a Latin Hypercube sample) and an excellent agreement was achieved. 5.2.1. Flutter analysis Several case studies were conducted in an attempt to determine the variability on the aeroelastic stability of the composite wing with respect to uncertain material properties applied to all the components of the wing box. The probability density function of the flutter speed was predicted using the Polynomial Chaos Expansion method described above, for variations in longitudinal Young’s modulus and the structural damping. The obtained PDFs are presented in Fig. 8(a) and (b), and it becomes apparent that the inclusion of the structural damping as independent variable has little effect upon the shape and the range of the derived distribution. Direct comparison of the PDFs obtained by simulation of the actual models and the corresponding approximated PDFs revealed the accuracy of developed PCE models. 5.2.2. Gust analysis In an effort to predict the uncertainty propagation of material properties on the gust response of the composite wing, several representative scenarios were performed and a typical ‘‘Interesting Quantity’’ of the wing structure was evaluated. The probability density function of the maximum and minimum root shear force was predicted using the Polynomial Chaos Expansion method, selecting as uncertain variables the longitudinal Young’s modulus and the structural damping. The derived PDFs are included in Fig. 9(a) and (b), and it is worth to mention that the shape as well as the dispersion of the values for the maximum and minimum root shear force presented significant differences. However, the PDFs approximations of the output responses were verified by the comparison with simulations of the actual models performed using the Latin Hypercube sampling technique and a very good agreement was achieved. 5.3. Probabilistic aeroelastic models with damage The effect of damage on the aeroelastic stability and response of the composite wing box was examined by inducing two damage modes, the transverse matrix cracking as well as the fiber breakage in a separate sequence at the outer layers of the bottom skin and at the middle layers of the leading spar. The severity of damage due to the transverse matrix cracking was investigated, considering the crack density as an uncertain Gaussian variable with mean value 0.3 and coefficient of variance 0.25. According to the self-consistent method, only the stiffness and compliance of composite laminates are changed due to the existence of open cracks, while the structural damping remains unaffected. However, previous studies (e.g. [34]) show that the damping in composite materials does tend to increase with the severity of damage. Thus, a combined case study was conducted in which the crack density and material properties, including longitudinal Young’s modulus and the structural damping were varied in a stochastic manner with mean value 136.52 GPa and coefficient of variance 0.03, and mean value 0.02 and coefficient of variance 0.2, respectively. Accurate predictions of the flutter speed and the maximum/minimum root shear force were achieved using the Polynomial Chaos Expansion method. A regression method was implemented in order to calculate the coefficients of the Hermite polynomials, using a small sample of observations (size 50) obtained by applying the Latin Hypercube sampling
Fig. 8. Probability density functions of the flutter speed: (a) for uncertain longitudinal Young’s modulus and (b) for uncertain longitudinal Young’s modulus and structural damping.
40
G. Georgiou et al. / Mechanical Systems and Signal Processing 32 (2012) 32–43
Fig. 9. Probability density functions of the wing root shear force for uncertain longitudinal Young’s modulus and structural damping: (a) maximum case and (b) minimum case.
Fig. 10. Probability density function of the flutter speed: (a) for uncertain crack density, (b) for uncertain longitudinal Young’s modulus, structural damping and crack density and (c) for uncertain longitudinal Young’s modulus, structural damping and inclusion of fiber breakages.
technique. The uncertainty quantification of the aeroelastic behavior for the composite wing was accomplished by performing 100,000 emulations of the developed PCE models. The statistical properties of the PCE models were verified against the relevant characteristics of the actual models by conducting 5000 simulations and through the evaluation of the probability density functions for the flutter speed, the maximum and minimum root shear force. 5.3.1. Flutter analysis Three discrete case scenarios were conducted in an effort to predict the uncertainty propagation of the damage metric (crack density) on the aeroelastic stability of the composite wing taking into account the variability of the crack density (Fig. 10(a)) as well as the Gaussian distribution of the longitudinal Young’s modulus, the structural damping and the crack density (Fig. 10(b)). It is worth to mention that the PDF of the flutter speed due to the uncertainty of the crack density presented a skewed distribution with small dispersion of the values. Consequently, the variability of the output response
G. Georgiou et al. / Mechanical Systems and Signal Processing 32 (2012) 32–43
41
Fig. 11. Probability density functions of the wing root shear force for uncertain crack density: (a) maximum case and (b) minimum case.
in the case where longitudinal Young’s modulus and the structural damping were considered as additional random variables (Fig. 10(b)) did not present significant differences in comparison with the corresponding results included in Fig. 8(b). A separate scenario was performed at which the longitudinal Young’s modulus and the structural damping were selected as uncertain parameters along with the inclusion of fiber breakages in pre-defined composite layers (Fig. 10(c)). The relevant shape of the derived PDF followed the normal distribution, although the mean value of the flutter speed was significantly reduced. In all the examined cases, the developed PCE models provided accurate PDFs predictions in comparison with the corresponding results obtained through the simulation of the actual models. Up to the edge of the stratosphere (11 km), the temperature decreases linearly (6.5 1C per km) and this has a major effect on the aerodynamic forces at a given altitude. The International Standard Atmosphere sets the sea level temperature at 15 1C, so it is of interest to note the effect of varying the effective sea level temperature between the possible extremes that might be encountered. Fig. A1 in Appendix 1 shows how the flutter speed varies over a range of about 4 m/s between the extreme ranges of 30 1C and þ40 1C. This range is much larger that might be expected in reality, so this variation in flutter speed can be considered to be of less consequence than the structural damage. This is an area that requires further detailed examination. 5.3.2. Gust analysis The severity of damage on the gust response of the composite wing was predicted through various case studies in which the probability density function of the maximum and minimum root shear force was evaluated using the Polynomial Chaos Expansion method. Initially, the crack density (Fig. 11(a) and (b)) was considered as the uncertain variable. It is worth to mention that the obtained distributions displayed significant dispersion of the values as well as negative and positive skewness for the maximum and minimum root shear force, respectively. This investigation was extended by incorporating as additional stochastic parameters the longitudinal Young’s modulus and the structural damping. The derived PDFs are presented in Fig. 12(a) and (b) in superposition with the corresponding results gathered for the uncracked material case (Fig. 9(a) and (b)). Direct comparison of the PDFs for the two discrete cases (uncracked and cracked material) revealed that the distribution for the maximum root shear force of the composite wing with cracked material was shifted towards higher values and the shape was changed considerably. However, the shape and dispersion of the values for the distribution of the minimum root shear force for the uncracked and cracked material cases did not present significant differences and the relevant PDF was only shifted towards lower values. Finally, similar behavior was observed in the scenario in which fiber breakages were included in the damage modeling of the composite layers (Fig. 12(c)). Particularly, the obtained PDF of the maximum root shear force was shifted towards higher values in comparison with the undamaged material case and the distribution presented a wider dispersion of the values. 6. Conclusions The Polynomial Chaos Expansion (PCE) method has been applied to predict the uncertainty propagation on the aeroelastic stability and gust response of a composite wing, due to uncertainties related to the material properties and the severity of damage, in an accurate and computationally efficient way. Initially, the material properties of the wing structure were used as stochastic variables with Gaussian distribution and characteristic output responses, such as the flutter speed, the maximum and minimum root shear force were considered as metrics for the aeroelastic behavior of the composite wing. The statistical properties of the output responses, in the form of probability density functions (PDF), were accurately determined by implementing PCE models, built from a small sample of observations. The PDF of the flutter speed was predicted using the PCE method, for variations in longitudinal Young’s modulus and the structural damping, where it became apparent that the inclusion of the structural damping as independent variable had little effect upon the
42
G. Georgiou et al. / Mechanical Systems and Signal Processing 32 (2012) 32–43
Fig. 12. Probability density functions of the wing root shear force for uncertain longitudinal Young’s modulus, structural damping and crack density: (a) maximum case, (b) minimum case and (c) maximum case with inclusion of fiber breakages.
shape and the range of the derived distribution. Concurrently, the PDFs of the maximum and minimum root shear force were predicted using the PCE method and it was found that the shape as well as the dispersion of the values for the interesting quantities presented significant variability due to the manufacturing uncertainties. The severity of damage in pre-defined locations of the composite layers was also considered as stochastic variable with Gaussian distribution and the effect of damage due to transverse matrix cracking was incorporated into the laminates using the self-consistent method. The extreme case of fiber fracture was also investigated. The PDF of the flutter speed was predicted using the PCE method, for variations in the material properties and taking into account the damage severity. The results did not present significant differences in comparison with the corresponding findings for the uncracked material, due mainly to the assumption for the damage locations. On the other hand, the mean value of the flutter speed was reduced when the fiber breakages were included in pre-defined composite layers. Different trends were observed for the maximum and minimum root shear force after the superposition of the results gathered for the uncracked and cracked material cases, where the shape and dispersion of the values for the distribution of the interesting quantities presented significant differences. All the derived PDFs were correlated against corresponding results obtained by performing an adequate number of simulations of the actual models, providing an excellent agreement and the future possibility to apply the developed methodology in order to efficiently investigate the effect of damage severity in robust design of composites wings.
Acknowledgments This work is supported through the Virtual Engineering Centre (VEC), which is a University of Liverpool initiative in partnership with the Northwest Aerospace Alliance, the Science and Technology Facilities Council (Daresbury Laboratory), BAE Systems, Morson Projects and Airbus (UK). The VEC is funded by the Northwest Regional Development Agency (NWDA) and European Regional Development Fund (ERDF) to provide a focal point for virtual engineering research, education and skills development, best practice demonstration, and knowledge transfer to the aerospace sector. Appendix 1 The effect of sea level temperature on the flutter speed for the composite wing model is depicted in Fig. A1.
G. Georgiou et al. / Mechanical Systems and Signal Processing 32 (2012) 32–43
43
Fig. A1. Flutter speed trend versus air temperature for the composite wing model.
References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28] [29] [30] [31] [32] [33] [34]
T.A. Weisshaar, Aeroelastic tailoring of forward swept composite wings, J. Aircr. 18 (8) (1981) 669–676. F.E. Eastep, V.A. Tischler, V.B. Venkayya, N.S. Khot, Aeroelastic tailoring of composite structures, J. Aircr. 36 (6) (1999) 1041–1047. T.A. Weisshaar, D.K. Duke, Induced drag reduction using aeroelastic tailoring with adaptive control surfaces, J. Aircr. 43 (1) (2006) 157–164. H. Arizono, K. Isogai, Application of genetic algorithm for aeroelastic tailoring of a cranked-arrow wing, J. Aircr. 42 (2) (2005) 493–499. M. Kameyama, H. Fukunaga, Optimum design of composite plate wings for aeroelastic characteristics using lamination parameters, Comput. Struct. 85 (2007) 213–224. S. Guo, Aeroelastic optimization of an aerobatic aircraft wing structure, Aerosp. Sci. Technol. 11 (2007) 396–404. C.L. Pettit, R.V. Grandhi, Optimization of a wing structure for gust response and aileron effectiveness, J. Aircr. 40 (6) (2003) 1185–1191. T.U. Kim, I.H. Hwang., Optimal design of composite wing subjected to gust loads, Comput. Struct. 83 (2005) 1546–1554. J.-M. Berthelot, Transverse cracking and delamination in cross-ply glass-fiber and carbon-fiber reinforced plastic laminates: static and fatigue loading, Appl. Mech. Rev. 56 (2003) 111–147. G.J. Dvora, N. Laws, M. Hejazi, Analysis of progressive matrix cracking in composite laminates. I. Thermoelastic properties of a ply with cracks, J. Compos. Mater. 19 (1985) 216–234. R. Talreja, Stiffness properties of composite laminates with matrix cracking and interior delamination, Eng. Fract. Mech. 25 (1986) 751–762. P. Soundappan, E. Nikolaidis, R.T. Haftka, R.V. Grandhi, R.A. Canfield, Comparison of evidence theory and Bayesian theory for uncertainty modeling, Reliab. Eng. Syst. Saf. 85 (2004) 295–311. K.L. Wood, K.N. Otto, E.K. Antonsson, Engineering design calculations with fuzzy parameters, Fuzzy Sets Syst. 52 (1992) 1–20. H. Bae, R.V. Grandhi, R.A. Canfield, Uncertainty quantification of structural response using evidence theory, AIAA J. 41 (10) (2003) 2026–2068. A. Hunter, S. Parsons, Applications of Uncertainty Formalisms, Springer, New York, 1998, pp. 8–32. G. Shafer, A Mathematical Theory of Evidence, Princeton University Press, Princeton, NJ, 1976. Y. Ben-Haim, Information-Gap Decision Theory: Decisions under Severe Uncertainty, Academic Press, 2001. A. Manan, J.E. Cooper, Design of composite wings including uncertainties: a probabilistic approach, J. Aircr. 46 (2) (2009) 601–607. A. Manan, J.E. Cooper, Robust Design of Composite Wings for Gust Response, Proceedings of the 50th AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics, and Materials Conference, Palm Springs, California, May 4–7, 2009. C.L. Pettit, Uncertainty quantification in aeroelasticity: recent results and research challenges, J. Aircr. 41 (5) (2004) 1217–1229. P.S. Beran, C.L. Pettit, D.R. Millman, Uncertainty quantification of limit-cycle oscillations, J. Comput. Phys. 217 (2006) 217–247. J. Kuttenkeuler, U. Ringertz, Aeroelastic tailoring considering uncertainties in material properties, Struct. Optim. 15 (1998) 157–162. S.C. Catravete, R.A. Ibrahim, Effect of stiffness uncertainties on the flutter of cantilever wing, AIAA J. 46 (4) (2008) 925–935. L.G. Crespo, D.P. Giesy, S.P. Kenny, Robustness analysis and robust design of uncertain systems, AIAA J. 46 (2) (2008) 388–396. F.S. Hover, M.S. Triantafyllou, Application of polynomial chaos in stability and control, Automatica 42 (2006) 789–795. L. Mathellin, M.Y. Hussaini, T.A. Zang, Stochastic approaches to uncertainty quantification in CFD simulations, Numer. Algorithms 38 (2005) 209–236. O.M. Knioa, O.P. Le Maˆıtre, Uncertainty propagation in CFD using polynomial chaos decomposition, J. Fluid Dyn. Res. 38 (2006) 616–640. S.K. Choi, R.A. Grandhi, R.A. Canfield, C.L. Pettit, Polynomial chaos expansion with Latin hypercube sampling for estimating response variability, AIAA J. 42 (6) (2004) 1191–1198. N. Wiener, The homogeneous chaos, Am. J. Math. 60 (4) (1938) 897–936. K. Ito, Multiple Wiener integrals, J. Math. Soc., Jpn 3 (1951) 157–169. R.G. Ghanem, P.D. Spanos, Stochastic Finite Elements: A Spectral Approach, Springer-Verlag Inc, New York, 1991. W.P. Rodden, E.H. Johnson, MSC.Nastran Aeroelastic Analysis User’s Guide v68, The MacNeal-Schwendler Corporation, 1994. J.R. Wright, J.E. Cooper, Introduction to Aircraft Aeroelasticity and Loads, John Wiley and Sons Ltd, 2007. Z. Zhang, G. Hartwig, Relation of damping and fatigue damage of unidirectional fibre composites, Int. J. Fatigue 24 (2002) 713–718.