NOMENCLATURE
- Amat
State matrix
- AC
Aerodynamic center
- b
Mid-chord
- Ca
Aerodynamic damping matrix
- Cs
Structural damping matrix
- C Lθ
Lift curve slope
- C gw
Center of gravity
- E
Elastic modulus
- EA
Elastic axis
- EI
Bending rigidity
- G
Shear modulus
- GJ
Torsional rigidity
- I
Wing cross-sectional moment of inertia
- I P
Wing moment of inertia about P
- J
Wing torsion constant
- Ka
Aeroelastic stiffness matrix
- Ks
Structural stiffness matrix
- L
Aerodynamic lift
- M
Aerodynamic moment
- Ma
Apparent mass matrix due to non-circulatory forces
- Ms
Inertia matrix
- M Q
Aerodynamic moment
- R
Reliability
- S
Flutter safety
- U F
Flutter speed
- U
Air stream velocity
- V
Volume of interference region
- V f
Failure volume
- V s
Safe volume
- k h
Flexural stiffness
- k θ
Torsional stiffness
- l
Wing length
- m
Typical section mass
- w
Wing bending deflection
- x θ
Chord-wise offset of the center of mass from the reference point
- χ
Possibility of each event
- γ
Modal damping
- λ
Eigenvalue
- λ 0
Induced flow velocity
$\rho$Air density
- $\theta$
Wing torsion deflection
- $\omega$
Modal frequency
- ζ
Uncertain fuzzy parameter
- $\boldsymbol{\Gamma}$
Modal damping vector
- $\boldsymbol{\Lambda}$
Eigenvalues vector
- $\boldsymbol{\Omega}$
Modal frequency vector
- $\boldsymbol{\Psi}$
State vector
- $\xi$
Generalised coordinate vector
- ζ
Vector of uncertain fuzzy parameters
Subscripts
- F
Flutter
- c
crisp value – nominal value
- f
failure region
- s
safe region
- min
minimum value
- max
maximum value
1.0 INTRODUCTION
Flutter is an undesirable phenomenon that may take place in an aeroelastic wing. One way to predict the flutter speed is via theoretical calculations based on experimentally obtained wing parameters(Reference Afonso, Vale, Oliveira, Lau and Suleman1–Reference Wright and Cooper5). Most of these parameters are physically uncertain due to manufacturing and operational conditions. In general, uncertainties in flight vehicles are divided into two major categories: internal and external sources. Structural and geometric uncertainties are examples of internal uncertainties(Reference Friswell and Mottershead6–Reference Massa, Mourier-Ruffin, Lallemand and Tison8) and aerodynamic and gust loads are examples of external uncertainties(Reference Fazelzadeh and Sadat-Hoseini9–Reference Sadat-Hoseini, Fazelzadeh, Rasti and Marzocca11).
Due to the lack of sufficient knowledge, the estimation of an appropriate model for the uncertain parameters is very important to estimate the flutter speed. Moreover, a reliability analysis is needed to determine the possibility of flutter failure; in this paper reliability is defined in terms of failure due to the flutter instability. Probabilistic and non-probabilistic methods are generally the two main approaches used for modeling uncertainties in structures. Probabilistic methods are based on generating a lot of data which leads to high cost calculations. Another problem with these methods is the lack of data that could be used to determine the statistical distribution of uncertain parameters. Non-probabilistic methods have been preferred in recent years due to their low-cost computations.
In this regard, the modeling of uncertain structures using a non-probabilistic interval method was conducted by Rao and Berke(Reference Rao and Berke12). A non-traditional uncertainty treatment for mechanical problems was investigated by Muhanna and Mullen(Reference Muhanna and Mullen13). They introduced uncertainties as bounded possible values. Moreover, they modeled uncertainty with interval arithmetic and applied this method to beams and trusses. The numerical algorithms of non-probabilistic convex models and an interval method for the static displacement of structures with uncertain parameters was presented by Qiu(Reference Qiu14). Qiu and Wang(Reference Qiu and Wang15) also investigated the non-probabilistic interval analysis method for the dynamical response of structures with interval parameters. Manson(Reference Manson16) used affine and interval arithmetic to solve a two-degrees-of-freedom eigenvalue problem. Qiu(Reference Qiu17) also used an interval analysis method to predict the effect of uncertain-but-bounded parameters on the buckling of composite structures. Moens and Vandepitte(Reference Moens and Vandepitte18) gave a review on the emerging non-probabilistic approaches for uncertainty treatment in finite element analysis. They discussed general theoretical and practical aspects of both interval and fuzzy finite element analysis. Muhanna et al.(Reference Muhanna, Zhang and Mullen19) combined the finite element method and interval analysis to analyze the system response subject to stiffness and loading uncertainties. The influence of uncertainty parameters using interval numbers on the flutter speed of a wing was conducted by Wang and Qiu(Reference Wang and Qiu20). They used a first-order Taylor series expansion to predict the lower and upper bounds of flutter speed. The problem of robust stability of a two-dimensional nonlinear aeroelastic system with uncertainties using the 
${\rm{\mu}}$-method was investigated by Yun and Hun(Reference Yun and Han21). The effect of parametric uncertainty on the stall flutter bifurcation behavior of a pitching airfoil was conducted by Sarkar et al.(Reference Sarkar, Witteveen, Loeven and Bijl22). Khodaparast et al.(Reference Khodaparast, Mottershead and Badcock23) studied the problem of linear flutter analysis in the presence of uncertainties. The use of eigenvalue stability to analyse very large dimension aeroelastic numerical models arising from the exploitation of computational fluid dynamics has been reviewed by Badcock et al.(Reference Badcock, Timme, Marques, Khodaparast, Prandina, Mottershead, Swift, Da Ronch and Woodgate24). Yang et al.(Reference Yang, Cai and Liu25) presented a new interval-based method for the analysis of uncertain structures using the Laplace transform. The upper and lower bounds of the natural frequencies of structures with uncertain but bounded parameters were evaluated by Sofi et al.(Reference Sofi, Muscolino and Elishakoff26). They also applied improved interval analysis via an extra unitary interval (EUI).
Wang et al.(Reference Wang, Xiong, Hu, Wang and Qiu27) developed a sequential multidisciplinary design optimisation and reliability analysis method under non-probabilistic theory to decouple the reliability analysis from the optimisation. They also used an improved dimension-wise method for multidisciplinary interval uncertainty analysis(Reference Wang, Xiong, Wang, Xu and Li28).
Mannini and Bartoli(Reference Mannini and Bartoli29) proposed a method to approach flutter instability and calculated the critical wind speed, starting from the probability distribution of the flutter derivatives. A probabilistic flutter analysis utilising a meta-modeling technique to evaluate the effect of parameter uncertainty on the flutter speed was conducted by Abbas and Morgenthal(Reference Abbas and Morgenthal30). Lokatt(Reference Lokatt31) approximated the aerodynamic model using a piece-wise continuous rational polynomial function. They proposed this method for efficient flutter analysis of aeroelastic systems including modelling uncertainties. Huan et al.(Reference Huan, Yuan and Chao32) used the polynomial chaos expansion (PCE) method for uncertainty quantification and showed that this method has less computational cost compared to Monte Carlo simulation.
Several other methods have been investigated to propagate uncertainty in mechanical structures. Friswell and Mottershead(Reference Friswell and Mottershead6) described various methods for parameter selection, error localisation, and sensitivity analysis and estimation, in mechanical structures. Also, many efforts were conducted on model updating in uncertain mechanical structures(Reference Boulkaibet, Marwala, Friswell, Khodaparast and Adhikari33–Reference Mottershead, Friswell, Ng and Brandon36). The fuzzy approach has been used for uncertainty modeling and propagation, and this non-probabilistic method is computationally low-cost compared to probabilistic methods(Reference Starczewski37). Chiang et al.(Reference Chiang, Dong and Wong38) modeled structures with fuzzy and random uncertainties. They also studied the response of structures with stiffness, damping and mass uncertainties. A fuzzy methodology to calculate the eigenvalues and eigenvectors of an uncertain mechanical structure was proposed by Massa et al.(Reference Massa, Lallemand, Tison and Level7). They described material and geometric parameters as imprecise fuzzy numbers. De Gersem et al.(Reference De Gersem, Moens, Desmet and Vandepitte39) proposed the fuzzy finite element and interval methods to carry out frequency response and eigenvalue analysis of structures with uncertain parameters. The flutter dynamic pressure of a semi-span super-sonic wind-tunnel model was predicted by Tartaruga et al.(Reference Tartaruga, Cooper, Georgiou and Khodaparast40). They used probabilistic and non-probabilistic approaches in their study. Khodaparast et al.(Reference Haddad Khodaparast, Govers, Dayyani, Adhikari, Link, Friswell, Mottershead and Sienz41) presented the application of fuzzy finite element model updating to the DLR AIRMOD structure. Rezaei et al.(Reference Rezaei, Fazelzadeh, Mazidi and Khodaparast42) investigated the flutter uncertainty analysis of an aircraft wing subjected to a thrust force using a fuzzy method. They modeled the uncertain parameters as triangle and trapezium membership functions. The eigenvalue problem with fuzzy input parameters was solved using the fuzzy Taylor expansion method, and a sensitivity analysis was performed. Also, the upper and lower bounds of the flutter region at different α-cuts were extracted. Wang et al.(Reference Wang, Xiong and Yang43) proposed a new reliability estimation model based on the level cut strategy and volume ratio theory.
Although many researchers have studied uncertainty propagation and identification in structures, there are limited works in the field of reliability. The reliability and free vibration of a cantilever composite beam under structural uncertainty was conducted by Oh and Librescu(Reference Oh and Librescu44). The structural uncertainty was propagated using Monte Carlo Simulation and the Stochastic Rayleigh-Ritz method to find the reliability of the beam at different frequencies. Di Sciuva and Lomario(Reference Di Sciuva and Lomario45) applied reliability methods to an isotropic beam and a laminated composite plate. Frangopoulos et al.(Reference Frangopol and Maute46) presented a review of life-cycle reliability optimisation with an emphasis on aerospace and civil structures. Wang and Qiu(Reference Wang and Qiu47) proposed a method to calculate the reliability of an aeroelastic wing using an interval approach. Increasing the reliability of aircraft structures subject to air loads was conducted by Bijl et al.(Reference Bijl, Lucor, Mishra and Schwab48). They also conducted some aeroelastic analysis and reliability studies to illustrate this key concept.
In this paper, uncertainty modeling is conducted using a possibility, rather than probability, approach. Possibility means something happens but with different quality, while probability means something happens or not. In conventional reliability, when the reliability of a system is 0.99, it means that if we have 100 systems, one of them may have failed on the defined criterion. In fuzzy reliability, when the reliability of a system is 0.99, it means that if we have 100 systems, none of them fail under the defined crierion but all of them have a 1% imperfection. Thus, fuzzy reliability speaks about the quality instead of quantity. The uncertain physical parameters are modeled as triangular fuzzy membership functions, although other membership functions, such as trapezoidal or Gaussian, may also be used. The uncertainty is then propagated through the aeroelastic wing model; a flutter analysis is conducted, and a fuzzy region of flutter is obtained instead of a deterministic flutter region. Furthermore, the reliability of the wing against flutter is determined from a pyramid based on the flutter speed and the air speed interference area.
To the best of the authors’ knowledge, in the pertinent literature, the reliability of an aeroelastic wing using this type of fuzzy approach has not yet been presented. This research intends to fill the gap in knowledge related to this problem. In this paper, the stability region is presented as a three-dimensional fuzzy pyramid shape. Furthermore, modal damping and frequency diagrams at different α-cuts are presented.
2.0 FLUTTER ANALYSIS OF A DETERMINISTIC WING MODEL
There are three general methods to estimate the wing flutter under unsteady aerodynamic loads, namely the K, PK and P methods(Reference Hodges and Pierce49). The P method, including the finite state unsteady aerodynamic loading of Peters et al.(Reference Peters, Karunamoorthy and Cao50), is most suited to the estimation of the flutter boundary(Reference Fazelzadeh, Mazidi and Kalantari51,Reference Fazelzadeh, Ghasemi and Mazidi52) and hence is used in this paper. In this method, the aeroelastic equations are converted to the state space form and the flutter boundary is obtained by solving an eigenvalue problem. The complex eigen solutions contain real and imaginary parts; by interpreting these values at different wind speeds, the flutter speed is determined from the stability of the eigenvalues.
The general discretised form of the wing aeroelastic governing equations can be expressed as:
\begin{equation}\left( {{{\bf{M}}_s} + {{\bf{M}}_a}} \right)\;{\ddot{{\boldsymbol{\rm{\xi}}}}}\left( t \right)+ \left( {{{\bf{C}}_s} + {{\bf{C}}_a}} \right)\;{\ddot{{\boldsymbol{\rm{\xi}}}}}\left( t \right){\bf{ + }}\left( {{{\bf{K}}_s} + {{\bf{K}}_a}} \right)\;{\boldsymbol{\rm{\xi}}}\left( t \right)= \left\{ {\bf{O}} \right\},\end{equation}where 
$\boldsymbol{\xi}$ is the generalised coordinates vector, Ms is the inertia matrix, Ma is the apparent mass matrix due to non-circulatory aerodynamic forces, Cs is the structural damping matrix, Ca is the aerodynamic damping matrix, Ks is the structural stiffness matrix and Ka is the aeroelastic stiffness matrix due to circularity forces(Reference Nguyen53,Reference Nguyen and Tuzcu54) . This second-order differential equation can be transformed into a set of first-order differential equations in state space:
\begin{equation}{\dot{\boldsymbol{\rm{\psi}}}}\left( {\bf{t}} \right) = {{\bf{A}}_{{\bf{mat}}}}{\boldsymbol{\rm{\psi}}}\left( {\bf{t}} \right),\end{equation}
where the state vector 
$\boldsymbol{\rm{\psi}}\bf{(t)}$ is defined as:
\begin{equation}\boldsymbol{\rm{\psi}}\bf{(t)} = \left[\begin{array}{ccccc}
\boldsymbol{\rm{\xi}}\!\left({\bf{t}}\right)\\[4pt]
{\dot{\boldsymbol{\rm{\xi}}}}\!\left({\bf{t}}\right)\end{array}\right].\end{equation}The system matrix Amat is obtained as:
\begin{equation}
{\bf{A}}_{\bf{mat}} = \left[\begin{array}{ccc}
{\bf{O}} & {\bf{I}}\\[5pt]
- \left({\bf{M}}_s + {\bf{M}}_a\right)^{-1}\left({\bf{K}}_s + {\bf{K}}_a \right) & - \left({\bf{M}}_s + {\bf{M}}_a\right)^{-1}
\left({\bf{C}}_s + {\bf{C}}_a\right)\end{array}\right]. \end{equation}After solving Equation (2), the eigenvalue vector is derived as
\begin{equation}{\boldsymbol{\Lambda}} = {\boldsymbol{\Gamma}} \pm {\bf{i}}{\boldsymbol{\Omega}},\end{equation}where 
${\boldsymbol{\Gamma}}$ is modal damping vector (with elements 
$\gamma_{j}$) and 
$\boldsymbol{\Omega}$ (with elements 
$\omega_{j}$) is the modal frequency vector. When 
$\gamma_{j} < 0$ for all j, any transient oscillations decay and the system is dynamically stable. As the wind speed increases, one component of the modal damping vector tends to zero and then becomes positive. The first airspeed at which this element become zero is the flutter speed, and the corresponding modal frequency is the flutter frequency.
3.0 FUZZY UNCERTAINTY APPROACH
The wing flutter speed is generally an uncertain parameter because it depends upon structural and aerodynamic parameters which are physically uncertain due to manufacturing and operational conditions. Airspeed is also not a certain parameter; when the speed is set at a specified value during the flight, the air speed may fluctuate around this value due to atmospheric conditions.
An eigenvalue problem for the P method including uncertain parameters, can be represented as(Reference Massa, Lallemand, Tison and Level7):
\begin{equation}{\dot{\boldsymbol{\rm{\psi}}}}\left({\bf{t}}\right) = {\bf{A}}_{\bf{mat}}\left({\tilde \zeta}_1,{\tilde \zeta}_2,...,{\tilde\zeta}_m \right)
{\boldsymbol{\rm{\psi}}}({\bf{t}}), \end{equation}where the 
$\widetilde{{\rm{\zeta}}_{i}}$ are the uncertain parameters and m is the number of uncertain parameters. The uncertain parameters are modeled as fuzzy numbers. The fuzzy parameter 
$\tilde{\zeta}$, which is shown in Fig. 1, is defined by a variation about a crisp value at each α-cut.
Figure 1. Triangular fuzzy membership function.
According to this figure, an α-cut is the set of all 
${\rm{\zeta}} $ such that 
$\mu \left( \zeta \right)$ is greater than or equal to α. The fuzzy vector 
${{\tilde {\rm{\zeta}} }}$ is defined by a crisp value ζc and the variation 
$\Delta {{\tilde {\rm{\zeta}} }}$ at a given α-cut as
\begin{equation}{\tilde{\boldsymbol{\rm{\zeta}}}}^\alpha = {\boldsymbol{\rm{\zeta}}}_c + \Delta{\tilde{\boldsymbol{\rm{\zeta}}}}^\alpha, \end{equation}or
\begin{equation}{\tilde{\rm{\zeta}}}^\alpha = \xi_{\bf{c}} + \left[{\underline{\Delta {{\rm{\zeta}}}^{\alpha}}};\overline{\Delta{{\rm{\zeta}}}^{\alpha}}\right],\end{equation}where 
$\overline{{\boldsymbol{\rm{\zeta}}}^\alpha} $ and 
$\underline {{{\boldsymbol{{\rm{\zeta}} }}^\alpha }} $ are the maximum and minimum of the fuzzy parameter vector 
${{\tilde{\boldsymbol{\rm{\zeta}}}}}$ for a given α-cut, respectively. The membership function is discretised by different intervals which are linked to α-cuts ranging from 0 to 1.
Many methods to determine output (response) intervals based on the input (parameter) intervals are available, many of which use the exact input-output functional relationship. Here, to demonstrate the proposed approach, the interval is divided into several intervals at each α-cut and a first order Taylor series expansion is used to determine the upper and lower bounds of the flutter frequency and modal damping for each interval.
The interval model can be used to describe nonlinear dynamic systems under uncertainty with low-order Taylor series expansions. However, the Taylor series-based interval method is only suitable for problems with small uncertainty levels(Reference Wu, Zhang, Chen and Luo55). Truncation errors exist in this linear model since the higher-order terms are neglected, but for highly nonlinear problems, the truncation error cannot be ignored(Reference Wu, Luo and Zhang56).
Using this method, the modal damping and frequency are expressed as
\begin{equation}{\gamma^\alpha}\left({\boldsymbol{\rm{\zeta}},U}\right) = {\gamma^\alpha}\left({\boldsymbol{\rm{\zeta}}}_{c},U\right) +
\sum\limits_{j = 1}^m \left|\frac{\partial\gamma^\alpha \left({\boldsymbol{\rm{\zeta}}}_{c},U \right)}{\partial{\zeta_i}^\alpha}\right| \Delta\zeta_{i}^{\alpha},\end{equation}
\begin{equation}{\omega^{\alpha}}\left({\boldsymbol{\rm{\zeta}},U}\right) = {\omega^\alpha}\left({\boldsymbol{\rm{\zeta}}}_{c},U\right) + \sum\limits_{j = 1}^{m}
\left|\frac{\partial{\omega^\alpha}\left({\boldsymbol{\rm{\zeta}}}_{c},U \right)}{\partial\zeta_{i}^\alpha} \right|\Delta \zeta_{i}^\alpha, \end{equation}After applying the interval operation, the lower and the upper bounds of the modal damping at each α-cut are defined, respectively, as:
\begin{equation}\underline{\gamma^{\alpha}} \left(\boldsymbol{\rm{\zeta}},U\right) = \gamma^{\alpha}
\left(\boldsymbol{\rm{\zeta}}_{c},U\right) + \sum\limits_{j = 1}^m \left|\frac{\partial\gamma^{\alpha}\left(\boldsymbol{\rm{\zeta}}_{c},U\right)}{\partial\zeta_{i}^\alpha}\right| \underline{\Delta\zeta_{i}^{\alpha}},\end{equation}
\begin{equation}\overline{\gamma^{\alpha}} \left(\boldsymbol{\rm{\zeta}},U\right) = \gamma^{\alpha}
\left(\boldsymbol{\rm{\zeta}}_{c},U\right) + \sum\limits_{j = 1}^m \left|\frac{\partial\gamma^{\alpha}\left(\boldsymbol{\rm{\zeta}}_{c},U\right)}{\partial\zeta_{i}^\alpha}\right| \overline{\Delta\zeta_{i}^{\alpha}},\end{equation}Also, the lower bound and the upper bound of the modal frequency at each α-cut are:
\begin{equation}\underline{\omega^{\alpha}} \left(\boldsymbol{\rm{\zeta}},U\right) = \omega^{\alpha}
\left(\boldsymbol{\rm{\zeta}}_{c},U\right) + \sum\limits_{j = 1}^m \left|\frac{\partial\omega^{\alpha}\left(\boldsymbol{\rm{\zeta}}_{c},U\right)}{\partial\zeta_{i}^\alpha}\right| \underline{\Delta\zeta_{i}^{\alpha}},\end{equation}
\begin{equation}\overline{\omega^{\alpha}} \left(\boldsymbol{\rm{\zeta}},U\right) = \omega^{\alpha}
\left(\boldsymbol{\rm{\zeta}}_{c},U\right) + \sum\limits_{j = 1}^m \left|\frac{\partial\omega^{\alpha}\left(\boldsymbol{\rm{\zeta}}_{c},U\right)}{\partial\zeta_{i}^\alpha}\right| \overline{\Delta\zeta_{i}^{\alpha}},\end{equation}Note that taking the absolute values of the sensitivities in Equations (11) to (14) may lead to conservative bounds. The flutter speed that corresponds to Equation (11) gives the upper bound of the flutter speed (
$\underline {U_F^\alpha } $), and the flutter speed that corresponds to Equation (12) provides the lower bound of the flutter speed (
$\overline {U_F^\alpha } $). The crisp values can be easily obtained by solving Equation (2). Analytical expressions of the partial derivatives 
${{\partial {\gamma ^\alpha }\left( {{{\boldsymbol{{\rm{\zeta}} }}_c},U} \right)} \mathord{\left/{\vphantom {{\partial {\gamma ^\alpha }\left( {{{\boldsymbol{{\rm{\zeta}} }}_c},U} \right)} {\partial {\zeta _i}^\alpha }}} \right.} {\partial {{\rm{\zeta}} _i}^\alpha }}$ and 
${{\partial {\omega ^\alpha }\left( {{{\boldsymbol{{\rm{\zeta}} }}_c},U} \right)} \mathord{\left/ {\vphantom {{\partial {\omega ^\alpha }\left( {{{\boldsymbol{{\rm{\zeta}} }}_c},U} \right)} {\partial {{\rm{\zeta}} _i}^\alpha }}} \right.} {\partial {\zeta _i}^\alpha }}$ cannot be obtained easily due to the complicated implicit functional relationship between 
$\gamma$ and 
$\omega$. One practical method to calculate these expressions is to use finite difference approximations(Reference Wang and Qiu47). After solving the above equations, the flutter speed bound at each α-cut is derived. These bounds are combined to achieve the example membership function shown in Fig. 2.
Figure 2. An example flutter speed membership function.
4.0 RELIABILITY OF THE WING FLUTTER SPEED
In this section, the reliability of the wing flutter speed is investigated. To this end, the flutter speed and air speed interference pyramid is obtained and then the fuzzy reliability approach is employed to determine the flutter reliability.
4.1 Airspeed and flutter speed interference
The method to obtain the flutter speed membership function was described in the previous section. The air speed can be modeled as another membership function 
$\widetilde{U_{w}}$. It is assumed that the crisp value of the air speed is 
$U_w^{\text{c}}$, and it will not be less than 
$U_w^{\min }$ nor greater than 
$U_w^{\max }$. If the triangular membership function is used for the air speed it means that the possibility of 
$U_w^{\text{c}}$ is α=1. For other interval points, the possibility reduces linearly from 1 to 0. When the air speed is equal to 
$U_w^{\min }$ or lower, or when the air speed is equal to 
$U_w^{\max }$or higher, the possibility of the air speed is zero. The possibility of the air speed membership function 
$\widetilde{U_{w}}$ can be written as:
\begin{equation}
\mu (U) = \left\{
\begin{array}{c@{\qquad}c}
0& {U \le U_w^{\min }}\\[5pt]
{\frac{{U - U_w^{\min }}}{{U_w^c - U_w^{\min }}}}&{U_w^{\min } < U < U_w^c}\\[5pt]
1&{U = U_w^c}\\[5pt]
{\frac{{U - U_w^{\max }}}{{U_w^{\text{c}} - U_w^{\max }}}}&{U_w^{\text{c}} < U
< U_w^{\max }}\\[5pt]
0&{U_w^{\max } \le U}
\end{array}
\right.\end{equation}The air speed and flutter speed membership functions and their two-dimensional interference are shown in Fig. 3.
Figure 3. Flutter speed and wind speed 2-D fuzzy interference.
If there is a region where this interference occurs, then there is a possibility of flutter, otherwise, the flutter possibility is zero. In this section, the flutter reliability of the wing is calculated based on the general non-probabilistic interval reliability model. Wang and Qiu(Reference Wang and Qiu47) used this method for flutter reliability analysis by means of rectangular membership functions for flutter speed and air speed. They represented the flutter speed and air speed in a plane, as shown in Fig. 4. The solid rectangle shows the region of variation of both 
$\widetilde{U_{w}}$ and 
$\widetilde{U_{F}}$. By crossing this region with the failure plane 
${U_{w}} = {U_{F}}$, the safe region and the failure region can be determined as shown in Fig. 4.
Figure 4. Space of variables and the occurrence of interference.
Here, due to the use of triangular membership functions for flutter speed and air speed, such a space of variables occurs at each α-cut. By assembling this space of variables at all α-cuts, a pyramid is created in the space of 
${U_{w}} - {U_{F}} - \alpha$, as shown in Fig. 5. The approach would also work for non-triangular membership functions although the shapes would be more complex and the computation more intensive.
Figure 5. The pyramid of variables in the case of triangular membership functions.
This pyramid shows the region of variations of 
$\widetilde{U_{w}}$, and 
$\widetilde{U_{F}}$ for each α-cut. By cutting this volume with the failure plane 
${U_{w}} = {U_{F}}$, the safe volume and the failure volume are defined as shown in Fig. 6. Using this approach, the effects of both air speed and flutter speed uncertainties are considered in the reliability analysis.
Figure 6. Safe region and failure region.
4.2 Fuzzy reliability of the wing flutter speed
In the design process, the air speed 
$\widetilde{U_{w}}$ is required to be smaller than the flutter speed 
$\widetilde{U_{F}}$. This implies that with respect to the crisp values (α = 1), the wing should be safe. As a result of the dispersion of the fuzzy areas, they may share the same numerical values, shown in Fig. 3 as the shaded region.
As shown before, the space of variables forms a pyramid. This volume is a function of 
${U_w}$ and 
${U_F}$, and defined as 
${F_V}\left( {{U_F},{U_w}} \right)$. By integrating this function, the total volume of the space of variables is obtained as:
$$V = \mathop{{\int\!\!\!\!\!\int\!\!\!\!\!\int}\mkern-31.2mu \bigodot}\limits_V {{F_V}\left( {{U_F},{U_w}} \right)} {\rm{ }}dV,$$The integration of Equation (16) is evaluated numerically. In this method the pyramid in the possible space is divided into cuboids based on the number of alpha cuts, as shown in Fig. 7 schematically. Then the volume of the pyramid is obtained by summing the volumes of all of the cuboids. To evaluate the failure volume, the cuboid for each α-cut is split into two polygon prisms (for example triangular prisms) defined by the failure plane. The volume of the failure region is obtained by summing the volume of all of the prisms within the failure region.
Figure 7. Numerical integration evaluation to obtain the volume of the possible space.
The safe volume 
$V_s$ is the part of this total volume V in which 
${U_w} \le {U_F}$, and the failure volume 
$V_{f}$ is the part of the total volume V in which 
${U_w} \ge {U_F}$ (see Fig. 6). To present the mathematical concept, the following function is defined.
\begin{equation} \Xi = \widetilde{U_F} - \widetilde{U_W},\end{equation}The safe volume is a volume in which Ξ > 0 and the failure volume is a volume in which Ξ < 0. For the case Ξ = 0 the failure plane is created. Assuming that χ is the possibility of each event, the possibility of the safe region is obtained as:
\begin{equation}\chi ( {\Xi
< 0} ) = \frac{{{V_s}}}{V},\end{equation}Similarly, the possibility of the failure region is:
\begin{equation}\chi ( {\Xi > 0} ) = \frac{{{V_f}}}{V},\end{equation}Then the reliability or safety probability of the flutter speed is formulated as:
\begin{equation}R = \chi ( {\Xi < 0} ) = \frac{{{V_s}}}{V} = 1 - \frac{{{V_f}}}{V},\end{equation}Figure 8 shows a case in which all parts of the space of variables are in the safe region. In this case, the flutter airspeed membership function lower bound is larger than the maximum wind speed and the reliability of the flutter occurrence is 1. This means that the flutter will not happen under any circumstances.
Figure 8. A case in which all parts of the interference volume are in the safe region.
If all parts of the space of variables are located in the region where 
$\Xi = \widetilde{U_{F}} - \widetilde{U_{w}} > 0$, as shown in Fig. 9, then the air speed is always larger than the maximum flutter airspeed. In this case, the reliability of flutter occurrence is 0, and hence flutter will definitely happen.
Figure 9. A case in which all parts of the interference volume are in the failure region.
The steps for this method are given as a flowchart in Fig. 10. First, uncertain parameters are determined as fuzzy membership functions. Then these fuzzy membership functions are divided to intervals, and the intervals are evaluated at each value of α. In the next step, the interval corresponding to each α-cut is propagated into the structural equations using a Taylor series expansion, and the upper and lower bounds of the flutter speed are obtained. By assembling the flutter speed at each α-cut the flutter speed membership function is obtained. The 3D interference between the obtained flutter speed and the airspeed membership functions forms a pyramid. By interpreting the volume of this pyramid by the mentioned cutting plane the flutter reliability is obtained.
Figure 10. Fuzzy method flowchart.
5.0 NUMERICAL EXAMPLES
Many theoretical and experimental studies on wing aeroelasticity have been performed in the literature using typical section models(Reference Ajaj, Friswell, Dettmer, Allegri and Isikveren57,Reference Fazelzadeh, Azadi and Azadi58) . These models can represent typical airfoil sections along a finite wing and are still widely used by researchers because of their simplicity. The most popular structural models used for wing aeroelasticity are beam models(Reference Fazelzadeh, Ghasemi and Mazidi52,Reference Goland59,Reference Mazidi, Kalantari and Fazelzadeh60) . These models consider the wing as a flexible beam that experiences, in general, bending in two orthogonal planes and torsion about the elastic axes. In this section, to illustrate the feasibility of the fuzzy reliability method, a typical section 2D wing model and a clean wing 3D model are used.
5.1 Example 1: Typical section 2D wing model
A typical section wing model(Reference Hodges and Pierce49) is shown in Fig. 11. This configuration could represent the case of a rigid, two-dimensional wind-tunnel model that is elastically mounted in a wind-tunnel test section, or could correspond to a typical airfoil section along a finite wing. In the latter case, the discrete springs would reflect the wing structural bending and torsional stiffness, and the reference point would represent the wing elastic axis.
Figure 11. A typical section wing model.
Suppose that the wing mass, moment of inertia and flexural and torsional stiffnesses are uncertain parameters. These parameters should be converted to fuzzy triangle membership functions. The membership functions are expressed as (
$\xi_{min}, \xi_{c}, \xi_{max}$) in which 
$\xi_{min}$ and 
$\xi_{max}$ are minimum and maximum values at 
$\alpha = 0$, and 
$\xi_{c}$ is the value at 
$\alpha = 1$. The fuzzy triangle membership functions of the uncertain parameters are assumed to be 
$\bar m = (0.95m,m,1.05m)$, 
$\overline {{I_P}} = (0.95{I_P},{I_P},1.05{I_P})$, 
$\overline {{k_h}} = (0.95{k_h},{k_h},1.05{k_h})$ and 
$\overline {{k_\theta }} = (0.95{k_\theta },{k_\theta },1.05{k_\theta })$.
The aeroelastic governing equations can be derived as:
\begin{equation}m\left( {\ddot h + b{x_\theta }\ddot \theta } \right) + {k_h}h = - L,\end{equation}
\begin{equation}{I_P}\ddot \theta + mb{x_\theta }\ddot h + {k_\theta }\theta = {M_Q} + b\left( {\frac{1}{2} + a} \right)L,\end{equation}where m is the typical section mass, 
${x_\theta }$ is the chord-wise offset of the center of mass from the reference point, 
$I_{P}$ is the moment of inertia about P, L is the aerodynamic lift, 
$M_{Q}$ is the aerodynamic moment and 
$k_{h}$ and 
$k_{\theta}$ are flexural and torsional stiffnesses, respectively.
Unsteady aerodynamic loads are simulated based on the model of Peters et al.(Reference Peters, Karunamoorthy and Cao50) as
\begin{equation}L = \pi \rho {b^2}\left[ {\ddot h + U\dot \theta \, - ba\ddot \theta } \right] + {C_{L\theta }}\rho Ub\,\,\left[ {\dot h + U\theta + b\left( {\frac{1}{2} - a} \right)\dot \theta - {\lambda _0}} \right],\end{equation}
\begin{equation}{M_Q} = - \pi \rho {b^3}\left( {\frac{1}{2}\ddot h + b\left( {\frac{1}{8} - \frac{a}{2}} \right)\ddot \theta + U\dot \theta } \right),\end{equation}where 
${\lambda _0} = \sum\limits_{n = 1}^N {{b_n}{\lambda _n}} $ is the induced flow velocity and is calculated through a system of N first order coupled differential equations(Reference Hodges and Pierce49).
To solve the above equations, the following dimensionless parameters are introduced.
\begin{equation}{r^2} = \frac{{{I_p}}}{{m{b^2}}}\,\,\,\,\,\sigma = \frac{{{\omega _h}}}{{{\omega _\theta }}}\,\,\,\,\,\mu = \frac{m}{{\pi \rho {b^2}}}\,\,\,\,\,V = \frac{U}{{b{\omega _\theta }}}.\end{equation}For validation, the modal damping versus dimensionless air speed is shown in Fig. 12 and compared with the results given in reference(Reference Hodges and Pierce49). This validation is performed to determine the accuracy of the current aeroelastic governing equations and the solution methodology.
Figure 12. Modal damping versus dimensionless airspeed for a = −0.2, e = 0.1, μ = 20, r 2 = 0.24 and Σ = 0.4.
Furthermore, for model validation, the deterministic typical section model is also compared with an Equivalent 2D Goland wing(Reference Goland59) in which the flexural and torsional stiffness are considered as 
${k_h} = {0.597^4}{\left( {{\pi \mathord{\left/ {\vphantom {\pi L}} \right. } L}} \right)^4}EI$ and 
${k_\theta } = {\left( {{\pi \mathord{\left/ {\vphantom {\pi {2L}}} \right. } {2L}}} \right)^2}GJ$. The parameters are given in Table 1. The flutter speed and frequency are obtained and compared in Table 2. The results show that this 2D model with finite state unsteady aerodynamic loading is in good agreement with the mentioned reference.
Table 1 Goland wing parameters(Reference Borello, Cestino and Frulla61)
Table 2 Deterministic flutter speed and frequency
To obtain the flutter fuzzy membership function the dominant flutter mode is considered. In this example, the bending mode experiences flutter and the other modes are stable.
Assuming that the relationship between the eigenvalues and the uncertain parameters is monotonic, and applying the fuzzy interval method, the modal damping for each α-cut is obtained as:
\begin{equation}\begin{array}{l}\overline {{\gamma ^\alpha }} \left( {m,I,{k_h},{k_\theta },U} \right) = {\gamma ^\alpha }\left( {{m_c},{I_c},{k_{hc}},{k_{\theta c}},U} \right) + \left| {\frac{{\partial {\gamma ^\alpha }\left( {{m_c},{I_c},{k_{hc}},{k_{\theta c}},U} \right)}}{{\partial {m^\alpha }}}} \right|\overline {\Delta m_{}^\alpha } \\[18pt]
+ \left| {\frac{{\partial {\gamma ^\alpha }\left( {{m_c},{I_c},{k_{hc}},{k_{\theta c}},U} \right)}}{{\partial {I^\alpha }}}} \right|\overline {\Delta I_{}^\alpha } + \left| {\frac{{\partial {\gamma ^\alpha }\left( {{m_c},{I_c},{k_{hc}},{k_{\theta c}},U} \right)}}{{\partial {k_h}^\alpha }}} \right|\overline {\Delta k_h^\alpha } \\[18pt]
+ \left| {\frac{{\partial {\gamma ^\alpha }\left( {{m_c},{I_c},{k_{hc}},{k_{\theta c}},,U} \right)}}{{\partial {k_\theta }^\alpha }}} \right|\overline {\Delta k_\theta ^\alpha } \end{array}\end{equation}
\begin{equation}\begin{array}{l}\underline {{\gamma ^\alpha }} \left( {m,I,{k_h},{k_\theta },U} \right) = {\gamma ^\alpha }\left( {{m_c},{I_c},{k_{hc}},{k_{\theta c}},U} \right) + \left| {\frac{{\partial {\gamma ^\alpha }\left( {{m_c},{I_c},{k_{hc}},{k_{\theta c}},U} \right)}}{{\partial {m^\alpha }}}} \right|\underline {\Delta m_{}^\alpha } \\[18pt]
+ \left| {\frac{{\partial {\gamma ^\alpha }\left( {{m_c},{I_c},{k_{hc}},{k_{\theta c}},U} \right)}}{{\partial {I^\alpha }}}} \right|\underline {\Delta I_{}^\alpha } + \left| {\frac{{\partial {\gamma ^\alpha }\left( {{m_c},{I_c},{k_{hc}},{k_{\theta c}},U} \right)}}{{\partial {k_h}^\alpha }}} \right|\underline {\Delta k_h^\alpha } \\[18pt]
+ \left| {\frac{{\partial {\gamma ^\alpha }\left( {{m_c},{I_c},{k_{hc}},{k_{\theta c}},U} \right)}}{{\partial {k_\theta }^\alpha }}} \right|\underline {\Delta k_\theta ^\alpha } \end{array}\end{equation}Using Equations (26) and (27) at each α-cut, the flutter airspeed membership function can be obtained, and is shown in Fig. 13. The flutter boundary range can be seen as a triangular fuzzy mountain shape. For each value of modal damping and at every α-cut, the upper and lower bounds of the flutter speed can be extracted. Values corresponding to α = 0 are the largest intervals, and the value corresponding to α = 1 is deterministic.
Figure 13. Three-dimensional plot of modal damping versus airspeed for different α-cuts.
In the following, the flutter safety of the wing for six different cases is examined. The 3D interference of air speed and flutter speed is shown in Fig. 14. A triangle membership function for airspeed is a reasonable assumption since the airspeed is not a crisp value but varies within intervals. The value which is most possible is the maximum value (α = 1) and the value which is least possible is the minimum value (α = 0). Furthermore, the triangle membership function is suitable for simplicity to show the concept. In this example the aim is to model airspeed as a fuzzy membership function, not using interval or probabilistic models, to calculate reliability. This method is general, and other arbitrary airspeed membership functions that are more realistic could be used.
Figure 14. Air speed and flutter speed 3D interference. (a) Case 1, (b) Case 2, (c) Case 3, (d) Case 4, (e) Case 5, (f) Case 6.
In Case 1, it is assumed that the fuzzy air speed is (115,120,125) m/s. This means that the possibility of the air speed for values less than 115 m/s and more than 125 m/s is zero, and the possibility that the air speed is 120 m/s equals one. Based on this assumption the 3D interface between the air speed and the flutter speed is shown in Fig. 14(a). Using Equation (20) the fuzzy flutter safety can be obtained as
\begin{equation}{S_F} = \frac{{{V_s}}}{V} = 1 - \frac{{{V_f}}}{V} = 1 - \frac{{{\text{0}}{\text{.41826}}}}{{{\text{123}}{\text{.5848}}}} = 99.66\% .\end{equation}In Case 2, by keeping the crisp value unchanged, the air speed interval is assumed to be larger. In this case, it is assumed that the fuzzy air speed is (100,120,140) m/s. As shown in Fig. 14(b), the failure region interference expands and the flutter safety decreases. The flutter safety value is decreased to 93.91%, in this case.
In Case 3, the air speed intervals are the same as in Case 1, but the crisp value is increased. In this case, it is assumed that the fuzzy air speed is (150, 155,160) m/s, so that the air speed region is larger than the flutter speed in this case. The flutter safety decreases significantly and the flutter reliability reduces to 1.654%. The interference for this case is shown in Fig. 14(c). In Case 4, it is assumed that the fuzzy air speed is (105, 110,115) m/s. Figure 14(d) shows that the reliability value increases to 100 due to the lack of interference between the air speed and the flutter speed regions. In Case 5, as can be seen in Fig. 14(e), an asymmetric air speed region is considered. It is assumed that the fuzzy wind speed is (100,120,125) m/s. In this case the flutter safety value is 99.87%. Finally, in Case 6, it is assumed that the fuzzy wind speed is (115, 120,140) m/s. Figure 14(f) shows that the difference between the minimum values and the crisp value is less than the difference between the maximum value and the crisp value. In this case, the flutter safety value reduces to 90.23%.
Finally, the impact of the number of α-cuts on the accuracy of the reliability is studied. Table 3 and Fig. 15 shows the reliability versus the number of α-cuts for the six cases defined earlier. The relative error is defined as
\begin{equation}Err = \frac{|R_e^{i} - R_e^{j}|}{R_e^{j}}, \left\{\begin{array}{cl}
i& :\text{number of}\ \alpha\text{-cut}\\
j& :\text{Maximum number of}\ \alpha\text{-cut}\end{array}\right.\end{equation}The results show that, when the number of α-cuts more than 100, the error is almost zero. In this study, the maximum number of α-cuts is set equal to 1000 to guarantee the accuracy of the simulation.
To verify the results the aforementioned method is also compared with Monte Carlo simulation to verify the integration accuracy. The results are shown in Table 4. The results show that the method has high accuracy relative to the Monte Carlo method with low computational cost.
Table 3 Number of α-cuts and its impact on the accuracy of the flutter reliability
Figure 15. Reliability versus number of α-cuts. (a) Case 1, (b) Case 2, (c) Case 3, (d) Case 4, (e) Case 5, (f) Case 6.
Table 4 The reliability of typical section wing flutter speed with 1000 α-cuts
5.2 Example 2: Clean wing
In this example the reliability of a clean wing, which is a famous benchmark model in wing aeroelasticity, is considered. This model is a cantilever beam with bending and torsional deflection as shown in Fig. 16.
Figure 16. (a) Schematic of the clean wing, and (b) the typical section wing.
Suppose that the wing bending and torsional rigidity, wing mass per unit length and moment of inertia per unit length are considered as uncertain fuzzy parameters. The fuzzy triangle membership functions of the uncertain parameters are assumed to be 
$\overline m = (.95m,m, 1.05m)$, 
$\overline {{I_P}} = (.95{I_P},{I_P},1.05{I_P})$, 
$\overline {EI} = (.95EI,EI,1.05EI)$ and 
$\overline {GJ} = (.95GJ,GJ,1.05GJ)$.
Using Hamilton’s principle, the wing equations of motion are obtained as
\begin{equation}\quad EIw'''' + m\ddot w + m{y_\theta }\ddot \theta = L\end{equation}
\begin{equation} - GJ\theta '' + {I_P}\ddot \theta + m{y_\theta }\ddot w = M\end{equation}where EI and GJ are the bending and torsional rigidity, respectively, m is the wing mass per unit length and I P is the wing moment of inertia per unit length. Also, L and M are aerodynamic lift and moment, respectively.
The aerodynamic lift and moment based on Peters unsteady aerodynamic theory(Reference Peters, Karunamoorthy and Cao50) are used in Equations (30) and (31). The aeroelastic governing equations are given by:
\begin{equation}\begin{array}{l}\left[ {{m_{\left( x \right)}} + \pi \rho {b^2}} \right]\ddot w + \left[ {{m_{\left( x \right)}}{y_\theta } + \pi \rho a{b^3}} \right]\ddot \theta \\[4pt]
\qquad+ U\left[ {{C_{L\theta }}\rho b} \right]\dot w + U\left[ { - \pi \rho {b^2} + {C_{L\theta }}\rho a{b^2}\, - {C_{L\theta }}\rho \,\,\frac{{{b^2}}}{2}\left( {\frac{{{C_{L\theta }}}}{\pi } - 1} \right)} \right]\dot \theta \\[4pt]
\qquad+ \left[ {EI} \right]w'''' + U\left[ {{C_{L\theta }}\rho b\,} \right]{\lambda _0}\left( t \right) - {U^2}\left[ {{C_{L\theta }}\rho b} \right]\theta = 0\end{array}\end{equation}
\begin{equation}\begin{array}{l}\left[ {{m_{\left( x \right)}}k_{EA}^2 + \pi \rho {b^4}\left(\frac{1}{8} + {a^2}\right)} \right]\ddot \theta + \left[ {{m_{\left( x \right)}}{y_\theta } + \pi \rho a{b^3}} \right]\ddot w + U\left[ \pi \rho \frac{{{b^3}}}{2}\left( {\frac{{{C_{L\theta }}}}{\pi } - 1} \right) - \pi \rho a{b^3} \right.\\[10pt] \qquad\left. +\, {C_{L\theta }}\rho a{b^3}\left( {a + \frac{1}{2}} \right) - {{C_{L\theta }}\rho \frac{{{b^3}}}{2}\left( {\frac{{{C_{L\theta }}}}{\pi } - 1} \right)\left( {a + \frac{1}{2}} \right)} \right]\dot \theta + U\left[ {{C_{L\theta }}\rho {b^2}\left( {a + \frac{1}{2}} \right)} \right]\dot w \\[10pt] \qquad -\, GJ\theta '' + U\left[ {{C_{L\theta }}\rho {b^2}\left( {a + \frac{1}{2}} \right)} \right]{\lambda _0}\left( t \right) - {U^2}\left[ {{C_{L\theta }}\rho {b^2}\left( {a + \frac{1}{2}} \right)} \right]\theta = 0\end{array}\end{equation}By discretizing the above equations with the Galerkin method, and solving the final eigenvalue problem, the flutter airspeed and frequency based on the Goland wing parameters are given in Table 2. The results show that this 3D model with finite state unsteady aerodynamic loading is in good agreement with the literature. Furthermore, as expected, the results are more accurate than the 2D equivalent Goland wing typical section model.
By applying the fuzzy interval method described previously, the modal damping for each α-cut is obtained as:
\begin{equation}\begin{array}{l}\overline {{\gamma ^\alpha }} \left[ {\left( {EI} \right),\left( {GJ} \right),m,{I_P},U} \right] = {\gamma ^\alpha }\left[ {{{\left( {EI} \right)}_c},{{\left( {GJ} \right)}_c},{m_c},{I_{Pc}},U} \right]\\[10pt]
+ \left| {\frac{{\partial {\gamma ^\alpha }\left[ {{{\left( {EI} \right)}_c},{{\left( {GJ} \right)}_c},{m_c},{I_{Pc}},U} \right]}}{{\partial {{\left( {EI} \right)}^\alpha }}}} \right|\overline {\Delta \left( {EI} \right)_{}^\alpha } + \left| {\frac{{\partial {\gamma ^\alpha }\left[ {{{\left( {EI} \right)}_c},{{\left( {GJ} \right)}_c},{m_c},{I_{Pc}},,U} \right]}}{{\partial {{\left( {GJ} \right)}^\alpha }}}} \right|\overline {\Delta \left( {GJ} \right)_{}^\alpha }, \\[18pt]
+ \left| {\frac{{\partial {\gamma ^\alpha }\left[ {{{\left( {EI} \right)}_c},{{\left( {GJ} \right)}_c},{m_c},{I_{Pc}},U} \right]}}{{\partial {m^\alpha }}}} \right|\overline {\Delta {m_{(x)}}_{}^\alpha } + \left| {\frac{{\partial {\gamma ^\alpha }\left[ {{{\left( {EI} \right)}_c},{{\left( {GJ} \right)}_c},{m_c},{I_{Pc}},,U} \right]}}{{\partial {I_P}^\alpha }}} \right|\overline {\Delta {I_{P(x)}}_{}^\alpha } \end{array}\end{equation}
\begin{equation}\begin{array}{l}\underline {{\gamma ^\alpha }} \left[ {\left( {EI} \right),\left( {GJ} \right),m,{I_P},U} \right] = {\gamma ^\alpha }\left[ {{{\left( {EI} \right)}_c},{{\left( {GJ} \right)}_c},{m_c},{I_{Pc}},U} \right]\\[10pt]
+ \left| {\frac{{\partial {\gamma ^\alpha }\left[ {{{\left( {EI} \right)}_c},{{\left( {GJ} \right)}_c},{m_c},{I_{Pc}},U} \right]}}{{\partial {{\left( {EI} \right)}^\alpha }}}} \right|\underline {\Delta \left( {EI} \right)_{}^\alpha } + \left| {\frac{{\partial {\gamma ^\alpha }\left[ {{{\left( {EI} \right)}_c},{{\left( {GJ} \right)}_c},{m_c},{I_{Pc}},,U} \right]}}{{\partial {{\left( {GJ} \right)}^\alpha }}}} \right|\underline {\Delta \left( {GJ} \right)_{}^\alpha } \,\\[18pt]
+ \left| {\frac{{\partial {\gamma ^\alpha }\left[ {{{\left( {EI} \right)}_c},{{\left( {GJ} \right)}_c},{m_c},{I_{Pc}},U} \right]}}{{\partial {m^\alpha }}}} \right|\underline {\Delta m_{}^\alpha } + \left| {\frac{{\partial {\gamma ^\alpha }\left[ {{{\left( {EI} \right)}_c},{{\left( {GJ} \right)}_c},{m_c},{I_{Pc}},,U} \right]}}{{\partial {I_P}^\alpha }}} \right|\underline {\Delta {I_P}_{}^\alpha } \end{array}\end{equation}To obtain the flutter fuzzy membership function the dominant mode which can cause flutter is considered. In this example the first bending mode is dominant and the other modes are stable.
Using Equations (31) and (32) at each α-cut, the flutter airspeed membership function can be obtained, as shown in Fig. 17. The flutter boundary range can be seen as a triangle fuzzy mountain shape. Indeed, for each value of the modal damping and at every α-cut, the upper and lower bounds of the flutter speed can be extracted from this figure.
The reliability analysis for all cases which were studied in Example 1 is also carried out for this example and the results are given in Table 5. As expected, the results are very similar to the results of the previous example because of the similar properties of the two models. Obviously, since a more accurate aeroelastic model was used in this example, the results are more accurate than the typical section model results.
Table 5 The fuzzy reliability of wing flutter speed
Figure 17. Three-dimensional plot of modal damping versus airspeed for different α-cuts.
6.0 CONCLUSION
In this paper, a new method for the flutter speed reliability analysis of aircraft wings using a fuzzy interval approach is investigated. Uncertain parameters are modeled as fuzzy membership functions and propagated through the wing aeroelastic model. The flutter region at each α-cut is obtained through the P method. By combining the flutter intervals at each α-cut the fuzzy flutter speed membership function is obtained. The interference between the fuzzy flutter area and the air speed area forms a pyramid. The reliability or flutter safety is then obtained from this interference volume.
The prominent advantage of this method is that only membership functions of uncertain parameters are required and the other statistical characteristics or the probabilistic distribution densities are not needed for the reliability analysis. In order to illustrate the feasibility of this fuzzy reliability method, a typical section 2D wing model and a clean wing 3D model were used. For both examples, the reliability analysis is performed for six wind speed conditions and the flutter reliability was determined for each condition. The results show that this approach is a suitable method to predict wing flutter safety.
