Random response analysis

Statistical measures

We define some measures of a variable that characterize it in a statistical sense.

The mean value of a random variable $x (t)$ is

E (x) = lim_{T \to \infty} \frac{1}{T} \int_{- \frac{T}{2}}^{\frac{T}{2}} x (t) d t .

Since the dynamic response is computed about a static equilibrium configuration, the mean value of any dynamic input or response variable will always be zero:

E (x) = 0.

The variance of a random variable measures the average square difference between the point value of the variable and its mean:

σ_{r}^{2} (x) = lim_{T \to \infty} \frac{1}{T} \int_{- \frac{T}{2}}^{\frac{T}{2}} {(x (t) - E (x))}^{2} d t .

Since $E (x) = 0$ for our applications, the variance is the same as the mean square value:

E (x^{2}) = lim_{T \to \infty} \frac{1}{T} \int_{- \frac{T}{2}}^{\frac{T}{2}} x^{2} (t) d t .

The units of variance are (amplitude)², so that—for example—the variance of a force has units of (force)². Generally we prefer to use the same units as the variable itself. Therefore, output variables in Abaqus/Standard are given as root mean square (“RMS”) values, $σ_{r} (x) = \sqrt{σ_{r}^{2}}$ .

Correlation

Correlation measures the similarity between two variables. Thus, the cross-correlation between two random functions of time, $x_{1} (t)$ and $x_{2} (t)$ , is the integration of the product of the two variables, with one of them shifted in time by some fixed value $τ$ to allow for the possibility that they are similar but shifted in time.

(Such a case would arise, for example, in studying a car driving along a rough road. If the separation of the axles is d and the car is moving at a steady speed v, the back axle sees the same road profile as the front axle, but delayed by a time $τ = d / v$ . Assume that the road profile moves each wheel which, in turn applies a force to the car frame through the suspension. If $F_{2}$ is the force applied to the rear axle (as a concentrated load in Abaqus) and $F_{1}$ is the force applied to the front axle, $F_{2} (t + d / v) = F_{1} (t)$ .)

Figure 1. Two random records to be correlated.

The cross-correlation function is, thus, defined as

R_{x_{1} x_{2}} (τ) = lim_{T \to \infty} \frac{1}{T} \int_{- \frac{T}{2}}^{\frac{T}{2}} x_{1} (t) x_{2} (t + τ) d t .

Since the mean value of any variable is zero, on average each variable has equal positive and negative content. If the variables are quite similar, their cross-correlation (for some values of $τ$ ) will be large; if they are not similar, the product $x_{1} x_{2}$ will sometimes be negative and sometimes positive so that the integral over all time will provide a much smaller value, regardless of the choice of $τ$ . A simple result is

\begin{matrix} R_{x_{1} x_{2}} (τ) & ​ = lim_{T \to \infty} \frac{1}{T} \int_{- \frac{T}{2}}^{\frac{T}{2}} x_{1} (t) x_{2} (t + τ) d t \\ ​ = lim_{T \to \infty} \frac{1}{T} \int_{- \frac{T}{2}}^{\frac{T}{2}} x_{1} (t - τ) x_{2} (t) d t \\ ​ = lim_{T \to \infty} \frac{1}{T} \int_{- \frac{T}{2}}^{\frac{T}{2}} x_{2} (t) x_{1} (t - τ) d t \\ ​ = R_{x_{2} x_{1}} (- τ) . \end{matrix}

For convenience the cross-correlation can be normalized to define the nondimensional normalized cross-correlation:

r_{x_{1} x_{2}} (τ) = \frac{R_{x_{1} x_{2}} (τ)}{\sqrt{σ_{r}^{2} (x_{1}) σ_{r}^{2} (x_{2})}} .

Thus, if $x_{1} = x_{2}$ , $r (0) = 1$ ; if $x_{1} = - x_{2}$ , $r (0) = - 1$ . If $x_{1}$ and $x_{2}$ are entirely dissimilar, $r \to 0$ . (When $r = 0$ the variables are said to be orthogonal.) Clearly $- 1 \leq r \leq 1$ , with small values of r indicating that $x_{1}$ and $x_{2}$ have quite different time histories.

Now consider the cross-correlation of a variable with itself: the autocorrelation. Intuitively we can see that if the variable is “very random,” its autocorrelation will be very small whenever $τ \neq 0$ : there will be no time shift that allows the variable to correlate with itself. However, if the variable is not so random—if it is just a vibration at a fixed frequency—the autocorrelation will be close to $\pm 1$ whenever $τ$ is chosen to be some integer multiple of half the period of the vibration. Thus, the autocorrelation provides a measure of how random a variable really is.

The autocorrelation of a variable $x (t)$ is, therefore,

R (τ) = lim_{T \to \infty} \frac{1}{T} \int_{- \frac{T}{2}}^{\frac{T}{2}} x (t) x (t + τ) d t .

Clearly, as $τ \to 0$ , $R (τ) \to σ_{r}^{2}$ : the autocorrelation equals the variance (the mean square value). We can, therefore, also use the normalized autocorrelation:

r (τ) = \frac{R (τ)}{σ_{r}^{2}} .

Obviously $R (τ)$ is symmetric about $τ = 0$ :

R (- τ) = R (τ);

and the value of $R (τ)$ never exceeds its value at $τ = 0$ :

R (τ) \leq R (0) = σ_{r}^{2} .

The autocorrelation function of records with very similar amplitude over the wide range of frequencies drops off rapidly as $τ$ increases. This kind of function is known as a “wide band” random function.

Figure 2. Wide band noise record and its autocorrelation.

The most extreme wide band random function would have an autocorrelation that is just a delta function:

R (τ) = δ (τ) = {\begin{matrix} σ_{r}^{2} & if τ = 0 \\ 0 & if τ \neq 0 . \end{matrix}

Such a function is called white noise. White noise has the same amplitude at all frequencies, and its autocorrelation is zero except at $τ = 0$ .

Figure 3. Autocorrelation of white noise.

Let us now consider the opposite case, known as a “narrow band” function. The extreme case of such a function is a simple sinusoidal vibration at a single frequency: $x = A \sin ω_{0} t$ . Then $R (τ)$ must also be periodic, since $R (τ)$ must attain the same value, $R (0) = σ_{r}^{2}$ , each time the shift, $τ$ , corresponds to the period of vibration. Performing the integration through time,

\begin{matrix} R (τ) & ​ = lim_{T \to \infty} \frac{1}{T} \int_{- \frac{T}{2}}^{\frac{T}{2}} A^{2} \sin ω_{0} t \sin ω_{0} (t + τ) d t \end{matrix}

\begin{matrix} ​ = A^{2} lim_{T \to \infty} \frac{1}{T} (\cos ω_{0} τ \int_{- \frac{T}{2}}^{\frac{T}{2}} \sin^{2} ω_{0} t d t + \sin ω_{0} τ \int_{- \frac{T}{2}}^{\frac{T}{2}} \sin ω_{0} t \cos ω_{0} t d t) \end{matrix}

\begin{matrix} ​ = A^{2} lim_{T \to \infty} \frac{1}{T} (\cos ω_{0} τ \int_{- \frac{T}{2}}^{\frac{T}{2}} \frac{1}{2} (1 - \cos 2 ω_{0} t) d t + \sin ω_{0} τ \int_{- \frac{T}{2}}^{\frac{T}{2}} \frac{1}{2} \sin 2 ω_{0} t d t) \\ = \frac{1}{2} A^{2} \cos ω_{0} τ = σ_{r}^{2} \cos ω_{0} τ . \end{matrix}

Figure 4. Sine wave and its autocorrelation.

The autocorrelation, $R (τ)$ , thus tells us about the nature of the random variable. If $R (τ)$ drops off rapidly as the time shift $τ$ moves away from $τ = 0$ , the variable has a broad frequency content; if it drops off more slowly and exhibits a cosine profile, the variable has a narrow frequency content centered around the frequency corresponding to the periodicity of $R (τ)$ .

Figure 5. Narrow band record and its autocorrelation.

We can extend this concept to detect the frequency content of a random variable by cross-correlating the variable with a sine wave: sweeping the wave over a range of frequencies and examining the cross-correlation tells us whether the random variable is dominated by oscillation at particular frequencies. We begin to see that the nature of stationary, ergodic random processes is best understood by examining them in the frequency domain.

As an illustration, consider a variable, $x (t)$ , which contains many discrete frequencies. We can write $x (t)$ in terms of a Fourier series expanded in N steps of a fundamental frequency $ω_{0}$ :

x (t) = \overset{N}{\sum_{n = 1}} [a_{n} \cos (i n ω_{0} t) + b_{n} \sin (i n ω_{0} t)] .

The series begins with the $n = 1$ term because the mean of the variable must be zero. We can write this series more compactly as a complex Fourier series (keeping in mind that we will be interested only in the real part):

x (t) = \overset{N}{\sum_{n = - N}} A_{n} \exp (i n ω_{0} t),

where $A_{n} = ℜ (A_{n}) + i ℑ (A_{n})$ is the complex amplitude of the nth term, $A_{n}^{*}$ is the complex conjugate of $A_{n}$ :

A_{n}^{*} = ℜ (A_{n}) - i ℑ (A_{n}),

and

A_{n} = {\begin{matrix} (a_{n} - i b_{n}) / 2 & for n > 0, \\ (a_{n} + i b_{n}) / 2 & for n < 0, \\ 0 & for n = 0 . \end{matrix}

The variance of $x (t)$ is

σ_{r}^{2} (x) = lim_{T \to \infty} \frac{1}{T} \int_{- \frac{T}{2}}^{\frac{T}{2}} \overset{N}{\sum_{n = - N}} A_{n}^{*} \exp (- i n ω_{0} t) A_{n} \exp (i n ω_{0} t) d t,

using the orthogonality of Fourier terms. Continuing,

\begin{matrix} σ_{r}^{2} (x) = ​ & lim_{T \to \infty} \frac{1}{T} \int_{- \frac{T}{2}}^{\frac{T}{2}} \overset{N}{\sum_{n = - N}} A_{n} A_{n}^{*} d t \\ = & \overset{N}{\sum_{n = - N}} A_{n} A_{n}^{*} = \overset{N}{\sum_{n = - N}} {| A_{n} |}^{2} = \overset{N}{\sum_{n = - N}} σ_{r}^{2} (x_{n}), \end{matrix}

where

x_{n} = ℜ (A_{n} \exp (i n ω_{0} t))

is the nth component of the Fourier series.

Thus, thanks to the orthogonality of Fourier terms, the variance (the mean square value) of the series is the sum of the variances (the mean square values) of its components. In particular, we see that $A_{n} A_{n}^{*}$ is the variance, or mean square value, of the variable at the frequency $n ω_{0}$ .

The contribution to the variance of x, $σ_{r}^{2} (x)$ , at the frequency $f_{n} = n ω_{0} / (2 π)$ , per unit frequency, is thus

S_{x} (f_{n}) = \frac{2 π}{ω_{0}} A_{n} A_{n}^{*},

since we are stepping up the frequency range in steps of $Δ f = ω_{0} / (2 π)$ . The variance can, therefore, be written as

σ_{r}^{2} (x) = \overset{N}{\sum_{n = - N}} S_{x} (f_{n}) Δ f .

As we examine x as a function of frequency, $S_{x} (f_{n})$ tells us the amount of “power” (in the sense of mean square value) contained in x, per unit frequency, at the frequency $f_{n}$ . As we consider smaller and smaller intervals, $Δ f \to 0$ , $S_{x}$ is the power spectral density (PSD) of the variable x:

σ_{r}^{2} (x) = \int_{- \infty}^{\infty} S_{x} (f) d f,

where f is the frequency in cycles per time (usually Hz).

Notice that $S_{x}$ has units of (variable)²/frequency, where (variable) is the unit of the variable (displacement, force, stress, etc.). In this case “frequency” is almost always given in Hz, although—since Abaqus does not have any built-in units—the frequency could be expressed in any other units of cycles per time. However, $S_{x}$ should not be given per circular frequency (radians per time): Abaqus assumes that $S_{x} (f)$ , not $S_{x} (ω)$ .

Fourier transforms

Since the variables of interest in random response analysis are characterized as functions of frequency, the Fourier transform plays a major role in converting from the time domain to the frequency domain and vice versa. The Fourier transform of $x (t)$ , which we write as $\bar{x} (f)$ , is defined by

x (t) = \int_{- \infty}^{\infty} \bar{x} (f) \exp (i 2 π f t) d f

or, in terms of the circular frequency $ω = 2 π f$ ,

x (t) = \frac{1}{2 π} \int_{- \infty}^{\infty} \bar{x} (ω) \exp (i ω t) d ω .

Simple manipulation provides

x (t) = \frac{1}{2 π} \int_{0}^{\infty} [(\bar{x} (ω) + \bar{x} (- ω)) \cos ω t + i (\bar{x} (ω) - \bar{x} (- ω)) \sin ω t] d ω,

so $(\bar{x} (ω) + \bar{x} (- ω)) / (2 π)$ is the amplitude of the cosine term in $x (t)$ at the circular frequency $ω$ , while $(\bar{x} (ω) - \bar{x} (- ω)) / (2 π)$ is the amplitude of the sine term. We, thus, see the physical meaning of the Fourier transform—it provides the complex magnitude (the amplitude and phase) of the content of $x (t)$ at a particular frequency.

If $x (t)$ is real only (which is the case for the variables we need to consider), this expression shows that we must have

ℑ (\bar{x} (ω)) = - ℑ (\bar{x} (- ω))

and

ℜ (\bar{x} (ω)) = ℜ (\bar{x} (- ω)),

which means that

\bar{x} (- ω) = {\bar{x}}^{*} (ω) .

For completeness we also note the inverse transformation,

\bar{x} (f) = \int_{- \infty}^{\infty} x (t) \exp (- i 2 π f t) d t,

which shows that the transformations between the time and frequency domains are rather symmetrical.

We now need Parseval's theorem:

\begin{matrix} \int_{- \infty}^{\infty} x_{1} (t) x_{2}^{*} (t) d t & ​ = \int_{- \infty}^{\infty} x_{2}^{*} (t) \int_{- \infty}^{\infty} {\bar{x}}_{1} (f) \exp (i 2 π f t) d f d t \\ ​ = \int_{- \infty}^{\infty} {\bar{x}}_{1} (f) \int_{- \infty}^{\infty} x_{2}^{*} (t) \exp (i 2 π f t) d t d f \\ ​ = \int_{- \infty}^{\infty} {\bar{x}}_{1} (f) {\bar{x}}_{2}^{*} (f) d f . \end{matrix}

Applying this theorem to the variance (the mean square value):

\begin{matrix} σ_{r}^{2} & ​ = lim_{T \to \infty} \frac{1}{T} \int_{- \frac{T}{2}}^{\frac{T}{2}} x^{2} (t) d t \\ ​ = \int_{- \infty}^{\infty} lim_{T \to \infty} \frac{1}{T} \bar{x} (f) {\bar{x}}^{*} (f) d f . \end{matrix}

Since $x (t)$ is real, we know that $\bar{x} (- f) = {\bar{x}}^{*} (f)$ , and so

\bar{x} (- f) {\bar{x}}^{*} (- f) = {\bar{x}}^{*} (f) \bar{x} (f) = \bar{x} (f) {\bar{x}}^{*} (f) .

We have already shown that we can write the variance in terms of the power spectral density as

σ_{r}^{2} (x) = \int_{- \infty}^{\infty} S_{x} (f) d f .

By comparison,

S_{x} (f) = lim_{T \to \infty} \frac{1}{T} \bar{x} (f) {\bar{x}}^{*} (f) .

We also see that

S_{x} (- f) = S_{x} (f) .

To avoid integration over negative frequencies, we write the variance as

\begin{matrix} σ_{r}^{2} (x) & ​ = \int_{0}^{\infty} [S_{x} (f) + S_{x} (- f)] d f \\ ​ = \int_{0}^{\infty} {\tilde{S}}_{x} (f) d f, \end{matrix}

where ${\tilde{S}}_{x} (f)$ is the single-sided PSD defined as

{\tilde{S}}_{x} (f) = S_{x} (f) + S_{x} (- f), for f \geq 0 .

Since $S_{x} (- f) = S_{x} (f)$ , we see that

{\tilde{S}}_{x} (f) = 2 S_{x} (f) .

Now consider the autocorrelation function:

\begin{matrix} R (τ) = lim_{T \to \infty} \frac{1}{T} \int_{- \frac{T}{2}}^{\frac{T}{2}} x (t) x (t + τ) d t \\ ​ = lim_{T \to \infty} \frac{1}{T} \int_{- \frac{T}{2}}^{\frac{T}{2}} x (t) [\int_{- \infty}^{\infty} \bar{x} (f) \exp (i 2 π f (t + τ)) d f] d t \\ ​ = \int_{- \infty}^{\infty} lim_{T \to \infty} \frac{1}{T} [\int_{- \infty}^{\infty} x (t) \exp (i 2 π f t) d t] \bar{x} (f) \exp (i 2 π f τ) d f \\ ​ = \int_{- \infty}^{\infty} (lim_{T \to \infty} \frac{1}{T} {\bar{x}}^{*} (f) \bar{x} (f)) \exp (i 2 π f τ) d f \\ ​ = \int_{- \infty}^{\infty} S_{x} (f) \exp (i 2 π f τ) d f, \end{matrix}

(using the result above). Thus, the power spectral density is the Fourier transform of the autocorrelation function. The inverse transform is

S_{x} (f) = \int_{- \infty}^{\infty} R (τ) \exp (- i 2 π f τ) d τ .

Since $R (τ)$ is symmetric about $τ = 0$ ( $R (- τ) = R (τ)$ ), we can also write this equation as

\begin{matrix} S_{x} (f) & ​ = \int_{0}^{\infty} R (τ) (\exp (- i 2 π f τ) + \exp (i 2 π f τ)) d τ \\ ​ = 2 \int_{0}^{\infty} R (τ) \cos (2 π f τ) d τ, \end{matrix}

so that

{\tilde{S}}_{x} (f) = 4 \int_{0}^{\infty} R (τ) \cos (2 π f τ) d τ .

Cross-spectral density

Following a similar argument to that used above to develop the idea of the power spectral density, we can define the cross-spectral density (CSD) function, $S_{x_{1} x_{2}} (f)$ , which gives the cross-correlation between two variables, $R_{x_{1} x_{2}} (τ)$ , as

R_{x_{1} x_{2}} (τ) = \int_{- \infty}^{\infty} S_{x_{1} x_{2}} (f) \exp (i 2 π f τ) d f,

with the inverse transformation

S_{x_{1} x_{2}} (f) = \int_{- \infty}^{\infty} R_{x_{1} x_{2}} (τ) \exp (- i 2 π f τ) d τ .

Transforming the original definition of $R_{x_{1} x_{2}} (τ)$ to the frequency domain provides

\begin{matrix} R_{x_{1} x_{2}} (τ) = lim_{T \to \infty} \frac{1}{T} \int_{- \frac{T}{2}}^{\frac{T}{2}} x_{1} (t) x_{2} (t + τ) d t \\ ​ = lim_{T \to \infty} \frac{1}{T} \int_{- \frac{T}{2}}^{\frac{T}{2}} x_{1} (t) \int_{- \infty}^{\infty} {\bar{x}}_{2} (f) \exp (i 2 π f (t + τ)) d f d t \\ ​ = \int_{- \infty}^{\infty} lim_{T \to \infty} \frac{1}{T} [\int_{- \infty}^{\infty} x_{1} (t) \exp (i 2 π f t) d t] {\bar{x}}_{2} (f) \exp (i 2 π f τ) d f \\ ​ = \int_{- \infty}^{\infty} (lim_{T \to \infty} \frac{1}{T} {\bar{x}}_{1}^{*} (f) {\bar{x}}_{2} (f)) \exp (i 2 π f τ) d f . \end{matrix}

By comparison,

S_{x_{1} x_{2}} (f) = lim_{T \to \infty} \frac{1}{T} {\bar{x}}_{1}^{*} (f) {\bar{x}}_{2} (f) .

We could also write

\begin{matrix} S_{x_{1} x_{2}} (f) & ​ = {[lim_{T \to \infty} \frac{1}{T} {\bar{x}}_{2}^{*} (f) {\bar{x}}_{1} (f)]}^{*} \\ = S_{x_{2} x_{1}}^{*} (f) = S_{x_{2} x_{1}} (- f) . \end{matrix}

Random response analysis

The general concept of random response analysis is now clear. A system is excited by some random loads or prescribed base motions, which are characterized in the frequency domain by a matrix of cross-spectral density functions, $S_{N M} (f)$ . Here we think of N and M as two of the degrees of freedom of the finite element model that are exposed to the random loads or prescribed base motions.

In typical applications the range of frequencies will be limited to those to which we know the structure will respond—we do not need to consider frequencies that are higher than the modes in which we expect the structure to respond.

The values of $S_{N M} (f)$ might be provided by Fourier transformation of the cross-correlation of time records or by the Fourier transformation of the autocorrelation of a single time record, together with known geometric data, as in the case of the car driving along a roughly grooved road, where the autocorrelation of the road surface profile, together with the speed of the car and the axle separation, allow $S_{12} (f)$ to be defined for the front (1) and rear (2) axles, as shown above. (If, in this case, the road profile seen by the wheels on the left side of the car is not similar to that seen by the wheels on the right side of the car, the cross-correlation of the left and right road surface profiles will also be required to define the excitation.)

The system will respond to this excitation. We are usually interested in looking at the power spectral densities of the usual response variables—stress, displacement, etc. The PSD history of any particular variable will tell us the frequencies at which the system is most excited by the random loading.

We might also compute the cross-spectral densities between variables. These are usually not of interest, and Abaqus/Standard does not provide them. (They might be needed if the analysis involves obtaining results that, in turn, will define the loading for some other system. For example, the response of a building to seismic loading might be used to obtain the motions of the attachment points for a piping system in the building so that the piping system can then be analyzed. The only option would be to model the entire system together.)

An overall picture is provided by looking at the variance (the mean square value) of any variable; the RMS value is provided for this purpose. The RMS value is used instead of variance because it has the same units as the variable itself. Abaqus/Standard computes it by integrating the single-sided power spectral density of the variable over the frequency range, since

σ_{r} = \sqrt{R (0)} = \sqrt{\int_{0}^{\infty} {\tilde{S}}_{x} (f) d f} .

This integration is performed numerically by using the trapezoidal rule over the range of frequencies specified for the random response step:

σ_{r} \approx \sqrt{\frac{1}{2} \overset{N}{\sum_{i = 2}} ({\tilde{S}}_{i} + {\tilde{S}}_{i - 1}) (f_{i} - f_{i - 1})},

where ${\tilde{S}}_{i}$ is the ${\tilde{S}}_{x}$ at the frequency $f_{i}$ and N is the number of points at which the response was calculated. N will depend on the number of eigenmodes used in the superposition and on the user-specified number of points between the eigenfrequencies.

The user must ensure that enough frequency points are specified so that this approximate integration will be sufficiently accurate.

The transformation of the problem into the frequency domain inherently assumes that the system under study is responding linearly: the random response procedure is considered as a linear perturbation analysis step.

What remains, then, is for us to consider how Abaqus/Standard finds the linear response to the random excitation.

The frequency response function

Random response is studied in the frequency domain. Therefore, we need the transformation from load to response as a function of frequency. Since the random response is treated as the integration of a series of sinusoidal vibrations, this transformation is based on the same steady-state response function used for a steady-state dynamic analysis and described in Steady-state linear dynamic analysis.

The discrete (finite element) linear dynamic system has the equilibrium equation

δ u^{N} (M^{N M} {\ddot{u}}^{M} + C^{N M} {\dot{u}}^{M} + K^{N M} u^{M}) = δ u^{N} F^{N},

where $M^{N M}$ is the mass matrix, $C^{N M}$ is the damping matrix, $K^{N M}$ is the stiffness matrix, $F^{N}$ are the external loads, $u^{N}$ is the value of degree of freedom N of the finite element model (usually a displacement or rotation component, or an acoustic pressure), and $δ u^{N}$ is an arbitrary virtual variation.

We project the problem onto the eigenmodes of the system. To do this the modes are first extracted from the undamped system:

(M^{N M} ω_{\bar{α}}^{2} - K^{N M}) ϕ_{\bar{α}}^{N} = 0.

(Here, and throughout the remainder of this section, repeated subscripts and superscripts are assumed to be summed over the appropriate range except when they are barred, like $\bar{α}$ above. Roman superscripts and subscripts indicate physical degrees of freedom; Greek superscripts and subscripts indicate modal variables.)

Typically the structural dynamic response is well represented by a small number of the lower modes of the model, so the number of modes is usually $O (10^{1})$ – $O (10^{2})$ , while the number of physical degrees of freedom might be $O (10^{3})$ – $O (10^{5})$ .

The eigenmodes are orthogonal across the mass and stiffness matrices:

\begin{matrix} ϕ_{α}^{N} M^{N M} ϕ_{β}^{M} & ​ = {\begin{matrix} m_{α} & if α = β \\ 0 & if α \neq β; \end{matrix} \\ ϕ_{α}^{N} K^{N M} ϕ_{β}^{M} & ​ = {\begin{matrix} k_{α} & if α = β \\ 0 & if α \neq β . \end{matrix} \end{matrix}

The eigenmodes can be displacement or mass normalized. For displacement normalization the largest entry in $ϕ_{α}^{N}$ is 1.0; for mass normalization the largest entry is $m_{α} = 1$ . In SIM architecture we normalize only with mass.

We assume that any damping is in the general form of “Rayleigh damping”:

C^{N M} = α M^{N M} + β K^{N M}

so that $C^{N M}$ will also project into a diagonal damping matrix $c_{α}$ .

The problem, thus, projects into a set of uncoupled modal response equations,

m_{\bar{α}} {\ddot{q}}_{\bar{α}} + c_{\bar{α}} {\dot{q}}_{\bar{α}} + k_{\bar{α}} q_{\bar{α}} = f_{α},

where

f_{α} = ϕ_{α}^{N} F^{N}

is the generalized load for mode $α$ and $q_{α}$ is the “generalized coordinate” (the modal amplitude) for mode $α$ .

Steady-state excitation is of the form

f_{α} = A \exp (i 2 π f t)

and creates response of similar form that we write as

q_{α} = H A \exp (i 2 π f t),

where $H (f)$ is the complex frequency response function defined in Steady-state linear dynamic analysis.

Response development

Random loading is defined by the cross-spectral density matrix $S_{N M}^{F} (f)$ , which links all loaded degrees of freedom (N and M). Projecting this matrix onto the modes provides the cross-spectral density function for the generalized (modal) loads:

S_{α β}^{f} (f) = ϕ_{α}^{N} S_{N M}^{F} (f) ϕ_{β}^{M} .

The complex frequency response function then defines the response of the generalized coordinates as

S_{α β}^{q} (f) = H_{\bar{α}} S_{\bar{α} \bar{β}}^{f} (f) H_{\bar{β}}^{*},

where $H_{α}^{*}$ is the complex conjugate of $H_{α}$ .

Finally, the response of the physical variables is recovered from the modal responses as

S_{N M}^{u} (f) = ϕ_{α}^{N} S_{α β}^{q} (f) ϕ_{β}^{M}

so that the power spectral density of degree of freedom $u^{N}$ is

S_{\bar{N} \bar{N}}^{u} (f) = ϕ_{α}^{\bar{N}} S_{α β}^{q} (f) ϕ_{β}^{\bar{N}} .

The PSDs for the velocity and acceleration of the same variable are

S_{\bar{N} \bar{N}}^{\dot{u}} (f) = {(2 π f)}^{2} S_{\bar{N} \bar{N}}^{u} (f)

and

S_{\bar{N} \bar{N}}^{\ddot{u}} (f) = {(2 π f)}^{4} S_{\bar{N} \bar{N}}^{u} (f) .

Recall that we may typically have $O (10^{1})$ – $O (10^{2})$ eigenmodes but many more ( $O (10^{3})$ – $O (10^{5})$ ) physical degrees of freedom. Therefore, if many of the physical degrees of freedom are loaded (as in the case of a shell structure exposed to random acoustic noise), it may be computationally expensive to perform operations such as

S_{α β}^{f} (f) = ϕ_{α}^{N} S_{N M}^{F} (f) ϕ_{β}^{M},

which involve products over all loaded physical degrees of freedom and must be done at each frequency in the range considered.

In contrast, operations such as

S_{α β}^{q} (f) = H_{\bar{α}} (f) S_{\bar{α} \bar{β}}^{f} (f) H_{\bar{β}}^{*} (f)

are relatively inexpensive since they are done independently for each combination of modes.

Forming

S_{\bar{N} \bar{N}}^{u} (f) = ϕ_{α}^{\bar{N}} S_{α β}^{q} (f) ϕ_{β}^{\bar{N}}

may be expensive if we choose to compute the results for a large selection of physical variables (displacements, velocities, accelerations, stresses, etc.). Abaqus/Standard will calculate the response only for the element and nodal variables requested. However, if a restart analysis is requested with the random response procedure, all variables are computed at the requested restart frequency, which can add substantially to the computational cost. The user is advised to write the restart file for the last increment only. To reduce the computational cost of random response analysis, Abaqus/Standard assumes that the cross-spectral density matrix for the loading can be separated into a frequency-dependent scalar function (containing the units of the CSD) and a set of coupling terms that are independent of frequency as follows:

\begin{matrix} S_{N M}^{F} (f) & ​ = \sum_{J} P^{J} (f) \sum_{I} Ψ_{N M}^{I J} for N < M, \\ S_{N M}^{F} (f) & ​ = {[S_{M N}^{F} (f)]}^{*} for N > M, \\ S_{N N}^{F} (f) & ​ = ℜ (\sum_{J} P^{J} (f) \sum_{I} Ψ_{N N}^{I J}) . \end{matrix}

Here J is the number of cross-correlation definitions included in the random response step. Each cross-correlation definition refers to an input complex frequency function, $P^{J} (f)$ . The spatial cross-correlations are then defined by the complex set of values $Ψ_{N M}^{I J}$ .

With this approach the cross-spectral density function for the generalized loads, $S_{α β}^{f} (f)$ , can be constructed as

\begin{matrix} S_{α β}^{f} (f) & ​ = \sum_{J} P^{J} (f) (ϕ_{α}^{N} ϕ_{β}^{M} \sum_{I} Ψ_{N M}^{I J}) for N < M, \\ S_{α β}^{f} (f) & ​ = {[S_{β α}^{f} (f)]}^{*} for N > M, \\ S_{α α}^{f} (f) & ​ = ℜ [\sum_{J} P^{J} (f) (ϕ_{α}^{N} ϕ_{α}^{M} \sum_{I} Ψ_{N M}^{I J})] . \end{matrix}

Since the $ϕ_{α}^{N} ϕ_{β}^{M} \sum_{I} Ψ_{N M}^{I J}$ are not functions of frequency, they can be computed once only, leaving the frequency-dependent operations to be done only in the space of the eigenmodes.

Although this procedure is not natural for typical correlated loadings (like road excitation or jet noise), loadings can always be defined this way by using enough cross-correlation definitions. The approach then reduces the computational cost for models with many loaded physical degrees of freedom. The approach works well for uncorrelated and fully correlated loadings, which are quite common cases.

Decibel conversion

Abaqus/Standard allows the user to provide an input PSD (say $P (f)$ ) in decibel units rather than units of power/frequency. There are various ways to convert from decibel units to units of power/frequency, depending on how the frequencies in one octave band are related to the frequencies in the next. A general formula relates the center (midband) frequencies $f_{c}$ between octaves as

\frac{f_{c}^{(i)}}{f_{c}^{(i - 1)}} = 2^{x},

where the superscript $(i)$ denotes $i th$ octave band and x is a chosen value. For example, $x = 1$ for full octave band conversion, and $x = 1 / 3$ for one-third octave band conversion. Abaqus/Standard uses full octave band conversion to convert from decibel units to units of power/frequency. For full octave band conversion, as shown by the above equation, the center frequency doubles from octave band to octave band.

Since decibel units are based upon log scales, the center frequency for an octave band bounded by lower frequency $f_{L}$ and upper frequency $f_{U}$ is given by

\log f_{c} = \frac{1}{2} \log f_{U} + \frac{1}{2} \log f_{L} .

Since $f_{U} / f_{L} = 2^{x}$ , we can easily show that

f_{U} = f_{c} 2^{0.5 x}

and

f_{L} = f_{c} 2^{- 0.5 x} .

Thus, the change in frequency within any given octave band is

Δ f = f_{U} - f_{L} = f_{c} [2^{0.5 x} - 2^{- 0.5 x}] .

To convert from one type of conversion formula to the next, we need the following general decibel to power/frequency conversion equation:

d b (f) = 10 \log \frac{P (f)}{P_{ref}^{(x)} / Δ f},

where $P_{ref}^{(x)}$ is a reference power value. The subscript $(x)$ means that the power reference is given for the type of conversion represented by x (e.g., full octave band conversion for $x = 1$ or one-third octave band conversion for $x = 1 / 3$ ). When $x = 1$ , we will simply use the notation $P_{ref}$ . Thus, since Abaqus/Standard uses a full octave band conversion, $Δ f = f_{c} / \sqrt{2}$ and

P (f) = \frac{\sqrt{2}}{f_{c}} P_{ref} 10^{d b (f) / 10} .

The PSD data can be given with respect to some other type of octave band frequency scale. In that case we can convert the PSD data at those frequencies coinciding with the full octave band scale by computing an equivalent full octave band reference power based on the following ratio:

\frac{P_{ref}}{P_{ref}^{(x)}} = \frac{1}{\sqrt{2} [2^{0.5 x} - 2^{- 0.5 x}]} .

For example, if we are given $P_{ref}^{(1 / 3)}$ (i.e., one-third octave band frequency scale), the equivalent full octave band reference power value would be

P_{ref} = 3.0536 P_{ref}^{(1 / 3)} .

This conversion would be valid only at the one-third octave band center frequencies that coincide with the full octave band center frequencies. Thus, only every third data point should be considered.

von Mises stress computation

The computation of PSD and RMS of von Mises stresses in Abaqus is based on work by Segalman, et al. (1998). Under this approach the PSD of von Mises stress at a node a is given by

S_{m i s e s}^{a} (f) = \overset{m}{\sum_{β = 1}} \overset{m}{\sum_{α = 1}} S_{α β} (f) T_{α β}^{a},

where m is the number of modes, $S_{α β} (f)$ are the elements of the PSD matrix of generalized displacements,

T_{α β}^{a} = {[ψ_{α}^{a}]}^{T} [A] [ψ_{β}^{a}],

$[ψ_{α}^{a}]$ are the modal stress components of the $α$ th mode at node a, and the constant matrix A is given by

(\begin{matrix} 1 & - 1 / 2 & - 1 / 2 & 0 & 0 & 0 \\ - 1 / 2 & 1 & - 1 / 2 & 0 & 0 & 0 \\ - 1 / 2 & - 1 / 2 & 1 & 0 & 0 & 0 \\ 0 & 0 & 0 & 3 & 0 & 0 \\ 0 & 0 & 0 & 0 & 3 & 0 \\ 0 & 0 & 0 & 0 & 0 & 3 \end{matrix}) .

Similarly, RMS of von Mises stress at a node a is computed as

R M I S E S^{a} (f) = \sqrt{\sum_{β = 1}^{m} \sum_{α = 1}^{m} V_{α β} (f) T_{α β}^{a}},

where $V_{α β} (f)$ are the elements of the variance matrix of generalized displacements.