Cylindrical Shell with Junctions: Uncertainty Quantification of Free Vibration and Frequency Response Analysis

Powered by TCPDF (www.tcpdf.org) This material is protected by copyright and other intellectual property rights, and duplication or sale of all or part of any of the repository collections is not permitted, except that material may be duplicated by you for your research use or educational purposes in electronic or print form. You must obtain permission for any other use. Electronic or print copies may not be offered, whether for sale or otherwise to anyone who is not an authorised user. Laaksonen, Mikael; Hakula, Harri; Kaarnioja, Vesa


Introduction
Numerical simulation of thin solids remains one of the challenges in computational mechanics. e advent of stochastic finite element methods has led to new possibilities in simulations, where it is possible to replace isotropy assumptions with statistical models of material parameters or incorporate manufacturing imperfections with parametrized computational domains and derive statistical quantities of interest such as expectation and variance of the solution. Here, we discuss two special classes of such uncertainty quantification problems: free vibration of shells of revolution and directly related frequency response analysis. What makes the specific shell configurations interesting from the application point of view is that they include junctions and shell-solid couplings. Free vibrations of cylinder-cone configurations have been studied by many authors before [1][2][3][4]. ere also exists an exhaustive literature review on shells with junctions [5].
In vibration problems, the sources of uncertainty relate to materials and geometry. In this paper, the stochastic vibration problems are solved using a nonintrusive approach, the stochastic collocation [6]. Collocation methods are particularly appealing in problems with parameterized (random) domains. e standard approach is to solve the problem at specified quadrature points and interpolate over the points. If the number of random variables is high, this leads to the curse of dimensionality which can be alleviated up to a point with special high-dimensional quadratures, the so-called sparse grids [7][8][9][10]. If the uncertainty is only in the material parameters, one can apply the intrusive approach such as stochastic Galerkin methods, see [11] and references therein.
Sparse grids are designed to satisfy given requirements for quadrature rules. It is, however, possible to construct similar interpolation operators even if the point set is limited to some subset of the parameter space (partial sparse grids), or, crucially, points are added to the sparse grid. In both cases, the construction is the same [12]. Every point in a sparse grid corresponds to some realization of the random field. If some real measurements become available and the data can be identified as points in the parameter space, they can be incorporated into the interpolation operator.
e two model problems studied here are a wind turbine tower resonance problem and the free vibration of a classical long cylinder-cone junction combined with T-junction via kinematic constraints. e wind turbine tower problem is inspired by an earlier study [13], where the tower is modeled as a solid due to solver limitations. e long cylinder case is in turn inspired by a photograph of a collapsed pipe ( [14], p. 82).
We have studied related stochastic shell eigenproblems in the previous work [15], but only for cylinders with constant thickness. In this study, in all cases one, of the random parameters is geometric which changes the character of the computational problems considerably. However, the underlying theoretical properties remain the same. One of the characteristic features of multiparametric eigenvalue problems is that the order of eigenpairs in the spectrum may vary over the parameter space, a phenomenon referred to as crossing of the modes.
In this paper, we show how the existing methods can be applied to shell problems that have characteristics that differ from those covered in the literature until now. We do observe crossing of the modes in the wind turbine tower problem which is satisfying since it highlights the fact that this feature is not an artificial mathematical concept but something that is present in standard engineering problems. e stochastic framework for eigenproblems is applied to frequency response analysis in a straightforward manner, yet the resulting method is new. We also demonstrate in connection with frequency response analysis that interpolation operators based on partial sparse grids can provide useful information.
All simulations have been computed using p-version of the finite element method [16]. e results derived using collocation on sparse grids have been verified using the Monte Carlo method. e rest of this paper is organized as follows. In Section 2, the two shell models used in this paper are introduced; in Section 3, the stochastic variant is given with the solution algorithms; Section 3.5 covers the multivariate interpolation construction; the frequency response analysis is briefly outlined in Section 4; the numerical experiments are covered in detail in Section 5; finally, conclusions are drawn in Section 6.

Shell Eigenproblem
Assuming a time harmonic displacement field, the free vibration problem for a general shell leads to the following abstract eigenvalue problem: find u ∈ R 3 and ω 2 ∈ R such that Su � ω 2 Mu, where u � u, v, w { } represents the shell displacement field, while ω 2 represents the square of the eigenfrequency. In the abstract setting, S and M are differential operators representing deformation energy and inertia, respectively. In the discrete setting, they refer to corresponding stiffness and mass matrices.
Simulation of thin solids using standard finite elements is difficult, since one (small thickness) dimension dominates the discretization. In the following, we shall derive two variants of problem (1) by employing different dimension reductions for shells of revolution. Instead of working in 3D, the elasticity equations are dimensionally reduced either in thickness or in the case of axisymmetric domains, using a suitable ansatz.

Shell Geometry.
In the following, we let a profile function y � f(x) revolve about the x-axis. Naturally, the profile function defines the local radius of the shell. Shells of revolution can formally be characterized as domains in R 3 of type where D is a (mid)surface of revolution, n(x) is the unit normal to D, and d(x) � d 1 (x) + d 0 (x) is the thickness of the shell measured in the direction of the normal. An illustration of such a configuration is given in Figure 1.

Constant (Dimensionless)
ickness. Let us have d(x) � d (constant) and define principal curvature coordinates, where only four parameters, the radii of principal curvature R 1 , R 2 , and the so-called Lamé parameters, A 1 , A 2 , which relate coordinates changes to arc lengths, are needed to specify the curvature and the metric on D.
e displacement vector field of the midsurface u � u, v, w { } can be interpreted as projections to directions where Ψ(x 1 , x 2 ) is a suitable parametrization of the surface of revolution, e 1 , e 2 are the unit tangent vectors along the principal curvature lines, and e 3 is the unit normal. In other words, Assuming that the shell profile is given by f(x 1 ), the geometry parameters are Both cylinders and cones are examples of parabolic shells according to the theory of surfaces by Gauss. is follows from the fact that, in both cases, f ″ (x 1 ) � 0. In particular the reciprocal function 1 For a cylinder with a constant radius, f(x 1 ) � 1, say, we get In the following, we shall keep the general form of the parameters in the equations.

Shock and Vibration
We can now define the dimensionless thickness t as t � d/R, where R ∼ R 2 is the unit length. In the following, we assume that R ∼ 1, and thus, use the two thickness concepts denoted by d and t interchangeably.

Reduction in ickness.
e free vibration problem for a general shell (1) in the case of shell of revolution with constant thickness t leads to the following eigenvalue problem: find u(t) and ω 2 (t) ∈ R such that where u(t) represents the shell displacement field, while ω 2 (t) represents the square of the eigenfrequency. e differential operators A m , A s , and A b account for membrane, shear, and bending potential energies, respectively, and are independent of t. Finally, M(t) is the inertia operator, which in this case can be split into the sum M(t) � tM l + t 3 M r , with M l (displacements) and M r (rotations) independent of t. Many well-known shell models fall into this framework.
Let us next consider the variational formulation of problem (8). Accordingly, we introduce the space V of admissible displacements and consider the problem: find where a m (·, ·), a s (·, ·), a b (·, ·), and m(t; ·, ·) are the bilinear forms associated with the operators A m , A s , A b , and M(t), respectively. Obviously, the space V and the three bilinear forms depend on the chosen shell model [17].

Two-Dimensional Model.
Our two-dimensional shell model is the so-called Reissner-Naghdi model [17], where the transverse deflections are approximated with low-order polynomials. e resulting vector field has five components u � (u, v, w, θ, ψ), where the first three are the standard displacements and the latter two are the rotations in the axial and angular directions, respectively. Here, we adopt the convention that the computational domain Ω t ⊂ R 2 is given by the surface parametrization, and the axial/angular coordinates are denoted by x and y: Deformation energy A(u, u) is divided into bending, membrane, and shear energies, denoted by subscripts B, M, and S, respectively.
Notice that we can safely cancel out one power of t. Bending, membrane, and shear energies are given as follows [18]: where v is the Poisson ratio (constant) and E(x, y) is Young's modulus with scaling 1/(12(1 − ] 2 )).
e symmetric 2D strains, bending (κ), membrane (β), and shear (ρ), are defined as In our experiments, the shell profiles are functions of x only, but not necessarily constant; hence, the strains are presented in a general form. e stiffness matrix S is obtained after integration and assembly. Here, the mass matrix M is the standard 2D one with density σ except for the scaling of the rotation components and surface differential A 1 (x, y)A 2 (x, y).

Reduction through Ansatz.
e 3D elasticity equations for shells of revolution can be reduced to two-dimensional ones using a suitable ansatz [19]. For shells of revolution, the eigenmodes u(x, y, z) or u(x, r, α) in cylindrical coordinates have either one of the forms: where K is the harmonic wavenumber. us, given a profile where the auxiliary functions d i (x) again simply indicate that the thickness d( Deformation energy A(u, u) is defined as where the symmetric 3D strains ϵ ij are e stiffness matrix S is obtained after integration and assembly. Here, the mass matrix M is the standard 3D one with density σ.

Boundary Layers.
e solutions of shell problems, static or dynamic, often include boundary layers and even internal layers, each of which has its own characteristic length scale [20]. Internal layers are generated by geometric features, such as cuts and holes, or changes in curvature, for instance, at shell junctions. e layers either decay exponentially toward the boundary or oscillate with exponentially decaying amplitude. For parabolic shells, it is known that the axial boundary layers decay exponentially with a characteristic length scale ∼ � t √ . If there are generators for internal layers, then the layers oscillate in the angular direction with characteristic length scale ∼ � t 4 √ . ere also exists an axial internal layer of a very long range ∼ 1/ � t √ . Identification of the boundary layer structure or layer resolution is a useful tool for assessing the quality of the discretized solution and providing guidance for setting up the discretization, for instance, in mesh generation. If the solution is not dominated by the layers, we simply observe the underlying smooth part of the solution.

Numerical Locking.
One of the numerical difficulties associated with thin structures is the so-called numerical locking, which means unavoidable loss of optimal convergence rate due to parameter-dependent error amplification [21]. We use the p-version of the finite element method [16] in order to alleviate the (possible) locking and ensure convergence.

Stochastic Shell Eigenproblem
e stochastic extension of the free vibration problem arises from introducing uncertainties in the material parameters and the geometry of the shell. In the context of this paper, we assume that these uncertainties may be parametrized using a finite number of uniformly distributed random variables. e approximate solution statistics are then computed employing a sparse grid stochastic collocation operator.

e Stochastic Eigenproblem.
We let the material and geometric uncertainties be represented by a vector of mutually independent random variables ξ � (ξ 1 , For uniformly distributed random variables, we may, without loss of generality, assume a scaling such that First, we assume a set of geometric parameters (ξ 1 , ξ 2 , . . . , ξ m′ ) and let the computational domain be given by a function D(ξ) � D(ξ 1 , ξ 2 , . . . , ξ m′ ). Second, we assume a set of material parameters (ξ m′+1 , ξ m′+2 , . . . , ξ M ), and let Young's modulus be a random field expressed in the form e parametrization (19) may, for instance, result from a Karhunen-Loéve expansion of the underlying random field E. We assume that the random field is strictly uniformly positive and uniformly bounded, i.e., there exists E min , E max > 0 such that for all ξ ∈ Γ. Let L 2 μ (Γ) denote a weighted L 2 -space, where μ is the uniform product probability measure associated with the probability distribution of ξ ∈ Γ. For functions in L 2 μ (Γ), we define the expected value Assuming stochastic models for the computational domain and Young's modulus, the eigenpairs of the free vibration problem now depend on ξ ∈ Γ.
e stochastic eigenvalue problem is obtained simply by replacing the differential operators in (8) with their stochastic counterparts and assuming that the resulting equation holds for every realization of ξ ∈ Γ. In the variational form, the problem reads find functions u(t) : where a m (ξ; ·, ·), a s (ξ; ·, ·), a b (ξ; ·, ·), and m(t, ξ; ·, ·) are stochastic equivalents of the deterministic bilinear forms in (9). We assume that the eigenvector u(t, ξ) is normalized with respect to the inner product m(t, ξ; ·, ·) for all ξ ∈ Γ.

Eigenvalue Crossings.
In the context of stochastic eigenvalue problems, special care must be taken in order to make sure that different realizations of the solution are in fact comparable. is is true for shell eigenvalue problems in particular, since the eigenvalues of thin cylindrical shells are typically tightly clustered. Double eigenvalues may appear due to symmetries and physical properties of the shell might be such that eigenvalues corresponding to different harmonic wavenumbers appear very close to each other. When the problems are brought to stochastic setting, we may face the issue of eigenvalue crossings. e ordering of the eigenvalues associated with each mode may be different for different realizations of the problem.
is issue has previously been illustrated in [11] and considered in the case of shell eigenvalue problems in [15].
In our examples, the eigenvalue crossings do not pose computational difficulties. When the problem is reduced through ansatz (14) or (15), the eigenmodes are separated by wavenumber K, and as a result for any fixed K, the eigenvalues of the reduced problem appear well-separated. However, one should bear in mind that the eigenvalues of the original 3D problem may still cross: In Section 5, we present numerical examples which illustrate that a kth smallest eigenvalue may be associated with eigenmodes with different wavenumbers for different values of the stochastic variables. Moreover, as revealed by the ansatz, each eigenvalue of the dimensionally reduced problem is actually associated with an eigenspace of dimension two.

e Spatially Discretized Eigenvalue Problem.
We employ standard high-order finite elements with polynomial order p ∈ N. For any fixed t > 0, the spatially discretized problem may be written as a parametric matrix eigenvalue problem: find λ p : Γ ⟶ R and y p : Γ ⟶ R n such that where n is the dimension of the discretization space. Here, M and S are the mass matrix and the stiffness matrix, which obviously depend on ξ ∈ Γ. For any fixed ξ ∈ Γ, problem (23) reduces to a positive-definite generalized matrix eigenvalue problem.
We let 〈·, ·〉 M(ξ) and ‖ · ‖ M(ξ) denote the inner product and norm induced by the mass matrix. We assume the eigenvector to be normalized according to ‖y p (ξ)‖ M(ξ) � 1 for all ξ ∈ Γ. Moreover, we fix its sign so that the inner product 〈y p (0), y p (ξ)〉 M(ξ) , where y p (0) is a suitable reference solution, is positive for every ξ ∈ Γ. In practice, we may disregard the possibility that the inner product is zero, Shock and Vibration 5 and therefore, as long as the discrete eigenvalues are simple, the solution is well defined for all ξ ∈ Γ.
For simplicity, we make the assumption that A is monotone in the following sense: whenever β ∈ N M 0 is such that β ≤ α for some α ∈ A, then β ∈ A. Let L p be the univariate Legendre polynomial of degree p. Denote by χ (p) k p k�0 the zeros of L p+1 . We define the one-dimensional Lagrange interpolation operators where ℓ (p) k p k�0 are the related Lagrange basis polynomials of degree p. Finally, we define the sparse collocation operator as e operator (25) may be rewritten in a computationally more convenient form where G α ≔ β ∈ N M 0 |α − 1 ≤ β ≤ α . We see that the complete grid of collocation points is now given by Statistics, such as the expected value and variance, of the solution may now be computed by applying the onedimensional Gauss-Legendre quadrature rules on the components of (26). e accuracy of the collocated approximation is ultimately determined by the smoothness of the solution as well as the choice of the multiindex set A ⊂ N M 0 , see [6] for a detailed analysis. Similar results in the context of source problems have been presented in [22][23][24][25].

Stochastic Collocation on Partial Sparse Grids.
In the framework of this paper, each collocation point in the stochastic space Γ corresponds to a single measurement configuration. In contrast to idealized mathematical models, the range of practical measurement settings may be limited to a small subset of Γ. In this section, we discuss one method to carry out the analysis of the response statistics of a random field in a subset of Γ.
We approach this problem from the viewpoint of polynomial interpolation. In the following, we assume for In addition, we denote We assume that the basis functions in B are ordered in degree lexicographic order. In particular, L 1 is constant.
Let β � (β 1 , . . . , β N ) denote a subsequence of 1, { 2, . . . , J} so that each element of β corresponds to one and only one element of this set. en for any such subsequence, we can construct the Vandermonde-like matrix V β � (L i (χ j )) i∈β,1≤j≤N . If V β is invertible, then for any input f : Γ ⟶ R we can set so that the Lagrange interpolating polynomial If β 1 � 1, then the expansion coefficients can be used to estimate the response statistics of f via Following [12], it turns out that the Vandermonde-like matrices related to the Smolyak interpolating polynomials are remarkably well conditioned. One may find a wellconditioned interpolating polynomial for a subset of a sparse grid by proceeding in the following way: where the coefficients c 1 , . . . , c N ∈ R are obtained as in (30) e use of the maximum volume submatrix ensures both that the Lagrange interpolation polynomial is well defined and that system (45) is sufficiently well conditioned.
In practical calculations, finding the analytical maximum volume submatrix is in general an NP-hard problem. However, there exist several algorithms in the literature which can be used to find the approximate maximum volume submatrix in a computationally efficient manner. We use the MaxVol algorithm proposed in [27]. A related algorithm has been discussed in [28].

Frequency Response Analysis
In frequency response analysis, the idea is to study the excitation or applied force in the frequency domain. erefore, the uncertainty in the eigenproblem is directly translated to frequency response. In the deterministic setting, the procedure is as follows [29]: starting from the equation of motion for the system, where M is the mass matrix, C is the viscous damping matrix, S is the stiffness matrix, f is the force vector, and x is the displacement vector. In our context, M and S are defined as mentioned above, and the viscous damping matrix is taken to be a constant diagonal matrix C � 2δI, where δ � 1/2000. In the case of harmonic excitation, a steady-state solution is sought, and the force and the corresponding response can be expressed as harmonic functions as x � x(ω)e iωt .

(35)
Taking the first and second derivatives of equation (34) and substituting using (35) leads to thus reducing to a linear system of equations: Any quantity of interest can then be derived from the solution x(ω). In the following, we consider maximal transverse deflections w max as the quantity of interest.

Frequency Response Analysis with Sparse Grids.
e collocation methods are invariant with respect to the quantity of interest, and therefore, the application to frequency response is straightforward.
Let us consider four-dimensional second-order Smolyak-Gauss-Legendre quadrature points X GL and the second-order Smolyak-Clenshaw-Curtis quadrature points X CC with 41 points. Let (P k ) ∞ k�0 denote the standard orthonormal univariate Legendre polynomials generated by the three-term recursion: where for α ∈ Z 4 + . With this convention, #B � #X GL � #X CC � 41 and the interpolation problem is well posed for the full sparse grids.
In order to make the problem well posed for some partial grids X GL ⊂ X GL and X CC ⊂ X CC , we proceed as in Section 3.5 and construct the rectangular Vandermonde-like matrix W � (L(ξ)) L∈B,ξ∈X and choose its maximum volume submatrix V using the MaxVol algorithm. We compute the expectation E[w max ] and upper confidence envelope E[w max ] + ��������� Var[w max ] using formulae (32) and (33) for each frequency, respectively.

Numerical Experiments
In the numerical experiments, at least one of the sources of uncertainty is related to some geometric feature, for instance, diameter of a cylinder or local thickness of a shell. In the finite element context, this means that in order to compute statistics over a set of solutions, it is necessary to introduce some nominal domain onto which every other realization of the discretized domain is mapped. In a general situation, this could be done via conformal mappings, and in specific situations as happens to be here, the meshes can be guaranteed to be topologically equivalent ensuring errors only in the immediate vicinity of the random parts of the domain.
Material constants adopted for all simulations are E � 2.069 × 10 11 MPa, ] � 1/3, and ρ � 7868 kg/m 3 , unless otherwise specified. Also, in the following examples, we employ the sparse collocation operator (25) with total degree multiindex sets where applicable. All computations were performed on an Intel Xeon(R) CPU E3-1230 v5 3.40 GHz (eight cores) desktop with 32 GB of RAM.

Tower Configuration: Free Vibration.
Our first example is an idealized wind turbine tower ( Figure 2). We can apply dimension reduction via ansatz as in Section 2.3 and arrive at a highly nontrivial shell-solid configuration. Notice that in this formulation, we model the shell as a solid which is computationally feasible in 2D, and the resulting systems can be solved for different harmonic wavenumbers as outlined above. e tower is assembled with four parts: a thin shell of thickness t � 1/100 and vertical length of 11 units (x ∈ [−1, 10]), a top ring of height C 1 � 2/5 and width 3/4, and two base rings, inner and outer ones, of height D 1 � 2/5 and random widths W 1 and W 2 , respectively. An additional source of randomness comes from varying thickness of the Shock and Vibration shell near the top ring, centered at x � 0. Here, the varying thickness can be interpreted as manufacturing imperfection or damage to the structure. e tower configuration is detailed also in Figure 2.
Let us define the computational domain Ω to be a union of five parameterized domains: t; a, b) , (43) e random variables and their ranges used in different experiments are (Figure 2) where a and b are parameters in the shape modeling the imperfection taken to be symmetric on the inner and outer surfaces: where s is the coordinate scaled to center the shape around x � 0. For a plot of the profiles, see Figure 3.
In each experiment, the dimension of the stochastic space is M � 4, and the random vector ξ ∈ Γ is obtained by scaling the geometric parameters to the interval [−1, 1]. e base of the tower is clamped, that is, fully kinematically constrained.
e p-version of FEM is used with p � 4 resulting in a linear system with 14247 d.o.f. e multiindex set is A L , with L � 3. Per collocation point, the time spent is approximately 40 seconds; 10 seconds for solution and 30 seconds for obtaining statistics.
In Tables 1 and 2, we have listed statistics for eight of the smallest eigenvalues of the tower configuration with the corresponding wavenumber K. Of interest is the pair of fourth and fifth smallest eigenvalues since the values are fairly close to each other; in other words, they form a cluster. Indeed, if we track the modes over the parameter space, we see that crossing occurs, i.e., different eigenmodes are associated with the fourth smallest eigenvalue for different realizations of the problem. Also notice that the actual squares of frequencies vary (at least) within range [4562.30, 4578.79] for K � 4, with the expected value � 4572.40. is means that in a multiparametric setting, the number of modes observed within a fixed range of frequencies is not necessarily constant over the whole parameter space. e statistics presented in Tables 1 and 2 were verified by a Monte Carlo simulation of S � 1000 samples. A comparison of the results obtained by Monte Carlo and collocation algorithms has been presented in Table 3. e central limit theorem guarantees (stochastic) error bounds for statistics computed via Monte Carlo: the error in the expected value of the quantity of interest is given by Table 3, we see that the errors between the two solutions fall well within the tolerances of the Monte Carlo results.  Shock and Vibration e effect of the manufacturing imperfection is illustrated for the eigenmodes associated with K � 1 in Figure 3. High-localized variance in the computed displacements means that the uncertainty in derived quantities of interest such as stresses is likely to be highest. In Figure 3, we see that the variance of the transverse deflection of the third eigenmode for K � 1 is localized exactly at the thinnest part indicating a highly likely location for failure of the tower.
is is of course what we would expect.

Cylinder with Junctions.
Our second example is a cylinder with two junctions: a T-junction with another cylinder and a cone extension at one end of the cylinder. Since we limit ourselves to shells of revolution, we remove the T-junction through clamped boundary conditions. Here, the sources of uncertainty are the random radius (diameter) of the circular T-junction, i.e., a circular hole in our modeling, D 1 , the slope of the cone, ζ, and Young's modulus. e radius of the cylindrical part is constant � 1, and the dimensionless thickness is constant t � 1/100. For the long cylinder, the distance of the center of the hole to ends is L 1 � 10π and the distance to the cone is L 2 � 40π/7, and for the short one, the lengths are simply scaled by 10. For a schematic of the construction, see Figure 4.    Now the computational domain is periodic in the angular direction: where ζ is the slope of the cone. Finally, the parameters ξ 3 and ξ 4 correspond Karhunen-Loéve basis functions, i.e., terms E 3 (x) � sin x and E 4 (x) � sin 2x in the expansion (19). Again, in each experiment the dimension of the stochastic space is M � 4, and the random vector ξ ∈ Γ is obtained by scaling the geometric parameters to the interval [−1, 1]. Both ends of the cylinder are clamped, as well as the T-junction interface, that is, fully kinematically constrained. e p-version of FEM is used with p � 3 resulting in a linear system with 41280 d.o.f. e multiindex set is A L with L � 2. Per collocation point, the time spent is approximately 60 seconds; 10 seconds for solution and 50 seconds for obtaining statistics.
We study both free vibration and frequency response for this configuration.

Free Vibration.
Let us first examine the lowest modes of the two structures. Interestingly, when comparing the two subfigures of Figure 5, the colours indicate that indeed in the long cylinder, the long-range internal layer of 1/ � t √ emerges in the lowest mode. However, none of the shorter characteristic length scales have strong amplitudes in the transverse deflection in this case.
In Table 4, we have listed statistics for four of the smallest eigenvalues of the short and long cylinder configurations. In this case, the eigenvalues appear well separated. Again, the statistics were verified by a Monte Carlo simulation of S � 1000 samples. A comparison of the results obtained by Monte Carlo and collocation algorithms has been presented in Table 5. As discussed before, we see from the figures in Table 5 that the errors between the two solutions fall well within the tolerances of the Monte Carlo results.
Statistics of the first eigenmode for the long cylinder with junctions have been presented in Figure 6. e long-range layer is clearly visible in the second, third, and fifth components. e effect of the cone junction is distinctly manifested in the statistics of the fourth component. e standard deviations of the two smallest eigenvalues are approximately 62 and 50 percent of the respective expected values.
e ratio of the expected values is approximately 0.21, and thus, the eigenvalues at least appear to be well separated.
Statistics of the first eigenmode for the short cylinder with junctions have been presented in Figure 7. In this case, we do not observe a separate long-range layer. e effect of the cone junction is never the less apparent in the statistics of the fourth component. e standard deviations of the two smallest eigenvalues are now approximately 31 and 49 percent of the respective expected values. e ratio of the expected values is approximately 0.24, and thus, the eigenvalues again appear to be well separated.
In Figure 8, the short-range effects of the T-junction have been highlighted. e statistics of the fifth component have been shown in vicinity of the junction for both cylinders. e relative strength of the long-range axial layer in the long cylinder can be observed, whereas in the short cylinder, the angular oscillatory layer is dominant. In both configurations, there is a hint of short-range axial layer emanating from the hole.

Frequency Response Analysis.
For the frequency response analysis of the short cylinder, the load is chosen to be F(x, y) � 1000 cos(y/2)N, and the angular frequency range is taken to be ω ∈ 2π 5, 10, . . . , 200 { }Hz. In our analysis, all modes are present; no attempt to choose, for instance, a subspace of lowest modes has been made.
We first carry out the solution of the frequency response statistics subject to the abovementioned second-order Smolyak-Gauss-Legendre quadrature points X GL and the second-order Smolyak-Clenshaw-Curtis quadrature points X CC and then consider their respective subsets with 25 points defined by setting where the choice ξ 1 � 0 is arbitrary. In the context of this numerical experiment, it means that the diameter of the hole is fixed at its expected value. Figure 4: Cylinder-cylinder-cone schematic. e cylinder diameter R 1 � 2, L 1 , L 2 , and C 1 are deterministic: L 1 + L 2 + C 1 � 20π; the diameter of the hole D 1 is random, and the cone diameter R 2 is random since it depends on the random slope of the cone. 10 Shock and Vibration Let us first consider sparse grids X GL and X CC . In Figure 9, we show two frequency response graphs where the expected values of the quantity of interest, the maximal transverse deflection w max , are shown along standard deviations. It is clear that the frequency response is sensitive to perturbations to the parameters. e two point sets are not hierarchic, and thus, one cannot expect to get exactly identical responses. Moreover, distributions of the values of the quantity of interest appear to be exponential.
is is illustrated in Figure 10. e maximal transverse deflection w max at any given frequency depends on how well the tail of the distribution is approximated.   Let us next move to the partial sparse grids X ∈ X CC , X GL . e results are displayed in Figure 11. e results should be compared with Figure 9; in particular, sampling the stochastic space with fixed first component already characterizes the locations of the maxima in the frequency response. However, the maximal amplitudes are different, again reflecting the fact that the points sets are not hierarchic.

Remark 1.
We emphasise that if it were possible to identify real measurements as points in the grid, with this approach, one could augment simulations with real data and arrive at more realistic interpolation operators.
is example underlines the fact that in multiparametric situations, it is always necessary to think in terms of distributions and relate the observed statistics to them.

Conclusions
Multiparametric vibration problems and especially those involving thin domains are and remain challenging. We have studied two vibration problems rooted in practice, yet artificial in terms of the parametrization of the random model. e efficacy of the stochastic collocation has been demonstrated, and its extension to frequency response analysis has been demonstrated. More specifically, in the context of thin shells, the boundary layers, including internal ones, are shown to have the predicted characteristic length scales, and their contribution to statistics such as variance is clearly illustrated in the results. In vibration-related problem, the frequency response analysis, the importance of   gaining complete statistical understanding of the results is underlined.

Data Availability
No data were used to support this study.