Hybrid Stochastic Finite Element Method for Mechanical Vibration Problems

We present and analyze a new hybrid stochastic finite element method for solving eigenmodes of structures with random geometry and random elastic modulus. The fundamental assumption is that the smallest eigenpair is well defined over the whole stochastic parameter space. The geometric uncertainty is resolved using collocation and random material models using Galerkin method at each collocation point. The response statistics, expectation and variance of the smallest eigenmode, are computed in numerical experiments. The hybrid approach is superior to alternatives in practical cases where the number of random parameters used to describe geometric uncertainty is much smaller than that of the material models.


Introduction
In standard engineering models many physical quantities such as material parameters are taken to be constant, even though their statistical nature is well understood.Similarly, assumptions of geometric constants, such as thickness of a structure, are not realistic due to manufacturing imperfections.In a detailed report on a state-of-the-art verification and validation process comparing modern simulations with the set of experiments performed in the Oak Ridge National Laboratory in the early 70s, Szabo and Muntges report discrepancies of over 20% in some quantities of interest [1].These discrepancies are attributed to machining imperfections not accounted for in the computations.Also, in important nonlinear problems such as buckling of a shell, it is known that variation between manufactured specimens has a profound effect in the actual performance [2].This suggests that a stochastic dimension should be added to the models.
The modern era of uncertainty quantification starts with the works of Babuska et al. [3,4] and the ETH-group led by Schwab and Gittelson (e.g., [5]) with provably faster convergence rates than the standard Monte Carlo methods.Although the so-called stochastic finite elements had been studied for a relative long time before (the classic reference is Ghanem and Spanos [6]), their application was thus limited to highly specific cases.
The solution methods can broadly speaking be divided into intrusive and nonintrusive ones.The same division applies to stochastic eigenvalue problems (SEVPs) as well.SEVPs have attracted a lot of attention recently and various algorithms have been suggested for computing approximate eigenpairs [7][8][9][10], specially the power iteration [11] and inverse iteration [12].
In this paper our focus is on effects of material models and manufacturing imperfections of geometric nature.It should be noted that in the context of this paper it is assumed that the problems are positive definite and the eigenpair of interest is the ground state, that is, the one with the smallest eigenvalue which, in theory, can be a double eigenvalue.Our experimental setup is nonsymmetric and, thus, a spectral gap exists and the inverse iteration converges to the desired eigenmode.
In stochastic eigenvalue problems one must address two central issues that do not arise in stochastic source problems: first, the eigenmodes are defined only up to a sign and, second, the eigenmodes must be normalized over the whole parameter space; that is, every realization must be normalized in the same way.Here our quantity of interest is the eigenvalue and, therefore, we do not necessarily have to fix the signs.The normalization is handled by solution of a nonlinear system of equations as in [12].The main result of this paper is the new hybrid algorithm which combines a nonintrusive outer loop (collocation) with an intrusive inner one (Galerkin).The randomness in geometry is resolved using collocation and that of materials with Galerkin at every collocation point.The model problem is an idealized engine bed vibration problem, where a 2D hollow square is clamped from one side (see Figure 1).The inner cavity is assumed to be of fixed shape but random orientation and Young's modulus is taken to be random in the normal direction into the body.The choice of the geometric uncertainty models the case of sand casting with moulds of fixed shape.The hybrid approach combines the geometric flexibility of collocation with much faster computation times of the Galerkin approach used in the inner loop and, thus, combines the best of both worlds.
The rest of the paper is organized as follows: first the abstract problem is introduced in Section 2. The representation of the random input is given in Section 3. Higher order finite elements are described in Section 4. Section 5 is central to our discussion; both solution methods, collocation and spectral inverse iteration, are presented.Numerical experiments are discussed in Section 6 and finally conclusions are drawn in Section 7.

Problem Statement
We start by presenting the deterministic system of Navier's equations of elasticity and the corresponding weak form.We then extend this to the stochastic setting in order to cover the case of uncertain domain and modulus of elasticity.We aim to compute statistics of the smallest eigenvalue of the resulting stochastic system.The geometry of the computational domain is illustrated in Figure 1.

Deterministic Formulation.
We use the Navier equations of elasticity to model the system of interest.Find the eigenvalue , the displacement field u = ( 1 ,  2 ), and the symmetric stress tensor  = (  ) 2  ,=1 , such that  = div (u) where  =   ∪   is a partitioned boundary of .The Lamé constants are with  and ] being Young's modulus and Poisson's ratio, respectively.Further, I is the identity tensor, n denotes the outward unit normal to   , and the strain tensor is The vector-valued tensor divergence is This formulation assumes a constitutive relation corresponding to linear isotropic elasticity with stresses and strains related by Hooke's generalized law: ] ] ] = D (, )  V . (5)

Weak Formulation.
Let us introduce a function space   = {V : V ∈  1 (), V|   = 0} and define the eigenproblem: find the eigenpair (, u) ∈ R × (  ×   ) such that where the bilinear form is the integrated tensor contraction and (⋅, ⋅)  denotes the mass matrix As usual, we define the kinematic relation and specify the constitutive matrix for the purpose of rewriting the bilinear form as 2.3.Stochastic Extension of the Eigenproblem.We introduce two probability spaces (Ω 1 , F 1 , P 1 ) and (Ω 2 , F 2 , P 2 ) to model randomness of the domain and the elastic modulus, respectively.Here Ω  denotes the set of outcomes, F  is a algebra of events, and P  is a probability measure on Ω  for each  = 1, 2. Furthermore we set to be the underlying product probability space.The geometry of the computational domain  is illustrated in Figure 1.We assume that the radius of the cavity  : Ω 1 → R and the angle  : Ω 1 → R are random variables satisfying  min ≤ ( 1 ) ≤  max and  min ≤ ( 1 ) ≤  max for some suitable constants  min ,  max ,  min ,  max > 0. Thus we may write  = ( 1 ) = (( 1 ), ( 1 )).In the numerical examples we set ( 1 ) and ( 1 ) to be uniformly distributed random variables.
In most cases we are mainly interested in computing statistical moments of the solution.Suppose that we have a random variable V( 1 ,  2 ).The expected value or mean of and its variance is In this paper our aim is to compute the expected value and variance of the smallest eigenvalue of system (15).

Spectral Representations
For purposes of numerical treatment we need to approximate the random Young modulus using only a finite number of random variables.A natural way of achieving this is to use the truncated Karhunen-Loève expansion.Once a finite-dimensional approximation has been obtained we may express our solution in a generalized polynomial chaos basis and apply solution schemes based on the stochastic Galerkin finite element method; see [6] for a thorough introduction to this approach.Instead of using the Karhunen-Loève expansion one could also fix a spatial basis and then estimate the polynomial chaos coefficients straight from measurement data; see, for instance, [13,14].

Karhunen-Loève Expansion.
The Karhunen-Loève expansion is a representation of a random field as a linear combination of the eigenfunctions of the associated covariance operator.Truncating the resulting series we obtain a finitedimensional approximation of the original random field.The Karhunen-Loève expansion is the optimal choice among linear expansions in the sense that it minimizes the mean square truncation error [6].
For a fixed  1 ∈ Ω 1 let ( 2 , x) = ( 1 ,  2 , x) be a random field on  = ( 1 ) with mean field and a covariance function that is symmetric and positive definite.The associated covariance operator is compact as an operator C  :  2 () →  2 () and therefore admits a countable number of positive eigenvalues {  } ∞ =1 and associated eigenfunctions {  } ∞ =1 .The eigenvalues {  } ∞ =1 tend to zero as  tends to infinity and the eigenfunctions {  } ∞ =1 form an orthonormal basis of  2 ().The eigenpairs may be computed numerically using, for example, fast multipole methods [15].
The Karhunen-Loève expansion of the random field ( 2 , x) is now given by where {  } ∞ =1 are centered and orthogonal but not necessarily independent random variables.
In numerical computations we replace the random field ( 2 , x) with an approximation obtained by truncating series (22) after  terms.The truncation error depends upon the decay of the Karhunen-Loève eigenvalues and eigenfunctions.See [16] for bounds on this decay.
Remark 1.It is essential for numerical algorithms that the random variables {  } ∞ =1 are mutually independent, even though this is not in general implied by their orthogonality.A solution to this issue has been proposed in [3], where suitable auxiliary density functions have been introduced.In the numerical examples we assume a representation of the random field with respect to a set of independent uniformly distributed random variables.Other types of random variables (e.g., Gaussian) could be considered as well.

Legendre Chaos.
We employ the generalized polynomial chaos framework which essentially means representing our solution on a basis of orthogonal polynomials.In our case we assume the input random variables to be uniformly distributed and the polynomial chaos basis is therefore given by tensorized Legendre polynomials.The use of orthogonal polynomials as a basis allows us to apply stochastic Galerkin based methods and ensures optimal convergence of these methods; see [17].In this paper we aim to solve the eigenmodes using a method of spectral inverse iteration introduced in [12].
Assume that y( 2 ) = ( 1 ( 2 ), . . .,   ( 2 )) is a vector of mutually independent random variables that are uniformly distributed in the interval [−1, 1].The expected value or mean of a random variable V(y) is now given by Here we have used a subscript to indicate that the expected value is taken over the sample space Ω 2 only.We associate each multi-index  = ( 1 , . . .,   ) with the multivariate Legendre polynomial where   () denotes the univariate Legendre polynomial of degree .We assume that the polynomials are normalized so that . This allows us to express any square integrable random variable V(y) in a polynomial chaos expansion where the coefficients are given by For numerical computations we truncate expansion (25).To this end we choose a finite set of multi-indices A ⊂ N  0 and approximate V(y) with the series Here the choice of the set A is key to obtaining accurate approximations and is currently a topic of active research.In our numerical examples we consider sets that for a fixed level  ∈ R are given by A , = { ∈ N  0 | ∑  =1   /√  ≤ }, where {  } ∞ =1 are the Karhunen-Loève eigenvalues.More sophisticated adaptive schemes for selecting the multi-index sets have been presented in [18][19][20].
The most important reason for using the polynomial chaos basis is that it allows fast computation of any expectations involved.This is because we may evaluate integrals over the stochastic domain analytically and are thus able to avoid numerical integration in high dimensions.Orthogonality of the Legendre polynomials yields for all multi-indices ,  ∈ N  0 .For notational convenience we set These coefficients are easy to evaluate analytically; see the appendix in [12].It is also worth noting that most of the coefficients are evaluated to zero.This fact can be exploited in order to speed up the solution process.

Higher Order Finite Elements
Here we give a short overview of the ℎ-FEM following closely the one in [21].In the -version of the FEM the polynomial degree  is used to control the accuracy of the solution while keeping the mesh fixed in contrast to the ℎ-version or standard finite element method, where the polynomial degree is constant but the mesh size varies.The ℎ-method simply combines the ℎand -refinements.
In our setting we can use topologically fixed meshes with high-order elements with  = 4 and  = 8 to ensure sufficient resolution of the deterministic part of the solution.
In the following one way to construct a -type quadrilateral element is given.The construction of triangles follows similar lines.First of all, the choice of shape functions is not unique.We use the so-called hierarchic integrated Legendre shape functions.
Legendre polynomials of degree  can be defined by a recursion formula where  0 () = 1 and  1 () = .
The derivatives can similarly be computed by using the recursion The integrated Legendre polynomials are defined for  ∈ [−1, 1] as and can be rewritten as linear combinations of Legendre polynomials The normalizing coefficients are chosen so that Using these polynomials we can now define the shape functions for a quadrilateral reference element over the domain [−1, 1] × [−1, 1].The shape functions are divided into three categories: nodal shape functions, side modes, and internal modes.
Note that some additional book-keeping is necessary.The Legendre polynomials have the property   (−) = (−1)    ().This means that every edge must be globally parameterized in the same way in both elements where it belongs.

Curved Boundaries.
Since we want to use fixed mesh topologies even with perturbed domains, it is important to represent curved boundary segments accurately.The linear blending function method of Gordon and Hall [22] is our choice for this purpose.
In the general case all sides of an element can be curved, but in our case only one side is-as in Figure 1.We assume that every side is parameterized: where  = 1, . . ., number of sides.Using capital letters as coordinates of the corner points (  ,   ), we can write the mapping for the global -coordinates of a quadrilateral as and symmetrically for the -coordinate.Note that if the side parameterizations represent straight edges, the mapping simplifies to the standard bilinear mapping of quadrilaterals.

The Solution Method
We will present a hybrid method of stochastic finite elements for solving the stochastic eigenvalue problem (15).More precisely, we apply stochastic collocation to resolve the dependency on the geometry and the stochastic Galerkin scheme for computing the effects of the random field.The method of stochastic collocation allows for an easy implementation, since the solution statistics are sampled from an ensemble of deterministic solutions.This makes the method attractive from the point of view of random geometry.On the other hand, if a sufficiently accurate parametric representation for the input random field exists, then the use of stochastic Galerkin based methods is well motivated by their high rate of convergence-especially if the number of parameters is large.
In this paper we apply a method of spectral inverse iteration introduced in [12] to compute the eigenpair of interest at each collocation point.

Discretization.
Let us consider the stochastic eigenproblem (15) for a fixed  1 ∈ Ω 1 .We replace the Young modulus  = ( 1 , ⋅, ⋅) with the expression obtained by truncating the Karhunen-Loève expansion (22) after  terms.Plugging (39) into (2) we may write the constitutive matrix (11) in the form We assume that {  }  =1 is a set of mutually independent and uniformly distributed random variables with supports scaled to the interval [−1, 1].The sample space Ω 2 is thus parameterized by the vector y = ( 1 , . . .,   ) ⊂ [−1, 1]  .With  1 ∈ Ω 1 being fixed we may now interpret the eigenvalue  as a function of y and the eigenfunction u as a function of y and x.
Next we address the spatial discretization of ( 15) on the fixed domain  = ( 1 ).Let {  }  =1 be a set of global basis functions for the approximation space of   ×   .Here we employ the -version of the finite element method so that the components of each   are some shape functions (34), (35), or (36) from Section 4. From (15) we obtain the deterministic Galerkin projection We define the stiffness matrix A () for each  = 0, . . .,  by setting and the mass matrix M via For a fixed  1 ∈ Ω 1 we have now reduced (15) to a stochastic matrix eigenvalue problem: find a random variable ( 1 , ⋅) :

Stochastic Collocation.
In our problem of interest the geometry of the computational domain depends on the random variables ( 1 ) and ( 1 ).Thus, these variables give us a parametrization of the sample space Ω 1 .We present an anisotropic collocation operator and apply it for computing the dependency of our solution on the parameters.Implementation of a sparse grid collocation operator (see [4,23]) would also be possible but since we only consider collocation in two dimensions we are restricted to full tensor collocation in the following.In our case we rely on the Gauss-Legendre quadrature, since the random variables ( 1 ) and ( 1 ) are uniformly distributed.Furthermore, for notational simplicity we assume these to be scaled so that ,  ∈ [−1, 1].Denote by { ()  }  =0 the abscissae of the Gauss-Legendre quadrature of order  and by { ()   }  =0 the associated quadrature weights.The quadrature points { ()   }  =0 are the zeros of the univariate Legendre polynomial of degree  + 1.The one-dimensional Lagrange interpolation operators I  with respect to the interpolation points { ()   }  =0 are defined via Here {ℓ ()  }  =0 are the Lagrange basis polynomials of degree .For information on orthogonal polynomials and computation of the quadrature weights see [24].
The full tensor Lagrange interpolation operator is defined as the tensor product of the univariate operators.In our twodimensional case this is where I ()  1 and I ()  2 denote the one-dimensional interpolation operators in the variables ( 1 ) and ( 1 ), respectively, and  = { 1 ,  2 } are the orders of quadrature.Applied on a random variable V(, ) the full tensor interpolation gives V ( The expected value of a random variable V(, ) is given by Here we have again used a subscript to designate that the expectation is taken over the sample space Ω 1 .From (47) we obtain V ( Hence we may compute the expected value of the solution once we have solved the underlying problem at each collocation point.

Spectral Inverse Iteration.
We choose to solve (44) at each collocation point using a method of spectral inverse iteration.This method has been introduced in [12] and there its convergence towards the eigenpair with the smallest eigenvalue has been verified in the case of an elliptic operator with random coefficients.However, assuming that the eigenpair of interest is of multiplicity one, the method is in fact applicable for a much wider range of problems.We may consider the spectral inverse iteration as a stochastic extension of the deterministic inverse iteration which has been given in Algorithm 2 for reference.For the generalized eigenvalue problem of a stiffness matrix A and a mass matrix M, the deterministic inverse iteration converges to the eigenpair for which the eigenvalue is closest to the given parameter κ.
The idea in the spectral inverse iteration is to interpret each equation involved in Algorithm 2 in a stochastic sense using Galerkin projections with respect to a given spectral basis.To this end we start by fixing a finite polynomial chaos basis {Λ  (y) |  ∈ A ⊂ N  0 } and proceed by formulating the respective Galerkin projections.Here we only recap the essentials of the process as it has already been thoroughly explained in [12].
For some κ ∈ R let us consider the equation that arises from the first step of Algorithm 2.Here we assume the matrix A to take the form as in (44).We replace the random vectors u(y) and k(y) with the truncated polynomial chaos expansions The Galerkin projection of (50) on the basis {Λ  } ∈A is now defined as Using the orthogonality relations (27) this reduces to the linear system where , , and Ψ are block vectors and matrices characterized by Let us next consider the normalization step in Algorithm 2. For a vector k(y) we aim to compute w(y) such that If we set (y) = ‖k(y)‖ then the respective Galerkin projection on the basis {Λ  } ∈A is defined as Shock and Vibration and reduces to This is a nonlinear system of equations for the coefficients {  } ∈A and can be solved using, for example, Newton's method with the initial guess   = ‖k  (y)‖.The Jacobian is given explicitly by The Galerkin projection of ( 56) is now defined as and reduces to the block linear system where For computing the eigenvalue we need to be able to evaluate the Rayleigh quotient for a vector u(y).As explained in [12] we may approximate (y) in the Galerkin sense to obtain the formula for the Galerkin coefficients {  } ∈A .Algorithm 3 now formulates the method of spectral inverse iteration.Similar to the deterministic case, the choice of the parameter κ affects the speed of convergence.However, if κ is chosen too close to the actual eigenvalue, the iteration might not converge as (50) becomes singular for some y ∈ [−1, 1]  .We refer to [12] for more information.In our numerical examples we set κ = 0 which ensures convergence to the eigenpair with the smallest eigenvalue as long as the matrix A(y) is positive definite for all y ∈ [−1, 1]  .Algorithm 3 (spectral inverse iteration for stochastic eigenvalue problems).Fix κ ∈ R and tol > 0. Let u (0) = {u (0)  } ∈A be a normalized initial guess for the eigenvector and set  (0) = {  (u (0) )} ∈A .For  = 1, 2, . . .do the following: (1) With u = {u (−1)  } ∈A solve system (54) for k = {k  } ∈A .

Response Statistics.
We may easily calculate statistical moments of a random variable V(, , y) once we know its polynomial chaos expansion at each collocation point.By using the orthogonality relations (27) we obtain The first and second moments of V(, , y) are now given by These can be computed by applying formula (49) on (66).From these we obtain approximations for the expected value and variance of V(, , y).

Numerical Experiments
Let us consider a problem with uncertain domain and uncertain field.The computational domain is In Figure 2 the effect of domain perturbation on the modes is illustrated.Notice that if one wanted to compute the statistics of the modes, it would be necessary to map every realisation onto the nominal domain via conformal mappings, for instance.
The (constant) material parameters are Poisson ratio ] = 1/3 and the mean field of the Young modulus  = 2.069 ⋅ 10 11 MPa.We have used the convention, however, where the reported results are normalized by the value of .
We will proceed in two steps by first letting only the geometry vary before adding uncertainty to the field as well.The cases are (A) uncertain domain with deterministic field, solved with collocation, (B) uncertain domain with uncertain field, solved with the new hybrid method.
In both cases the nominal domain is the same (see Figure 1).In the absence of exact solutions overkill solutions have been used in convergence graphs unless otherwise specified.We have also verified the results by comparing to a full Monte Carlo simulation of 620000 draws.

Case A.
The values of the smallest eigenvalue over the parameter space are shown in Figure 3 using both contour and surface plots.We consider a collocation sequence, where the collocation points are taken to span all 16 cases from 1 × 1 to 4 × 4 points over the  × -parameter space.The relative convergence of the smallest eigenvalue over the whole collocation sequence is shown in Figure 4. We have used constant polynomial order  = 8 (1728 degrees of freedom) throughout and  = 16 (6528 degrees of freedom) for the reference computation.
The dashed line shows the convergence in the case when the deterministic discretization of the reference solution is used also in collocation.The exponential convergence in expectation is due to convergence in parameter space only.The poorly performing grids in Figure 4 are the anisotropic ones with more points for .This is due to relative effects of  and  not being exactly balanced (see Figure 3).The solid line indicates that already at rule (2, 2) the spatial discretization used is not sufficient when the reference has been computed with much higher resolution.The convergence stalls completely since the deterministic error dominates.Convergence in variance is slower, as suggested by theory.

Case B.
We proceed to consider a fixed collocation grid with 4 × 4 collocation points for the parameters  and .Assuming normalization of the mean field ( = 1) we set the coefficients in the Karhunen-Loève expansion of the Young modulus (39) to be   (x) = cos((x)), where (x) is the distance from x to the inner boundary and   = ( + 1) −3 .Thus the series can be considered to converge at an algebraic

Figure 1 :
Figure 1: Two realisations of the domain with meshes.The nominal domain is the reference.The rotation and scaling of the cavity as well as the direction of the uncertainty in the Young's modulus is indicated.

Figure 2 :
Figure 2: Case A: effect of perturbation on eigenmodes.Contour lines of displacement.