- Email: [email protected]
Finite element method for random response of structures due to stochastic excitation
COMPUTER METHODS IN APPLIED IviECHANICS @ NORTH-HOLLAND PUBLISHING COMPANY
AND ENGINEERING
20 (1979) 173-194
FINITE ELEMENT METHOD FOR RANDOM RESPONSE OF STRUCTURES DUE TO STOCHASTIC EXCITATION S.S. DEY* Institut fiir Statik und Dynamik.der
Luft-und Raumfahrthonstruktionen,
Universitiit Stuttgart, Germany
Received 31 January 1978 Revised 20 November
1978
This paper describes the application of the finite element method to analyse the response of multidegree-linear-elastic structures subjected to stationary random stochastic loading. Two different methods, namely (a) complex matrix inversion method and, (b) normal mode method have been used to compute the responses. The basic theory for these two approaches is presented. In the first approach the global stiffness and mass matrices are used as input data, while in the other approach the eigenvalues and eigenvectors are used as input for random response calculations. For test problems the potentiality of ASKA I and ASKA II packages with their standard elements has been made use of to form the global stiffness and mass matrices, and subsequently the eigenvalues and eigenvectors have been computed. To validate the numerical procedures, further verification was obtained for various problems by comparison with the analytic solution.
Notation A a 3(m) c D Tit,
matrix K + id - 02M diagonal matrix of areas area allocated at ith nodal point admittance matrix damping matrix damping matrix in diagonalized form Fourier transform of f vector of generalized forces
lim L T fdtlfdt) dt Wlf2(0> time average of product of f,(t) and f2(t), i.e. ~402T I -T
!z h(t)
J
K
structural damping coefficient real part of admittance matrix unit impulse response function imaginary part of admittance matrix stiffness matrix
*Present address: Civil Engineering
Department,
Indian Institute of Technology,
Kharagpur,
W.B., India.
174
M n
Pi Q
40) %I S, t W X xk a1 a2 T!(t) A hk 7
S.S. Dey, Finite element method
for random response of structures
mass matrix number of degrees of freedom pressure applied over area Uj Fourier transform of q vector of generalized displacements cross-power spectral density for generalized forces cross-power spectral density for generalized displacements time cross-spectral density matrix matrix of eigenvectors k th eigenvector t- 8, t-&+7 vector of generalized displacements in diagonalized system diagonal matrix of eigenvalues k th eigenvalue time lag circular frequency cut-off frequency natural frequency for nth mode of plate
1. Introduction
In recent years considerable interest has been focussed on the response of complex structures subjected to random excitations. Examples that immediately draw one’s attention are the following: (1) aircraft flying in gusty weather condition, (2) high-rise buildings, suspension bridges, suspended roofs due to wind and earthquake loads, (3) aerospace structures and missiles subjected to jet noise excitations, (4) automotive vehicle running over a rough road, (5) ships and off-shore structures subjected to the hostile environment of the sea, (6) aircraft landing on an imperfect runway. The common feature of the design of such varieties of structures is that some of the loadings are random in nature and cannot be assessed in a deterministic sense. A great deal of effort has been spent in the past to define accurately a deterministic input of the random phenomena so that the existing knowledge of structural vibrations can be employed to determine the structural response. Deterministic analysis based on well-defined external loads is hypothetical or academic in nature. In recent years research workers have begun to realize that the problem is essentially stochastic in nature and should be described as random processes rather than as deterministic functions. In most cases either it is possible to predict the probability distribution of external loads or one can compute the cross-correlation functions of the loads based on experimental measurement.
S.S. Dey, Finite element method for random response of structures
175
It is generally found economical and convenient to carry out the random response analysis in the frequency domain by utilizing the power spectral density function associated with random excitations and the frequency response characteristic of the excited structures. The power spectral densities of the response are calculated for a set of discrete values of the frequency. These power spectral densities are then integrated numerically over the frequency to yield the mean square amplitude for response.
2. Theory of random response 2.1. Method
1 - complex
matrix inversion
The matrix theory and finite element techniques developed by Argyris [l]-[4] for static analysis of complex structure have been extended further for dynamic analysis [S]. The equation of motion for a linear structural system with n degrees of freedom can be written in matrix notation as Mtj +
cij + Kq
= f(t),
(2.1.1)
where M, C, K are the mass, damping and stiffness matrices, respectively, and f(t) is the time-dependent random force vector acting upon certain degrees of freedom. Moreover, q, 4 and ij represent vectors of generalized displacement, velocity and acceleration, respectively. Let us assume that q(t) and f(t) are stationary random processes with zero mean value. Let q(t) and f(t) be so truncated that they become zero outside the interval (-T, T). Taking the Fourier transform of eq. (2.1.1) gives (K + ioC - w*M)Q = F,
(2.1.2)
wherew is the frequency of vibration, and Q and F are the truncated Fourier transforms of q and f, respectively. Eq. (2.1.2) can be written as AQ=F, where A is the complex matrix of order n x it and can be separated parts as follows:
(2.1.3) into real and imaginary
A = A’ + iA”,
(2.1.4)
A’=!ReA=K-w*M
(2.1.5)
A” = 3mA = WC.
(2.1.6)
where and
Eq. (2.1.3) may be inverted to give Q = A-‘F = BF.
(2.1.7)
176
S.S. Dey. Finite element method for random
The matrix B is called the admittance follows:
response of structures
matrix, and it can be split up into real and imaginary
B = B’ + iB”,
parts as
(2.1 .X)
where and zrnB = B” = -(A”+
A’A” ‘A’) ‘_
(21.10)
provided C is nonsingular. Expressing eq. (2.1.7) in index notation, Q, =
2 BijF,.
(2.1.11)
/=I
Multiplying Qi by the complex conjugate of Q,,,, dividing by 2T and taking the limit as T tends to density for the generalized displacement for the ith infinity give Sq,,, which is the cross-spectral and mth degrees of freedom, i.e. (2.1.12)
where the bar denotes the conjugate of a complex number. The cross-spectral density for the generalized random forces points j and 1 is given by
f, and fi acting
at the nodal
1 s,, = Fz 2T F,fi. Substituting
(2.1.13)
eq. (2.1.13) into eq. (2.1.12) gives (2.1.14)
The cross-spectral
sfd=
density
matrix
is generally
Hermitian
[ 131 in nature:
s,,.
The cross-spectral density aginary parts as follows:
(2.1.15) for the generalized
displacement
is separated
into
real
and im-
SA,,,= 2 2 (BijBA/ + B’;jBLl)Sh, - (B’:jBA,- BijBL,)S);il, j=l I=1
(2.1.16)
SG,, = 2 2 (BijBh, + B’;jBG,)S’;,, - (B’:jBAr - BijB&,)S;,,. j=l I=1
(2.1.17)
S.S. Dey, Finite element method for random response of structures
177
The cross-spectral density matrix for the total structure is obtained by taking all the possible combinations for i and m and can be combined as S, = Si+ iSi.
(2.1.18)
The joint moment of the generalized displacements integration of the cross-spectral density matrix
for the structure can be obtained from the
(2.1.19)
w =2
The diagonal elements Wii of the above matrix represent the mean square value (q;) for the generalized displacements, and the off-diagonal terms Wij are the time average of the products of the pairs of generalized displacements (qiqj) at different node points. The approximate integration of the matrix is performed up to a preassigned cut-off frequency by using Simpson’s rule as follows: W=~AO{S,,+~(S,+S,,+~~)+~(S,+S,,+~~)+S,,},
(2.1.20)
where n is odd and do is the frequency step length. The procedure described above is well-suited for an arbitrary nonsingular damping, but the computation time and the cost of the analysis would be high for a reasonably sized problem since it requires matrix inversion at each step of the frequency. An alternative approach, which avoids matrix inversion and uses a specific form of damping which uncouples the equation of motion, is presented in the next section. 2.2. Method 2- normal mode The normal mode approach provides an elegant and powerful method for dealing with the random response of complex structures. It does not suffer from any of the disadvantages mentioned in the previous method. The admittance matrix B(o) in eq. (2.1.9) can be determined without matrix inversion if the equation of motion is uncoupled when the displacements are expressed in terms of normal mode shapes. The mode shapes are determined from the undamped, homogenous equation of motion Mti+Kq=o.
(2.2.1)
The solution to this equation has the form (2.2.2)
q(t) = Xj eh’,
where Xj is the jth mode shape, and Oj is the jth natural frequency. Eq. (2.2.1) can be reduced to the form (K - ~~M)Xjeb’ =
O.
(2.2.3)
178
S.S. Dey, Finite element method for random response of structures
For a nontrivial solution the determinant (2.2.3) must be equal to zero:
of the square matrix on the left-hand
side of eq.
det (K - w:M) = 0.
(2.2.4)
For nondefective systems it is possible to find it eigenfrequencies Oj and the associated eigenvectors Xj. By use of the coordinate transformation the generalized displacement can be written in the form (2.25)
4=x%
where each column of X is a normal mode Xjui,and q(t) is a column matrix of coordinates called the principal coordinates. Since X merely performs a coordinate transformation, it is independent of time t. Hence, lj=X7j,
(2.2.6)
4=X$
(2.2.7)
Substituting
eqs. (2.2.5)-(2.2.7) in eq. (2.1.1) and premultiplying
by X’, one gets
+ X’CXTj + x’KXr_l = X’f(t).
X’MXij
The inertia and stiffness matrices are orthogonal
(2.2.8) with respect to normal modes. Hence,
X’MX = A,
(2.2.9)
X’KX = I,
(2.2.10)
where A = [A,] is the diagonal matrix of eigenvalues The structure damping matrix C of eq. (2.1.1) is not formulation of the normal mode approach a form C way (through the normal coordinate transformation)
of the undamped system, and oj = A;“‘. known in an explicit form. In the usual is selected which uncouples in the same as the stiffness and mass matrix. Hence.
X’CX = D.
(2.2.11)
Finally, eq. (2.2.8) takes the form
nij + Dlj The jth decoupled Aj;i,
Taking
+
DjTjj
the Fourier
+ I7J= X’f(t).
(2.2.12)
equation of motion may be written as +
T)j
=
(2.2.13)
X:f(t).
transform
of the above equation
for the unit impulse
applied
for all
179
S.S. Dey, Finite element method for random response of struchues
degrees of freedom in succession, the following result is obtained: 1 rsi(w) = -w2Aj + hDj + 1 $7
(2.2.14)
*j(O) = Bj(W)XJ,
(2.2.15)
or
where jj(o) is the Fourier transform of qj(t), and pi(w) is the admittance Using eq. (2.2.5), one obtains Qj(0.J)
=
of the jth mode.
XjBj(W)X;.
(2.2.16)
The right-hand side may be arranged in a matrix form if the same procedure succession for all degrees of freedom of the structure, i.e. Q(w) = X
[B(o)]
is applied in
(2.2.17)
X'.
The admittance function is defined as the Fourier transform of the output of the system to a unit impulse. Hence, B(w) = Q(w) = X [B(m)] X'.
(2.2.18)
The diagonal matrix [B] is an array of generalized admittances Bj(w) =
defined by
1 --ORAL+ iaDj + 1’
Once the free vibration analysis is complete, the admittance matrix inversion.
(2.2.19) matrix can be determined
without
2.3. Damping The damping in a structure, in general, will be a combination of structural (hysteretic) and viscous damping. Since damping is small in most structures, it can be represented adequately by considering it to be completely structural or completely viscous. In this report only structural damping has been considered. The structural damping matrix can be written as (see [6]) C = igK,
(2.3.1)
where g is a small constant known as the structural damping factor. For structural damping the jth admittance is given by Bj(w) =
1 1 - (W/Oj)’ +
ig ’
(2.3.2)
180
S.S. Dey, Finite element method for random response of structures
The above expression can be separated Bj(O)=
into real and imaginary parts:
iJl(O),
Hj(O)-
(2.3.3)
where (2.3.4)
qco>=h [1 -
(O.J/OJj)*]*
+
g*
(2.3.5)
’
2.4. Random pressure loads In the finite element procedure, loads applied at locations other than the nodal points have to be appropriately converted to the equivalent nodal loads. A similar concept [7] has been used in the case of distributed random loads with a uniform spatial correlation. Random nodal loads can be obtained by multiplying the random forces by the area allocated at the nodal points (fig. 2.1). That is to say, the same random force acts on the entire area with full correlation.
Fig. 2.1. Equivalent
random
nodal point load.
The response at the point i of a structure due to a pressure pj applied over an area aj at the nodal point j is the convolution of the unit impulse response function h(t) and the pressure:
q{(t) = a, cLh{(t- O,)pj(O,)dO,. I -I
(2.4.1)
Similarly, the response at the point r due to pressure applied at the point k is (2.4.2)
If the loads and the response are assumed to be ergodic and stationary random processes, the
S.S. Dey, Finite element method for random response of structures
cross-correlation
181
of the two responses is (2.4.3)
for a time lag T. Substituting
eqs. (2.4.1) and (2.4.2) in the above equation, eq. (2.4.3) becomes (2.4.4)
where Rpj, is the cross-correlation of the pressure at the points i and k, and (Y, and a2 are the appropriate time transformations. By taking the Fourier transform of eq. (2.4.4) the crossspectral density of the response can be derived: S&r(W) =
UjUkSpjk
(W)Bi(o)BtCo)7
(2.4.5)
where Spjk is the cross-power spectral density of the pressure at the points i and k, and @ is the complex conjugate of Bj. The generalized displacement at the point i when all the points of the structure are loaded is the sum of the components resulting from each load:
qi(l)
=
2
(2.4.6)
&(f)7
j=l
where II is the number of load points. Similarly, the cross-correlation when all the points are loaded is Rq,(w) = i: j=l
2 Rq’qW
of two displacements
(2.4.7)
k=l
Thus, the cross-power spectral density for displacements
at i and r is (2.4.8)
which can be cast in matrix form as
SqhJ)= R4 bl SPW bl w4,
(2.4.9)
where Sq(o) and Q(o) are the cross-spectral density matrices for generalized displacement and pressure, respectively, and [a] is a diagonal matrix of areas associated with each node. Eq. (2.4.9) can be simplified to s&o) =
B(f+&)w4
(2.4.10)
S.S. Dey, Finite element method for rundom
182
response of struc‘turrs
where
S&J) = [aI m4r4. Substituting
(2.3.11)
eq. (2.2.1X) in eq. (2.4.10) the resulting equation may be written in the form
Sq(0) =
x [B(o)] xts,(o)x
[B(o)] X’.
(2.3.12)
The joint moment of the various generalized displacements for the total structure can be obtained from the integration of the cross-spectra1 density matrix as shown in eq. (21.19).
3. Examples and discussion The theories outlined in sections 2.1 and 2.2 have been incorporated into two separate computer programs to evaluate the response of structures to stochastic excitations. A few structures have been analysed by the proposed procedures and the numerical results are compared. Example
1. Simply supported beam Elastic data EI = 1 .O
Geometric data
Length = 1.0 Cross-sectional Density = 1.0
Structural damping g = 0.1
area A = 1.O
The beam is excited by clipped white noise; for the assumed excitation the generalised spectra1 density S,, in eq. (2.1.14) becomes S,, =
z&I__a (F(LJj,T)F(&
T +
power
s)>e-““’ ds
cz
=-I7i_% 1 1C sin s0 s e _1~3ds =
01
where o, is the cut-off frequency. The beam is idealised by the finite element method using the ASKA system with BECOS elements [S] (fig. 3.1). The beam is idealised with 4, 8 and 16 elements. At each idealization the global stiffness and mass matrices are formed using ASKA I and ASKA II, which are used subsequently as input data in the random response programs. In the normal mode approach the eigenvalues and eigenvectors of the undamped system are used in the response analysis. At each idealization the power spectral densities for displacements and slopes are computed by both procedures and compared with the analytic solution. The comparison is given in tables 3.la and 3.lb. This shows that the results tend to converge to the analytic solution as the number of elements in the structure is increased. It may be mentioned here that in the normal mode approach the
S.S. Dey, Finite element method for random
Noise BECOS
Element
BECOS
Element
Noise
Propaaation
c
Propaaation
@@@@@Q@ 2
3
4
Noise BECOS
response of structures
5
Nodal Point Numbering
6
Propaflat
ion
Element
Fig. 3.1. Simply supported beam under clipped white noise.
32
Elements
72 Elements(25
(9
Simply
Supported
External
IeSrees
of Freedom)
External
Dedrees
of Freedom)
ldealisation
128 Elements{49 Fig.
3.2. Simply supported
External
Debrees
of
Freedom)
plate under clipped white noise.
183
184
S.S. Dey, Finite element method
for random
Table 3.la. Simply supported beam (w = 0) Power spectral density for displacement at the mid-point Numerical Complex matrix inverse method 4 elements with 8 degrees of freedom 8 elements with 16 degrees of freedom 16 elements with 32 degrees of freedom First peak value of power for W = 9.9
1.5150x
solution Normal mode method
1.5150 x 1V4
I .6369 x 1W4
1.6369 x 1W
1.6682 x lO-“
1.6673x lW1
density
for displacement
Numerical
4 elements with 8 degrees of freedom 8 elements with 16 degrees of freedom 16 elements with 32 degrees of freedom
of the beam
1W4
spectral
response of structures
Analytic solution [12]
1.69x I or3
at the mid-point
of the beam
solution
Complex matrix inverse method
Normal mode method
1.53211 x IO-
I.53232 x lo-z
1.65266 x 1o-Z
1.65860 x lo-?
I .69108X lWZ
1.69074x lo-”
Analytic solution [ 121
I .708S x 1(1-’
-Mean square
displacement
at the centre
of the beam up to o, = X3.0 Numerical
4 elements with 8 degrees of freedom 8 elements with 16 degrees of freedom 16 elements with 32 degrees of freedom
solution
Complex matrix inverse method
Normal mode method
0.04648787
0.04713212
(3.05027Y38
0.050347Y2
0.05 I WA)6
0.0508210 1
Analytic solution [ 121
o.oso3
power spectral densities are computed by taking the contributions of all the natural modes of the system. The root mean square displacement at the mid-point of the beam and the root mean square slope at the support points are obtained by numerical integration using Simpson’s rule up to the cut-off frequency w, = 83. These values are compared with the di~erential equation solution. The results are given in tables 3.la and 3.lb. These results can be further improved by taking the smaller value of step length in the numerical integration.
S.S. Dey, Finite element method for random response of structures Table 3.lb. Simply supported beam (w = 0) Power spectral density for slope at the support point of the beam Numerical method
4 elements with 8 degrees of freedom 8 elements with 16 degrees of freedom 16 elements with 32 degrees of freedom
Complex matrix inverse method
Normal mode method
1.5107 x 1o-3
1.5107 x 1o-3
1.6656 x 1O-3
1.6656 x 1O-3
1.7055 x 1o-3
1.7046 x 1o-3
Analytic solution [ 121
1.69 x 1O-3
First peak value of power spectral density for slope at the support point of the beam for w = 9.9 Numerical method Complex matrix inverse method
Normal mode method
0.151220
0.151241
0.163131
0.163717
0.166926
0.166896
4 elements with 8 degrees of freedom 8 elements with 16 degrees of freedom 16 elements with 32 degrees of freedom
Analytic solution [12]
0.1686
Mean square slope at the support point of the beam up to oc = 83.0 Numerical method Complex matrix inverse method
Normal mode method
0.45882334
0.4654052
0.49643168
0.49710781
0.50692389
0.50180664
4 elements with 8 degrees of freedom 8 elements with 16 degrees of freedom 16 elements with 32 degrees of freedom
Analytic solution [ 121
0.523559
Example 2. Simply supported square plate subjected to clipped white noise Geometric data a = 1.0 b= 1.0
Length Breadth Thickness Density
t= p=
1.0 1.0
Elastic data
Structural damping
D=l.O
g = 0.1
Poisson’s ratio v = 0.0
185
The plate is excited by clipped white noise whose power spectral density is the same as that af example 1. The plate is idealized by TUBA-3 elements as shown in fig. 3.2. A condensation technique is applied, and the degrees of freedom for transverse deflections are retained in the analysis. The comparison of the results with increasing number of eiements is given in table 3.2.
Table 3.2. Simply supported square plate (w = 0) Power spectral density for displacement at the centre Complex matrix method
Normal mode method
0.00001355
0.00001354
0.00001502
0.0#01502
~.~~~~9
0.~15~~
32 elements 9 external degrees of freedom 72 elements 25 external degrees of freedom 128 elements 49 external degrees of freedom
The first peak value of PSD for displacement
displacement
at the centre
Normal mode method
O.OOt35262
0.0013525
0.00148406
0.00 14839
O.OOlS393
0.0015393
at the centre
0.00001714
Analytic solution [12]
0.001731
of the plate up to or = 100
Complex matrix method 32 elements 9 external degrees of freedom 72 elements 925 external degrees of freedom 128 elements 49 external degrees of freedom
Analytic solution [ 121
of the plate for w = 20
Complex matrix method 32 elements 9 external degrees of freedom 72 eiements 25 external degrees of freedom 128 elements 49 external degrees of freedom Mean square
of the plate
Normal mode method
0.0085548 1
0.0085548
0.009%6 1
fCWG95556
0.00996389
0.0099638 1
Analytic solution [ 121
0.0107
It is evident from the table that the results tend to converge to the differential equation solution as the number of elements in the structure is increased. The root mean square values are computed by using Simpson’s rule with step length equal to unity. In fact, these results can be further improved by using a finer step length for the frequency in the numerical integration.
187
S.S. Dey, Finite element method for random response of structures
Example 3. Five-bay continuous beam on equidistant simple supports Geometric data
l=l.O p= 1.0 A= 1.0
Length of each span Density Cross-sectional area Nondimensional speed of sound
Elastic data
Structural damping
Flexural stiffness = 1.O
g = 0.2
Co= 6.
The random response of this type of beam has drawn the attention of several investigators [9], [lo]. This simple model represents the panel-rib-stringer configuration of aircraft structures. The geometric and elastic data of the continuous beam is given below. Numerical calculations were carried out for the random excitation known as plane wave propagation of clipped white noise with the direction of propagation being along the span of the beam. This type of excitation is of practical interest since it is representative of the pressure field induced downstream by jet engine exhausts. The power spectral density of a acoustic noise is given in [lo] as for 10~15o,, for If.01> o,, where pi o, Co 5
is is is is
the the the the
mean square pressure, cut-off frequency, nondimensional speed of sound, separation distance between the points under consideration.
The continuous pressure field along the span of the beam is approximated loads at the nodal points of the beam. Jet r
2
Element
3
4
Noise
Propaaation BECOS
Number
‘5
6
7
as discrete lateral
Element7
8
4
i-
Nodal
Point
Numberin
Fig. 3.3. Five-span continuous
beam under jet noise.
Each bay of the beam is idealized by four identical BECOS elements. The ASKA program was used to build up the uncondensed and condensed global stiffness and mass matrices. Eigenvalues and eigenvectors of the system were computed using the ASKA II Dynamics package. The first fifteen eigenvalues for the structure are presented in table 3.3a. It is clear from the
188
S.S. Dey, Finite element method
Table beam
3.3~~
Mode
Five-bay
Circular
continuous
frequency
1 2 3 4 5
9.87261540 10.9539594 i 3.7008657 17.2636354 20.73568 1x
6 7
39.74~7 42.069792 1 47.3931741 53.8498151 59.9188X1 1
8 9 10 11 12 13 14 1s
for random response of structures
94.0306014 97.7029126 105.935288 114.996815 122.059915
table that the natural frequencies appear to be clustered in distinct bands. The number of frequencies in each band is equal to the number of bays. This makes it different from the first two structures, whose natural frequencies are distinct and widely separated. The uncondensed and condensed stiffness and mass matrices and their associated eigenvalues and eigenvectors were used in the random response calculations. The results of both methods were compared. The power spectral densities for displacements and slopes at w = 0 are given in tables 3.3b and 3.3c, which indicate that both methods yield identical results. The power spectra for displacements and slopes for a few nodal points of the beam where the response is maximum are plotted and shown in figs. 3.4 and 3.5. Fig. 3.4 gives the deflection spectra corresponding to nodal points 18 to 21. It is seen that all the response curves reveal distinct bands with individual peaks within each band. This was to be expected since the natural frequencies were grouped in bands of five each. This indicates that the response of such structures cannot be represented with any single mode. On the contrary, many modes must be included. It may be noted that the first three modes respond very strongly, whereas the fourth and fifth ones are virtually absent, showing up only small peaks in the spectra curve. The next seven modes show distinct peaks, but the amplitudes of these peaks are small relative to the first three. The distribution of the root mean square displacements and slopes along the length of the beam are shown in fig. 3.6. It is seen that the maximum responses of both displacement and slope occur in the first and last bays. ‘Ibis is probably due to reduced boundary restraints at each end as compared to the inner supports where the beam is continuous. It is also to be noted that the results are not symmetric with respect to the eentre of the beam. This may be explained as a result of orientation of the noise field.
S.S. Dey, Finite element method for random response of structures
Table 3.3b. Five-bay continuous beam (o = 0) Power spectra density for displacement Condensed Complex matrix method 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21
system Normal mode method
Uncondensed Complex matrix method
system Normal mode method
O.oooooooO
O.oooooooO
O.oooooooO
O.oooooooO
0.00002434 0.00003845 0.00001150 O.oooooooO O.OOOOOO25 o.OOOOO248 o.OOOOOO79 o.oOOOOOO0 0.00000342 O.ooooO972 O.OOOOO342 O.oooooooO o.OOOOOO79 0.0004M248 O.OOOOOO25 O.oooooooO 0.00001150 0.00003345 0.00002434 O.oooooooO
0.00902435 0.00003847 0.00001151 O.oooooooO O.OOOOOO25 O.OOOOO248 o.OOOOOO79 o.oOOOOOO0 O.OOOOO342 0.00000972 0.00000342 o.oOOOOOO0 o.OOoOOO79 o.OOoOO248 O.OOOOOO25 O.oooooooO o.OObO1151 0.00003847 0.00002435 O.oooooooO
0.00002434 0.00003845 0.00001150 O.oooooooO O.OOOOOO25 0.00000248 o.OOOOOO79 O.oooooooO o.OOOOO342 O.OOOOO972 O.OOOOO342 O.oooooooO o.oooooo79 O.OOOOO248 O.OOOOOO25 O.oooooooO o.oOOO115o 0.00003845 0.00002434 o.oOOOOOO0
0.00002435 0.00003847 0.00001151 O.OOOOO#O O.OOOOOO25 O.ooooO248 o.OOOOOO79 O.oooooooO O.OOOOO342 O.OOOOO972 O.OOOOO342 O.oooooooO o.OOOOOO79 o.OOooO248 o.OOOOOO25 O.oooooooO 0.00001151 0.00003847 0.00002435 O.oooooooO
i
Nodal
Point-
18 19 20
-2
-------.-.-.-
-4 ;; 4 -0 -6
-8
Circular
Frequency
-
Fig. 3.4. Power spectral density for displacement
at various nodal points.
189
190
S.S. Dey, Finite element method for random response of structures
Table 3.3~. Five-bay continuous beam (o = 0) Power spectral density for slope Uncondensed Complex matrix method 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21
0.00051124 0.00019537 0.00001690 0.00025706 0.00003803 0.00003615 0.OOOOO106 0.00003023 o.OOOOO423 0.00007812 O.oooooooO 0.00007812 O.OOOOO423 0.00003023 O.OOOOOlO6 0.00003615 0.00003803 0.00025706 0.00001690 o.c!0019537 0.00051124
system Normal mode method 0.00051144 0.00019545 0.00001691 0.00025716 0.00003804 0.00003616 0.OOOOO106 0.00003025 O.ooooo423 0.00007815 O.oooooooO 0.00007815 O.OOWO423 0.00003025 0.ooooO106 0.00003616 0.00003804 0.00025716 0.00001691 0.00019545 0.00051144
Nodal
Point
18 1g_____----3-J -.-2,
-8
-
-
-~-..-.-
-
-101 0
Circular
Fig. 3.5. Power
spectral
75
50
25
density
Frequency
100
L
for slope at various
nodal points.
191
S.S. Dey, Finite element method for random response of structures
C
1
0
2 LONGITUDINAL
3
4
5
X/L
4
5
X/L
POSITIONXIL
I
2’2
s u-l
uil
-
zi
& 8 1
OS 0
I
2
Fig. 3.6. Distribution
3
of root mean square displacement
and slope.
Example 4. Two-dimensional framework Geometrical data
Elastic data
L=8m h=lm Cross-section
Density Poisson’s ratio, Young’s modulus
A = lop7 m2
The framework in [ll] by
Structural damping
p = 8 x lo7 kg/m3 v = 0.3 E = 2.1 x 1O1lN/m2
g = 0.02
is excited by a random point loading whose power spectral density is given
S(w) = co2exp (-0.08~)
Fig. 3.7. Framework under concentrated
random loading.
S.S. Dey, Finite element method for random response of structures
192
The framework is idealized by the finite element method using the ASKA system with FLA2 elements. The first fourteen natural frequencies of the structure are given in table 3.4a. The power spectral density for horizontal and vertical deflections are computed using both the procedures. The value of the power spectral densities for all the nodal points close to the first natural frequency of the structure is given in table 3.4b. These tables indicate that both methods yield identical results. The root mean square displacements for different nodes are calculated by using Simpson’s rule of integration. A frequency step length of 0.5 is used in this example. The numerical integration is continued up to frequency 60, which is sufficient to cover all the natural frequencies of the structure. The mean square values are presented in table 3.4~.
Table 3.4a. framework Mode 1 2 3 4 5 6 7 8 9 10 11 12 13 14
Two-dimensional
Circular
frequency
4.5871440 8.5284681 11.779205 19.415803 21.27268 1 22.679357 28.587067 32.888721 35.313580 37.593613 38.847446 40.210280 48.534597 49.325798
Table 3.4b. Two-dimensional
framework
(w = 4.5) Power spectral density for vertical deflection
Power spectral density for horizontal deflection Nodal point 2 3 4 s 6 7 8
Complex matrix inverse method 0.00007824 0.00000358 0.00001381 o.oOOOOOO0 0.00001381 0.00000358 O.oooO7824
Normal mode method 0.00007825 0.OOOOO358 0.00001381 O.(H_) 0.00001381 0.oOOOO358 0.00007825
Nodal point
Complex matrix inverse method
Normal mode method
2 3 4 5 6 7 8
0.00013552 0.00054206 0.00092176 0.00114813 0.00092176 0.00054206 0.00013552
0.000135ss 0.00054216 0.00092193 0.00114834 0.00092193 0.00054216 0.00013555
193
S.S. Dey, Finite element method for random response of structures
Table 3.4c. Two-dimensional
framework (w, = 60) Mean square vertical deflection
Mean square horizontal deflection Nodal point
Complex matrix inverse method
Normal mode method
Nodal point
Complex matrix inverse method
Normal mode method
2 3 4 5 6 7 8
0.00026138 0.00007777 0.00004788 O.oooooooO 0.00004788 0.00007777 0.00026138
0.00026139 0.00007777 0.00004788 O.oooooooO 0.00004788 0.00007777 0.00026139
2 3 4 5 6 7 8
0.00034987 0.00089565 0.00139746 0.00185153 0.00139746 0.00089565 0.00034987
0.00034989 0.00089578 0.00139769 0.00185181 0.00139769 0.00089578 0.00034989
5. Conclusion This report has described the general theories for the frequency domain approach of multidegree linear elastic structures subjected to random excitation forces. Various examples are solved and the results are presented. The size of the problems solved was limited to in-core capacity of the computer. In actual practice the method can be easily extended to the analysis of more complex structures. The rich element library of ASKA and its modelling capability can be fully exploited along with the ASKA II Dynamics package for random response calculations.
Acknowledgment The author wishes to express his sincere gratitude hospitality and support.
to Professor
J.H. Argyris for his
References [l] J.H. Argyris, Energy theorems and structural analysis, Aircraft Eng. 26 (1954) 347-356,383-387, 394; 27 (1955) 42-58, 80-84, 125-134, 145-158. Also (Butterworths, London, 1960, 1973). [2] J.H. Argyris, Continua and d&continua, Opening address to International Conference on Matrix Methods of Structural Mechanics, Dayton, Ohio, Wright-Patterson USAF Base, Oct. 1%5, (1%7) l-198. [3] J.H. Argyris, Die Matrizentheorie der Statik, Ing.-Arch. 25 (1957) 174-192. [4] J.H. Argyris, Recent advances in matrix methods of structural analysis, in: Progress in Aeronautical Science (Pergamon, London, 1964). [5] J.H. Argyris, Some results on the free-free oscillations of aircraft type structures, IUTAM Symposium on Recent Advances in the Mechanics of Linear Vibrations, Apr. 1%5, Rev. Francaise Met. 15 (1%5) 59. [6] W.W. Soroka, Note on the relations between viscous and structural damping coefficients. J. Aeron. Sci. 7 (1949) 409-411. [7] W.T. Thomson, Continuous structures excited by correlated random forces, Int. J. Mech. Sci. 4 (1962) 109-114. [8] ASKA I user’s reference manual, ISD Rept. 78, rev. B (Inst. Stat. Dyn., Stuttgart, 1973).
194
S.S. Dey, Finite element method for random response of structures
[9] C.A. Mercer. Response of a multi-supported beam to a random pressure field. J. Sound Vib. (1965) 29%X)6. [lo] M.D. Olson, A consistent finite element method for random response problems. Camps. Structs. 2 (1972) 163-180. [ I11 I.B. Alpay and S. Utku, On response of initially stressed structures to random excitations. Camps. Structs. 3 (1973) 297-314. [12] S.S. Dey. Finite element method for random response of structures due to stochastic excitation, ISD Rept. 221 (Inst. Stat. Dyn., Stuttgart. 1976). [13] J.S. Bendat and A. G.Piersol. Measurement and analysis of random data (Wiley, 19%).
AND ENGINEERING
20 (1979) 173-194
FINITE ELEMENT METHOD FOR RANDOM RESPONSE OF STRUCTURES DUE TO STOCHASTIC EXCITATION S.S. DEY* Institut fiir Statik und Dynamik.der
Luft-und Raumfahrthonstruktionen,
Universitiit Stuttgart, Germany
Received 31 January 1978 Revised 20 November
1978
This paper describes the application of the finite element method to analyse the response of multidegree-linear-elastic structures subjected to stationary random stochastic loading. Two different methods, namely (a) complex matrix inversion method and, (b) normal mode method have been used to compute the responses. The basic theory for these two approaches is presented. In the first approach the global stiffness and mass matrices are used as input data, while in the other approach the eigenvalues and eigenvectors are used as input for random response calculations. For test problems the potentiality of ASKA I and ASKA II packages with their standard elements has been made use of to form the global stiffness and mass matrices, and subsequently the eigenvalues and eigenvectors have been computed. To validate the numerical procedures, further verification was obtained for various problems by comparison with the analytic solution.
Notation A a 3(m) c D Tit,
matrix K + id - 02M diagonal matrix of areas area allocated at ith nodal point admittance matrix damping matrix damping matrix in diagonalized form Fourier transform of f vector of generalized forces
lim L T fdtlfdt) dt Wlf2(0> time average of product of f,(t) and f2(t), i.e. ~402T I -T
!z h(t)
J
K
structural damping coefficient real part of admittance matrix unit impulse response function imaginary part of admittance matrix stiffness matrix
*Present address: Civil Engineering
Department,
Indian Institute of Technology,
Kharagpur,
W.B., India.
174
M n
Pi Q
40) %I S, t W X xk a1 a2 T!(t) A hk 7
S.S. Dey, Finite element method
for random response of structures
mass matrix number of degrees of freedom pressure applied over area Uj Fourier transform of q vector of generalized displacements cross-power spectral density for generalized forces cross-power spectral density for generalized displacements time cross-spectral density matrix matrix of eigenvectors k th eigenvector t- 8, t-&+7 vector of generalized displacements in diagonalized system diagonal matrix of eigenvalues k th eigenvalue time lag circular frequency cut-off frequency natural frequency for nth mode of plate
1. Introduction
In recent years considerable interest has been focussed on the response of complex structures subjected to random excitations. Examples that immediately draw one’s attention are the following: (1) aircraft flying in gusty weather condition, (2) high-rise buildings, suspension bridges, suspended roofs due to wind and earthquake loads, (3) aerospace structures and missiles subjected to jet noise excitations, (4) automotive vehicle running over a rough road, (5) ships and off-shore structures subjected to the hostile environment of the sea, (6) aircraft landing on an imperfect runway. The common feature of the design of such varieties of structures is that some of the loadings are random in nature and cannot be assessed in a deterministic sense. A great deal of effort has been spent in the past to define accurately a deterministic input of the random phenomena so that the existing knowledge of structural vibrations can be employed to determine the structural response. Deterministic analysis based on well-defined external loads is hypothetical or academic in nature. In recent years research workers have begun to realize that the problem is essentially stochastic in nature and should be described as random processes rather than as deterministic functions. In most cases either it is possible to predict the probability distribution of external loads or one can compute the cross-correlation functions of the loads based on experimental measurement.
S.S. Dey, Finite element method for random response of structures
175
It is generally found economical and convenient to carry out the random response analysis in the frequency domain by utilizing the power spectral density function associated with random excitations and the frequency response characteristic of the excited structures. The power spectral densities of the response are calculated for a set of discrete values of the frequency. These power spectral densities are then integrated numerically over the frequency to yield the mean square amplitude for response.
2. Theory of random response 2.1. Method
1 - complex
matrix inversion
The matrix theory and finite element techniques developed by Argyris [l]-[4] for static analysis of complex structure have been extended further for dynamic analysis [S]. The equation of motion for a linear structural system with n degrees of freedom can be written in matrix notation as Mtj +
cij + Kq
= f(t),
(2.1.1)
where M, C, K are the mass, damping and stiffness matrices, respectively, and f(t) is the time-dependent random force vector acting upon certain degrees of freedom. Moreover, q, 4 and ij represent vectors of generalized displacement, velocity and acceleration, respectively. Let us assume that q(t) and f(t) are stationary random processes with zero mean value. Let q(t) and f(t) be so truncated that they become zero outside the interval (-T, T). Taking the Fourier transform of eq. (2.1.1) gives (K + ioC - w*M)Q = F,
(2.1.2)
wherew is the frequency of vibration, and Q and F are the truncated Fourier transforms of q and f, respectively. Eq. (2.1.2) can be written as AQ=F, where A is the complex matrix of order n x it and can be separated parts as follows:
(2.1.3) into real and imaginary
A = A’ + iA”,
(2.1.4)
A’=!ReA=K-w*M
(2.1.5)
A” = 3mA = WC.
(2.1.6)
where and
Eq. (2.1.3) may be inverted to give Q = A-‘F = BF.
(2.1.7)
176
S.S. Dey. Finite element method for random
The matrix B is called the admittance follows:
response of structures
matrix, and it can be split up into real and imaginary
B = B’ + iB”,
parts as
(2.1 .X)
where and zrnB = B” = -(A”+
A’A” ‘A’) ‘_
(21.10)
provided C is nonsingular. Expressing eq. (2.1.7) in index notation, Q, =
2 BijF,.
(2.1.11)
/=I
Multiplying Qi by the complex conjugate of Q,,,, dividing by 2T and taking the limit as T tends to density for the generalized displacement for the ith infinity give Sq,,, which is the cross-spectral and mth degrees of freedom, i.e. (2.1.12)
where the bar denotes the conjugate of a complex number. The cross-spectral density for the generalized random forces points j and 1 is given by
f, and fi acting
at the nodal
1 s,, = Fz 2T F,fi. Substituting
(2.1.13)
eq. (2.1.13) into eq. (2.1.12) gives (2.1.14)
The cross-spectral
sfd=
density
matrix
is generally
Hermitian
[ 131 in nature:
s,,.
The cross-spectral density aginary parts as follows:
(2.1.15) for the generalized
displacement
is separated
into
real
and im-
SA,,,= 2 2 (BijBA/ + B’;jBLl)Sh, - (B’:jBA,- BijBL,)S);il, j=l I=1
(2.1.16)
SG,, = 2 2 (BijBh, + B’;jBG,)S’;,, - (B’:jBAr - BijB&,)S;,,. j=l I=1
(2.1.17)
S.S. Dey, Finite element method for random response of structures
177
The cross-spectral density matrix for the total structure is obtained by taking all the possible combinations for i and m and can be combined as S, = Si+ iSi.
(2.1.18)
The joint moment of the generalized displacements integration of the cross-spectral density matrix
for the structure can be obtained from the
(2.1.19)
w =2
The diagonal elements Wii of the above matrix represent the mean square value (q;) for the generalized displacements, and the off-diagonal terms Wij are the time average of the products of the pairs of generalized displacements (qiqj) at different node points. The approximate integration of the matrix is performed up to a preassigned cut-off frequency by using Simpson’s rule as follows: W=~AO{S,,+~(S,+S,,+~~)+~(S,+S,,+~~)+S,,},
(2.1.20)
where n is odd and do is the frequency step length. The procedure described above is well-suited for an arbitrary nonsingular damping, but the computation time and the cost of the analysis would be high for a reasonably sized problem since it requires matrix inversion at each step of the frequency. An alternative approach, which avoids matrix inversion and uses a specific form of damping which uncouples the equation of motion, is presented in the next section. 2.2. Method 2- normal mode The normal mode approach provides an elegant and powerful method for dealing with the random response of complex structures. It does not suffer from any of the disadvantages mentioned in the previous method. The admittance matrix B(o) in eq. (2.1.9) can be determined without matrix inversion if the equation of motion is uncoupled when the displacements are expressed in terms of normal mode shapes. The mode shapes are determined from the undamped, homogenous equation of motion Mti+Kq=o.
(2.2.1)
The solution to this equation has the form (2.2.2)
q(t) = Xj eh’,
where Xj is the jth mode shape, and Oj is the jth natural frequency. Eq. (2.2.1) can be reduced to the form (K - ~~M)Xjeb’ =
O.
(2.2.3)
178
S.S. Dey, Finite element method for random response of structures
For a nontrivial solution the determinant (2.2.3) must be equal to zero:
of the square matrix on the left-hand
side of eq.
det (K - w:M) = 0.
(2.2.4)
For nondefective systems it is possible to find it eigenfrequencies Oj and the associated eigenvectors Xj. By use of the coordinate transformation the generalized displacement can be written in the form (2.25)
4=x%
where each column of X is a normal mode Xjui,and q(t) is a column matrix of coordinates called the principal coordinates. Since X merely performs a coordinate transformation, it is independent of time t. Hence, lj=X7j,
(2.2.6)
4=X$
(2.2.7)
Substituting
eqs. (2.2.5)-(2.2.7) in eq. (2.1.1) and premultiplying
by X’, one gets
+ X’CXTj + x’KXr_l = X’f(t).
X’MXij
The inertia and stiffness matrices are orthogonal
(2.2.8) with respect to normal modes. Hence,
X’MX = A,
(2.2.9)
X’KX = I,
(2.2.10)
where A = [A,] is the diagonal matrix of eigenvalues The structure damping matrix C of eq. (2.1.1) is not formulation of the normal mode approach a form C way (through the normal coordinate transformation)
of the undamped system, and oj = A;“‘. known in an explicit form. In the usual is selected which uncouples in the same as the stiffness and mass matrix. Hence.
X’CX = D.
(2.2.11)
Finally, eq. (2.2.8) takes the form
nij + Dlj The jth decoupled Aj;i,
Taking
+
DjTjj
the Fourier
+ I7J= X’f(t).
(2.2.12)
equation of motion may be written as +
T)j
=
(2.2.13)
X:f(t).
transform
of the above equation
for the unit impulse
applied
for all
179
S.S. Dey, Finite element method for random response of struchues
degrees of freedom in succession, the following result is obtained: 1 rsi(w) = -w2Aj + hDj + 1 $7
(2.2.14)
*j(O) = Bj(W)XJ,
(2.2.15)
or
where jj(o) is the Fourier transform of qj(t), and pi(w) is the admittance Using eq. (2.2.5), one obtains Qj(0.J)
=
of the jth mode.
XjBj(W)X;.
(2.2.16)
The right-hand side may be arranged in a matrix form if the same procedure succession for all degrees of freedom of the structure, i.e. Q(w) = X
[B(o)]
is applied in
(2.2.17)
X'.
The admittance function is defined as the Fourier transform of the output of the system to a unit impulse. Hence, B(w) = Q(w) = X [B(m)] X'.
(2.2.18)
The diagonal matrix [B] is an array of generalized admittances Bj(w) =
defined by
1 --ORAL+ iaDj + 1’
Once the free vibration analysis is complete, the admittance matrix inversion.
(2.2.19) matrix can be determined
without
2.3. Damping The damping in a structure, in general, will be a combination of structural (hysteretic) and viscous damping. Since damping is small in most structures, it can be represented adequately by considering it to be completely structural or completely viscous. In this report only structural damping has been considered. The structural damping matrix can be written as (see [6]) C = igK,
(2.3.1)
where g is a small constant known as the structural damping factor. For structural damping the jth admittance is given by Bj(w) =
1 1 - (W/Oj)’ +
ig ’
(2.3.2)
180
S.S. Dey, Finite element method for random response of structures
The above expression can be separated Bj(O)=
into real and imaginary parts:
iJl(O),
Hj(O)-
(2.3.3)
where (2.3.4)
qco>=h [1 -
(O.J/OJj)*]*
+
g*
(2.3.5)
’
2.4. Random pressure loads In the finite element procedure, loads applied at locations other than the nodal points have to be appropriately converted to the equivalent nodal loads. A similar concept [7] has been used in the case of distributed random loads with a uniform spatial correlation. Random nodal loads can be obtained by multiplying the random forces by the area allocated at the nodal points (fig. 2.1). That is to say, the same random force acts on the entire area with full correlation.
Fig. 2.1. Equivalent
random
nodal point load.
The response at the point i of a structure due to a pressure pj applied over an area aj at the nodal point j is the convolution of the unit impulse response function h(t) and the pressure:
q{(t) = a, cLh{(t- O,)pj(O,)dO,. I -I
(2.4.1)
Similarly, the response at the point r due to pressure applied at the point k is (2.4.2)
If the loads and the response are assumed to be ergodic and stationary random processes, the
S.S. Dey, Finite element method for random response of structures
cross-correlation
181
of the two responses is (2.4.3)
for a time lag T. Substituting
eqs. (2.4.1) and (2.4.2) in the above equation, eq. (2.4.3) becomes (2.4.4)
where Rpj, is the cross-correlation of the pressure at the points i and k, and (Y, and a2 are the appropriate time transformations. By taking the Fourier transform of eq. (2.4.4) the crossspectral density of the response can be derived: S&r(W) =
UjUkSpjk
(W)Bi(o)BtCo)7
(2.4.5)
where Spjk is the cross-power spectral density of the pressure at the points i and k, and @ is the complex conjugate of Bj. The generalized displacement at the point i when all the points of the structure are loaded is the sum of the components resulting from each load:
qi(l)
=
2
(2.4.6)
&(f)7
j=l
where II is the number of load points. Similarly, the cross-correlation when all the points are loaded is Rq,(w) = i: j=l
2 Rq’qW
of two displacements
(2.4.7)
k=l
Thus, the cross-power spectral density for displacements
at i and r is (2.4.8)
which can be cast in matrix form as
SqhJ)= R4 bl SPW bl w4,
(2.4.9)
where Sq(o) and Q(o) are the cross-spectral density matrices for generalized displacement and pressure, respectively, and [a] is a diagonal matrix of areas associated with each node. Eq. (2.4.9) can be simplified to s&o) =
B(f+&)w4
(2.4.10)
S.S. Dey, Finite element method for rundom
182
response of struc‘turrs
where
S&J) = [aI m4r4. Substituting
(2.3.11)
eq. (2.2.1X) in eq. (2.4.10) the resulting equation may be written in the form
Sq(0) =
x [B(o)] xts,(o)x
[B(o)] X’.
(2.3.12)
The joint moment of the various generalized displacements for the total structure can be obtained from the integration of the cross-spectra1 density matrix as shown in eq. (21.19).
3. Examples and discussion The theories outlined in sections 2.1 and 2.2 have been incorporated into two separate computer programs to evaluate the response of structures to stochastic excitations. A few structures have been analysed by the proposed procedures and the numerical results are compared. Example
1. Simply supported beam Elastic data EI = 1 .O
Geometric data
Length = 1.0 Cross-sectional Density = 1.0
Structural damping g = 0.1
area A = 1.O
The beam is excited by clipped white noise; for the assumed excitation the generalised spectra1 density S,, in eq. (2.1.14) becomes S,, =
z&I__a (F(LJj,T)F(&
T +
power
s)>e-““’ ds
cz
=-I7i_% 1 1C sin s0 s e _1~3ds =
01
where o, is the cut-off frequency. The beam is idealised by the finite element method using the ASKA system with BECOS elements [S] (fig. 3.1). The beam is idealised with 4, 8 and 16 elements. At each idealization the global stiffness and mass matrices are formed using ASKA I and ASKA II, which are used subsequently as input data in the random response programs. In the normal mode approach the eigenvalues and eigenvectors of the undamped system are used in the response analysis. At each idealization the power spectral densities for displacements and slopes are computed by both procedures and compared with the analytic solution. The comparison is given in tables 3.la and 3.lb. This shows that the results tend to converge to the analytic solution as the number of elements in the structure is increased. It may be mentioned here that in the normal mode approach the
S.S. Dey, Finite element method for random
Noise BECOS
Element
BECOS
Element
Noise
Propaaation
c
Propaaation
@@@@@Q@ 2
3
4
Noise BECOS
response of structures
5
Nodal Point Numbering
6
Propaflat
ion
Element
Fig. 3.1. Simply supported beam under clipped white noise.
32
Elements
72 Elements(25
(9
Simply
Supported
External
IeSrees
of Freedom)
External
Dedrees
of Freedom)
ldealisation
128 Elements{49 Fig.
3.2. Simply supported
External
Debrees
of
Freedom)
plate under clipped white noise.
183
184
S.S. Dey, Finite element method
for random
Table 3.la. Simply supported beam (w = 0) Power spectral density for displacement at the mid-point Numerical Complex matrix inverse method 4 elements with 8 degrees of freedom 8 elements with 16 degrees of freedom 16 elements with 32 degrees of freedom First peak value of power for W = 9.9
1.5150x
solution Normal mode method
1.5150 x 1V4
I .6369 x 1W4
1.6369 x 1W
1.6682 x lO-“
1.6673x lW1
density
for displacement
Numerical
4 elements with 8 degrees of freedom 8 elements with 16 degrees of freedom 16 elements with 32 degrees of freedom
of the beam
1W4
spectral
response of structures
Analytic solution [12]
1.69x I or3
at the mid-point
of the beam
solution
Complex matrix inverse method
Normal mode method
1.53211 x IO-
I.53232 x lo-z
1.65266 x 1o-Z
1.65860 x lo-?
I .69108X lWZ
1.69074x lo-”
Analytic solution [ 121
I .708S x 1(1-’
-Mean square
displacement
at the centre
of the beam up to o, = X3.0 Numerical
4 elements with 8 degrees of freedom 8 elements with 16 degrees of freedom 16 elements with 32 degrees of freedom
solution
Complex matrix inverse method
Normal mode method
0.04648787
0.04713212
(3.05027Y38
0.050347Y2
0.05 I WA)6
0.0508210 1
Analytic solution [ 121
o.oso3
power spectral densities are computed by taking the contributions of all the natural modes of the system. The root mean square displacement at the mid-point of the beam and the root mean square slope at the support points are obtained by numerical integration using Simpson’s rule up to the cut-off frequency w, = 83. These values are compared with the di~erential equation solution. The results are given in tables 3.la and 3.lb. These results can be further improved by taking the smaller value of step length in the numerical integration.
S.S. Dey, Finite element method for random response of structures Table 3.lb. Simply supported beam (w = 0) Power spectral density for slope at the support point of the beam Numerical method
4 elements with 8 degrees of freedom 8 elements with 16 degrees of freedom 16 elements with 32 degrees of freedom
Complex matrix inverse method
Normal mode method
1.5107 x 1o-3
1.5107 x 1o-3
1.6656 x 1O-3
1.6656 x 1O-3
1.7055 x 1o-3
1.7046 x 1o-3
Analytic solution [ 121
1.69 x 1O-3
First peak value of power spectral density for slope at the support point of the beam for w = 9.9 Numerical method Complex matrix inverse method
Normal mode method
0.151220
0.151241
0.163131
0.163717
0.166926
0.166896
4 elements with 8 degrees of freedom 8 elements with 16 degrees of freedom 16 elements with 32 degrees of freedom
Analytic solution [12]
0.1686
Mean square slope at the support point of the beam up to oc = 83.0 Numerical method Complex matrix inverse method
Normal mode method
0.45882334
0.4654052
0.49643168
0.49710781
0.50692389
0.50180664
4 elements with 8 degrees of freedom 8 elements with 16 degrees of freedom 16 elements with 32 degrees of freedom
Analytic solution [ 121
0.523559
Example 2. Simply supported square plate subjected to clipped white noise Geometric data a = 1.0 b= 1.0
Length Breadth Thickness Density
t= p=
1.0 1.0
Elastic data
Structural damping
D=l.O
g = 0.1
Poisson’s ratio v = 0.0
185
The plate is excited by clipped white noise whose power spectral density is the same as that af example 1. The plate is idealized by TUBA-3 elements as shown in fig. 3.2. A condensation technique is applied, and the degrees of freedom for transverse deflections are retained in the analysis. The comparison of the results with increasing number of eiements is given in table 3.2.
Table 3.2. Simply supported square plate (w = 0) Power spectral density for displacement at the centre Complex matrix method
Normal mode method
0.00001355
0.00001354
0.00001502
0.0#01502
~.~~~~9
0.~15~~
32 elements 9 external degrees of freedom 72 elements 25 external degrees of freedom 128 elements 49 external degrees of freedom
The first peak value of PSD for displacement
displacement
at the centre
Normal mode method
O.OOt35262
0.0013525
0.00148406
0.00 14839
O.OOlS393
0.0015393
at the centre
0.00001714
Analytic solution [12]
0.001731
of the plate up to or = 100
Complex matrix method 32 elements 9 external degrees of freedom 72 elements 925 external degrees of freedom 128 elements 49 external degrees of freedom
Analytic solution [ 121
of the plate for w = 20
Complex matrix method 32 elements 9 external degrees of freedom 72 eiements 25 external degrees of freedom 128 elements 49 external degrees of freedom Mean square
of the plate
Normal mode method
0.0085548 1
0.0085548
0.009%6 1
fCWG95556
0.00996389
0.0099638 1
Analytic solution [ 121
0.0107
It is evident from the table that the results tend to converge to the differential equation solution as the number of elements in the structure is increased. The root mean square values are computed by using Simpson’s rule with step length equal to unity. In fact, these results can be further improved by using a finer step length for the frequency in the numerical integration.
187
S.S. Dey, Finite element method for random response of structures
Example 3. Five-bay continuous beam on equidistant simple supports Geometric data
l=l.O p= 1.0 A= 1.0
Length of each span Density Cross-sectional area Nondimensional speed of sound
Elastic data
Structural damping
Flexural stiffness = 1.O
g = 0.2
Co= 6.
The random response of this type of beam has drawn the attention of several investigators [9], [lo]. This simple model represents the panel-rib-stringer configuration of aircraft structures. The geometric and elastic data of the continuous beam is given below. Numerical calculations were carried out for the random excitation known as plane wave propagation of clipped white noise with the direction of propagation being along the span of the beam. This type of excitation is of practical interest since it is representative of the pressure field induced downstream by jet engine exhausts. The power spectral density of a acoustic noise is given in [lo] as for 10~15o,, for If.01> o,, where pi o, Co 5
is is is is
the the the the
mean square pressure, cut-off frequency, nondimensional speed of sound, separation distance between the points under consideration.
The continuous pressure field along the span of the beam is approximated loads at the nodal points of the beam. Jet r
2
Element
3
4
Noise
Propaaation BECOS
Number
‘5
6
7
as discrete lateral
Element7
8
4
i-
Nodal
Point
Numberin
Fig. 3.3. Five-span continuous
beam under jet noise.
Each bay of the beam is idealized by four identical BECOS elements. The ASKA program was used to build up the uncondensed and condensed global stiffness and mass matrices. Eigenvalues and eigenvectors of the system were computed using the ASKA II Dynamics package. The first fifteen eigenvalues for the structure are presented in table 3.3a. It is clear from the
188
S.S. Dey, Finite element method
Table beam
3.3~~
Mode
Five-bay
Circular
continuous
frequency
1 2 3 4 5
9.87261540 10.9539594 i 3.7008657 17.2636354 20.73568 1x
6 7
39.74~7 42.069792 1 47.3931741 53.8498151 59.9188X1 1
8 9 10 11 12 13 14 1s
for random response of structures
94.0306014 97.7029126 105.935288 114.996815 122.059915
table that the natural frequencies appear to be clustered in distinct bands. The number of frequencies in each band is equal to the number of bays. This makes it different from the first two structures, whose natural frequencies are distinct and widely separated. The uncondensed and condensed stiffness and mass matrices and their associated eigenvalues and eigenvectors were used in the random response calculations. The results of both methods were compared. The power spectral densities for displacements and slopes at w = 0 are given in tables 3.3b and 3.3c, which indicate that both methods yield identical results. The power spectra for displacements and slopes for a few nodal points of the beam where the response is maximum are plotted and shown in figs. 3.4 and 3.5. Fig. 3.4 gives the deflection spectra corresponding to nodal points 18 to 21. It is seen that all the response curves reveal distinct bands with individual peaks within each band. This was to be expected since the natural frequencies were grouped in bands of five each. This indicates that the response of such structures cannot be represented with any single mode. On the contrary, many modes must be included. It may be noted that the first three modes respond very strongly, whereas the fourth and fifth ones are virtually absent, showing up only small peaks in the spectra curve. The next seven modes show distinct peaks, but the amplitudes of these peaks are small relative to the first three. The distribution of the root mean square displacements and slopes along the length of the beam are shown in fig. 3.6. It is seen that the maximum responses of both displacement and slope occur in the first and last bays. ‘Ibis is probably due to reduced boundary restraints at each end as compared to the inner supports where the beam is continuous. It is also to be noted that the results are not symmetric with respect to the eentre of the beam. This may be explained as a result of orientation of the noise field.
S.S. Dey, Finite element method for random response of structures
Table 3.3b. Five-bay continuous beam (o = 0) Power spectra density for displacement Condensed Complex matrix method 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21
system Normal mode method
Uncondensed Complex matrix method
system Normal mode method
O.oooooooO
O.oooooooO
O.oooooooO
O.oooooooO
0.00002434 0.00003845 0.00001150 O.oooooooO O.OOOOOO25 o.OOOOO248 o.OOOOOO79 o.oOOOOOO0 0.00000342 O.ooooO972 O.OOOOO342 O.oooooooO o.OOOOOO79 0.0004M248 O.OOOOOO25 O.oooooooO 0.00001150 0.00003345 0.00002434 O.oooooooO
0.00902435 0.00003847 0.00001151 O.oooooooO O.OOOOOO25 O.OOOOO248 o.OOOOOO79 o.oOOOOOO0 O.OOOOO342 0.00000972 0.00000342 o.oOOOOOO0 o.OOoOOO79 o.OOoOO248 O.OOOOOO25 O.oooooooO o.OObO1151 0.00003847 0.00002435 O.oooooooO
0.00002434 0.00003845 0.00001150 O.oooooooO O.OOOOOO25 0.00000248 o.OOOOOO79 O.oooooooO o.OOOOO342 O.OOOOO972 O.OOOOO342 O.oooooooO o.oooooo79 O.OOOOO248 O.OOOOOO25 O.oooooooO o.oOOO115o 0.00003845 0.00002434 o.oOOOOOO0
0.00002435 0.00003847 0.00001151 O.OOOOO#O O.OOOOOO25 O.ooooO248 o.OOOOOO79 O.oooooooO O.OOOOO342 O.OOOOO972 O.OOOOO342 O.oooooooO o.OOOOOO79 o.OOooO248 o.OOOOOO25 O.oooooooO 0.00001151 0.00003847 0.00002435 O.oooooooO
i
Nodal
Point-
18 19 20
-2
-------.-.-.-
-4 ;; 4 -0 -6
-8
Circular
Frequency
-
Fig. 3.4. Power spectral density for displacement
at various nodal points.
189
190
S.S. Dey, Finite element method for random response of structures
Table 3.3~. Five-bay continuous beam (o = 0) Power spectral density for slope Uncondensed Complex matrix method 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21
0.00051124 0.00019537 0.00001690 0.00025706 0.00003803 0.00003615 0.OOOOO106 0.00003023 o.OOOOO423 0.00007812 O.oooooooO 0.00007812 O.OOOOO423 0.00003023 O.OOOOOlO6 0.00003615 0.00003803 0.00025706 0.00001690 o.c!0019537 0.00051124
system Normal mode method 0.00051144 0.00019545 0.00001691 0.00025716 0.00003804 0.00003616 0.OOOOO106 0.00003025 O.ooooo423 0.00007815 O.oooooooO 0.00007815 O.OOWO423 0.00003025 0.ooooO106 0.00003616 0.00003804 0.00025716 0.00001691 0.00019545 0.00051144
Nodal
Point
18 1g_____----3-J -.-2,
-8
-
-
-~-..-.-
-
-101 0
Circular
Fig. 3.5. Power
spectral
75
50
25
density
Frequency
100
L
for slope at various
nodal points.
191
S.S. Dey, Finite element method for random response of structures
C
1
0
2 LONGITUDINAL
3
4
5
X/L
4
5
X/L
POSITIONXIL
I
2’2
s u-l
uil
-
zi
& 8 1
OS 0
I
2
Fig. 3.6. Distribution
3
of root mean square displacement
and slope.
Example 4. Two-dimensional framework Geometrical data
Elastic data
L=8m h=lm Cross-section
Density Poisson’s ratio, Young’s modulus
A = lop7 m2
The framework in [ll] by
Structural damping
p = 8 x lo7 kg/m3 v = 0.3 E = 2.1 x 1O1lN/m2
g = 0.02
is excited by a random point loading whose power spectral density is given
S(w) = co2exp (-0.08~)
Fig. 3.7. Framework under concentrated
random loading.
S.S. Dey, Finite element method for random response of structures
192
The framework is idealized by the finite element method using the ASKA system with FLA2 elements. The first fourteen natural frequencies of the structure are given in table 3.4a. The power spectral density for horizontal and vertical deflections are computed using both the procedures. The value of the power spectral densities for all the nodal points close to the first natural frequency of the structure is given in table 3.4b. These tables indicate that both methods yield identical results. The root mean square displacements for different nodes are calculated by using Simpson’s rule of integration. A frequency step length of 0.5 is used in this example. The numerical integration is continued up to frequency 60, which is sufficient to cover all the natural frequencies of the structure. The mean square values are presented in table 3.4~.
Table 3.4a. framework Mode 1 2 3 4 5 6 7 8 9 10 11 12 13 14
Two-dimensional
Circular
frequency
4.5871440 8.5284681 11.779205 19.415803 21.27268 1 22.679357 28.587067 32.888721 35.313580 37.593613 38.847446 40.210280 48.534597 49.325798
Table 3.4b. Two-dimensional
framework
(w = 4.5) Power spectral density for vertical deflection
Power spectral density for horizontal deflection Nodal point 2 3 4 s 6 7 8
Complex matrix inverse method 0.00007824 0.00000358 0.00001381 o.oOOOOOO0 0.00001381 0.00000358 O.oooO7824
Normal mode method 0.00007825 0.OOOOO358 0.00001381 O.(H_) 0.00001381 0.oOOOO358 0.00007825
Nodal point
Complex matrix inverse method
Normal mode method
2 3 4 5 6 7 8
0.00013552 0.00054206 0.00092176 0.00114813 0.00092176 0.00054206 0.00013552
0.000135ss 0.00054216 0.00092193 0.00114834 0.00092193 0.00054216 0.00013555
193
S.S. Dey, Finite element method for random response of structures
Table 3.4c. Two-dimensional
framework (w, = 60) Mean square vertical deflection
Mean square horizontal deflection Nodal point
Complex matrix inverse method
Normal mode method
Nodal point
Complex matrix inverse method
Normal mode method
2 3 4 5 6 7 8
0.00026138 0.00007777 0.00004788 O.oooooooO 0.00004788 0.00007777 0.00026138
0.00026139 0.00007777 0.00004788 O.oooooooO 0.00004788 0.00007777 0.00026139
2 3 4 5 6 7 8
0.00034987 0.00089565 0.00139746 0.00185153 0.00139746 0.00089565 0.00034987
0.00034989 0.00089578 0.00139769 0.00185181 0.00139769 0.00089578 0.00034989
5. Conclusion This report has described the general theories for the frequency domain approach of multidegree linear elastic structures subjected to random excitation forces. Various examples are solved and the results are presented. The size of the problems solved was limited to in-core capacity of the computer. In actual practice the method can be easily extended to the analysis of more complex structures. The rich element library of ASKA and its modelling capability can be fully exploited along with the ASKA II Dynamics package for random response calculations.
Acknowledgment The author wishes to express his sincere gratitude hospitality and support.
to Professor
J.H. Argyris for his
References [l] J.H. Argyris, Energy theorems and structural analysis, Aircraft Eng. 26 (1954) 347-356,383-387, 394; 27 (1955) 42-58, 80-84, 125-134, 145-158. Also (Butterworths, London, 1960, 1973). [2] J.H. Argyris, Continua and d&continua, Opening address to International Conference on Matrix Methods of Structural Mechanics, Dayton, Ohio, Wright-Patterson USAF Base, Oct. 1%5, (1%7) l-198. [3] J.H. Argyris, Die Matrizentheorie der Statik, Ing.-Arch. 25 (1957) 174-192. [4] J.H. Argyris, Recent advances in matrix methods of structural analysis, in: Progress in Aeronautical Science (Pergamon, London, 1964). [5] J.H. Argyris, Some results on the free-free oscillations of aircraft type structures, IUTAM Symposium on Recent Advances in the Mechanics of Linear Vibrations, Apr. 1%5, Rev. Francaise Met. 15 (1%5) 59. [6] W.W. Soroka, Note on the relations between viscous and structural damping coefficients. J. Aeron. Sci. 7 (1949) 409-411. [7] W.T. Thomson, Continuous structures excited by correlated random forces, Int. J. Mech. Sci. 4 (1962) 109-114. [8] ASKA I user’s reference manual, ISD Rept. 78, rev. B (Inst. Stat. Dyn., Stuttgart, 1973).
194
S.S. Dey, Finite element method for random response of structures
[9] C.A. Mercer. Response of a multi-supported beam to a random pressure field. J. Sound Vib. (1965) 29%X)6. [lo] M.D. Olson, A consistent finite element method for random response problems. Camps. Structs. 2 (1972) 163-180. [ I11 I.B. Alpay and S. Utku, On response of initially stressed structures to random excitations. Camps. Structs. 3 (1973) 297-314. [12] S.S. Dey. Finite element method for random response of structures due to stochastic excitation, ISD Rept. 221 (Inst. Stat. Dyn., Stuttgart. 1976). [13] J.S. Bendat and A. G.Piersol. Measurement and analysis of random data (Wiley, 19%).