On the analytic representation of the correlation function of linear random vibration systems

J. Gruner  J. vom Scheidt  R. Wunderlich

Abstract

This paper is devoted to the computation of statistical characteristics of the response of discrete vibration systems with a random external excitation. The excitation can act at multiple points and is modeled by a time-shifted random process and its derivatives up to the second order. Statistical characteristics of the response are given by expansions as to the correlation length of a weakly correlated random process which is used in the excitation model. As the main result analytic expressions of some integrals involved in the expansion terms are derived.

1 Introduction

Mathematical modeling of real-world vibration systems (e. g. vehicles moving on a rough guideway, rotating generator shafts excited by random fluctuations of the generator torque) results in a system of ordinary differential equations (ODE) containing random parameters. For a linear system with n degrees of freedom (DOF) and a random external excitation we get

Ap̈ + Bṗ + Cp = f(t, ω)  (1)

with the initial conditions p(t_{0}) = p_{0} , ṗ(t_{0}) = p_{1} . Thereby the n-dimensional vector p contains the response variables, A, B, C are n × n matrices and f(., ω) is a n-dimensional random vector function representing the random excitation terms.

In the case of vehicle vibrations we have to model time-shifted excitations by the road irregularities at multiple points of the vehicle suspension (spring and damper combinations). This can be done assuming the right hand side of (1) to have the representation

f(t, ω) = ∑_{l=0}^{2} P_{l}f̂^{(l)}(t, ω)  (2)
with  f̂(t, ω) = (g(t + τ_{1}, ω), …, g(t + τ_{r}, ω), 0, …, 0)^{τ},

where P_{l}, l = 0, 1, 2, are real n × n matrices, τ_{p}, p = 1, …, r, r ≤ n, are the time shifts and g(t, ω) denotes the time dependent random road profile.

Since the right hand side of (1) contains random terms (which are explained in Section 2) the response p is a random function, too. Section 3 deals with the problem of computation of statistical characteristics (such as mean, variance, correlation and spectral density function) of the random response p = p(t, ω) to given characteristics of the random excitation term f = f(t, ω). For the solution of that problem we use an expansion method presented in [5, 6] from where the ideas and the notations are taken. For a modified expansion method see [7]. Section 4 investigates the analytic representation of some integrals involved in the expansion coefficients. With the help of such representations (computer) time expensive numerical integration can be avoided. They allow the treatment of random vibration systems with larger numbers of DOF.

For numerical results for linear and nonlinear vehicle vibrations we refer to [6, 10, 9] and for linear torsional vibrations of a generator shaft to [11]. The case of an excitation of vehicles by two correlated random tracks is considered in [1, 8]

2 Random Excitation

While the structure of the excitation term f on the right hand side of (1) is given by (2) now the statistical properties of the random profile g(t, ω) are investigated.

In many applications weakly stationary random excitations (e.g. road profiles) are modeled by random processes with a correlation function

R_{gg}(t) = E {g(s)g(s + t)} = σ^{2}e^{−γ|t|} ,  γ > 0,  (3)

where E {.} denotes the expectation, and a corresponding spectral density

S(α) = \frac{σ^{2}}{π} \frac{γ}{γ^{2} + α^{2}}

(see [2], [3], [4]). The correlation function (3) is not differentiable at t = 0, and therefore these processes need not to be differentiable. Hence, processes of this type are not suitable as excitations in (1) if derivatives of the excitations are included.

To obtain differentiable processes of a similar type the process g(t, ω) is expressed by a integral functional of a weakly correlated and weakly stationary process h_{ε}(t, ω) (see [5],[6])

g(t, ω) = ∫_{−∞}^{t} Q(t − u)h_{ε}(u, ω) du  (4)

where Q(t) is a twice continuously differentiable function with

Q(0) = 0 ,  Q′(0) = 0 ,  Q′′(0) = 0.

Further, Q(.) and its first and second derivative are assumed to be square-integrable. Hence, the process g(t, ω) is twice continuously differentiable with an uniform representation of the process and its derivatives,

g^{(l)}(t, ω) = ∫_{−∞}^{t} Q^{(l)}(t − u)h_{ε}(u, ω) du, l = 0, 1, 2.

The essential property of a weakly correlated function h_{ε}(s, ω) is that the influence of the values of the function does not reach far. Values of the function at two points do not correlate (are independent) if the distance of the two points exceeds the correlation length ε. The quantity ε is assumed to be sufficiently small to apply expansions as to ε for approximations of the distribution of linear functionals as (4) (see [5]). For the correlation function we get the expansion

R_{gg}(t, s) = E {g(t)g(s)}
= ∫_{−∞}^{t} ∫_{−∞}^{s} Q(t − u)Q(s − v)E {h_{ε}(u)h_{ε}(v)} du dv
= ε · a ∫_{−∞}^{min(t,s)} Q(t − u)Q(s − u) du + o(ε).  (5)

Here, a is the so-called intensity defined by

a = lim_{ε↓0} \frac{1}{ε} ∫_{−ε}^{ε} E \{h_{ε}(t)h_{ε}(t + u)\} du.

If the non-random kernel function in (4) is chosen as

Q(t) = Q_{0}(t)e^{−γt}, γ > 0,

where for δ > 0

Q_{0}(t) = \{ \begin{matrix}0 & for & t < 0 \\ 6 (\frac{t}{δ} )^{5} − 15 (\frac{t}{δ} )^{4} + 10 (\frac{t}{δ} )^{3} & for & 0 ≤ t ≤ δ \\ 1 & for & t > δ, \end{matrix}  (6)

we get an excitation which is centered and wide-sense stationary with an autocorrelation function for which the expansion (5) gives

R_{gg}(τ) = E {g(s)g(s + t)}
= ε · (\frac{a}{2γ} e^{−γ|t|} + O(δ)) + o(ε)
≈ ε · \frac{a}{2γ} e^{−γ|t|} for ε, δ small.

This approximation of the correlation function corresponds to the model (3) if we put σ^{2} = \frac{aε}{2γ}. Therefore, this method can also be interpreted as a smoothing of the correlation function (3) in a δ-neighbourhood of t = 0. The parameters of Q(t) can be determined by comparison with results from measurements (see [6]).

As can be seen in the next Section, the response variables can be represented as sums of integral functionals of weakly correlated functions as above in (4) with appropriate kernel functions. Hence, the correlation function (as well as other characteristics) of the response can be expanded as to ε in a similar way.

3 Linear Random Vibration Systems

Now we consider the solution of the linear system of second-order ODEs (1). By forming the state vector z = (ṗ, p)^{τ} ∈ R^{2n} the system (1) is transformed into a system of 2n first-order ODEs

ż(t) = R z(t) + D \binom{f(t, ω)}{0}, z(t_{0}) = z_{0},  (7)

with 2n × 2n matrices

R = − ( \begin{matrix}A & O \\ O & I_{n} \end{matrix} )^{−1} ( \begin{matrix}B & C \\ −I_{n} & O \end{matrix} ) and D = ( \begin{matrix}A & O \\ O & I_{n} \end{matrix} )^{−1} ,

where A is assumed to be regular, I_{n} denotes the n × n identity matrix and z_{0} = (p_{1}, p_{0})^{τ}.

We assume stability of the system (7), i.e. the eigenvalues of A have negative real parts only. Furthermore, the eigenvalues are assumed to be simple. Then there exists the decomposition

R = VΛV^{−1} with Λ = diag(λ_{1}, …, λ_{2n}), Reλ_{i} < 0, λ_{i} ≠ λ_{j},

for i, j = 1, …, 2n, i ≠ j, and a 2n × 2n-matrix V containing the eigenvectors of R. By the substitution z = Vx the system (7) is transformed into a decoupled complex system of first-order ODE

ẋ(t) = Λ x(t) + f̃(t, ω), x(t_{0}) = x_{0} = V^{−1}z_{0},  (8)

with f̃(t, ω) = V^{−1}D \binom{f(t, ω)}{0}.

To study the statistical properties of the system response x(t, ω) in terms of stationary random vector processes we consider the so-called stationary solution of (8), that is the solution for t_{0} → −∞. It is found to be

x(t, ω) = ∫_{−∞}^{t} e^{Λ(t−s)}f̃(s, ω)ds.  (9)

Equation (9) can be used for the stochastic analysis of the response. For the response mean we get E {x(t)} = 0. It is equal to the solution of the so-called averaged problem, which is obtained from (8) by replacing the random input term f̃(t, ω) by its (non-random) expectation E {f̃(t)} = 0.

For the evaluation of the variances and correlation functions we have to consider the fluctuation x(t, ω) − E {x(t)} = x(t, ω) − 0. Substituting our excitation models (2) and (4) we find for the i-th response variable and its first derivative x_{i}^{(k)}(t, ω), i = 1, … 2n, k = 0, 1, (see [6])

x_{i}^{(k)}(t, ω) = ∑_{q=1}^{r} ∫_{−∞}^{t+τ_{q}}G_{iq}^{(k)}(t − u + τ_{q})h_{ε}(u, ω)du,  (10)

with the matrix function

G^{(k)}(t) = ∑_{l=0}^{2} ∫_{0}^{t} e^{Λ(t−v)} P̃_{l} Q^{(l+k)}(v)dv  (11)

where  P̃_{l} = V^{−1}D ( \begin{matrix}P_{l} & O \\ O & O \end{matrix} ) , l = 0, 1, 2.  (12)

The right hand side of (10) is a sum of integral functionals of h_{ε}(s, ω) for which we can apply the expansion (5) to find the following approximation of the correlation function

R_{x_{j_{1}}^{(k_{1})}x_{j_{2}}^{(k_{2})}}(t_{1}, t_{2}) = E \{x_{j_{1}}^{(k_{1})} \overline{x}_{j_{2}}^{(k_{2})}\} = ε · R_{x_{j_{1}}^{(k_{1})}x_{j_{2}}^{(k_{2})}}^{1}(t_{1}, t_{2}) + o(ε)  (13)

where the first expansion term is found to be

R_{x_{j_{1}}^{(k_{1})}x_{j_{2}}^{(k_{2})}}^{1}(t_{1}, t_{2}) = a ∑_{q_{1}=1}^{r} ∑_{q_{2}=1}^{r} ∫_{−∞}^{min(t_{1}+τ_{q_{1}},t_{2}+τ_{q_{2}})}G_{j_{1}q_{1}}^{(k_{1})}(t_{1} − u + τ_{q_{1}})\overline{G}_{j_{2}q_{2}}^{(k_{2})}(t_{2} − u + τ_{q_{2}})du

for k_{1}, k_{2} = 0, 1 und j_{1}, j_{2} = 1, ..., 2n and \overline{(.)} denotes complex conjugate.

For t_{1} + τ_{q_{1}} ≤ t_{2} + τ_{q_{2}} we find for the integral

∫_{−∞}^{t_{1}+τ_{q_{1}}}G_{j_{1}q_{1}}^{(k_{1})}(t_{1}−u+τ_{q_{1}})\overline{G_{j_{2}q_{2}}^{(k_{2})}}(t_{2}−u+τ_{q_{2}})du = ∫_{0}^{∞}G_{j_{1}q_{1}}^{(k_{1})}(u)\overline{G_{j_{2}q_{2}}^{(k_{2})}}(t_{2}+τ_{q_{2}} −t_{1}−τ_{q_{1}} +u)du.

For t_{1} + τ_{q_{1}} ≥ t_{2} + τ_{q_{2}} follows

∫_{−∞}^{t_{2}+τ_{q_{2}}}G_{j_{1}q_{1}}^{(k_{1})}(t_{1}−u+τ_{q_{1}})\overline{G_{j_{2}q_{2}}^{(k_{2})}}(t_{2}−u+τ_{q_{2}})du = ∫_{0}^{∞}G_{j_{1}q_{1}}^{(k_{1})}(t_{1}+τ_{q_{1}} −t_{2}−τ_{q_{2}} +u)\overline{G_{j_{2}q_{2}}^{(k_{2})}}(u)du.

From the last two equations it can be seen that the first expansion term in (13) is a function of t = t_{2} − t_{1} only and can be represented by

R_{x_{j_{1}}^{(k_{1})}x_{j_{2}}^{(k_{2})}}^{1}(t) = a ∑_{q_{1}=1}^{r} ∑_{q_{2}=1}^{r} T_{j_{1},j_{2} q_{1},q_{2}}^{k_{1},k_{2}}(t + τ_{q_{2}} − τ_{q_{1}}).

with

T_{j_{1},j_{2} q_{1},q_{2}}^{k_{1},k_{2}}(t) = \{ \begin{matrix}∫_{0}^{∞}G_{j_{1}q_{1}}^{(k_{1})}(u − t)\overline{G_{j_{2}q_{2}}^{(k_{2})}}(u)du & for t ≤ 0 \\ ∫_{0}^{∞}G_{j_{1}q_{1}}^{(k_{1})}(u)\overline{G_{j_{2}q_{2}}^{(k_{2})}}(u + t)du & for t ≥ 0. \end{matrix}  (14)

The corresponding correlation function of the original system (7) can be found by applying the inverse transformation. Using matrix notation we get

R_{zz}(t) = E {z(s)z^{τ}(s + t)} =_{z=Vx} V R_{xx}(t) V^{∗}.

R_{zz}(τ) contains the correlation functions for the response variables p_{i}(t, ω) and its first derivatives ṗ_{i}(t, ω), i = 1, …, n. Spectral density functions can be computed applying the Fourier tranformation.

4 Computation of the matrix function G^{(k)}(.)

For the approximation of the response correlation function using the expansion method presented in the last Section the integrals in (14) have to be evaluated. These integrals contain the product of two elements of the matrix function G^{(k)}(.), k = 0, 1. Especially in the case of large systems numerical integation of these integral as well as of the integral involved in G^{(k)}(.) itself (see (11)) can be very time consuming.

Since the integrand of G^{(k)}(.) contains a product of an exponential function and a polynomial it is possible to derive analytic expressions for these integrals. Subsequently, for the case of the integral involved in G^{(k)} this derivation is presented.

From (11) results

G^{(k)}(t) = ∑_{l=0}^{2} ∫_{0}^{t} e^{Λ(t−v)} Q^{(l+k)}(v) dv P̃_{l}

= ∑_{l=0}^{2} F_{(l+k)}(t)P̃_{l}  (15)

where

F_{(l)}(t) = ∫_{0}^{t} e^{Λ(t−v)} Q^{(l)}(v) dv, l = 0, … 3.

The matrices P̃_{l} are given in (12) and

Q(t) = Q_{0}(t)e^{−γt},

with Q_{0}(t) given in (6) which is a polynomial of degree 5 in [0, δ], δ > 0. To evaluate F_{(l)}(t), l = 0, …, 3 we consider the following decomposition

F_{(l)}(t) = ∫_{0}^{t} e^{Λ(t−v)} Q^{(l)}(v) dv

= η_{δ}(t) ^{2}F_{(l)}(t, t) + (1 − η_{δ}(t)) (^{2}F_{(l)}(t, δ) +^{1}F_{(l)}(t)) .  (16)

Here, the function η_{δ}(t) is defined by

η_{δ}(t) := \{ \begin{matrix}1 & for δ ≥ t \\ 0 & for δ < t \end{matrix} ,  (17)

hence

(1 − η_{δ}(t)) = \{ \begin{matrix}0 & for δ ≥ t \\ 1 & for δ < t \end{matrix} .

Further (16) contains the complex diagonal matrix functions

1F_{(l)}(t) = ∫_{δ}^{t} e^{Λ(t−v)} Q^{(l)}(v) dv  (18)

2F_{(l)}(t, β) = ∫_{0}^{β} e^{Λ(t−v)} Q^{(l)}(v) dv,  (19)

with β = min(t, δ).

For the diagonal entries ^{1}F_{j(l)}(t), j = 1, .., 2n, l = 0, …, 3, of the first function we find

1F_{j(l)}(t) = ∫_{δ}^{t} e^{λ_{j}(t−v)} Q^{(l)}(v) dv

= ∫_{δ}^{t} e^{λ_{j}(t−v)} (−γ)^{l} e^{−γv} dv

= (−γ)^{l} e^{λ_{j}t} ∫_{δ}^{t} e^{−(λ_{j}+γ) v} dv

= \{ \begin{matrix}(−γ)^{l} e^{−γ t} (λ_{j} + γ)^{−1} (e^{(λ_{j}+γ) (t−δ)} − 1) & for λ_{j} ≠ −γ \\ (−γ)^{l} e^{λ_{j}t} (t − δ) & for λ_{j} = −γ. \end{matrix} (20)

For the matrix function ^{2}F_{(l)}(t, β), l = 0, …, 3, β = min(t, δ) we derive by l-fold partial integration and using of Q^{(l)}(0) = 0, l = 0, 1, 2,

2F_{(l)}(t, β) = ∫_{0}^{β} e^{Λ(t−v)} Q^{(l)}(v) dv

= e^{Λt} ∫_{0}^{β} e^{−Λv} Q^{(l)}(v) dv

= e^{Λt} (e^{−Λv}Q^{(l−1)}(v)|_{0}^{β} − ∫_{0}^{β} −Λ e^{−Λv} Q^{(l−1)}(v) dv)

= e^{Λt} (e^{−Λβ}Q^{(l−1)}(β) + Λ ∫_{0}^{β} e^{−Λv} Q^{(l−1)}(v) dv)

⋮

= e^{Λt} (e^{−Λβ} ∑_{k=1}^{l} Λ^{k−1} Q^{(l−k)}(β) + Λ^{l} ∫_{0}^{β} e^{−Λv} Q(v) dv) . (21)

Since Q(t) is a product of an exponetial function and a polynomial in [0, δ] we derive for the remaining integral in (21) for λ_{j} ≠ −γ by 6-fold partial integration

∫_{0}^{β} e^{−λ_{j}v} Q(v) dv = ∫_{0}^{β} e^{−(λ_{j}+γ)v} Q_{0}(v) dv

= − (λ_{j} + γ)^{−1} e^{−(λ_{j}+γ)v} Q_{0}(v)|_{0}^{β} − ∫_{0}^{β} − (λ_{j} + γ)^{−1} e^{−(λ_{j}+γ)v} Q_{0}^{(1)}(v) dv

⋮

= ∑_{k=0}^{5} − (λ_{j} + γ)^{−k−1} e^{−(λ_{j}+γ)v} Q_{0}^{(k)}(v)|_{0}^{β} + (λ_{j} + γ)^{−6} ∫_{0}^{β} e^{−(λ_{j}+γ)v} Q_{0}^{(6)}(v) dv

= −e^{−(λ_{j}+γ)β} (λ_{j} + γ)^{−1} ∑_{k=0}^{5} (λ_{j} + γ)^{−k} Q_{0}^{(k)}(β)

+ (λ_{j} + γ)^{−1} ∑_{k=3}^{5} (λ_{j} + γ)^{−k} Q_{0}^{(k)}(0),  (22)

since Q_{0}^{(6)}(v) ≡ 0 and Q_{0}^{(p)}(0) = 0 for p = 0, 1, 2. In the case λ_{j} = −γ we have

∫_{0}^{β} e^{−λ_{j}v} Q(v) dv = ∫_{0}^{β} e^{−(λ_{j}+γ)v} Q_{0}(v) dv

= ∫_{0}^{β} 6 (\frac{v}{δ} )^{5} − 15 (\frac{v}{δ} )^{4} + 10 (\frac{v}{δ} )^{3} dv

= \frac{β^{6}}{δ^{5}} − \frac{3β^{5}}{δ^{4}} + \frac{10β^{4}}{4δ^{3}} .  (23)

Substitution of (21), (22), (23) into (19) and (20) into (18) leads to the analytic representation of the matrix function F_{(l)}(t) in (16) and it follows

Theorem 1 For the matrix function G^{(k)}, k = 0, 1, defined in (11) it holds

G^{(k)}(t) = ∑_{l=0}^{2} F_{(l+k)}(t)P̃_{l}

with P̃_{l} given in (12) and the diagonal matrix function

F_{(l)}(t) = η_{δ}(t) ^{2}F_{(l)}(t, t) + (1 − η_{δ}(t)) (^{2}F_{(l)}(t, δ) +^{1}F_{(l)}(t))

where the diagonal entries for j = 1, …, 2n are given by

1F_{j(l)}(t) = \{ \begin{matrix}(−γ)^{l} e^{−γ t} (λ_{j} + γ)^{−1} (e^{(λ_{j}+γ) (t−δ)} − 1) & for λ_{j} ≠ −γ \\ (−γ)^{l} e^{λ_{j}t} (t − δ) & for λ_{j} = −γ \end{matrix}

and

2F_{j(l)}(t, β) = \{ \begin{matrix}e^{λ_{j}t} e^{−λ_{j}β} ∑_{k=1}^{l} λ_{j}^{k−1} Q^{(l−k)}(β) \\ −e^{λ_{j}t} λ_{j}^{l} e^{−(λ_{j}+γ)β} ∑_{k=0}^{5} (λ_{j} + γ)^{−k−1} Q_{0}^{(k)}(β) \\ +e^{λ_{j}t} λ_{j}^{l} ∑_{k=3}^{5} (λ_{j} + γ)^{−k−1} Q_{0}^{(k)}(0) & for λ_{j} ≠ −γ \\ e^{λ_{j}t} e^{−λ_{j}β} ∑_{k=1}^{l} λ_{j}^{k−1} Q^{(l−k)}(β) \\ +e^{λ_{j}t} λ_{j}^{l} (\frac{β^{6}}{δ^{5}} − \frac{3β^{5}}{δ^{4}} + \frac{10β^{4}}{4δ^{3}} ) & for λ_{j} = −γ \end{matrix}

and η_{δ}(t) is defined in (17).

By means of this theorem an analytic expression of the stationary response variables and its first derivatives given in (10) is available.

5 Conclusions

Expansions of statistical characteristics of the response of discrete vibration systems with an external random excitation have been considered. A random process which generates the time-shifted excitation at multiple points has been modeled by an integral functional of a weakly correlated random process with a correlation length ε. The expansions are given with respect to this correlation length for ε ↓ 0.

To avoid time expensive numerical integration the system of second order ODEs has been transformed into complex first order system with a diagonal sytem matrix by modal decompostion. Futher, analytic expressions of some integrals involved in the expansion terms have been found. They allow the analytic representation of the stationary response.

References

[1] J. Gruner, J. vom Scheidt, and R. Wunderlich. Numerical treatment of nonlinear random vibrations of vehicles. In Proceedings of 5th Mini Conf. on Vehicle System Dynamics, Identification and Anomalies, Budapest, Nov. 11-13, 1996, 9 pages, to appear.

[2] M. Mitschke. Dynamik der Kraftfahrzeuge, Band B: Schwingungen. Springer, Berlin, 1997.

[3] K. Popp and W. Schiehlen. Fahrzeugdynamik. Teubner, Stuttgart, 1993.

[4] T.T. Soong and M. Grigoriu. Random vibration of mechanical and structural systems. Prentice Hall, Englewood Cliffs NJ, 1993.

[5] J. vom Scheidt. Stochastic Equations of Mathematical Physics. Akademie– Verlag, Berlin, 1990.

[6] J. vom Scheidt, B. Fellenberg, and U. Wöhrl. Analyse und Simulation stochastischer Schwingungssysteme. Leitfäden der angewandten Mathematik und Mechanik, Band 71. B. G. Teubner, Stuttgart, 1994.

[7] J. vom Scheidt, H.-J. Starkloff, and R. Wunderlich. Asymptotic expansions of second-order moments of integral functionals of weakly correlated random functions. Preprint 97-17, TU Chemnitz-Zwickau, 1997.

[8] J. vom Scheidt and U. Wöhrl. Nonlinear vibration systems with two parallel random excitations. Journal for Analysis and its Applications, 16(1):217–228, 1997.

[9] J. vom Scheidt and R. Wunderlich. Nonlinear random vibrations of vehicles. Z. Angew. Math. Mech., 76-S5:447–448, 1996.

[10] J. vom Scheidt and R. Wunderlich. Nonlinear stochastic vibrations of vehicles. In H. Neunzert, editor, Progress in Industrial Mathematics at ECMI 1994, pages 521 – 528, New York, Stuttgart, 1996. Wiley, Teubner.

[11] R. Wunderlich, J. Gruner, and J. vom Scheidt. Model reduction of random vibration systems. In Proceedings of the 9th Conference of ECMI, Copenhagen, June 25-29, 1996, 8 pages to appear.