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

Authors: Harri Hakula; Vesa Kaarnioja; Mikael Laaksonen
Journal: Shock and Vibration
 (2018)

---

## Abstract

Numerical simulation of thin solids remains one of the challenges in computational mechanics. The 3D elasticity problems of shells of revolution are dimensionally reduced in different ways depending on the symmetries of the configurations resulting in corresponding 2D models. In this paper, we solve the multiparametric free vibration of complex shell configurations under uncertainty using stochastic collocation with thep-version of finite element method and apply the collocation approach to frequency response analysis. In numerical examples, the sources of uncertainty are related to material parameters and geometry representing manufacturing imperfections. All stochastic collocation results have been verified with Monte Carlo methods.


---

## Body

## 1. Introduction

Numerical simulation of thin solids remains one of the challenges in computational mechanics. The 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–4]. There 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 anonintrusive approach, the stochastic collocation [6]. Collocation methods are particularly appealing in problems with parameterized (random) domains. The 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–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.The 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. The wind turbine tower problem is inspired by an earlier study [13], where the tower is modeled as a solid due to solver limitations. The 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. The 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 usingp-version of the finite element method [16]. The results derived using collocation on sparse grids have been verified using the Monte Carlo method.The rest of this paper is organized as follows. In Section2, 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.


## 2. Shell Eigenproblem

Assuming a time harmonic displacement field, the free vibration problem for a general shell leads to the following abstracteigenvalue problem: find u∈ℝ3 and ω2∈ℝ such that(1)Su=ω2Mu,+boundary  conditions,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.


### 2.1. Shell Geometry

In the following, we let a profile functiony=fx 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 R3 of type(2)Ω=x+ηxnxx∈D,−d0x<ηx<d1x,where D is a (mid)surface of revolution, nx is the unit normal to D, and dx=d1x+d0x is the thickness of the shell measured in the direction of the normal. An illustration of such a configuration is given in Figure 1.Figure 1
Shell of revolution with variable axial thickness; midsurface indicated with (a) red colour and (b) dashed line. (a) Shell of revolution. (b) Cross section; midsurfacef(x) = 1 (dashed), lower surface f(x) − d0(x) and upper surface f(x) + d1(x) (solid).
(a)
(b)


#### 2.1.1. Constant (Dimensionless) Thickness

Let us havedx=d (constant) and define principal curvature coordinates, where only four parameters, the radii of principal curvature R1, R2, and the so-called Lamé parameters, A1, A2, which relate coordinates changes to arc lengths, are needed to specify the curvature and the metric on D. The displacement vector field of the midsurface u=u,v,w can be interpreted as projections to directions(3)e1=1A1∂Ψ∂x1,e2=1A2∂Ψ∂x2,e3=e1×e2,where Ψx1,x2 is a suitable parametrization of the surface of revolution, e1,e2 are the unit tangent vectors along the principal curvature lines, and e3 is the unit normal. In other words,(4)u=ue1+ve2+we3.Assuming that the shell profile is given byfx1, the geometry parameters are(5)A1x1=1+f′x12,A2x1=fx1,R1x1=−A1x13f″x1,R2x1=A1x1A2x1.Both cylinders and cones are examples of parabolic shells according to the theory of surfaces by Gauss. This follows from the fact that, in both cases,f″x1=0. In particular the reciprocal function(6)1R1x1=0.For a cylinder with a constant radius,fx1=1, say, we get(7)A1x1=A2x1=R2x1=1,1R1x1=0.In the following, we shall keep the general form of the parameters in the equations.We can now define thedimensionless thickness t as t=d/R, where R∼R2 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.


### 2.2. Reduction in Thickness

The 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 ut and ω2t∈ℝ such that(8)tAmut+tAsut+t3Abut=ω2tMtut,+boundary  conditions.where ut represents the shell displacement field, while ω2t represents the square of the eigenfrequency. The differential operators Am, As, and Ab account for membrane, shear, and bending potential energies, respectively, and are independent of t. Finally, Mt is the inertia operator, which in this case can be split into the sum Mt=tMl+t3Mr, with Ml (displacements) and Mr (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 ut,ω2t∈V×R such that(9)tamut,v+tasut,v+t3abut,v=ω2tmt;ut,v∀v∈V,where am⋅,⋅, as⋅,⋅, ab⋅,⋅, and mt;⋅,⋅ are the bilinear forms associated with the operators Am, As, Ab, and Mt, respectively. Obviously, the space V and the three bilinear forms depend on the chosen shell model [17].


#### 2.2.1. 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. The 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⊂R2 is given by the surface parametrization, and the axial/angular coordinates are denoted by x and y:(10)Ω^t=x,ya≤x≤b,0≤y<2π.Deformation energyAu,u is divided into bending, membrane, and shear energies, denoted by subscripts B, M, and S, respectively.(11)Au,u=t2ABu,u+AMu,u+ASu,u.Notice that we can safely cancel out one power oft.Bending, membrane, and shear energies are given as follows [18]:(12)t2ABu,u=t2∫Ω^tEx,yνκ11u+κ22u2+1−ν∑i,j=12κiju2A1x,yA2x,ydxdy,AMu,u=12∫Ω^tEx,yνβ11u+β22u2+1−ν∑i,j=12βiju2A1x,yA2x,ydxdy,ASu,u=61−ν∫Ω^tEx,yρ1u2+ρ2u2×A1x,yA2x,ydxdy,where v is the Poisson ratio (constant) and Ex,y is Young’s modulus with scaling 1/121−ν2.The symmetric 2D strains, bending (κ), membrane (β), and shear (ρ), are defined as(13)κ11=1A1∂θ∂x+ψA1A2∂A2∂y,κ22=1A2∂ψ∂y+θA1A2∂A2∂x,κ12=121A1∂ψ∂x+1A2∂θ∂y−θA1A2∂A2∂y−ψA1A2∂A2∂x−1R11A2∂u∂y−vA1A2∂A2∂x−1R21A1∂v∂x−uA1A2∂A1∂y,β11=1A1∂u∂x+vA1A2∂A1∂y+wR1,β22=1A2∂v∂y+uA1A2∂A2∂x+wR2,β12=121A1∂v∂x+1A2∂u∂y−uA1A2∂A1∂y−vA1A2∂A2∂x,ρ1=1A1∂w∂x−uR1−θ,ρ2=1A2∂w∂y−vR2−ψ.In our experiments, the shell profiles are functions ofx only, but not necessarily constant; hence, the strains are presented in a general form.The stiffness matrixS 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 A1x,yA2x,y.


### 2.3. Reduction through Ansatz

The 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 ux,y,z or ux,r,α in cylindrical coordinates have either one of the forms:(14)u1x,y,z=u1x,r,α=ux,rcosKαvx,rsinKαwx,rcosKα,(15)u2x,y,z=u2x,r,α=ux,rsinKαvx,rcosKαwx,rsinKα,where K is the harmonic wavenumber. Thus, given a profile function r=y=fx, with x∈a,b, the computational domain Ω^ is 2D only:(16)Ω^=x,r+ηxa≤x≤b,−d0x<ηx<d1x,where the auxiliary functions dix again simply indicate that the thickness dx=d1x+d0x need not be constant.Deformation energyAu,u is defined as(17)Au,u=r3∫Ω^Ex,yν1+ν1−2νtrϵ2+11+ν∑i,j=13ϵij2dxdη,where the symmetric 3D strains ϵij are(18)ϵ11=∂u∂x,ϵ12=12−Kr+ηu+∂v∂x,ϵ13=12∂u∂η+∂w∂x,ϵ22=1r+ηKv+w,ϵ23=12∂v∂η−1r+ηv+Kw,ϵ33=∂w∂η.The stiffness matrixS is obtained after integration and assembly. Here, the mass matrix M is the standard 3D one with density σ.


### 2.4. Special Features of Shell Problems

#### 2.4.1. Boundary Layers

The 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. The 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 ∼t4. There also exists an axial internal layer of a very long range ∼1/t.Identification of the boundary layer structure orlayer 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.


#### 2.4.2. Numerical Locking

One of the numerical difficulties associated with thin structures is the so-callednumerical 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.


## 2.1. Shell Geometry

In the following, we let a profile functiony=fx 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 R3 of type(2)Ω=x+ηxnxx∈D,−d0x<ηx<d1x,where D is a (mid)surface of revolution, nx is the unit normal to D, and dx=d1x+d0x is the thickness of the shell measured in the direction of the normal. An illustration of such a configuration is given in Figure 1.Figure 1
Shell of revolution with variable axial thickness; midsurface indicated with (a) red colour and (b) dashed line. (a) Shell of revolution. (b) Cross section; midsurfacef(x) = 1 (dashed), lower surface f(x) − d0(x) and upper surface f(x) + d1(x) (solid).
(a)
(b)


### 2.1.1. Constant (Dimensionless) Thickness

Let us havedx=d (constant) and define principal curvature coordinates, where only four parameters, the radii of principal curvature R1, R2, and the so-called Lamé parameters, A1, A2, which relate coordinates changes to arc lengths, are needed to specify the curvature and the metric on D. The displacement vector field of the midsurface u=u,v,w can be interpreted as projections to directions(3)e1=1A1∂Ψ∂x1,e2=1A2∂Ψ∂x2,e3=e1×e2,where Ψx1,x2 is a suitable parametrization of the surface of revolution, e1,e2 are the unit tangent vectors along the principal curvature lines, and e3 is the unit normal. In other words,(4)u=ue1+ve2+we3.Assuming that the shell profile is given byfx1, the geometry parameters are(5)A1x1=1+f′x12,A2x1=fx1,R1x1=−A1x13f″x1,R2x1=A1x1A2x1.Both cylinders and cones are examples of parabolic shells according to the theory of surfaces by Gauss. This follows from the fact that, in both cases,f″x1=0. In particular the reciprocal function(6)1R1x1=0.For a cylinder with a constant radius,fx1=1, say, we get(7)A1x1=A2x1=R2x1=1,1R1x1=0.In the following, we shall keep the general form of the parameters in the equations.We can now define thedimensionless thickness t as t=d/R, where R∼R2 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.


## 2.1.1. Constant (Dimensionless) Thickness

Let us havedx=d (constant) and define principal curvature coordinates, where only four parameters, the radii of principal curvature R1, R2, and the so-called Lamé parameters, A1, A2, which relate coordinates changes to arc lengths, are needed to specify the curvature and the metric on D. The displacement vector field of the midsurface u=u,v,w can be interpreted as projections to directions(3)e1=1A1∂Ψ∂x1,e2=1A2∂Ψ∂x2,e3=e1×e2,where Ψx1,x2 is a suitable parametrization of the surface of revolution, e1,e2 are the unit tangent vectors along the principal curvature lines, and e3 is the unit normal. In other words,(4)u=ue1+ve2+we3.Assuming that the shell profile is given byfx1, the geometry parameters are(5)A1x1=1+f′x12,A2x1=fx1,R1x1=−A1x13f″x1,R2x1=A1x1A2x1.Both cylinders and cones are examples of parabolic shells according to the theory of surfaces by Gauss. This follows from the fact that, in both cases,f″x1=0. In particular the reciprocal function(6)1R1x1=0.For a cylinder with a constant radius,fx1=1, say, we get(7)A1x1=A2x1=R2x1=1,1R1x1=0.In the following, we shall keep the general form of the parameters in the equations.We can now define thedimensionless thickness t as t=d/R, where R∼R2 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.


## 2.2. Reduction in Thickness

The 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 ut and ω2t∈ℝ such that(8)tAmut+tAsut+t3Abut=ω2tMtut,+boundary  conditions.where ut represents the shell displacement field, while ω2t represents the square of the eigenfrequency. The differential operators Am, As, and Ab account for membrane, shear, and bending potential energies, respectively, and are independent of t. Finally, Mt is the inertia operator, which in this case can be split into the sum Mt=tMl+t3Mr, with Ml (displacements) and Mr (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 ut,ω2t∈V×R such that(9)tamut,v+tasut,v+t3abut,v=ω2tmt;ut,v∀v∈V,where am⋅,⋅, as⋅,⋅, ab⋅,⋅, and mt;⋅,⋅ are the bilinear forms associated with the operators Am, As, Ab, and Mt, respectively. Obviously, the space V and the three bilinear forms depend on the chosen shell model [17].


### 2.2.1. 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. The 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⊂R2 is given by the surface parametrization, and the axial/angular coordinates are denoted by x and y:(10)Ω^t=x,ya≤x≤b,0≤y<2π.Deformation energyAu,u is divided into bending, membrane, and shear energies, denoted by subscripts B, M, and S, respectively.(11)Au,u=t2ABu,u+AMu,u+ASu,u.Notice that we can safely cancel out one power oft.Bending, membrane, and shear energies are given as follows [18]:(12)t2ABu,u=t2∫Ω^tEx,yνκ11u+κ22u2+1−ν∑i,j=12κiju2A1x,yA2x,ydxdy,AMu,u=12∫Ω^tEx,yνβ11u+β22u2+1−ν∑i,j=12βiju2A1x,yA2x,ydxdy,ASu,u=61−ν∫Ω^tEx,yρ1u2+ρ2u2×A1x,yA2x,ydxdy,where v is the Poisson ratio (constant) and Ex,y is Young’s modulus with scaling 1/121−ν2.The symmetric 2D strains, bending (κ), membrane (β), and shear (ρ), are defined as(13)κ11=1A1∂θ∂x+ψA1A2∂A2∂y,κ22=1A2∂ψ∂y+θA1A2∂A2∂x,κ12=121A1∂ψ∂x+1A2∂θ∂y−θA1A2∂A2∂y−ψA1A2∂A2∂x−1R11A2∂u∂y−vA1A2∂A2∂x−1R21A1∂v∂x−uA1A2∂A1∂y,β11=1A1∂u∂x+vA1A2∂A1∂y+wR1,β22=1A2∂v∂y+uA1A2∂A2∂x+wR2,β12=121A1∂v∂x+1A2∂u∂y−uA1A2∂A1∂y−vA1A2∂A2∂x,ρ1=1A1∂w∂x−uR1−θ,ρ2=1A2∂w∂y−vR2−ψ.In our experiments, the shell profiles are functions ofx only, but not necessarily constant; hence, the strains are presented in a general form.The stiffness matrixS 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 A1x,yA2x,y.


## 2.2.1. 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. The 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⊂R2 is given by the surface parametrization, and the axial/angular coordinates are denoted by x and y:(10)Ω^t=x,ya≤x≤b,0≤y<2π.Deformation energyAu,u is divided into bending, membrane, and shear energies, denoted by subscripts B, M, and S, respectively.(11)Au,u=t2ABu,u+AMu,u+ASu,u.Notice that we can safely cancel out one power oft.Bending, membrane, and shear energies are given as follows [18]:(12)t2ABu,u=t2∫Ω^tEx,yνκ11u+κ22u2+1−ν∑i,j=12κiju2A1x,yA2x,ydxdy,AMu,u=12∫Ω^tEx,yνβ11u+β22u2+1−ν∑i,j=12βiju2A1x,yA2x,ydxdy,ASu,u=61−ν∫Ω^tEx,yρ1u2+ρ2u2×A1x,yA2x,ydxdy,where v is the Poisson ratio (constant) and Ex,y is Young’s modulus with scaling 1/121−ν2.The symmetric 2D strains, bending (κ), membrane (β), and shear (ρ), are defined as(13)κ11=1A1∂θ∂x+ψA1A2∂A2∂y,κ22=1A2∂ψ∂y+θA1A2∂A2∂x,κ12=121A1∂ψ∂x+1A2∂θ∂y−θA1A2∂A2∂y−ψA1A2∂A2∂x−1R11A2∂u∂y−vA1A2∂A2∂x−1R21A1∂v∂x−uA1A2∂A1∂y,β11=1A1∂u∂x+vA1A2∂A1∂y+wR1,β22=1A2∂v∂y+uA1A2∂A2∂x+wR2,β12=121A1∂v∂x+1A2∂u∂y−uA1A2∂A1∂y−vA1A2∂A2∂x,ρ1=1A1∂w∂x−uR1−θ,ρ2=1A2∂w∂y−vR2−ψ.In our experiments, the shell profiles are functions ofx only, but not necessarily constant; hence, the strains are presented in a general form.The stiffness matrixS 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 A1x,yA2x,y.


## 2.3. Reduction through Ansatz

The 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 ux,y,z or ux,r,α in cylindrical coordinates have either one of the forms:(14)u1x,y,z=u1x,r,α=ux,rcosKαvx,rsinKαwx,rcosKα,(15)u2x,y,z=u2x,r,α=ux,rsinKαvx,rcosKαwx,rsinKα,where K is the harmonic wavenumber. Thus, given a profile function r=y=fx, with x∈a,b, the computational domain Ω^ is 2D only:(16)Ω^=x,r+ηxa≤x≤b,−d0x<ηx<d1x,where the auxiliary functions dix again simply indicate that the thickness dx=d1x+d0x need not be constant.Deformation energyAu,u is defined as(17)Au,u=r3∫Ω^Ex,yν1+ν1−2νtrϵ2+11+ν∑i,j=13ϵij2dxdη,where the symmetric 3D strains ϵij are(18)ϵ11=∂u∂x,ϵ12=12−Kr+ηu+∂v∂x,ϵ13=12∂u∂η+∂w∂x,ϵ22=1r+ηKv+w,ϵ23=12∂v∂η−1r+ηv+Kw,ϵ33=∂w∂η.The stiffness matrixS is obtained after integration and assembly. Here, the mass matrix M is the standard 3D one with density σ.


## 2.4. Special Features of Shell Problems

### 2.4.1. Boundary Layers

The 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. The 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 ∼t4. There also exists an axial internal layer of a very long range ∼1/t.Identification of the boundary layer structure orlayer 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.


### 2.4.2. Numerical Locking

One of the numerical difficulties associated with thin structures is the so-callednumerical 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.


## 2.4.1. Boundary Layers

The 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. The 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 ∼t4. There also exists an axial internal layer of a very long range ∼1/t.Identification of the boundary layer structure orlayer 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.


## 2.4.2. Numerical Locking

One of the numerical difficulties associated with thin structures is the so-callednumerical 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.


## 3. Stochastic Shell Eigenproblem

The 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. The approximate solution statistics are then computed employing a sparse grid stochastic collocation operator.


### 3.1. The Stochastic Eigenproblem

We let the material and geometric uncertainties be represented by a vector of mutually independent random variablesξ=ξ1,ξ2,…,ξM taking values in a suitable domain Γ⊂RM for some M∈N. For uniformly distributed random variables, we may, without loss of generality, assume a scaling such that Γ≔−1,1M.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(19)Ex,ξ=E0x+∑m=m′+1MEmxξm,x∈Dξ.The 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 Emin,Emax>0 such that(20)Emin≤ess  infx∈DξEx,ξ≤ess  supx∈DξEx,ξ≤Emax,for all ξ∈Γ.LetLμ2Γ denote a weighted L2-space, where μ is the uniform product probability measure associated with the probability distribution of ξ∈Γ. For functions in Lμ2Γ, we define the expected value(21)Ev=∫Γvξdμξ,and variance Varv=Ev−Ev2.Assuming stochastic models for the computational domain and Young’s modulus, the eigenpairs of the free vibration problem now depend onξ∈Γ. The 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 ut:Γ⟶V and ω2t:Γ⟶R such that for all ξ∈Γ, we have(22)tamξ;ut,ξ,v+tasξ;ut,ξ,v+t3abξ;ut,ξ,v=ω2t,ξmt,ξ;ut,ξ,v,∀v∈V,where amξ;⋅,⋅, asξ;⋅,⋅, abξ;⋅,⋅, and mt,ξ;⋅,⋅ are stochastic equivalents of the deterministic bilinear forms in (9). We assume that the eigenvector ut,ξ is normalized with respect to the inner product mt,ξ;⋅,⋅ for all ξ∈Γ.


### 3.2. 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. This 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. The ordering of the eigenvalues associated with each mode may be different for different realizations of the problem. This 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.


### 3.3. The Spatially Discretized Eigenvalue Problem

We employ standard high-order finite elements with polynomial orderp∈N. For any fixed t>0, the spatially discretized problem may be written as a parametric matrix eigenvalue problem: find λp:Γ⟶R and yp:Γ⟶Rn such that(23)Sξypξ=λpξMξypξ,∀ξ∈Γ,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 ypξMξ=1 for all ξ∈Γ. Moreover, we fix its sign so that the inner product yp0,ypξMξ, where yp0 is a suitable reference solution, is positive for every ξ∈Γ. In practice, we may disregard the possibility that the inner product is zero, and therefore, as long as the discrete eigenvalues are simple, the solution is well defined for all ξ∈Γ.


### 3.4. Stochastic Collocation on Sparse Grids

We introduce an anisotropic Smolyak-type sparse grid collocation operator ([6, 22, 23]) for resolving the effects of the random variables ξ∈Γ.Assume a finite multiindex setA⊂N0M. For α,β∈A we write α≤β if αm≤βm for all m≥1. For simplicity, we make the assumption that A is monotone in the following sense: whenever β∈N0M is such that β≤α for some α∈A, then β∈A. Let Lp be the univariate Legendre polynomial of degree p. Denote by χkpk=0p the zeros of Lp+1. We define the one-dimensional Lagrange interpolation operators ℐpm via(24)ℐpmvξm=∑k=0pvχkpℓkpξm,where ℓkpk=0p are the related Lagrange basis polynomials of degree p. Finally, we define the sparse collocation operator as(25)ℐA≔∑α∈A⊗m=1Mℐαmm−ℐαm−1m.The operator (25) may be rewritten in a computationally more convenient form(26)ℐA=∑α∈A∑β∈Gα−1α−β1⊗m=1Mℐβmm,where Gα≔β∈N0Mα−1≤β≤α. We see that the complete grid of collocation points is now given by(27)XA≔∪α∈A∪β∈Gα∏m=1Mχkβmk=0βm=∪α∈A∏m=1Mχkαmk=0αm.Statistics, such as the expected value and variance, of the solution may now be computed by applying the one-dimensional Gauss–Legendre quadrature rules on the components of (26). The accuracy of the collocated approximation is ultimately determined by the smoothness of the solution as well as the choice of the multiindex set A⊂N0M, see [6] for a detailed analysis. Similar results in the context of source problems have been presented in [22–25].


### 3.5. 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 simplicity thatΓ=−1,1M. Let X=χ1,…,χN⊂Γ denote a set of collocation points, and let ℬ=ℒ1,…,ℒJ, J≥N, be a set of orthogonal polynomial basis functions with respect to the measure dμ satisfying(28)∫Γℒiξℒjξdμξ=0wheneveri≠j.In addition, we denote(29)γk=∫Γℒkξ2dμξ1/2fork∈1,…,J.We assume that the basis functions inℬ are ordered in degree lexicographic order. In particular, ℒ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. Then for any such subsequence, we can construct the Vandermonde-like matrix Vβ=ℒiχji∈β,1≤j≤N. If Vβ is invertible, then for any input f:Γ⟶R we can set(30)c1⋯cN=fχ1⋯fχNVℐ−1,so that the Lagrange interpolating polynomial(31)Pξ=∑i=1Nciℒβiξ,satisfies Pχi=fχi for all i∈1,…,N. If β1=1, then the expansion coefficients can be used to estimate the response statistics of f via(32)Ef≈γ1−1c1,(33)Varf≈γβ2−1c22+⋯+γβN−1cN2.Following [12], it turns out that the Vandermonde-like matrices related to the Smolyak interpolating polynomials are remarkably well conditioned. One may find a well-conditioned interpolating polynomial for a subset of a sparse grid by proceeding in the following way:(i)
LetX=χ1,…,χN be a subset of a sparse grid X=χ1,…,χJ(ii)
Letℬ=ℒ1,…,ℒJ be the Smolyak polynomial basis functions (cf. [26])(iii)
Letβ=β1,…,βN denote a subsequence corresponding to the row indices of the maximum volume submatrix of W=ℒiχj1≤i≤J,1≤j≤N(iv)
Form the Lagrange interpolating polynomialPfξ=∑i=1Nciℒβiξ, where the coefficients c1,…,cN∈R are obtained as in (30)The 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].


## 3.1. The Stochastic Eigenproblem

We let the material and geometric uncertainties be represented by a vector of mutually independent random variablesξ=ξ1,ξ2,…,ξM taking values in a suitable domain Γ⊂RM for some M∈N. For uniformly distributed random variables, we may, without loss of generality, assume a scaling such that Γ≔−1,1M.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(19)Ex,ξ=E0x+∑m=m′+1MEmxξm,x∈Dξ.The 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 Emin,Emax>0 such that(20)Emin≤ess  infx∈DξEx,ξ≤ess  supx∈DξEx,ξ≤Emax,for all ξ∈Γ.LetLμ2Γ denote a weighted L2-space, where μ is the uniform product probability measure associated with the probability distribution of ξ∈Γ. For functions in Lμ2Γ, we define the expected value(21)Ev=∫Γvξdμξ,and variance Varv=Ev−Ev2.Assuming stochastic models for the computational domain and Young’s modulus, the eigenpairs of the free vibration problem now depend onξ∈Γ. The 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 ut:Γ⟶V and ω2t:Γ⟶R such that for all ξ∈Γ, we have(22)tamξ;ut,ξ,v+tasξ;ut,ξ,v+t3abξ;ut,ξ,v=ω2t,ξmt,ξ;ut,ξ,v,∀v∈V,where amξ;⋅,⋅, asξ;⋅,⋅, abξ;⋅,⋅, and mt,ξ;⋅,⋅ are stochastic equivalents of the deterministic bilinear forms in (9). We assume that the eigenvector ut,ξ is normalized with respect to the inner product mt,ξ;⋅,⋅ for all ξ∈Γ.


## 3.2. 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. This 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. The ordering of the eigenvalues associated with each mode may be different for different realizations of the problem. This 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.


## 3.3. The Spatially Discretized Eigenvalue Problem

We employ standard high-order finite elements with polynomial orderp∈N. For any fixed t>0, the spatially discretized problem may be written as a parametric matrix eigenvalue problem: find λp:Γ⟶R and yp:Γ⟶Rn such that(23)Sξypξ=λpξMξypξ,∀ξ∈Γ,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 ypξMξ=1 for all ξ∈Γ. Moreover, we fix its sign so that the inner product yp0,ypξMξ, where yp0 is a suitable reference solution, is positive for every ξ∈Γ. In practice, we may disregard the possibility that the inner product is zero, and therefore, as long as the discrete eigenvalues are simple, the solution is well defined for all ξ∈Γ.


## 3.4. Stochastic Collocation on Sparse Grids

We introduce an anisotropic Smolyak-type sparse grid collocation operator ([6, 22, 23]) for resolving the effects of the random variables ξ∈Γ.Assume a finite multiindex setA⊂N0M. For α,β∈A we write α≤β if αm≤βm for all m≥1. For simplicity, we make the assumption that A is monotone in the following sense: whenever β∈N0M is such that β≤α for some α∈A, then β∈A. Let Lp be the univariate Legendre polynomial of degree p. Denote by χkpk=0p the zeros of Lp+1. We define the one-dimensional Lagrange interpolation operators ℐpm via(24)ℐpmvξm=∑k=0pvχkpℓkpξm,where ℓkpk=0p are the related Lagrange basis polynomials of degree p. Finally, we define the sparse collocation operator as(25)ℐA≔∑α∈A⊗m=1Mℐαmm−ℐαm−1m.The operator (25) may be rewritten in a computationally more convenient form(26)ℐA=∑α∈A∑β∈Gα−1α−β1⊗m=1Mℐβmm,where Gα≔β∈N0Mα−1≤β≤α. We see that the complete grid of collocation points is now given by(27)XA≔∪α∈A∪β∈Gα∏m=1Mχkβmk=0βm=∪α∈A∏m=1Mχkαmk=0αm.Statistics, such as the expected value and variance, of the solution may now be computed by applying the one-dimensional Gauss–Legendre quadrature rules on the components of (26). The accuracy of the collocated approximation is ultimately determined by the smoothness of the solution as well as the choice of the multiindex set A⊂N0M, see [6] for a detailed analysis. Similar results in the context of source problems have been presented in [22–25].


## 3.5. 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 simplicity thatΓ=−1,1M. Let X=χ1,…,χN⊂Γ denote a set of collocation points, and let ℬ=ℒ1,…,ℒJ, J≥N, be a set of orthogonal polynomial basis functions with respect to the measure dμ satisfying(28)∫Γℒiξℒjξdμξ=0wheneveri≠j.In addition, we denote(29)γk=∫Γℒkξ2dμξ1/2fork∈1,…,J.We assume that the basis functions inℬ are ordered in degree lexicographic order. In particular, ℒ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. Then for any such subsequence, we can construct the Vandermonde-like matrix Vβ=ℒiχji∈β,1≤j≤N. If Vβ is invertible, then for any input f:Γ⟶R we can set(30)c1⋯cN=fχ1⋯fχNVℐ−1,so that the Lagrange interpolating polynomial(31)Pξ=∑i=1Nciℒβiξ,satisfies Pχi=fχi for all i∈1,…,N. If β1=1, then the expansion coefficients can be used to estimate the response statistics of f via(32)Ef≈γ1−1c1,(33)Varf≈γβ2−1c22+⋯+γβN−1cN2.Following [12], it turns out that the Vandermonde-like matrices related to the Smolyak interpolating polynomials are remarkably well conditioned. One may find a well-conditioned interpolating polynomial for a subset of a sparse grid by proceeding in the following way:(i)
LetX=χ1,…,χN be a subset of a sparse grid X=χ1,…,χJ(ii)
Letℬ=ℒ1,…,ℒJ be the Smolyak polynomial basis functions (cf. [26])(iii)
Letβ=β1,…,βN denote a subsequence corresponding to the row indices of the maximum volume submatrix of W=ℒiχj1≤i≤J,1≤j≤N(iv)
Form the Lagrange interpolating polynomialPfξ=∑i=1Nciℒβiξ, where the coefficients c1,…,cN∈R are obtained as in (30)The 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].


## 4. Frequency Response Analysis

In frequency response analysis, the idea is to study the excitation or applied force in the frequency domain. Therefore, 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,(34)Mx¨+Cx˙+Sx=f,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(35)f=f^ωeiωt,x=x^ωeiωt.Taking the first and second derivatives of equation (34) and substituting using (35) leads to(36)−ω2Mx^ωeiωt+iωCx^ωeiωt+Sx^ωeiωt=f^ωeiωt,thus reducing to a linear system of equations:(37)−ω2M+iωC+Sx^ω=f^ω.Any quantity of interest can then be derived from the solutionx^ω. In the following, we consider maximal transverse deflections wmax as the quantity of interest.


### 4.1. Frequency Response Analysis with Sparse Grids

The 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 pointsXGL and the second-order Smolyak-Clenshaw-Curtis quadrature points XCC with 41 points. Let Pkk=0∞ denote the standard orthonormal univariate Legendre polynomials generated by the three-term recursion:(38)P0x=1,P1x=x,k+1Pk+1x=2k+1xPkx−kPk−1x,k≥1,and define the sequence m0=0, m1=1, and mk=2k−1+1 for k>1. The standard Smolyak polynomial basis functions for the dimension 4, second-level Smolyak rule are then given by [26](39)ℬ=∪α∈Z+44≤α1≤4+2Bα,where(40)Bα=Pi1−1ξ1⋯Pi4−1ξ4;mαj−1+1≤ij≤mαj,1≤j≤4for α∈Z+4. With this convention, #ℬ=#XGL=#XCC=41 and the interpolation problem is well posed for the full sparse grids.In order to make the problem well posed for some partial gridsXGL⊂XGL and XCC⊂XCC, we proceed as in Section 3.5 and construct the rectangular Vandermonde-like matrix W=ℒξℒ∈ℬ,ξ∈X and choose its maximum volume submatrix V using the MaxVol algorithm. We compute the expectation Ewmax and upper confidence envelope Ewmax+Varwmax using formulae (32) and (33) for each frequency, respectively.


## 4.1. Frequency Response Analysis with Sparse Grids

The 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 pointsXGL and the second-order Smolyak-Clenshaw-Curtis quadrature points XCC with 41 points. Let Pkk=0∞ denote the standard orthonormal univariate Legendre polynomials generated by the three-term recursion:(38)P0x=1,P1x=x,k+1Pk+1x=2k+1xPkx−kPk−1x,k≥1,and define the sequence m0=0, m1=1, and mk=2k−1+1 for k>1. The standard Smolyak polynomial basis functions for the dimension 4, second-level Smolyak rule are then given by [26](39)ℬ=∪α∈Z+44≤α1≤4+2Bα,where(40)Bα=Pi1−1ξ1⋯Pi4−1ξ4;mαj−1+1≤ij≤mαj,1≤j≤4for α∈Z+4. With this convention, #ℬ=#XGL=#XCC=41 and the interpolation problem is well posed for the full sparse grids.In order to make the problem well posed for some partial gridsXGL⊂XGL and XCC⊂XCC, we proceed as in Section 3.5 and construct the rectangular Vandermonde-like matrix W=ℒξℒ∈ℬ,ξ∈X and choose its maximum volume submatrix V using the MaxVol algorithm. We compute the expectation Ewmax and upper confidence envelope Ewmax+Varwmax using formulae (32) and (33) for each frequency, respectively.


## 5. 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 areE=2.069×1011 MPa, ν=1/3, and ρ=7868 kg/m3, unless otherwise specified. Also, in the following examples, we employ the sparse collocation operator (25) with total degree multiindex sets(41)A=AL≔α∈N0M∑m=1αm≤L,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.


### 5.1. Tower Configuration: Free Vibration

Our first example is an idealized wind turbine tower (Figure2). 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.Figure 2
Tower configuration: Computational domain;B indicates the location of the assumed manufacturing imperfections at x=0. (a) Shell-solid model; diameter of the base is random. (b) Computational domain. (c) Lowest mode; normalized transverse detection w; rainbow colour scale.
(a)
(b)
(c)The tower is assembled with four parts: a thin shell of thicknesst=1/100 and vertical length of 11 units (x∈−1,10), a top ring of height C1=2/5 and width 3/4, and two base rings, inner and outer ones, of height D1=2/5 and random widths W1 and W2, respectively. An additional source of randomness comes from varying thickness of the shell near the top ring, centered at x=0. Here, the varying thickness can be interpreted as manufacturing imperfection or damage to the structure. The tower configuration is detailed also in Figure 2.Let us define the computational domainΩ^ to be a union of five parameterized domains:(42)Ω^=∪i=15Ω^i,where(43)Ω^1=x,r+η−1≤x≤−1−C1,−3/4−t/2≤η≤t/2,Ω^2=x,r+η−1−C1≤x≤−b,−t/2≤η≤t/2,Ω^3=x,r+η−b≤x≤b,−ϕbx,t;a,b≤η≤ϕbx,t;a,b,Ω^4=x,r+ηb≤x≤10−D1,−t/2≤η≤t/2,Ω^5=x,r+η10−D1≤x≤10,−t/2−W1≤η≤t/2+W2.The random variables and their ranges used in different experiments are (Figure2)(44)TowerW11/4,3/4W21/4,5/4a0,1/2b3/40,1/8,where a and b are parameters in the shape modeling the imperfection taken to be symmetric on the inner and outer surfaces:(45)ϕs,t;a,b=ta−1b4s4−2a−1b4s2+ab42b4,where s is the coordinate scaled to center the shape around x=0. For a plot of the profiles, see Figure 3.Figure 3
Variances of the tower configuration near the manufacturing defect: components of the first (left) and third (right) eigenmodes forK=1. Solution computed with p=4 and L=3 (165 collocation points). (a) Component μ of the first eigenmode. (b) Component μ of the third eigenmode. (c) Component v of the first eigenmode. (d) Component v of the third eigenmode. (e) Component w of the first eigenmode. (f) Component w of the third eigenmode.
(a)
(b)
(c)
(d)
(e)
(f)In each experiment, the dimension of the stochastic space isM=4, and the random vector ξ∈Γ is obtained by scaling the geometric parameters to the interval −1,1. The base of the tower is clamped, that is, fully kinematically constrained. The p-version of FEM is used with p=4 resulting in a linear system with 14247 d.o.f. The multiindex set is AL, with L=3. Per collocation point, the time spent is approximately 40 seconds; 10 seconds for solution and 30 seconds for obtaining statistics.In Tables1 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. This 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.Table 1
Statistics and wavenumbers of the eight smallest (in expected value) eigenvalues of the tower configuration as defined by the wavenumberK.
E
ω
2
Var
ω
2
Var
ω
2
K
79.7358
0.00175638
0.0419092
1
2055.46
218.179
14.7709
3
3449.66
2194.93
46.8501
2
4548.42
989.703
31.4595
0
4572.40
30.0598
5.48268
4
5958.11
2550.51
50.5026
0
6664.07
446.766
21.1368
4
7296.49
3535.53
59.4603
3Table 2
The eigenmodes evaluated at two different points inΓ (solution computed with p=4 and L=3 (165 collocation points)).
ξ
=
−
0.9
,
0,0,0
ξ
=
0.9
,
0,0,0
ω
2
K
ω
2
K
79.7681
1
79.6528
1
2072.60
3
2028.05
3
3503.49
2
3361.97
2
4578.79
4
4485.45
0
4581.93
0
4562.30
4
6013.69
0
5861.08
0
6688.77
4
6625.28
4
7365.13
3
7185.36
3The statistics presented in Tables1 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. The 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 Var⋅/S. From Table 3, we see that the errors between the two solutions fall well within the tolerances of the Monte Carlo results.Table 3
The errors between eigenvalues computed by collocation (ωC2) and Monte Carlo (ωMC2) methods. The eight smallest (in expected value) eigenvalues of the tower configuration.
ε
mean
ε
stdev
ε
MC
0.000453485
0.000637457
0.00130513
0.0120954
0.0580925
0.465259
0.0444178
0.205332
1.47504
0.0921339
0.163266
0.989675
0.00269920
0.0209490
0.172715
0.0590391
0.293009
1.58777
0.00789916
0.0910774
0.665526
0.0626177
0.263587
1.87197
The error measures are as follows: error in expected valueεmean=EωMC2−EωC2, error in standard deviation εstdev=VarωMC2−VarωC2, and Monte Carlo standard error εMC=VarωMC2/S.The effect of the manufacturing imperfection is illustrated for the eigenmodes associated withK=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. This is of course what we would expect.


### 5.2. 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,D1, the slope of the cone, ζ, and Young’s modulus. The 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 L1=10π and the distance to the cone is L2=40π/7, and for the short one, the lengths are simply scaled by 10. For a schematic of the construction, see Figure 4.Figure 4
Cylinder-cylinder-cone schematic. The cylinder diameterR1=2, L1, L2, and C1 are deterministic: L1+L2+C1=20π; the diameter of the hole D1 is random, and the cone diameter R2 is random since it depends on the random slope of the cone.Now the computational domain is periodic in the angular direction:(46)Ω^t=−10π,10π×0,2πDisk0,0,D1/2.Notice that the change in the slope does not affect the computational domain.The random variables and their ranges used in different experiments are(47)Long  ShortD1/2π/36,7π/36π/84,π/12ζ0,tan1/100,tan1/10ξ3−1,1−1,1ξ4−1,1−1,1,where ζ is the slope of the cone. Finally, the parameters ξ3 and ξ4 correspond Karhunen–Loéve basis functions, i.e., terms E3x=sinx and E4x=sin2x 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. The p-version of FEM is used with p=3 resulting in a linear system with 41280 d.o.f. The multiindex set is AL 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.


#### 5.2.1. Free Vibration

Let us first examine the lowest modes of the two structures. Interestingly, when comparing the two subfigures of Figure5, 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.Figure 5
Shell configuration:t=1/100; lowest mode; normalized transverse deflection w; rainbow colour scale. (a) Long cylinder with T-junction and cone junction; frequency = 1 Hz; x∈−10π,10π; cone: 40π/7,10π. (b) Short cylinder with T-junction and cone junction; frequency = 10 Hz; x∈−π,π; cone: 4π/7,π.
(a)
(b)In Table4, 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.Table 4
Statistics of the four smallest eigenvalues of the long (a) and short (b) cylinder with junctions (solution computed withp=3 and L=2 (45 collocation points)).
E
ω
2
Var
ω
2
Var
ω
2
(a) Long cylinder
3.12173
3.70583
1.92505
15.1519
57.9719
7.61393
25.7390
161.805
12.7203
88.1692
239.444
15.4740
(b) Short cylinder
58.6626
335.750
18.3235
240.370
13908.9
117.936
655.337
4200.64
64.8123
1137.43
10215.8
101.073Table 5
The errors between eigenvalues computed by collocation (ωC2) and Monte Carlo (ωMC2) methods. The four smallest eigenvalues of the long (a) and short (b) cylinders.
ε
mean
ε
stdev
ε
MC
(a) Long cylinder
0.0652920
0.0150364
0.0604000
0.287854
0.187262
0.234852
0.430142
0.00845208
0.401983
0.503874
0.199000
0.483037
(b) Short cylinder
0.0134338
0.473668
0.564461
0.163070
1.13330
3.76530
0.422822
1.92808
1.98858
1.12535
9.27454
3.48951
The error measures are as follows: error in expected valueεmean=EωMC2−EωC2, error in standard deviation εstdev=VarωMC2−VarωC2, and Monte Carlo standard error εMC=VarωMC2/S.Statistics of the first eigenmode for the long cylinder with junctions have been presented in Figure6. The long-range layer is clearly visible in the second, third, and fifth components. The effect of the cone junction is distinctly manifested in the statistics of the fourth component. The standard deviations of the two smallest eigenvalues are approximately 62 and 50 percent of the respective expected values. The ratio of the expected values is approximately 0.21, and thus, the eigenvalues at least appear to be well separated.Figure 6
Long cylinder with junctions: expected value and variance of the first eigenmode. Solution computed withp=3 and L=2 (45 collocation points). (a) Expected value (A) and variance (B) of the component μ. (b) Expected value (A) and variance (B) of the component v. (c) Expected value (A) and variance (B) of the component w. (d) Expected value (A) and variance (B) of the component θ. (e) Expected value (A) and variance (B) of the component ψ.
(a)
(b)
(c)
(d)
(e)Statistics of the first eigenmode for the short cylinder with junctions have been presented in Figure7. In this case, we do not observe a separate long-range layer. The effect of the cone junction is never the less apparent in the statistics of the fourth component. The standard deviations of the two smallest eigenvalues are now approximately 31 and 49 percent of the respective expected values. The ratio of the expected values is approximately 0.24, and thus, the eigenvalues again appear to be well separated.Figure 7
Short cylinder with junctions: expected value and variance of the first eigenmode. Solution computed withp=3 and L=2 (45 collocation points). (a) Expected value (A) and variance (B) of the component μ. (b) Expected value (A) and variance (B) of the component v. (c) Expected value (A) and variance (B) of the component w. (d) Expected value (A) and variance (B) of the component θ. (e) Expected value (A) and variance (B) of the component ψ.
(a)
(b)
(c)
(d)
(e)In Figure8, the short-range effects of the T-junction have been highlighted. The statistics of the fifth component have been shown in vicinity of the junction for both cylinders. The 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.Figure 8
Statistics of the fifth field component ψ near the T-junction. (a) Long cylinder: expected value (A) and variance (B). (b) Short cylinder: expected value (A) and variance (B).
(a)
(b)


#### 5.2.2. Frequency Response Analysis

For the frequency response analysis of the short cylinder, the load is chosen to beFx,y=1000cosy/2N, and the angular frequency range is taken to be ω∈2π5,10,…,200Hz. 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 pointsXGL and the second-order Smolyak–Clenshaw–Curtis quadrature points XCC and then consider their respective subsets with 25 points defined by setting(48)XGL=ξ1,…,ξ4∈XGLξ1=0,XCC=ξ1,…,ξ4∈XCCξ1=0,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.Let us first consider sparse gridsXGL and XCC. In Figure 9, we show two frequency response graphs where the expected values of the quantity of interest, the maximal transverse deflection wmax, are shown along standard deviations. It is clear that the frequency response is sensitive to perturbations to the parameters. The 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. This is illustrated in Figure 10. The maximal transverse deflection wmax at any given frequency depends on how well the tail of the distribution is approximated.Figure 9
Frequency response: comparison of sparse grids:XGL vs XCC; maximal transverse deflection wmax vs frequency; dashed line represents the added standard deviation. (a) Gauss XGL p = 2, 41 points. (b) Clenshaw–Curtis XCC p = 2, 41 points. (c) Combined.
(a)
(b)
(c)Figure 10
Frequency response distribution: Transverse deflections over the combined sparse gridsXGL and XCC. (a) Frequency = 65 Hz. (b) Frequency = 85 Hz.
(a)
(b)Let us next move to the partial sparse gridsX∈XCC,XGL. The results are displayed in Figure 11. The 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.Figure 11
Frequency response: comparison of partial sparse grids:XGL vs XCC; maximal transverse deflection wmax vs frequency; dashed line represents the added standard deviation. (a) Gauss XGL p = 2, 25 points. (b) Clenshaw–Curtis XCC p = 2, 25 points. (c) Combined.
(a)
(b)
(c)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.
This 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.


## 5.1. Tower Configuration: Free Vibration

Our first example is an idealized wind turbine tower (Figure2). 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.Figure 2
Tower configuration: Computational domain;B indicates the location of the assumed manufacturing imperfections at x=0. (a) Shell-solid model; diameter of the base is random. (b) Computational domain. (c) Lowest mode; normalized transverse detection w; rainbow colour scale.
(a)
(b)
(c)The tower is assembled with four parts: a thin shell of thicknesst=1/100 and vertical length of 11 units (x∈−1,10), a top ring of height C1=2/5 and width 3/4, and two base rings, inner and outer ones, of height D1=2/5 and random widths W1 and W2, respectively. An additional source of randomness comes from varying thickness of the shell near the top ring, centered at x=0. Here, the varying thickness can be interpreted as manufacturing imperfection or damage to the structure. The tower configuration is detailed also in Figure 2.Let us define the computational domainΩ^ to be a union of five parameterized domains:(42)Ω^=∪i=15Ω^i,where(43)Ω^1=x,r+η−1≤x≤−1−C1,−3/4−t/2≤η≤t/2,Ω^2=x,r+η−1−C1≤x≤−b,−t/2≤η≤t/2,Ω^3=x,r+η−b≤x≤b,−ϕbx,t;a,b≤η≤ϕbx,t;a,b,Ω^4=x,r+ηb≤x≤10−D1,−t/2≤η≤t/2,Ω^5=x,r+η10−D1≤x≤10,−t/2−W1≤η≤t/2+W2.The random variables and their ranges used in different experiments are (Figure2)(44)TowerW11/4,3/4W21/4,5/4a0,1/2b3/40,1/8,where a and b are parameters in the shape modeling the imperfection taken to be symmetric on the inner and outer surfaces:(45)ϕs,t;a,b=ta−1b4s4−2a−1b4s2+ab42b4,where s is the coordinate scaled to center the shape around x=0. For a plot of the profiles, see Figure 3.Figure 3
Variances of the tower configuration near the manufacturing defect: components of the first (left) and third (right) eigenmodes forK=1. Solution computed with p=4 and L=3 (165 collocation points). (a) Component μ of the first eigenmode. (b) Component μ of the third eigenmode. (c) Component v of the first eigenmode. (d) Component v of the third eigenmode. (e) Component w of the first eigenmode. (f) Component w of the third eigenmode.
(a)
(b)
(c)
(d)
(e)
(f)In each experiment, the dimension of the stochastic space isM=4, and the random vector ξ∈Γ is obtained by scaling the geometric parameters to the interval −1,1. The base of the tower is clamped, that is, fully kinematically constrained. The p-version of FEM is used with p=4 resulting in a linear system with 14247 d.o.f. The multiindex set is AL, with L=3. Per collocation point, the time spent is approximately 40 seconds; 10 seconds for solution and 30 seconds for obtaining statistics.In Tables1 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. This 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.Table 1
Statistics and wavenumbers of the eight smallest (in expected value) eigenvalues of the tower configuration as defined by the wavenumberK.
E
ω
2
Var
ω
2
Var
ω
2
K
79.7358
0.00175638
0.0419092
1
2055.46
218.179
14.7709
3
3449.66
2194.93
46.8501
2
4548.42
989.703
31.4595
0
4572.40
30.0598
5.48268
4
5958.11
2550.51
50.5026
0
6664.07
446.766
21.1368
4
7296.49
3535.53
59.4603
3Table 2
The eigenmodes evaluated at two different points inΓ (solution computed with p=4 and L=3 (165 collocation points)).
ξ
=
−
0.9
,
0,0,0
ξ
=
0.9
,
0,0,0
ω
2
K
ω
2
K
79.7681
1
79.6528
1
2072.60
3
2028.05
3
3503.49
2
3361.97
2
4578.79
4
4485.45
0
4581.93
0
4562.30
4
6013.69
0
5861.08
0
6688.77
4
6625.28
4
7365.13
3
7185.36
3The statistics presented in Tables1 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. The 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 Var⋅/S. From Table 3, we see that the errors between the two solutions fall well within the tolerances of the Monte Carlo results.Table 3
The errors between eigenvalues computed by collocation (ωC2) and Monte Carlo (ωMC2) methods. The eight smallest (in expected value) eigenvalues of the tower configuration.
ε
mean
ε
stdev
ε
MC
0.000453485
0.000637457
0.00130513
0.0120954
0.0580925
0.465259
0.0444178
0.205332
1.47504
0.0921339
0.163266
0.989675
0.00269920
0.0209490
0.172715
0.0590391
0.293009
1.58777
0.00789916
0.0910774
0.665526
0.0626177
0.263587
1.87197
The error measures are as follows: error in expected valueεmean=EωMC2−EωC2, error in standard deviation εstdev=VarωMC2−VarωC2, and Monte Carlo standard error εMC=VarωMC2/S.The effect of the manufacturing imperfection is illustrated for the eigenmodes associated withK=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. This is of course what we would expect.


## 5.2. 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,D1, the slope of the cone, ζ, and Young’s modulus. The 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 L1=10π and the distance to the cone is L2=40π/7, and for the short one, the lengths are simply scaled by 10. For a schematic of the construction, see Figure 4.Figure 4
Cylinder-cylinder-cone schematic. The cylinder diameterR1=2, L1, L2, and C1 are deterministic: L1+L2+C1=20π; the diameter of the hole D1 is random, and the cone diameter R2 is random since it depends on the random slope of the cone.Now the computational domain is periodic in the angular direction:(46)Ω^t=−10π,10π×0,2πDisk0,0,D1/2.Notice that the change in the slope does not affect the computational domain.The random variables and their ranges used in different experiments are(47)Long  ShortD1/2π/36,7π/36π/84,π/12ζ0,tan1/100,tan1/10ξ3−1,1−1,1ξ4−1,1−1,1,where ζ is the slope of the cone. Finally, the parameters ξ3 and ξ4 correspond Karhunen–Loéve basis functions, i.e., terms E3x=sinx and E4x=sin2x 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. The p-version of FEM is used with p=3 resulting in a linear system with 41280 d.o.f. The multiindex set is AL 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.


### 5.2.1. Free Vibration

Let us first examine the lowest modes of the two structures. Interestingly, when comparing the two subfigures of Figure5, 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.Figure 5
Shell configuration:t=1/100; lowest mode; normalized transverse deflection w; rainbow colour scale. (a) Long cylinder with T-junction and cone junction; frequency = 1 Hz; x∈−10π,10π; cone: 40π/7,10π. (b) Short cylinder with T-junction and cone junction; frequency = 10 Hz; x∈−π,π; cone: 4π/7,π.
(a)
(b)In Table4, 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.Table 4
Statistics of the four smallest eigenvalues of the long (a) and short (b) cylinder with junctions (solution computed withp=3 and L=2 (45 collocation points)).
E
ω
2
Var
ω
2
Var
ω
2
(a) Long cylinder
3.12173
3.70583
1.92505
15.1519
57.9719
7.61393
25.7390
161.805
12.7203
88.1692
239.444
15.4740
(b) Short cylinder
58.6626
335.750
18.3235
240.370
13908.9
117.936
655.337
4200.64
64.8123
1137.43
10215.8
101.073Table 5
The errors between eigenvalues computed by collocation (ωC2) and Monte Carlo (ωMC2) methods. The four smallest eigenvalues of the long (a) and short (b) cylinders.
ε
mean
ε
stdev
ε
MC
(a) Long cylinder
0.0652920
0.0150364
0.0604000
0.287854
0.187262
0.234852
0.430142
0.00845208
0.401983
0.503874
0.199000
0.483037
(b) Short cylinder
0.0134338
0.473668
0.564461
0.163070
1.13330
3.76530
0.422822
1.92808
1.98858
1.12535
9.27454
3.48951
The error measures are as follows: error in expected valueεmean=EωMC2−EωC2, error in standard deviation εstdev=VarωMC2−VarωC2, and Monte Carlo standard error εMC=VarωMC2/S.Statistics of the first eigenmode for the long cylinder with junctions have been presented in Figure6. The long-range layer is clearly visible in the second, third, and fifth components. The effect of the cone junction is distinctly manifested in the statistics of the fourth component. The standard deviations of the two smallest eigenvalues are approximately 62 and 50 percent of the respective expected values. The ratio of the expected values is approximately 0.21, and thus, the eigenvalues at least appear to be well separated.Figure 6
Long cylinder with junctions: expected value and variance of the first eigenmode. Solution computed withp=3 and L=2 (45 collocation points). (a) Expected value (A) and variance (B) of the component μ. (b) Expected value (A) and variance (B) of the component v. (c) Expected value (A) and variance (B) of the component w. (d) Expected value (A) and variance (B) of the component θ. (e) Expected value (A) and variance (B) of the component ψ.
(a)
(b)
(c)
(d)
(e)Statistics of the first eigenmode for the short cylinder with junctions have been presented in Figure7. In this case, we do not observe a separate long-range layer. The effect of the cone junction is never the less apparent in the statistics of the fourth component. The standard deviations of the two smallest eigenvalues are now approximately 31 and 49 percent of the respective expected values. The ratio of the expected values is approximately 0.24, and thus, the eigenvalues again appear to be well separated.Figure 7
Short cylinder with junctions: expected value and variance of the first eigenmode. Solution computed withp=3 and L=2 (45 collocation points). (a) Expected value (A) and variance (B) of the component μ. (b) Expected value (A) and variance (B) of the component v. (c) Expected value (A) and variance (B) of the component w. (d) Expected value (A) and variance (B) of the component θ. (e) Expected value (A) and variance (B) of the component ψ.
(a)
(b)
(c)
(d)
(e)In Figure8, the short-range effects of the T-junction have been highlighted. The statistics of the fifth component have been shown in vicinity of the junction for both cylinders. The 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.Figure 8
Statistics of the fifth field component ψ near the T-junction. (a) Long cylinder: expected value (A) and variance (B). (b) Short cylinder: expected value (A) and variance (B).
(a)
(b)


### 5.2.2. Frequency Response Analysis

For the frequency response analysis of the short cylinder, the load is chosen to beFx,y=1000cosy/2N, and the angular frequency range is taken to be ω∈2π5,10,…,200Hz. 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 pointsXGL and the second-order Smolyak–Clenshaw–Curtis quadrature points XCC and then consider their respective subsets with 25 points defined by setting(48)XGL=ξ1,…,ξ4∈XGLξ1=0,XCC=ξ1,…,ξ4∈XCCξ1=0,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.Let us first consider sparse gridsXGL and XCC. In Figure 9, we show two frequency response graphs where the expected values of the quantity of interest, the maximal transverse deflection wmax, are shown along standard deviations. It is clear that the frequency response is sensitive to perturbations to the parameters. The 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. This is illustrated in Figure 10. The maximal transverse deflection wmax at any given frequency depends on how well the tail of the distribution is approximated.Figure 9
Frequency response: comparison of sparse grids:XGL vs XCC; maximal transverse deflection wmax vs frequency; dashed line represents the added standard deviation. (a) Gauss XGL p = 2, 41 points. (b) Clenshaw–Curtis XCC p = 2, 41 points. (c) Combined.
(a)
(b)
(c)Figure 10
Frequency response distribution: Transverse deflections over the combined sparse gridsXGL and XCC. (a) Frequency = 65 Hz. (b) Frequency = 85 Hz.
(a)
(b)Let us next move to the partial sparse gridsX∈XCC,XGL. The results are displayed in Figure 11. The 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.Figure 11
Frequency response: comparison of partial sparse grids:XGL vs XCC; maximal transverse deflection wmax vs frequency; dashed line represents the added standard deviation. (a) Gauss XGL p = 2, 25 points. (b) Clenshaw–Curtis XCC p = 2, 25 points. (c) Combined.
(a)
(b)
(c)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.
This 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.


## 5.2.1. Free Vibration

Let us first examine the lowest modes of the two structures. Interestingly, when comparing the two subfigures of Figure5, 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.Figure 5
Shell configuration:t=1/100; lowest mode; normalized transverse deflection w; rainbow colour scale. (a) Long cylinder with T-junction and cone junction; frequency = 1 Hz; x∈−10π,10π; cone: 40π/7,10π. (b) Short cylinder with T-junction and cone junction; frequency = 10 Hz; x∈−π,π; cone: 4π/7,π.
(a)
(b)In Table4, 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.Table 4
Statistics of the four smallest eigenvalues of the long (a) and short (b) cylinder with junctions (solution computed withp=3 and L=2 (45 collocation points)).
E
ω
2
Var
ω
2
Var
ω
2
(a) Long cylinder
3.12173
3.70583
1.92505
15.1519
57.9719
7.61393
25.7390
161.805
12.7203
88.1692
239.444
15.4740
(b) Short cylinder
58.6626
335.750
18.3235
240.370
13908.9
117.936
655.337
4200.64
64.8123
1137.43
10215.8
101.073Table 5
The errors between eigenvalues computed by collocation (ωC2) and Monte Carlo (ωMC2) methods. The four smallest eigenvalues of the long (a) and short (b) cylinders.
ε
mean
ε
stdev
ε
MC
(a) Long cylinder
0.0652920
0.0150364
0.0604000
0.287854
0.187262
0.234852
0.430142
0.00845208
0.401983
0.503874
0.199000
0.483037
(b) Short cylinder
0.0134338
0.473668
0.564461
0.163070
1.13330
3.76530
0.422822
1.92808
1.98858
1.12535
9.27454
3.48951
The error measures are as follows: error in expected valueεmean=EωMC2−EωC2, error in standard deviation εstdev=VarωMC2−VarωC2, and Monte Carlo standard error εMC=VarωMC2/S.Statistics of the first eigenmode for the long cylinder with junctions have been presented in Figure6. The long-range layer is clearly visible in the second, third, and fifth components. The effect of the cone junction is distinctly manifested in the statistics of the fourth component. The standard deviations of the two smallest eigenvalues are approximately 62 and 50 percent of the respective expected values. The ratio of the expected values is approximately 0.21, and thus, the eigenvalues at least appear to be well separated.Figure 6
Long cylinder with junctions: expected value and variance of the first eigenmode. Solution computed withp=3 and L=2 (45 collocation points). (a) Expected value (A) and variance (B) of the component μ. (b) Expected value (A) and variance (B) of the component v. (c) Expected value (A) and variance (B) of the component w. (d) Expected value (A) and variance (B) of the component θ. (e) Expected value (A) and variance (B) of the component ψ.
(a)
(b)
(c)
(d)
(e)Statistics of the first eigenmode for the short cylinder with junctions have been presented in Figure7. In this case, we do not observe a separate long-range layer. The effect of the cone junction is never the less apparent in the statistics of the fourth component. The standard deviations of the two smallest eigenvalues are now approximately 31 and 49 percent of the respective expected values. The ratio of the expected values is approximately 0.24, and thus, the eigenvalues again appear to be well separated.Figure 7
Short cylinder with junctions: expected value and variance of the first eigenmode. Solution computed withp=3 and L=2 (45 collocation points). (a) Expected value (A) and variance (B) of the component μ. (b) Expected value (A) and variance (B) of the component v. (c) Expected value (A) and variance (B) of the component w. (d) Expected value (A) and variance (B) of the component θ. (e) Expected value (A) and variance (B) of the component ψ.
(a)
(b)
(c)
(d)
(e)In Figure8, the short-range effects of the T-junction have been highlighted. The statistics of the fifth component have been shown in vicinity of the junction for both cylinders. The 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.Figure 8
Statistics of the fifth field component ψ near the T-junction. (a) Long cylinder: expected value (A) and variance (B). (b) Short cylinder: expected value (A) and variance (B).
(a)
(b)


## 5.2.2. Frequency Response Analysis

For the frequency response analysis of the short cylinder, the load is chosen to beFx,y=1000cosy/2N, and the angular frequency range is taken to be ω∈2π5,10,…,200Hz. 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 pointsXGL and the second-order Smolyak–Clenshaw–Curtis quadrature points XCC and then consider their respective subsets with 25 points defined by setting(48)XGL=ξ1,…,ξ4∈XGLξ1=0,XCC=ξ1,…,ξ4∈XCCξ1=0,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.Let us first consider sparse gridsXGL and XCC. In Figure 9, we show two frequency response graphs where the expected values of the quantity of interest, the maximal transverse deflection wmax, are shown along standard deviations. It is clear that the frequency response is sensitive to perturbations to the parameters. The 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. This is illustrated in Figure 10. The maximal transverse deflection wmax at any given frequency depends on how well the tail of the distribution is approximated.Figure 9
Frequency response: comparison of sparse grids:XGL vs XCC; maximal transverse deflection wmax vs frequency; dashed line represents the added standard deviation. (a) Gauss XGL p = 2, 41 points. (b) Clenshaw–Curtis XCC p = 2, 41 points. (c) Combined.
(a)
(b)
(c)Figure 10
Frequency response distribution: Transverse deflections over the combined sparse gridsXGL and XCC. (a) Frequency = 65 Hz. (b) Frequency = 85 Hz.
(a)
(b)Let us next move to the partial sparse gridsX∈XCC,XGL. The results are displayed in Figure 11. The 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.Figure 11
Frequency response: comparison of partial sparse grids:XGL vs XCC; maximal transverse deflection wmax vs frequency; dashed line represents the added standard deviation. (a) Gauss XGL p = 2, 25 points. (b) Clenshaw–Curtis XCC p = 2, 25 points. (c) Combined.
(a)
(b)
(c)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.
This 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.


## 6. 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. The 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.



---

*Source: 5817940-2018-12-02.xml*