Delft Aerospace Computational Science 0000 0000 0000 0000 0000 0000 0000 0000 0000 1111 1111 1111 1111 1111 1111 1111 1111 1111
000000000
000000000
000000000
000000000
000000000
111111111
111111111
111111111
111111111
111111111
ENGINEERING MECHANICSREPORT
March 2006Stochastic Flutter Analysis
C.V. Verhoosel, M.A. Guti´errez and S.J. Hulshoff
c.v.verhoosel@tudelft.nl
All rights reserved. This work may not be translated or copied in whole or in part without the written permission of the author.
TUD/LR/EM
Kluyverweg 1, 2629HS Delft, The Netherlands
P.O. Box 5058, 2600 GB Delft, The Netherlands
Phone: +31 15 278 5460
Fax: +31 15 261 1465
C.V. Verhoosel, M.A. Guti´errez and S.J. Hulshoff
Delft University of Technology, Faculty of Aerospace Engineering P.O. Box 5058, 2600 GB Delft, The Netherlands
ABSTRACT
The field of fluid-structure interaction is combined with the field of stochastics to perform a stochastic flut-ter analysis. Various methods to directly incorporate the effects of uncertainties in the flutflut-ter analysis are investigated. The panel problem with a supersonic fluid flowing over it is considered as a testcase.
The stochastic moments (mean, standard deviation, etc.) of the flutter point are computed by an uncer-tainty analysis. Sensitivity-based methods are used to determine the stochastic moments of the flutter point. This is done by implicit differentiation of the flutter requirement. The moments can also be determined using the spectral method, which can be considered as a projection method. An iterative solution to the general random eigenvalue problem is proposed for determining the spectral expansion of the flutter point. It turns out that the asymptotic method, which is a sensitivity-based method, is the most efficient method for approximating the moments of the flutter point. The success of this method can be explained by the fact that the relation between the random Mach number and the random field of elastic properties is close to linear.
The probability of the occurrence of flutter below a specified Mach number can be computed using a reliability analysis. Fully sampling-based techniques are not applicable, since the required sample size would be to large. The probability of failure can be approximated without the use of sampling by the use of the first-and second-order reliability method. The results of these methods are inaccurate however. The importance sampling method combines the speed of the reliability methods with the accuracy of the sampling techniques to obtain an efficient approximation of the probability of flutter.
Keywords and Phrases:stochastic finite element methods, fluid-structure interaction, uncertainty analysis, reliability analysis, random eigenvalue problem.
1. Introduction
Modern aircraft have to be both lightweight and slender in order to be competitive. This requirement makes them very flexible in general. As a consequence, the interaction between the aerodynamic forces acting on an aircraft and the deformations of that aircraft plays an important role in modern aircraft design. The field of problems associated with this interaction is in general referred to by the term aeroelasticity.
One of the most important aeroelastic phenomena is flutter. Flutter is the problem that arises when the aircraft deformation and the aerodynamic forces resonate. This can lead to undesirable dynamic response, or worse, in-flight structural failure. The main task of a flutter analysis is to determine the boundary between the domain at which an oscillation is damped and the domain for which this oscillation can grow without bound. This boundary is called the flutter boundary and a single point on this boundary is referred to as a flutter point.
Several methods to analyse flutter problems have been developed and are commonly used. None of these methods directly incorporates the effect of uncertainties. These uncertainties can be divided into two categories:
• Uncertainties in the structure
Safety factors are used to incorporate the effects of these uncertainties. Since the influence of these uncertainties is unknown, conservative choices for the safety factors are necessary. This leads to over-dimensioned structures, increased masses and as a consequence decreased economical efficiency.
A superior approach would be to directly incorporate the uncertainties in the computation, which is common practice in the field of structural mechanics. There are two approaches to incorporate the uncertainties in the computations. The first is to formulate the stochastic differential equations and solve them numerically. The second approach is to determine stochastic properties by making use of the deterministic problem as far as possible. The advantage of the second approach is that existing numerical techniques can be used. Stochastic finite element methods (SFEM) [13] are an example of the second approach to including uncertainties in the computation. In this report only the SFEM will be considered.
In general there are two types of stochastic analysis. First, the influence of the uncertainties on the stochastic moments (mean, standard deviation, etc.) of the flutter point can be investigated. This type of analysis is often referred to as uncertainty analysis. Second, the probability of the occurrence of some rare event (e.g. flutter) can be studied. This type of analysis is called reliability analysis.
Although commonly used in the field of structural mechanics, these stochastic analyses have not yet been applied to the flutter problem. In this report, the possibilities for applying stochastic methods to the flutter problem will be investigated. The goal is to get an overview of the methods that are available for stochastic flutter analysis. The methods will be compared based on their robustness, accuracy, speed and efficiency.
The 2-D panel problem at supersonic Mach number with a random field of elastic properties is used to investigate stochastic flutter. This problem is selected because on the one hand it contains most of the aspects that are encountered when performing a flutter analysis of an aircraft. On the other hand the deterministic problem remains relatively simple, such that one can focus on the stochastic analysis. An additional advantage of the 2-D panel problem is that it has been extensively investigated [22, 8], such that the deterministic results can easily be verified.
In Chapter 2 the physical problem of finding the flutter point will be transformed in a mathematical formulation. In Chapter 3 the deterministic solution of the flutter problem will be considered, since it is an essential part of most of the stochastic methods. After that the methods for uncertainty analysis will be considered. In Chapter 4 the sensitivity-based methods are explored. In Chapter 5, the spectral method will be explained. In this chapter also a new method to find the spectral solution of the general random eigenvalue problem is proposed. In Chapter 6 the various methods for reliability analysis will be explored. To be able to compare the various methods, numerical simulations will be performed in Chapter 7. Finally conclusions will be drawn in Chapter 8.
The computation of the statistical moments of the flutter point using a sensitivity-based method, as described in parts of Chapter 4 has been published in the conference proceedigs:
• M.A. Guti´errez, C.V. Verhoosel, and S.J. Hulshoff. Fluid-structure interaction analysis of a plate with a rondom field of elastic properties. In A. Augusti et al., editor, 9th International Conference on Structural Safety and Reliability, Rotterdam, 2005. Millpress
Furthermore, the iterative method for the computation of the spectral expansion of random eigenvalues has been published in:
• C.V. Verhoosel, M.A. Guti´errez, and S.J. Hulshoff. Iterative solution of the random eigenvalue problem with application to spectral stochastic finite element systems. International Journal for Numerical Methods in Engineering, 2006. in press DOI: 10.1002/nme.1712
2. Problem statement
M∞, a∞, ρ∞, p∞ h p∞ L y x w(x)
Figure 1: Schematic representation of the 2-D panel problem.
discused. Finally the mathematical formulation of the flutter problem is presented and an outline for the various available methods to solve this problem is given.
2.1 Panel problem
The panel problem consists of a plate with thickness h and length L that is fully clamped at both ends (Figure 1). The density of the panel is ρp, the modulus of elasticity is Epand the Poisson ratio
is indicated by νp. Due to the fact that a two-dimensional problem is considered, the modulus of
elasticity and the Poisson ratio can be combined in a single variable as E = Ep
1−ν2
p. A one dimensional
horizontal fluid is flowing over the plate at supersonic Mach number. The freestream Mach number of the fluid is M∞, the freestream speed of sound a∞, the freastream presssure p∞and the freestream
density ρ∞.
The modulus of elasticity of the panel is assumed to be uncertain. Therefore it is assumed not only to be a function of the spatial coordinate x, but also of a set of m random variables. The dependence of the modulus of elasticity on these variables will be discussed in the next section.
2.2 Discretisation of the random field of elastic properties
Material properties such as the modulus of elasticity are described using random fields. Where a random variable describes the random distribution of a scalar variable, a random field describes the random distribution of a (scalar) function. In the case of the panel problem, the modulus of elasticity could be defined by a random variable at infinitely many points. This would not be realistic, since there is a spatial restriction on the randomness. For example, suppose that the panel is considerably more stiff at a certain location, then the points that are in the neighbourhood of this point will also be stiffer. Due to this spatial correlation, the random modulus of elasticity of the panel can accurately be modelled using a finite set of random variables.
To model this spatial correlation, the modulus of elasticity is discretised using the Karhunen-Loeve expansion [10, ch. 2]. The K-L expansion is a normal random field that approximates an exact spatial correlation function. The random field for a panel with local mean modulus of elasticity µE and local
standard deviation σE and spatial correlation length Lc can be approximated using
E(x, ˜z) = µE+ σE ∞ X i=1 ˜ zigi(x) ≈ µE+ σE m X i=1 ˜ zigi(x), (2.1)
where {˜zi} (also written as ˜z) is a set of standard normal random variables. The ˜• is used to emphasise
that a variable is random. The exact spatial covariance function is taken to be Cov [x1, x2] = e−
|x1−x2|
Lc , (2.2)
which has extensively been used in a variety of fields. In this covariance function, Lcis the correlation
correlation length means that imperfections only influence the close neighbourhood of the imperfection. The functions {gi(x)} can be obtained by solving the eigenvalue problem (Appendix F)
Z L 2
−L2
e−|x1−x2|Lc fi(x2)dx2 = λifi(x1), for i = 1 . . . m. (2.3)
The functions {gi(x)} are then given by
gi(x) =
p
λifi(x), for i = 1 . . . m. (2.4)
In Ref. [17] it is proposed that in the case that the eigenfunctions gi can be computed exactly, the
K-L expansion is the most efficient method to discretise a random field. This means that given a certain accuracy, the random field is discretised using the minimum number of random variables. In this particular case the K-L discretisation is therefore the best choice. In the case that no exact solution for the eigenvalue problem would exist other discretisation methods like the Optimal Linear Estimation Method (OLEM) [17] are more efficient.
2.3 Equations of motion
In this section the discrete dynamical system will be derived. First Hamilton’s principle will be used to relate the motion of the panel to its energy. Then the aerodynamic operator will be discussed and a finite element (FE) discretisation will be used to obtain generalised coordinates. Finally the Euler-Lagrange equation will be used to obtain a discrete dynamical system.
2.3.1 Hamilton’s principle
To study the stability of the panel, the equations of motion have to be derived first. According to Hamilton’s principle, the motion of the panel in the time interval t ∈ R[t1, t2] is defined by the
variational problem
Z t2
t1
(δT − δU + δW ) dt = 0, (2.5)
where T is the kinetic energy of the plate, U is the potential energy and δW is the virtual work performed by the aerodynamic forces.
The potential energy per unit length of the panel can be written as U (t, ˜z) = Z L 0 1 2 h3 12E(x, ˜z) ∂2w ∂x2 2 dx (2.6)
and the kinetic energy per unit length is equal to T (t) = Z L 0 1 2ρph ∂w ∂t 2 dx. (2.7)
The work done by the aerodynamic forces can be written as δW (t) = − Z L 0 ∆pδwdx = − Z L 0 (pu− pl) δwdx, (2.8)
where the pressure difference ∆p(x, t) has to be described by some aerodynamic operator. 2.3.2 Aerodynamic operator
The aerodynamic forces in the case of supersonic flow can be approximated using piston theory. This theory is valid for high supersonic Mach number flows, typically in the range from 2 to 5. The pressure variation at the upper side of the panel is approximated using
x
1 2 3 Ne− 2 Ne− 1
Figure 2: Schematic representation of the nodes and elements.
On the lower side of the panel, the pressure remains equal to the freestream pressure. The local pressure difference between the upper and lower side of the panel then equals
∆p(x, t) = pu− pl= ρ∞a∞ ∂w ∂t + v∞ ∂w ∂x . (2.10)
Since the focus in this report is not on the aerodynamic operator, the choice for piston theory is justified. A significant improvement of the aerodynamic model can be obtained by considering the higher-order piston theory [9, 22, 8]. The lower bound for the Mach number at which this theory can be applied is approximately 1.6. The pressure variation at the upper side of the panel using the higher-order piston theory is predicted as
p′ u(x, t) = ρ∞a∞ M∞ p M2 ∞− 1 M2 ∞− 2 M2 ∞− 1 ∂w ∂t + v∞ ∂w ∂x , (2.11)
which reduces to first-order piston theory if M∞≫ 1. In this study first-order piston theory is used.
Despite the non-linearity in higher-order piston theory, all the considered methods can be adjusted to incorporate this aerodynamic operator.
2.3.3 Spatial finite element discretisation
The integrals (2.6), (2.7) and (2.8) can be approximated using a finite element discretisation. The panel is divided in Neelements with length Le= NLe (Figure 2). Since the displacement and rotation
of the clamps is known, only the internal nodes have to be considered. Since second-order derivatives with respect to x need to be evaluated, it is required that w(x) satisfies
w(x) ∈ H2[0, L]. (2.12)
To satisfy this requirement, the element shape functions should also be in H2[0, L]. In 1-D the
H2 requirement is equivalent with the requirement that the first derivative should be continuously
differentiable. The Hermitian shape functions are therefore used (Figure 3):
φ1 = 1 − 3η2+ 2η3; (2.13a) φ2 = Le η − 2η2+ η3; (2.13b) φ3 = 3η2− 2η3; (2.13c) φ4 = Le −η2+ η3 , (2.13d) with η = x¯ Le = x−x1
Le . The vertical displacement for a given element is then approximated as
1 2 x Le ¯ x w(x) x1 x2 w1 w2 ∂w ∂x ˛ ˛1 ∂w ∂x ˛ ˛ 2
Figure 3: Definition of the elements.
Substituting the approximate solution in the equations (2.6), (2.7) and (2.8) gives the energies for a single element as Ue = 1 24h 3Z Le 0 E(x, ˜z) 4 X i=1 d2φ i dx2 ui !2 d¯x; (2.15a) Te = 1 2ρph Z Le 0 4 X i=1 φi dui dt !2 d¯x; (2.15b) δWe = −ρ∞a∞ Z Le 0 4 X i=1 φi dui dt + v∞ dφi dxui ! 4 X i=1 φiδui ! d¯x. (2.15c)
According to Hamilton’s principle, the energies should satisfy the Euler-Lagrange equation d dt ∂Te ∂ ˙uj −∂T∂ue j +∂Ue ∂uj = δWe δuj , for j = 1 . . . 4. (2.16)
The element matrices can then easily be evaluated with Me(j, i) = ρph Z Le 0 φiφjd¯x; (2.18a) De(j, i) = ρ∞a∞ Z Le 0 φiφjd¯x; (2.18b) ˜ Ke(j, i) = 1 12h 3Z Le 0 E(x, ˜z)d 2φ i dx2 d2φ j dx2 d¯x | {z } ˜ Kelas e (j,i) + ρ∞a2∞M∞ Z L 0 φjdφi dxd¯x | {z } Kaero e (j,i) . (2.18c)
When all the element matrices are evaluated, these matrices can be assembled to yield
M¨u+ D ˙u + K (˜z) u = 0, (2.19) in which u=h w1 ∂w∂x 1 w2 ∂w ∂x 2 · · · wNe−1 ∂w ∂x Ne−1 iT . (2.20)
The number of degrees of freedom obviously equals N = 2(Ne− 1).
As can be seen in equation (2.18) the damping matrix can be written as
D= 2ξM (2.21) with ξ = 1 2 a∞ h ρ∞ ρp . (2.22)
Furthermore it can be seen that the stiffness matrix K has a symmetric part Kelas depending on the
modulus of elasticity E(x, ˜z) and a non-symmetric part Kaero which depends on the Mach number. It is this non-symmetric part that allows the system to become undamped.
2.4 Stability analysis
In this section, the stability of the panel will be investigated. First, asymptotic stability analysis in the stochastic case will be explained. After that, attention will be paid to the transient motion of the panel.
2.4.1 Asymptotic stability
To be able to determine the asymptotic stability of the panel, a harmonic solution is assumed u(x, t, ˜z) = N X i=1 ˆ ui(˜z)epi(˜z)t. (2.23)
The functions pi(˜z) can be found by substituting this solution in equation (2.19), which defines the
spectrum P as P:=p(˜z)| detp˜2M+ ˜pD + K(˜z)= 0 , (2.24) with K(˜z) = K0(M∞, µE) + m X i=1 ˜ ziKi(σE, Lc). (2.25)
Since D = 2ξM, this can be interpreted as a random eigenvalue problem. Define
S(M∞, ˜z) = M−1K(M∞, ˜z) , (2.26)
to write
In the case that an arbitrary initial state is considered, the motion of the plate will only be asymp-totically stable if
pr(M∞, ˜z) ≤ 0, ∀˜p ∈ P (M∞, ˜z) . (2.28)
with pr = Re [p]. Obviously the stability of the plate depends on the random variables. From this
stability requirement it is obvious that the flutter point will also be a function of the random variables. Another possibility is to formulate the stability requirement in terms of the spectrum of S, which is defined as
L(M∞, ˜z) :=nλ(˜z)| dethS(M∞, ˜z) − ˜λIi= 0o. (2.29) The stability criterion (2.28) can be written in terms of the eigenvalues of S as
λr(M∞, ˜z) ≥
λˆı(M∞, ˜z)2
4ξ2 , ∀˜λ ∈ L . (2.30)
with λr= Re [λ] and λˆı= Im [λ].
2.4.2 Transient stability
In the previous section the stability criterion is based on the asymptotic stability of the panel. In the case that the governing operator would be normal, which would be equivalent with a symmetric matrix S(since the basis of eigenvectors is then orthogonal), this criterion would imply absolute stability.
Since the matrix S as used in the previous section is non-symmetric and non-normal, the transient behaviour of the panel can be irregular. If the solution is asymptotically stable, the solution might initially grow. The problem then is that failure can occur before the solution reaches its asymptote. Checking the transient behaviour is therefore necessary to ensure that it is realistic to considerer asymptotic stability.
The transient stability of the panel for any initial state with finite energy can be investigated. This is done by rewriting the second-order system (2.19) as a twice as large first-order system
d dt u ˙u = D M I 0 −1 −K 0 0 I | {z } A u ˙u | {z } ¯ u . (2.31)
The solution of this system is
¯
u= etAu¯0, (2.32)
where the subscript 0 indicates the initial state. To look at the stability of the panel, consider the energy norm of the panel as defined by the inner product
h¯v, ¯wiE = 1 2 v ˙v T Kelas 0 0 M w ˙ w . (2.33)
The energy norm is then defined as
k¯ukE =
q
h¯u, ¯uiE. (2.34)
The vector space U0 is defined as the space containing all initial vectors which are bounded with
respect to this norm. The norm of the first-order system (2.32) equals k¯ukE =
etAu¯0 E, ∀¯u0∈ U0, (2.35)
of which the upper bound can be approximated as [20, p. 408] k¯ukE ≤
etA Ek¯u0kE, ∀¯u0∈ U0. (2.36)
If the initial energy of the panel is not equal to 0. The upper bound of the energy of the panel relative to the initial energy can be computed as
ρ = sup t>0 k¯ukE k¯u0kE ≤ sup t>0 etA E, ∀¯u0∈ U0. (2.37)
2.5 Flutter point
The flutter point is defined as the Mach number for which the oscillation of the panel becomes unbounded in time. The asymptotic stability is in general an appropriate measure to determine if the motion of the panel is damped. To be sure that the flutter point is physically relevant, the transient behaviour of the panel should be checked.
2.5.1 Random flutter problem
The flutter problem can be defined in terms of the reduced frequencies as Find the smallest Mf l(˜z)|∃p ∈ P (Mf l, ˜z) ,
satisfying pr(Mf l, ˜z) = 0, ∀˜z ∈ Rm.
(2.38)
with
P :=p(˜z)| detp˜2M+ ˜pD + K= 0 . (2.39) Where the notion smallest refers to the smallest expected value. Another option is to define it in terms of the spectrum of the matrix ˜S, which yields
Find the smallest Mf l(˜z)|∃λ ∈ L (Mf l, ˜z) ,
satisfying λr(Mf l, ˜z) = λi(Mf l, ˜ z)2 4ξ2 , ∀˜z ∈ R m. (2.40) with L :=nλ(˜z)| dethS(M∞, ˜z) − ˜λIi= 0o. (2.41) Both definitions are equivalent, but sometimes the first definition is most convenient (e.g. deterministic approach) and in other cases the second definition is more convenient (e.g. spectral approach).
The solution of the stochastic flutter problem essentially consists of two parts: Solving the random general eigenvalue problem (2.41) and applying the stochastic flutter requirement (2.40). The problem is that the eigenvalues and eigenvectors are now a function of the random variables. These functions can in general not be found analytically, so approximations of these functions need to be made. The application of the flutter requirement also becomes more complicated, since functions instead of scalars are involved.
2.5.2 Solutions of the flutter problem
The stochastic flutter problem has now been formulated. Before the stochastic solution methods for this problem are described, the deterministic problem is given some attention. The reason for this is that the deterministic problem is a crucial part of most of the stochastic methods that are considered in this study.
When the deterministic solution procedure is described, the stochastic methods are considered. Since analytic solutions are in general not available, approximate solutions need to be found. The two different types of approximation will now be discussed.
Uncertainty analysis The purpose of uncertainty analysis is to obtain the stochastic moments of the flutter point. This type of analysis can therefore be considered as an approximation of the centre of the distribution (around ˜z≈ 0). As mentioned above, the random eigenvalue problem needs to be solved as a part of the solution procedure for the flutter problem. Therefore various methods to solve this random eigenvalue problem will now briefly be discussed. It will be determined whether these methods are also applicable to the flutter problem. An overview of stochastic methods that are so far available for the solution of the discrete random eigenvalue problem can be found in Ref. [24, 2].
Hessian of the eigenvalues. The connection with the flutter problem is made by application of the chain rule. Application of gradient-based methods will be explained in Chapter 4.
In Ref. [24] it is suggested that in some cases, the eigenvalues and eigenvectors can be considered as being independent. This assumption would simplify the problem considerably. But unfortunately it does not hold in the case of the problem considered in this report.
In Ref. [24] methods are described to construct the probability density functions of the random eigenvalues. The most important method is the crossing-theory method. Although this method would be applicable for constructing the probability density function of the eigenvalues, the application for the solution of the random flutter problem is not possible. The problem is that the out-crossing rate of the eigenvalues can be computed, but that no relation with the out-crossing rate of the flutter point can be found. Application of this method in the context of reliability analysis is theoretically possible. But since the method focuses on the centre of the distribution instead of on the tails, this application seems to be irrelevant.
A method that is not discussed in the above-mentioned references is the spectral method, which is basically a projection on a polynomial space. Attempts have been made to solve the random eigenvalue problem by using the spectral method [11], but with limited success. In this report a solution using the spectral method is iteratively obtained. This method is described in Chapter 5.
Reliability analysis The purpose of reliability analysis is to approximate the probability of occurrence of some rare event as accurately as possible. Reliability analysis can be considered as an approximation of the tail of the distribution (|˜z| ≫ 1). A good overview of the methods available for reliability analysis is given in Ref. [19]. The methods described in this reference are applicable to the flutter problem.
Basically three types of methods are available. The first type is the fully sampling-based Monte-Carlo method. The second type include the first- and second-order reliability method which do not require sampling at all. Finally importance sampling is considered as a hybrid method. Reliability methods are used to increase the efficiency of the Monte-Carlo method. The methods available for reliability analysis will be described in Chapter 6.
3. Deterministic problem
Before the stochastic flutter problem is considered, the solution procedure for the deterministic case is studied. This is necessary since the deterministic solution procedure is a part of most of the stochastic solution procedures (sensitivity-based methods, Monte-Carlo simulation). Furthermore the deterministic solution will be used to investigate the spectra of the eigenvalues, to get more insight in the relation between the random eigenvalue problem and the flutter problem.
3.1 Deterministic flutter problem
In the case that there is no uncertainty in the problem, the modulus of elasticity of the panel is simply a positive constant (E ∈ R+). The flutter problem (2.38) can then be rewritten as
Find the smallest Mf l∈ R+|∃p ∈ P (Mf l) ,
satisfying pr(Mf l) = 0.
(3.1)
Where the spectrum P is defined as
P(M∞) :=p ∈ C ∃u ∈ CN Su= − p2+ 2ξpu . (3.2) The problem basicly consists of two parts: The deterministic general eigenvalue problem (3.2) and the flutter requirement (3.1). The solution procedure for the deterministic flutter point is schematically shown in Figure 4.
To work properly, first the mode that causes flutter has to be determined. For physically relevant settings, flutter is caused by the mode corresponding to the smallest eigenvalue. Typically the solution of the iteration process will converge in less than 5 steps.
In Figure 5 the flutter diagram of the first three modes of the setting as used in Ref. [22] is shown. As can be seen, flutter is indeed caused by the lowest mode. The first four mode shapes are shown in Figure 6.
The transient stability is investigated using equation (2.37). The upper boundary of the energy of the panel for the settings as used in Ref. [22] is shown in Figure 7. It turns out that ρ = 21.6, which means that for the worst possible initial state, the energy is increased 26.1 times. This number puts a restriction on the (initial) disturbances of the panel that can be considered to produce a physically relevant flutter point.
3.2 Eigenvalue spectrum
From the general eigenvalue problem is it obvious that the eigenvalues of ˜Sare dependent on ˜z. Before trying to find approximations of these functions, the spectra can be generated by doing the eigenvalue computation for several values of ˜z. If a certain ˜z is selected, the problem can be considered to be deterministic.
The root-locus plot [21, ch. 6] for the setting used above for the first 4 random variables is shown in Figure 8. The variation coefficient of the modulus of elasticity is 0.1 and the correlation length is 20 percent of the panel length. The results for other physically relevant settings are similar. The figure shows the eigenvalue spectrum for the first 2 modes. The higher-order modes behave the same as the second mode, but the corresponding eigenvalue is larger. The root-locus plots show the variation of the spectrum on ˜zfor the range −6 to 6, which basically means that the probability of an eigenvalue outside the plotted range in the corresponding stochastic case is approximately one in a billion.
Since flutter can only occur if the eigenvalue is complex, it is observed that flutter is always caused by the first mode. Furthermore it can be seen that the eigenvalues are most sensitive to the first random variable and that the sensitivity decreases with increasing random variable number. Finally the root-locus plots also indicate that the eigenvalues are distinct, which is a necessary requirement for most of the stochastic methods that will be considered in the following chapters.
It should be noted that this root-locus plot is basically a special case of the pseudospectrum. In the case of a pseudospectrum the sensitivity of a matrix to an arbitrary disturbance is shown. In this case the disturbance is random but not arbitrary, since the disturbance matrices have a prescribed form.
4. Uncertainty analysis using sensitivity-based methods
The influence of uncertainties in the input parameters can be investigated by computing the corre-sponding sensitivities. Stochastic moments can be obtained using these sensitivities. In this chapter two methods that are based on the sensitivities will be discussed. First the perturbation method will be investigated. After that, the asymptotic method will be explained.
4.1 Sensitivity computation
Since all methods discussed in this chapter are based on the sensitivities of the flutter point, the computation of these sensitivities will first be discussed. Since for the case considered flutter is in general caused by the first mode, the eigenvalue and frequency corresponding to the first mode flutter point are respectively indicated by λ and p. The subscript i is dropped for convenience.
4.1.1 Gradient of the flutter point
The gradient operator for a scalar function of the variables {xi}mi=1 is defined as
Input parameters. Output: Mf l. True False Secant iteration to find Mf l.
Construct flutter diagram
Check flutter requirement: pr≥ 0.
Compute {pi}.
Solve the eigenval. problem (3.2). Select a value for Mf l.
2 2.5 3 3.5 4 4.5 5 0 200 400 600 800 1000 1200 2 2.5 3 3.5 4 4.5 5 −200 −100 0 100 200 M∞ pr |pi | Mode 1 Mode 2 Mode 3
Figure 5: Example of a deterministic flutter diagram with Ne= 4.
0 0 0 0 x L L L L M o d e 1 M o d e 2 M o d e 3 M o d e 4
0 1 2 3 4 5 6 7 8 9 10 10−1 100 101 102 103 t k ¯u kE / k ¯u0 kE Lower bound Upper bound
Figure 7: Example of the transient stability at the flutter point with Ne= 4.
0 2 4 6 8 10 12 14 16 18 x 105 −2 0 2x 10 5 0 2 4 6 8 10 12 14 16 18 x 105 −2 0 2x 10 5 0 2 4 6 8 10 12 14 16 18 x 105 −2 0 2x 10 5 0 2 4 6 8 10 12 14 16 18 x 105 −2 0 2x 10 5 M o d e 1 M o d e 2 M o d e 3 M o d e 4 λˆı λˆı λˆı λˆı λr
The Mach number at which flutter occurs is a function of the input parameters. Since the gradient with respect to the Gaussian random variables is required the other parameters are considered to be fixed. The flutter point can then be considered as only depending on these variables
Mf l: Rm7→ R+. (4.2)
This function is implicitly known by the flutter requirement
pr(Mf l, ˜z1, · · · , ˜zm) = 0. (4.3)
The fact that the flutter point is only known implicitly is not a problem for the computation of the gradient, since the chain rule can be applied. Taking the gradient of equation (4.3) gives
∂pr
∂Mf l∇
zMf l+ ∇zpr= 0, (4.4)
which can be rewritten to yield
∇zMf l= − ∂pr
∂Mf l
−1
∇zpr. (4.5)
The partial derivatives of pr with respect to a real variable θ can be found using
∂pr ∂θ = Re ∂p ∂θ . (4.6)
This derivative can be evaluated using
p2+ 2pξ + λ = 0, (4.7) with ξ = 1 2 a∞ h ρ∞ ρm . (4.8)
Application of the chain rule yields 2p∂p ∂θ+ 2ξ ∂p ∂θ + 2p ∂ξ ∂θ+ ∂λ ∂θ = 0, (4.9)
which can be rewritten as
∂p ∂θ = − 1 2 (p + ξ) 2p∂ξ ∂θ+ ∂λ ∂θ . (4.10)
In the case that θ = Mf l this gives
∂p ∂Mf l = − 1 2 (p + ξ) ∂λ ∂Mf l . (4.11) And if θ = zi it results in ∂p ∂zi = − 1 2 (p + ξ) ∂λ ∂zi. (4.12)
The gradient of Mf lcan be determined when the eigenvalue sensitivities are known. The computation
4.1.2 Hessian of the flutter point
The Hessian operator for a scalar function of the variables {xi}mi=1 is defined as
∇(2)x (•) := ∂2(•) ∂x2 1 · · · ∂2(•) ∂x1∂xm .. . . .. ... ∂2(•) ∂xm∂x1 · · · ∂2(•) ∂x2 m . (4.13)
The Hessian of Mf lcan then be computed by taking the gradient of (4.4) to yield
∇Tz ∂pr ∂Mf l ∇zMf l+∂ 2p r ∂M2 f l ∇TzMf l∇zMf l+∂pr Mf l∇ (2) z Mf l+∇(2)z pr+∇TzMf l ∂ ∂Mf l(∇ zpr) = 0. (4.14)
This can be rewritten as ∇(2)z Mf l= ∂pr Mf l −1 ∇Tz ∂pr ∂Mf l ∇zMf l+ ∂2p r ∂M2 f l ∇TzMf l∇zMf l+ +∇(2) z pr+ ∇TzMf l ∂ ∂Mf l(∇ zpr) . (4.15)
The second-order derivatives of p in this expression can be computed by differentiation of equation (4.9) 2∂p ∂φ ∂p ∂θ+ 2p ∂2p ∂θ∂φ+ 2 ∂ξ ∂φ ∂p ∂θ + 2ξ ∂2p ∂θ∂ξ + 2 ∂p ∂φ ∂ξ ∂θ+ 2p ∂2ξ ∂θ∂φ+ ∂2λ ∂θ∂φ. (4.16)
This expression can be rewritten to yield ∂2p ∂θ∂φ = − 1 2(p + ξ) 2∂p ∂φ ∂p ∂θ+ 2 ∂ξ ∂φ ∂p ∂θ+ 2 ∂p ∂φ ∂ξ ∂θ+ 2p ∂2ξ ∂θ∂φ+ ∂2λ ∂θ∂φ . (4.17)
In the case that both θ = Mf land φ = Mf lthis gives
∂2p ∂2M f l = − 1 2(p + ξ) 2 ∂p ∂Mf l 2 + ∂ 2λ ∂2M f l ! . (4.18)
In the case that θ = Mf l and φ = zi it is found that
∂2p ∂Mf l∂zi = − 1 2(p + ξ) 2∂p ∂zi ∂p ∂Mf l+ ∂2λ ∂Mf l∂zi . (4.19)
And finally in the case that θ = zi and φ = zj it is found that
∂2p ∂zi∂zj = − 1 2(p + ξ) 2∂p ∂zj ∂p ∂zi + ∂ 2λ ∂zi∂zj . (4.20)
The Hessian of the flutter point can then be evaluated as soon as the second-order derivatives of the eigenvalues are known. The computation of these derivatives is explained in Appendix C.
4.2 Mean-centred perturbation method
Let ˜Mf lbe the random flutter Mach number that depends on the set of uncorrelated normal random
variables {˜zi}m1 with mean 0 and standard deviation 1. The function Mf l: Rm7→ R can be expanded
around the point 0 as
Mf l(˜z) = Mf l(0) + ∇zMf l(0)˜z+
1 2˜z
T
The second-order approximation for the expectation of ˜Mf lis then obtained by taking the expectation
and neglecting higher-order terms µMf l= E h ˜ Mf l i ≈ Mf l(0) + 1 2tr ∇(2)z Mf l(0) . (4.22)
And the first-order approximation of the variance is equal to σM2f l ≈ ∇zMf l(0)
T
∇zMf l(0). (4.23)
In the case of Gaussian random variables, the theory of quadratic forms can be used to compute arbitrary moments [16, ch. 29] µ1Mf l = Mf l(0) + 1 2tr ∇(2)z Mf l(0) , (4.24a) µrMf l = E h ˜ Mf l− µMf l ri = r! 2∇zMf l(0) h ∇(2)z Mf l(0) ir−2 ∇zMf l(0)T +(r − 1)! 2 tr h ∇(2)z Mf l(0) ir , r ≥ 2.(4.24b) 4.3 α-centred perturbation method
The expansion of the function ˜Mf laround the deterministic value is not necessarily the most accurate
one. An expansion around z = α that is optimal in a certain sense can give better results [2] Mf l(˜z) = Mf l(α) + [∇zMf l(α)] (˜z− α) + 1 2(˜z− α) Th ∇(2)z Mf l(α) i (˜z− α) + O(˜z− α)3. (4.25) The vector α should be chosen such that the mean of ˜Mf lis computed most accurately. To determine
this point, consider the integral expression for the mean µMf l= Z Rm ˜ Mf lpzdz = Z Rm e−h(z)dz, (4.26) with h(z) = ln ˜Mf l− 1 2z Tz. (4.27)
The largest constribution to the integral expression is obtained when h is at its global minimum. Expansion of the ˜Mf laround that point gives the best approximation for the first moment. The point
αcan then be found by
∂h(α) ∂zi
= 0, for i = 1 . . . m. (4.28)
The system of m equations and m unknowns can then be solved iteratively using αk+1i = 1 ˜ Mf l(αk) ∂ ˜Mf l(αk) ∂zi , for i = 1 . . . m. (4.29)
Which results in µ1Mf l = a + 1 2tr ∇(2)z Mf l(α) , (4.32a) µrMf l = E h ˜ Mf l− µMf l ri = r! 2b h ∇(2)z Mf l(α) ir−2 bT +(r − 1)! 2 tr h ∇(2)z Mf l(α) ir , for r ≥ 2. (4.32b) 4.4 Asymptotic method
The asymptotic integral approximation can be used to compute the stochastic moments of the flutter point [2]. The probability density function of z is given by
pz(z) = 1 (2π)m/2e −1 2z Tz . (4.33)
The mean value of the flutter point can then be obtained by solving the integral E [Mf l] = 1 (2π)m/2 Z Rm Mf l(z)e− 1 2z Tz dz, (4.34)
which can be rewritten as
E [Mf l] = 1 (2π)m/2 Z Rm e−(12z Tz−ln (M f l(z)))dz. (4.35) Now let h(z) =1 2z Tz− ln (M f l(z)), (4.36)
then the integral can be approximated by (Appendix D) µMf l= E [Mf l] ≈ e
−h(θ)|∇(2) z h(θ)|−
1
2, (4.37)
in which θ represents the global minimum that can be found using ∇zh = θ −
1 Mf l(θ)∇
zMf l(θ) = 0. (4.38)
This can be rewritten as
θk+1= 1
Mf l(θk)
∇zMf l(θk), (4.39)
which can be solved iteratively. The Hessian of h(z) can be computed as ∇(2)z h = I − 1 Mf l(z)∇ (2) z Mf l+ 1 Mf l(z)2∇ T zMf l∇zMf l. (4.40)
The higher-order moments can best be computed by first computing the raw moments and then applying a binomial transform to find the corresponding central moments. The reason for this is that the assumption of small higher-order terms in the Taylor expansion of h is then satisfied.
The higher-order raw moments are given by the integral EMf lr = 1 (2π)m/2 Z Rm Mf l(z)re− 1 2z Tz dz. (4.41)
The moments can be found by computation of the Hessian of h(z) =1
2z
Tz
which is equal to ∇(2)z h = I − r Mf l(z)∇ (2) z Mf l+ r Mf l(z)2∇ T zMf l∇zMf l. (4.43)
This Hessian matrix needs to be evaluated at the global minimum of h, which can be found by solving
θk+1= r
Mf l(θk)
∇zMf l(θk), (4.44)
iteratively. The raw moments can then be obtained using EMf lr
≈ e−h(θ)|∇(2)z h(θ)|−
1
2. (4.45)
Once the raw moments are known, the central moments can easily be computed using the binomial transformation E Mf l− µMf l r = r X k=0 r k (−1)r−kEMf lk µr−kMf l. (4.46)
5. Uncertainty analysis using the spectral method
The influence of uncertainties on the flutter point can also be modelled using a projection method. This method is often referred to as the spectral method. In this chapter the spectral solution of the random eigenvalue problem will first be considered. After that, the calculation of the stochastic flutter point using the spectral method will be discussed.
5.1 Random eigenvalue problem
Solving the random eigenvalue directly using the spectral method is not possible, therefore effort is made to solve the problem iteratively. Recently a new method was proposed for solving the random eigenvalue problem iteratively using polynomial chaos expansion [11]. This method will be discussed briefly and after that a more efficient method will be suggested.
5.1.1 Non-linear system iteration
The idea proposed in Ref. [11] is to apply the spectral method to the random eigenvalue problem ˜
Su˜ = ˜λ˜u, (5.1)
with ˜Sbeing a random symmetric matrix, in exactly the same way as is done for a linear system of equations. The result is a non-linear system of N ×(n+1) independent equations and (N +1)×(n+1) unknowns, with n being the number of chaos polynomials in the expansion. In order to solve this system, n + 1 additional equations need to be obtained. These equations are obtained by using the fact that the norm of the eigenvectors is constant. The result is a non-linear system of (N +1)×(n+1) equations and (N + 1) × (n + 1) unknowns. This non-linear system can be solved using an iterative solver, e.g. Newton-Raphson. In order to ensure convergence, an apppropriate initial estimate for the spectral expansion of both the eigenvalue and eigenvector is required. In general, sampling is demanded to obtain these initial estimates.
5.1.2 Stochastic inverse power method
In Ref. [11] it is chosen to solve the random eigenvalue problem with the spectral method by con-structing a large system of non-linear equations and solving this system iteratively. Deterministic eigenvalues can also be obtained iteratively [6, ch. 9]. Such an iterative method can be adjusted to the stochastic case.
Since the deterministic eigenvalues are known and variances are typically moderate, a good initial guess of the random eigenvalues is known. In the deterministic case, a good estimate of the eigenvalues would mean that the exact eigenvalues can be found efficiently by application of the inverse power method. Therefore it makes sense to adjust this deterministic inverse power method such that it also works for random eigenvalue problems.
Inverse power method for the deterministic general eigenvalue problem In Appendix E it is demon-strated that the inverse power method is able to efficiently determine an eigenvalue and corresponding eigenvector of a general matrix if:
• There is a large separation between the eigenvalues. • There is a good approximation of the desired eigenvalue.
Since both conditions are satisfied for the problem considered in this report, the inverse power method can be used to determine the eigenvalue and eigenvector that cause flutter. To explain how this method works, consider the eigenvalue problem
Su= λu. (5.2)
To be able to find a eigenvalue near q, the modified eigenvalue problem
[S − qI] v = φv (5.3)
can be solved. This modified eigenvalue problem has the same eigenvectors as the original problem (5.2) and has eigenvalues
φ = λ − q. (5.4)
Since the matrix S is non-symmetric, the eigenvalues and eigenvectors are in general complex. For convenience the problem is split up into a real and an imaginary part as
[S − (qr+ qˆıˆı) I] (vr+ vˆıˆı) = (φr+ φˆıˆı) (vr+ vˆıˆı) , (5.5)
which can be rewritten as S− qrI qˆıI −qˆıI S− qrI vr vˆı = φrI −φˆıI φˆıI φrI vr vˆı . (5.6)
To find the eigenvalue φ near 0 and the corresponding eigenvector v, this equation can be solved iteratively. If a guess for a eigenvalue and corresponding eigenvector is given by
φk = φkr+ φˆkıˆı; (5.7a)
vk = vk
r+ vˆıkˆı, (5.7b)
the eigenvector can be updated using vnew r vnew ˆ ı = S − qrI qˆıI −qˆıI S− qrI −1 φk rI −φˆkıI φk ˆ ıI φkrI vk r vk ˆı . (5.8)
Since an eigenvector can be multiplied with an arbitrary complex constant, the angle in the complex plane of each of the eigenvector components is not fixed. From a convergence point of view, it turns out to be favourable if this angle is fixed. This is done by multiplication of one of the non-zero entries (index α) with its complex conjugate
vnew→ vnew
α vnew. (5.9)
For convenience, the updated eigenvector can be normalised using vnew→ v
new
kvnewk. (5.10)
Using the updated eigenvector vnew, the corresponding eigenvalue φnew can be computed by
φnew= v
k+1HA vk+1
Then a new initial guess can be made for φ and v using
vk+1 = (1 − κ)vk+ κvnew; (5.12a)
φk+1 = φnew, (5.12b)
with
κ ∈ R(0, 1]. (5.13)
The parameter κ can be chosen such that the method converges. If κ is chosen to be equal to 1, the eigenvector is fully updated. Often this choice for κ leads to divergence of the method. Convergence can be obtained by choosing a lower value of κ. In that case the eigenvector is only partially updated. That the method still converges to the correct solution can be proven by substituting equation (5.8) in equation (5.12a), which gives
vk+1 r vk+1 ˆ ı = (1 − κ) vrk vˆık + κ S− qrI qˆıI −qˆıI S− qrI −1 φk rI −φˆkıI φk ˆıI φkrI vkr vkˆı . (5.14)
Once the method has converged it holds that vk+1 r vk+1 ˆı = vkr vk ˆ ı , (5.15) which leads to −κ vr vˆı + κ S− qrI qˆıI −qˆıI S− qrI −1 φrI −φˆıI φˆıI φrI vr vˆı = 0. (5.16)
This can be rewritten as κ S− qrI qˆıI −qˆıI S− qrI vr vˆı = κ φrI −φˆıI φˆıI φrI vr vˆı , (5.17)
which is the same as the exact problem (5.6) if κ is not equal to 0.
Once the eigenvector v and corresponding eigenvalue φ have converged, equation (5.4) can be used to compute the eigenvalue of the original eigenvalue problem (5.2) by
λ = φ + q. (5.18)
In Appendix E it is demonstrated that the inverse power method converges to the desired eigenvalue and eigenvector.
Inverse power method for the stochastic general eigenvalue problem The stochastic general eigenvalue problem can be written in the same form as (5.6)
As in the deterministic case, an initial guess for the eigenvector (˜vk) and the eigenvalue ( ˜φk) has to be made. For the first step the deterministic results plus a small deviation can be used. The eigenvector can then be updated by solving
˜
G˜vnew= ˜Λk˜vk. (5.21)
This equation can be solved using the spectral method, by using the polynomial chaos expansion to write ˜ G = n X i=0 Giψ˜i; (5.22a) ˜ vk = n X i=0 vkiψ˜i; (5.22b) ˜ Λk = n X i=0 Λikψ˜i, (5.22c)
where the polynomials are orthogonal in the sense that [13] Ehψ˜iψ˜j
i
= ciδij. (5.23)
Substitution of these expansions in equation (5.19) yields
n X i=0 n X j=0 Gjψ˜iψ˜j vnewi = n X i=0 n X j=0 Λkivkjψ˜iψ˜j. (5.24)
Multiplication of this equation with ˜ψrand taking the expectation then results in n X i=0 n X j=0 GjEhψ˜iψ˜jψ˜ri vnewi = n X i=0 n X j=0 ΛkivkjEhψ˜iψ˜jψ˜ri, for r = 0 . . . n. (5.25) Now define Ωri = n X j=0 GjEhψ˜iψ˜jψ˜ri; (5.26a) ξkr = n X i=0 n X j=0 ΛkivkjE h ˜ ψiψ˜jψ˜r i , (5.26b)
such that these equation can be assembled as Ω00 · · · Ω1n .. . . .. ... Ωn1 · · · Ωnn vnew 0 .. . vnnew = ξk0 .. . ξkn . (5.27)
The updated eigenvector ˜vnew that corresponds to the eigenvalue ˜φk can then be obtained by solving this linear system. As in the deterministic case, the complex angle of the updated eigenvector needs to be fixed and the updated eigenvector needs to be normalised.
The stochastic eigenvalue corresponding to the better approximation of the eigenvector can be found using the Rayleigh quotient
˜
φnew= (˜v
new)HA˜(˜vnew)
Substitution of the chaos expansions then gives
n
X
i=0
φnewi Γnewri = ηrnew, for r = 0 . . . n, (5.29)
with Γnewri = n X j=0 n X p=0 vnew j H vnew p Ehψ˜iψ˜jψ˜rψ˜p i ; (5.30a) ηnewr = n X i=0 n X j=0 n X p=0 vnewj HAi vnewp Ehψ˜iψ˜jψ˜rψ˜pi. (5.30b)
The updated eigenvalues can then be found by solving the linear system Γnew 00 · · · Γnew1n .. . . .. ... Γnew n1 · · · Γnewnn φnew 0 .. . φnew n = ηnew 0 .. . ηnew n . (5.31)
As in the deterministic case, the new initial guess for the eigenvector and eigenvalue is made using ˜
vk+1 = (1 − κ)˜vk+ κ˜vnew, with κ ∈ R(0, 1]; (5.32a) ˜
φk+1 = φ˜new, (5.32b)
where κ is again selected such that the method converges.
In Ref. [11] it is metioned that the spectral expansion converges to the correct solution if the polynomial chaos basis is appropriate. For the randomness considered in this study, it turns out that the flutter point as a function of the random variables can be approximated using linear polynomials in the spectral expansion. In the case that large variations are considered, it is possible that linear polynomials are no longer capable of modelling the behaviour of the flutter point. The polynomial chaos basis is in that case not appropriate, which leads to convergence problems. These problems can be overcome by including higher-order polynomials in the expansion.
5.2 Computation of the stochastic flutter point
The problem of solving the random general eigenvalue problem using the spectral method has now been explained. The goal however is to find the spectral representation of the flutter point. The flutter problem (2.40) can best be considered. The flutter requirement can be written as
f (˜λ) = ˜λr− ˜ λ2 ˆ ı 4ξ2 = 0. (5.33)
This equation can be solved by using an iterative method (e.g. Broyden’s quasi-Newton method). In that case the random eigenvalue ˜λ has to be evaluated each iteration step. The stochastic inverse power method can be used to do this. To be able to apply this method, the random general matrix ˜
Sneeds to be written in the form
SM˜f l, ˜z=
n
X
i=0
Siψ˜i. (5.34)
In the case that the linear piston theory is used, this is straight forward. In the case that the higher-order piston theory would be applied, linearisation of the Mach number function is required. This can for example be done by making a Taylor expansion around the deterministic Mach number.
The convergence of this method can e.g. be monitored by checking the L2-norm of the explicit
Input parameters, incl. initial guess for ˜M∞. Convergence: f˜ ≤ ǫ ? Solve the random eigenvalue problem.
Compute the resid. ˜ f using (5.33). Output: ˜ λ and ˜Mf l. True False First step? Yes No Broyden Update the Jacobian. Compute the Jacobian. Update ˜M∞ using requirement 5.33.
5.2.1 Algorithmic aspects
As can be seen in Figure 9 an initial guess for the stochastic flutter point has to be made. It turns out that when the deterministic flutter point is used, the method has problems to reach the correct stochastic distribution. To much weight is given to the higher-order K-L terms. This problem can be avoided by adding a deviation of e.g. 10 percent to the deterministic flutter point. The convergence of the method in that case is improved considerably. For a correlation length of 20 percent and a coefficient of variation of the modulus of elasticity of 15 percent, typically 15 to 20 iteration steps are required. For problem with a smaller uncertainty, convergence is reached considerably faster.
Since the computational effort of the spectral method strongly depends on the number of poly-nomials in the expansion, the second-order spectral method (SOSM) is computationally much more intensive then the first-order spectral method (FOSM). In the case that a second-order spectral ap-proximation is required, it is therefore recommended first to obtain the first-order spectral expansion of the flutter point using the FOSM and then using this expansion as a starting vector for the SOSM to rapidly find the second-order expansion. It turns out that the robustness of the SOSM is increased, while the computational effort is considerably decreased. Typically only a few iterations are required to obtain convergence.
In Ref. [13, p. 672] the occurrence of singularities in the spectral formulation is mentioned to be a problem. The theoretical occurrence of negative stiffnesses due to the randomness is the cause of an appearing singularity. In the case of the random eigenvalue problem, this kind of singularities does not occur. The reason for this is that the relation between the random field of elastic properties and the random flutterpoint is close to linear, as will be shown in Chapter 7.
5.2.2 Validation
The spectral method as explained in this section can be validated using the Monte-Carlo method. Not only the moments can be checked but also the coefficients of the spectral expansion [23]. This is done for the spectral expansion of the flutter point as well as for the corresponding spectral expansion of the eigenvalue ˜λ.
The accuracy of the spectral method will be investigated in the chapter on numerical simulations, here a concise numerical validation will be performed. The setting as used in Ref. [22] is used to perform the validation. To be able to use the second order expansion, only two Karhunen-Loeve terms are used. The correlation length is taken to be 20% of the total length and the coefficient of variation is taken to be 5% of the mean modulus of elasticity.
The validation results for the moments are assembled in Table 1. The Monte-Carlo results have a confidence level of at least (for the third moment) 99%. Obviously the first order method is not capable of approximating the third moment, the second order method solves this problem. The results for the approximations of the moments can be said to be accurate.
The coefficients of the spectral expansion of the flutterpoint as well as of the corresponding eigen-value are respectively assembled in Table 2 and Table 3. The higher-order components are difficult to be validated. The reason for this is that a large number of samples is required in order to reach a sufficient confidence level. Besides that it turns out that the higher-order coefficients are sensitive to the skewness of the random set that is used. As a consequence it turns out that it is very hard to have a converged Monte-Carlo solution. The first-order coefficients are easier to determine, as can be seen in the tables, the coefficients can be predicted reasonably accurately. It can be seen that the lower the contribution of a coefficient, the lower the accuracy. This is a good thing when the computation of the moments is considered, since the coefficients with the smaller contribution will be of less importance and so the associated error will not be a problem.
1st 2nd 3rd 4th
MC (15K) 2.530 − 4.912 · 10−3 − −1.10 · 10−6 − 7.322 · 10−5 −
FOSM 2.530 0.00% 4.908 · 10−3 −0.08% − − 7.226 · 10−5 −1.31%
SOSM 2.530 0.00% 4.915 · 10−3 +0.06% −1.21 · 10−6 −10.0% 7.248 · 10−5 −1.01%
Table 1: Validation of the first- and second-order spectral method based on the stochastic moments of the flutter point.
MC FOSM SOSM 0 2.530 2.530 0.00% 2.530 0.00% 1 7.009 · 10−2 7.006 · 10−2 −0.04% 7.009 · 10−2 0.00% 2 −2.196 · 10−4 −2.127 · 10−4 +3.14% 2.131 · 10−4 −2.96% 3 − − − −4.064 · 10−5 − 4 − − − −1.525 · 10−6 − 5 − − − −1.057 · 10−3 −
Table 2: Validation of the coefficients of the first- and second-order spectral expansion of the flutter point.
predicting the moments.
6. Reliability analysis
In this chapter the computation of the probability of the occurrence of flutter will be discussed. In most of the applications in structural mechanics, the probability of failure pf is computed. Although
flutter is not a failure mechanism, the probability of the occurence of flutter will in this chapter be considered as a probability of failure.
First the most natural approach to compute the probability of failure, the crude Monte-Carlo method, will be investigated. After that methods that do not require sampling are discussed and finally a combination of both methods, the importance sampling method, will be considered.
6.1 Crude Monte-Carlo simulation
The probability of failure pf can be computed by defining an indicator function I (˜z), that has value
0 if z is in the safe domain (Sm) and has value 1 if ˜zis in the failure domain (Fm) as shown in Figure 10. The probability of failure can then be computed by
pf = Z RmI (z) p zdz ≈ 1 N N X i=1 I (zi) , (6.1) MC FOSM SOSM 0 2.12 · 105− 0.151 · 105ˆı 2.12 · 105− 0.151 · 105ˆı 0.00% 2.12 · 105− 0.151 · 105ˆı 0.00% 1 5.57 · 103− 0.198 · 103ˆı 5.56 · 103− 0.198 · 103ˆı +0.18% 5.56 · 103− 0.198 · 103ˆı +0.18% 2 5.48 · 101− 0.199 · 101ˆı 5.56 · 101− 0.190 · 101ˆı +1.47% 5.48 · 101− 0.198 · 101ˆı +0.018% 3 − − − −2.40 + 0.138ˆı − 4 − − − −6.15 + 0.218ˆı − 5 − − − −4.77 · 101+ 0.170 · 101ˆı −
where the vectors zi are distributed corresponding to the probability density function pz. Estimates
of the required sample size can be found in [19, ch. 3]. For example, consider a confidence level C with confidence intervalµMf l− κ, µMf l+ κ
, with κ chosen such that PµMf l− κ ≤ Mf l≤ µMf l+ κ
= C, (6.2)
where C is typically close to 1. A first estimate for the required samples size can be obtained using Ns> − ln(1 − C)
pf
. (6.3)
This approximate sample size should only be used as an indication. The required sample size can best be determined by checking the convergence of the Monte-Carlo method. Although the quantitative meaning of this equation is doubtfull, it gives some qualitative information on the required sample size:
• The required sample size is inversely proportional to the probability of failure.
• An increasing confidence level will also increase the required sample size, but this influence is weaker than the influence of the probability of failure.
Since the probability of failure is in practice required to be very small, the sample size is required to be very large. In practice the crude Monte-Carlo method is therefore inefficient to determine the probability of failure.
6.2 Reliability methods
The Monte-Carlo method is a fully sampling based method, which makes it inefficient for reliability analysis. On the contrary the first- and second-order reliability methods (FORM/SORM) do not require sampling at all. The idea is to make a linear or quadratic approximation of the limit state function (Figure 10), such that the integral expression
pf =
Z
RmI (z) p
zdz (6.4)
can be evaluated analytically. 6.2.1 β-point computation
The approximation of the limit state function Mlimit(z) can best be done by using a Taylor expansion
around the point that is closest to the origin since that point has the largest contribution to the integral in equation (6.4). This point is indicated by z∗ (with kz∗k = β) and is often referred to as
the β-point or the design point.
In the case of the 2-D panel problem with realistic parameters, the limit state function is smooth and close to linear. The β-point can be found by solving the optimisation problem
minimise ||z||;
subject to Mf l(z) = Mlimit. (6.5)
Due to the smoothness of Mf l(z), a gradient based solver can be used to find this optimum. In
Ref. [12] this method is referred to as the HL-RF algorithm. The idea is that given a point zk, an improved estimate zk+1 can be found by making use of a first order Taylor expansion around that
point
Mf l(zk+1) = Mf l(zk) + ∇z=zkMf l zk+1− zk
, (6.6)
where Mf l(zk+1) should satisfy
Mf l(zk+1) = Mlimit; (6.7a)
The first requirement makes sure that the better approximation is closer to the limit state, while the second requirement assures that the normal vector of the hyperplane contains the origin, which is equivalent with the minimisation requirement. Combination of these expressions then yields
zk+1= ∇z=zkMf lzk− Mf l(zk) [∇z=zkMf l]T k∇z=zkMf lk2 . (6.8)
The computation of the gradient of Mf l has already been discussed in Chapter 4. The design point
z∗can therefore be found by using the above equation iteratively until convergence is obtained. 6.2.2 First-order reliability method (FORM)
Once the β-point has been computed, the limit state function can be approximated by rotation of the random variables {zi} as shown in Figure 10. This rotation is performed by
η= Rz, (6.9)
such that the integral equation (6.1) can be rewritten as pf = Z RmI (η) p η R−1 dη = Z RmI (η) p η 1 |R|dη. (6.10)
Since R is a rotation matrix, the determinant of R equals 1, which gives pf=
Z
RmI (η) pη
dη. (6.11)
Due to the linear approximation of the limit state function, the indicator function is only a function of η1, which has value 1 if η1< −β and has value 0 if η1≥ −β. The integral is then reduced to
pf = Z RI (η 1) pη1dη1= Z −β −∞ I (η 1) pη1dη1= Φ (−β) . (6.12)
Where Φ is the cumulative density function of the standard normal random variable. The probability of failure can then easily be computed.
6.2.3 Second-order reliability method (SORM)
The first-order result can be improved by application of the Breitung [5] correction factor. This is called the second-order reliability method because the curvature of the limit state function is also taken into account. The probability of failure can then be computed by
pfquadratic = Φ (−β) m−1 Y i=1 1 √ 1 + βκi , (6.13)
where {κi} are the main curvatures of the limit state function at the design point. The procedure
to compute these main curvatures is described in e.g. Ref. [18]. To obtain the main curvatures, the Hessian of the Mach number at the design point is required. This computation has been explained in Chapter 4.
Although fast, the FORM and SORM are not accurate when compared to the crude Monte-Carlo method. The accuracy of the crude Monte-Carlo simulation can be combined with the speed of the reliability methods to obtain an efficient solution for the reliability problem.
6.3 Monte-Carlo simulation with variance reduction
z1 z2 β Mlimit Mlinear Mquadratic Fm, I = 1 Sm, I = 0 η1 η2
Figure 10: Schematic representation of the limit state function.
6.3.1 Variance reduction
Consider the random variable ˜g with mean value µg that can be approximated by
µg≈ N X i=1 gi Ns . (6.14)
Using a central limit theorem (e.g. Lindeberg-L´evy [3, p. 103]) it can be proven that in the case that {gi} are independent and identically distributed, µgwill follow a normal distribution with expectation
E [µg] and variance Var [µg]. Since {gi} are indeed idependent and identically distributed, it follows
that
µg(˜z) = E [µg] + Var [µg]
1
2z.˜ (6.15)
Using the fact that E [µg] = µg it is obtained that
µg(˜z) = µg+ Var [µg]
1
2z.˜ (6.16)
This equation states that the outcome of the sampling can be written as a deterministic mean plus some random error, where the error can be defined as
ǫ(˜z) = Var [µg]
1
2z.˜ (6.17)
It is then obvious that the error of the sampled mean is proportional to the standard deviation (or square root of the variance) of that mean. The variance of the mean value can therefore be considered as a good measure for the accuracy of the sampling method. This variance can be rewritten as
Var [µg] ≈ Ns X i=1 Var [gi] N2 s + Ns X i=1 Ns X j=1,j6=i Cov [gi, gj] N2 s . (6.18)
Since {gi} are independent and identically distributed, this can be rewritten as
Var [µg] ≈
Var [˜g] Ns
From this equation it is concluded that the variance of the mean value is proportional to the variance of ˜g itself. This implies that the accuracy of the Monte-Carlo method can be improved by reduction of the variance of the observed values {gi} of the random variable ˜g.
6.3.2 Variance of the Crude Monte-Carlo simulation
The crude Monte-Carlo method can be written in the form of equation (6.14) by taking ˜
g = I (˜z) , (6.20)
where the vector ˜zfollows a uncorrelated normal distribution. The variance of g can then be computed by Var [˜g] = Z Rm g2pzdz − E [g]2= E g2− E [g]2= pf− p2f. (6.21)
In general the probability of failure is very small (pf ≪ 1) which results in
Var [˜g] = O (pf) . (6.22)
Since the {gi} are independent and identically distributed, the order of the variance of pf can be
determined using equation (6.19) as
Var [pf] = Var [µg] = O pf Ns . (6.23) 6.3.3 Importance sampling
A method to reduce the variance of the observed values is importance sampling. The idea is not to sample around the mean, but around the design point. To demonstrate that the variance is indeed reduced in that case, rewrite the integral expression for the probability of failure as
pf = Z RmI (z) pz hz hzdz, (6.24)
with hz(z) being a probability density function. It then follows that
pf ≈ 1 Ns Ns X i=1 I (zi) pz(zi) hz(zi) , (6.25)
where the vectors {zi} are generated according to the probability density function hz(z). This can be
rewritten in the form of equation (6.19) by defining ˜
g = I (˜z)phz(˜z)
z(˜z)
(6.26) Since the {gi} are independent and identically distributed, the variance of pf can be approximated
using equation (6.19) as
Var [pf] ≈
Var [˜g] Ns
. (6.27)
The variance of the probability of failure can then be reduced significantly by selecting hz(˜z) such
that the mean of ˜zcoincides with the design point. A simple choice is to shift the probability density function pz(˜z) which results in
hz(˜z) = 1 (2π)m2 e −1 2(˜z−z ∗)T(˜z−z∗) . (6.28)
In that case ˜g is given by
˜
g = I (˜z) e−1 2(2˜z−z
∗)Tz∗
