Projection methods for stochastic structural dynamics

A set of novel hybrid projection approaches are proposed for approximating the response of stochastic partial differential equations which describe structural dynamic systems. An optimal basis for the response of a stochastic system has been computed from the eigen modes of the parametrized structural dynamic system. The hybrid projection methods are obtained by applying appropriate approximations and by reducing the modal basis. These methods have been further improved by an implementation of a sample based Galerkin error minimization approach. In total four methods are presented and compared for numerical accuracy and efficiency by analysing the bending of a Euler-Bernoulli cantilever beam.


Introduction
Propagation of uncertainties in complex engineering dynamical systems is receiving increasing attention. When uncertainties are taken in to account, the equations of motion of discretised dynamical systems can be expressed by coupled ordinary differential equations with stochastic coefficients. The computational cost for the solution of such system mainly depends on the number of degrees of freedom and number of random variables. Among various numerical methods developed for such systems, the polynomial chaos based Galerkin projection approach [1] shows significant promise because it is more accurate compared to the classical expansion based methods [2,3] and computationally more efficient compared to the Monte Carlo simulation based methods [4]. However, the computational cost increases significantly with the number of random variables and the results tend to become less accurate for longer length of time. In this study, a set of novel projection approaches are developed to address these issues. The rationale behind proposing a set of methods is to analyse the effect of altering the nature of the coefficients and vectors associated with projection methods. Altering the nature of the coefficients and the vectors could result in great computational savings, however, it could also induce additional error.

The stochastic finite element method
The stochastic parameter associated with structures which are governed by stochastic hyperbolic partial differential equations can be characterised by the random parameter a(x, θ) on a bounded domain and a probability space (Θ, F , P). Through the use of the well established stochastic finite element method, a set of discretised equations can be obtained to represent and approximate the response of the governing stochastic equations.

Discretisation of the random field
In order to proceed with the stochastic finite element method, it is necessary to discretise the random field that is associated with the governing stochastic equations. In this study, we consider the stochastic parameter a(x, θ) to be a Gaussian random field with a autocovariance kernel C a : D × D → R defined on the domain D. The covariance function is positive definite, symmetric and square bounded. The random field a(x, θ) can be expressed by a truncated Karhunen-Loève series expansion. This expansion is achieved by performing a spectral decomposition of the covariance function of the field, thus resulting in where a 0 (x) is the mean function of the random field and � ξ i are uncorrelated standard random variables. As the random field under consideration is Gaussian, the random variables are deemed as uncorrelated standard normal random variables with zero mean and an unit variance. � λ i and � φ i (x) correspond to the eigenvalues and eigenvectors which satisfy the following Fredholm integral equation of the second kind � If the eigenvalues rapidly decay, the value of M can be kept relatively small in order to obtain an accurate depiction of the Gaussian random field.

The equation of motion of the stochastic system
The system matrices inherit the randomness of the input stochastic parameters and hence the stochastic linear system for structural dynamics takes the form of where M(θ), C(θ) and K(θ) are the random mass, damping and stiffness matrices respectively, u(t; θ) and its dotted variants are the system response vector and its time derivatives respectively, and f 0 (t) is the deterministic forcing vector. Following from the discretized spectral representation of the random field in (1), the system matrices can be expanded in terms of the functions of the input random variables as Here the mass and stiffness matrices have been expressed in terms of their deterministic components (M 0 and K 0 ) and the corresponding random contributions (M i and K i ) obtained from discretizing the mass and stiffness parameters with finite number of random variables (µ i (θ) and ν i (θ)). The total number of random variables utilized to represent the stochastic system matrices is M = p 1 + p 2 . For most cases the damping parameter is expressed as linear combination of the mass matrix and the system stiffness matrix which is the 'proportional damping model' and this has been adopted in the present work.

Formulating dynamic systems in the frequency domain
The methods for obtaining the discretised random form of the governing partial differential equations are well-established in stochastic finite element literature. By utilising the stochastic finite element method, the response a multiple-degrees-of-freedom structural vibration problem in the frequency domain can be expressed as where D 0 (ω) ∈ C N×N represents the complex deterministic part of the system and D j (ω) ∈ R N×N the random components. The total number of random variables, M, corresponds to the number of random variables associated with the Karhunen-Loève series expansion. For the systems under consideration, the expressions for D 0 and D j are as follows where M 0 , K 0 , C 0 are deterministic mass, stiffness and damping matrices, and M i and K i are symmetric matrices which contribute towards the random components of the mass and stiffness matrices.ũ(ω, θ) andf 0 (ω) are the dynamic response and the forcing in the frequency domain. For this study, ξ(θ) = {ξ i (θ)} for i = 1, . . . , M is considered to be a set of uncorrelated standard Gaussian random variables. The first p 1 random variables are used to model the random mass matrices, and the remaining random variables used to model the random stiffness matrices. The damping framework utilised in this study is a constant modal damping model [5]. Therefore, the matrix ζ is a diagonal matrix which contains modal damping factors, ζ = diag[ζ 1 , ζ 2 , . . . ζ N ] ∈ R N×N . It is assumed that the all the diagonal entries are equal, therefore ζ 1 = ζ 2 = · · · = ζ N . In the subsequent sections, different projection methods for efficiently approximating the response of Equation (5) are proposed and compared. The methods are subsequently applied for all θ ∈ Θ and for every frequency value ω ∈ Ω.

Derivation of the projection methods
We will initially consider two projection methods which have different characteristics. In order to compare the accuracy and effectiveness methods, a benchmark solution can be obtained by implementing a direct Monte Carlo approach [DMCS] for each frequency and realisation. The rationale behind proposing different methods is to analyse the effect of the nature of the coefficients and their associated vectors. The first two methods under consideration have the following characteristics: • Projecting onto a stochastic basis with stochastic coefficients [SS] • Projecting onto a deterministic basis with deterministic coefficients [DD] For both methods, we aim to keep the basis vector independent of the frequency. This is done in an attempt to reduce the computational effort if more than one frequency value were to be analysed. We initially consider the case which incorporates the whole stochastic nature of Equation (5).

Projecting onto a stochastic basis with stochastic coefficients [SS]
For such an approach, we initially consider the generalised undamped eigenvalue problem where λ k (θ) and φ k (θ) are the kth undamped random eigenvalue and eigenvector. For convenience, matrices that contain the whole set of random eigenvalues and eigenvectors are defined as follows The eigenvalues are arranged in ascending order so λ 1 (θ) < λ 2 (θ) < . . . < λ n (θ) and their corresponding eigenvectors are mass normalised and arranged in the same order. Therefore, MATEC Web of Conferences 211, 01003 (2018) https://doi.org/10.1051/matecconf/201821101003 VETOMAC XIV it is apparent that Φ T (θ)M(θ)Φ(θ) = I and Φ T (θ)K(θ)Φ(θ) = Ω 2 (θ). As the undamped eigenvectors from a complete basis, it is possible to obtain the response of Equation (5) through projecting on the undamped eigenvectors. By combining the stated identities with the stochastic finite representations introduced in the previous section, an expression for the response of a stochastically parametrised system can be expressed as where N corresponds to the number of degrees of freedom associated with the dynamic structure. Therefore, the response of the dynamic stochastic system under consideration has been represented in the required form.

Projecting onto a deterministic basis with deterministic coefficients [DD]
As a worst-case scenario, we consider the case of projecting deterministic coefficients onto a deterministic basis. The method is initiated by considering the deterministic equivalent of Equation (11), therefore, we deem all the eigensolution to be deterministic. For this case, the response vector takes the following form where γ 0 j (ω) ∈ C and c j ∈ C N are deterministic scalars and vectors respectively. If both the undamped random eigenvalues and eigenvectors seen in Equation (11) are exchanged for their deterministic counterpart, the response vector takes the following form where λ 0 j and φ 0 j denote the jth undamped deterministic eigenvalue and eigenvector respectively. Due to all the terms in Equation (13) being deterministic, the stochastic nature of the system is not at all incorporated into the response vector. It can be deduced that Equation (13) provides the deterministic solution and therefore can be established as a worst case scenario. However, if the coefficient of variation associated with the stochastic process is low, this method could provide an adequate approximation of the mean of the true solution.

Computational reduction
At present, the computational time associated with both the proposed methods can be considered quite high, especially for a high degree of freedom finite element system. This is due: • A large number of terms are present in the summations given by Equations (11) and (13).
At present, the number of terms in the series correspond to the number of degrees of freedom.
• Calculating the random eigensolutions is computationally expensive. This only applies to the first method.
Combining these with the need to simulate the methods for a large number of sample points result in a large computational effort. As a result, these issues have been addressed in the following sections.

Approximating the undamped eigensolutions
Calculating the exact undamped random eigensolutions can be extremely expensive, especially if the number of degrees of freedom is large. Thus a sensitivity approach to approximate MATEC Web of Conferences 211, 01003 (2018) https://doi.org/10.1051/matecconf/201821101003 VETOMAC XIV the eigensolutions could computationally be a better option. Due to the low computational effort associated with the first order perturbation method, this method has been used to approximate both the undamped random eigenvalues and eigenvectors. Therefore the random eigenvalues and eigenvectors can be approximated by the following expansions where λ 0 j and φ 0 j are the jth deterministic undamped eigenvalue and eigenvector respectively, and dξ k (θ) is set of random variables. Expressions for the derivatives of the undamped random eigenvalues and eigenvectors with respect to ξ k are given by The full algebraic details for obtaining the derivatives have been omitted, however for a comprehensive discussion about their derivations can be found in [6]. In the instance of Equation (15), the values of both ∂M ∂ξ k and ∂K ∂ξ k are where M k and K k − p 1 correspond to the random components of the mass and stiffness matrices.

Modal basis reduction
At present, all the methods described in Section 3 require the calculation and summation of N terms. However, we aim to lower the computational effort of the proposed methods by applying suitable truncations. By revisiting the ordering of the eigenvalues seen in Equation (10), it can be deduced that where λ j (θ) corresponds to the jth eigenvalue. From the scalar terms α j (ω, θ) and γ 0 j (ω) seen Equations (11) and (13), it can be observed that the eigenvalues appear in the denominator. The scalar α j (ω, θ) is shown for illustration For the values of j satisfying λ j (θ) + 2i � λ j (θ)ωζ > ω 2 , it is apparent that the value of the denominator increases as the value of j increases. Therefore, it is established that the value of α j (ω, θ) generally decreases as the value of j increases. Consequently the upper limits of the summations seen in Equations (11) and (13) can be lowered. In turn, these equations can be expressed asũ respectively, where n r < N. The value of n r can either be predefined or it can defined by applying an iterative framework.

Sample based Galerkin error minimisation
Two different projection methods have been proposed in Section 3. The first projects random scalars onto a stochastic basis whilst the third method projects deterministic scalars onto a deterministic basis. We have shown that it is possible to approximate the random eigensolutions that arise in the proposed methods in order to lower the computational effort. However these approximations, in addition to the modal reduction introduced in Section 4.2 introduces error into the calculation. This has motivated an error minimisation technique through applying a sample based Galerkin approach [7]. As a result, in addition to the methods introduced in Section 3, the following two projection methods are proposed: • A sample based Galerkin approach with projecting onto a stochastic basis with stochastic coefficients [SSG] • A sample based Galerkin approach with projecting onto a deterministic basis with deterministic coefficients [DDG] The case of incorporating a sample based Galerkin approach with projecting onto a stochastic basis with stochastic coefficients is initially considered.

A sample based Galerkin approach with projecting onto a stochastic basis with stochastic coefficients [SSG]
The response vector for the given case is modified to take the following series representationũ Here α j (ω, θ) and φ j (θ) correspond to the random scalars and random eigenvectors seen in Equation (11), whilst c j (ω, θ) ∈ C are constants which need to be obtained for each realisation. This can be done by applying a sample based Galerkin approach. We initially consider the following residual where ξ 0 = 1 is used in order to simplify the summation. D i (ω), ξ i (θ) andf 0 (ω) correspond to the terms arising in Equations (5). By making the residual orthogonal to a basis function, the unknown c j (ω, θ) can be computed. As Equation (21) can be viewed as a projection onto a stochastic basis, the residual is made orthogonal to the undamped random eigenvectors As a sample based Galerkin approach is considered, applying the orthogonality condition results in By defining the vector c 1 (ω, θ) = � c 1 (ω, θ) c 2 (ω, θ) . . . c n r (ω, θ) � T , Equation (25) can be rewritten as . . n r and y 1 (ω, θ) = φ T k (θ)f 0 (ω). The number of equations that need to be solved in order to calculate the unknown vector c(ω, θ) corresponds to the value of n r . By increasing the number of terms from n r to n r +1, the number of terms in Z 1 (ω, θ) increases by 2n +1. Therefore the lower the dimension of the reduced system, the fewer the number of equations that need to be solved. This is of importance as the given procedure needs to be repeated for every realisation and for every frequency under consideration.

A sample based Galerkin approach with projecting onto a deterministic basis with deterministic coefficients [DDG]
A Galerkin approach can also be considered for the case that contains undamped deterministic eigenvalues and eigenvectors. For this case, the response vector is defined as follows where γ 0 j (ω) and φ 0 j correspond to the deterministic scalars and the undamped deterministic eigenvectors introduced in Equation (13), and c j (ω, θ) ∈ C are unknown constants which need to obtained for each realisation of each frequency. By applying a sample based Galerkin approach, the following set of equations can be derived to obtain the unknown constants Z 2 (ω, θ)c 2 (ω, θ) = y 2 (ω, θ) j, k = 1, 2, . . . , n r where . . n r and c 2 (ω, θ) is the vector that contains the unknown constants c j (ω, θ) Similarly to the SSG method, Equation (28) needs to be solved for every realisation in each considered frequency. The computational effort associated with this method is considerably lower than the SSG method as the scalars γ 0 j only need to be calculated once for each given frequency. The aim of this method is to incorporate the whole stochastic nature of system within the unknown scalars c j (ω, θ). However it is not known if the Galerkin approach can substantially lower the error as all the eigensolutions used are deterministic. This method is of significant interest as it is known that the behaviour of deterministic and stochastic systems can differ substantially, especially if the coefficient of variation is significantly large.

Numerical example
The four methods proposed and discussed in Sections 3 and 5 are applied to analyse a classical Euler-Bernoulli cantilever beam. As the beam is clamped at one end, the displacement and rotational degrees of freedom at that particular end are zero. The beam is of length L, width W and thickness H and is discretised into 100 elements by applying the stochastic finite element method. This results in the cantilever beam having 200 degrees of freedom. The values of the parameters and the properties of the cantilever beam are given in Table 1  the cantilever beam, its bending rigidity, EI, is assumed to be a stochastic parameter. This stochastic parameter is assumed to be a stationary Gaussian random field of the following form where x corresponds to the coordinate along the length of the beam, a(x, θ) is a zero mean stationary Gaussian random field, and E 0 I 0 is the mean of the bending rigidity. The autocovariance kernel of the random field is given by where µ a corresponds to the correlation length. The values of the correlation length and input standard deviation are set to µ a = L 2 and σ a = 0.20 respectively, whilst four terms have been retained in the Karhunen-Loève expansion. The case of an unit amplitude harmonic point load acting on the free tip of the beam is considered for the frequency range 0 − 450 Hz at an interval of 2 Hz. This corresponds to considering 226 frequency values. For the given beam, this allows for the first eight resequence frequencies to be studied. 10, 000 Monte Carlo simulation samples are considered for each frequency step. It has been numerically verified that using 10, 000 Monte Carlo samples gives a satisfactory level of convergence for the first two moments of the quantities of interest. By taking the DMCS method to be a benchmark, the proposed methods, SS, DD, SSG and DDG have been compared through the use of relative error indicators. The approximate L 2 relative error of the mean of the response vector at each frequency step is defined as followŝ where µ DMCS denotes the mean of the response vector obtained by using the DMCS method, and µ CM denotes the mean of the response vector obtained by using a comparable method. This method ensures that the error arising from each of the projection methods can be characterised by a single value for each ω ∈ Ω. Figure 1 depicts the log of the approximate L 2 relative error of the mean of the response vector for different values of n r and for each frequency step. It can be easily deduced that the DD method introduces considerably more error than the other methods at the resonance frequencies. The visible troughs seen in the relative error arising from the SS, SSG and DDG methods correspond to the resonance frequencies. For a given value of n r the trend of the approximate relative error increases with the frequency. This is to be expected as the higher order terms in the summations become more important as the frequency increases. The relative errors induced by both the sample based Galerkin methods are almost identical. This suggests that computing the stochastic eigensolutions is non-essential if a sample based Galerkin method is used. The expression for the approximate L 2 relative error of the standard deviation takes a similar form where σ DMCS denotes the standard deviation of the response vector obtained by using the DMCS method, and σ CM denotes the standard deviation of the response vector obtained by using a comparable method. Figure 2 depicts the log of the approximate L 2 relative error of the standard deviation of the response vector for different values of n r at each frequency step. It is apparent that the DD method does not capture the standard deviation of the DMCS very well, however this was to be expected. The SS method does not capture the standard deviation of the DMCS very well either, however, all the sample based Galerkin methods capture the necessary standard deviation very well, especially when the frequency step corresponds to a resonance value. Similarly to the case of the approximate L 2 relative error of the mean of the response vector, the approximate L 2 relative error of the standard deviation of the response vector are extremely similar for both the Galerkin methods.

Conclusion
In this study, a set of hybrid projection methods has been proposed to analyse the effect of altering the nature of the coefficients and vectors which are associated with projection MATEC Web of Conferences 211, 01003 (2018) https://doi.org/10.1051/matecconf/201821101003 VETOMAC XIV methods. Following the implementation of a stochastic finite element method, two projection methods have been developed by utilising the random eigenvalue problem. The first method utilises both random eigenvalues and eigenvectors, and the second uses deterministic eigensolutions. In order to reduce the computational effort associated with the methods, the random eigensolutions have been approximated by a first order perturbation and only the dominant projection terms have been retained. Due to the approximations and the reduced modal basis, two additional projection methods have been proposed. These methods utilise a sample based Galerkin error minimisation approach in an attempt to lower the error. In the example provided, it is apparent that if the Galerkin error minimisation approach is applied, calculating the stochastic eigensolutions is unwarranted. The Galerkin method compensates for the error induced while utilising the baseline eigenmodes. Therefore, a Monte Carlo type sample-based method to evaluate the coefficients and their associate basis can be avoided. The application of the Galerkin error minimisation approach in conjunction with projecting onto a deterministic basis with deterministic coefficients [DDG] produces a level of accuracy comparable to any of the other proposed methods. Our study leads us to suggest that this simple approach has significant potential for analysing stochastic structural systems.